基于中心差分格式的冻结过程模拟

基于中心差分格式的冻结过程模拟

季节性冷冻土壤是一种特殊的土壤水系统。土壤冻结过程非常复杂,伴随物理、化学、力学现象和过程,包括水分、热量迁移、水分的相变和盐分的积聚。国内对冻土水热耦合迁移的数值模拟问题研究较晚,起步于20世纪80年代末期,但发展较快。雷志栋等采用Harlan模型模拟了水平、垂直土柱的冻结过程,并定性分析了土壤初始含水率对土壤冻胀量的影响;尚松浩等模拟了室内土柱冻结试验和田间的水热状况,改进了求解方法使计算速度加快;郑秀清对山西省汾河灌区土壤冻融过程水热动态进行仿真分析;李扬基于非饱和土多孔介质理论的水热迁移耦合模型,采用有限元法求解,模拟了长春市季节性冻土水分迁移,经试验验证,效果较好。鉴此,本文采用Matlab7.1软件编制计算程序,应用中心差分格式并用线性化迭代法求解冻土一维垂直水热耦合迁移数学模型,对沈阳农业大学水利学院综合试验基地土壤在冻结过程中的水、热动态进行了模拟,为揭示季节性冻土水热耦合迁移状况提供依据。1数据采集和参数测量1.1温度测区温度和含水率测量2007年11月26日～12月8日在沈阳农业大学水利学院综合试验基地采集野外数据。测区为一个由不透水保温板围成的3m×3m×26m(长×宽×高)试验区。测区埋设TRM-ZS1气象及生态环境监测系统探头,可收集气象数据及地表至地下1.5m处土壤各剖面的温度数据;1.5m深的冻土器可测量冻深;在固定位置安放中子仪测管和TDR测管测量总含水率和未冻水含水率;测区土壤为潮棕壤土,0～20cm的土样干容重为1.39g/cm3,饱和含水率0.48cm3/cm3;20～150cm的土样干容重为1.47g/cm3,饱和含水率0.45cm3/cm3。1.2土壤水吸力与d的关系土壤非饱和扩散率D(θ)(单位:cm2/min)通过室内水平土柱吸渗法确定。当土层深度为0～20cm时:{D(θ)=0.0984e7.482θθ≤0.38D(θ)=0.0005e22.144θθ＞0.38(1)当土层深度为20～150cm时:D(θ)=0.0006e18.422θ(2)土壤非饱和导水率K(θ)(单位:cm/min)与D(θ)土壤水吸力s(θ)的关系为:Κ(θ)=-D(θ)dθds(3)s(θ)(单位:cm)由压力仪法测定。当土层深度为0～20cm时:s(θ)={5×108e-41.379θθ＞0.336166374θ2-17403θ+4554θ≤0.336(4)当土层深度为20～150cm时:s(θ)={2×107e-35.119θθ＞0.30419721θ2-5749.3θ+2102.9θ≤0.304(5)1.3水土流失量i冻结区的扩散率和导水率很难测定,本文采用阻抗系数法。分别用未冻结区的扩散率和导水率除以阻抗系数I表示,I值取决于土壤中的含冰率:I=1010θI(6)DI(θ)=D(θ)/I(7)KI(θ)=K(θ)/I(8)式中,θ1为含水率,cm3/cm3。1.4土壤负温对未冻水含量的影响冻土中含有未冻水是由土壤自身性质决定。未冻水在冻土中可自由移动。影响其含量的因素很多,其中影响最大、起决定性作用的是土壤的负温。未冻水含量随土壤负温的增大而减小,每种土壤均有较固定的未冻水含量与负温的关系曲线。表1为根据现场试验结果拟合的未冻水含量与负温关系,拟合方程为:θ=aT-b(9)式中,θ为未冻水含水率,cm3/cm3;T为土层负温绝对值,℃;a、b为经验常数。1.5w+a/a的热容比估算土壤热导率Kh十分困难,Kh不仅取决于土壤中各组成物的比例,还与各组成物的形状等有关。本文采用的计算公式为:Κh=1σ/Κ1+(1-σ)/Κ2(10)当θ+θI<θs时:{Κ1=1X/λs+θ/λw+θΙ/λΙ+A/λaΚ2=Xλs+θΙλΙ+θλw+A/λa当θ+θI≥θs时:{Κ1=1X/λs+θ/λw+θΙ/λΙΚ2=Xλs+θΙλΙ+θλw式中,K1、K2分别为串、并联土层的热导率;σ为串、并联土层的体积比,一般取0.5;θs为饱和含水率;X、θ、θI、A分别为土颗粒、未冻水、冰、空气所占的体积;λs、λw、λI、λa分别为土颗粒、水、冰、空气的热传导系数,W/(m·K),数值分别为3.2445、0.5620、2.1800、0.0250。土壤体积比热容具有可加性,是土壤中各组成物热容量之和:cv=cvsxs+cvwxw+cvlxvl+cvaxva(11)式中,xs、xw、xvI、xva分别为单位体积土壤矿物、水、冰、空气所占体积;cvs、cvw、cvI、cva分别为相应组成物的体积比热容,J/(cm3·K),数值分别为2.223、4.187、1.730、0.0012。2数学模型2.1水热迁移的机制①假设土壤骨架无膨胀和收缩变形、均质且各向同性,水热迁移主要发生在垂向,可近似为一维问题。②土壤内部满足局部热力学平衡条件。③冻土中水分迁移以液态形式运动,认为冰固定不动,忽略气相的影响。2.2冰的密度和水率冻结条件下,水热耦合运移问题必须联立求解水分运动和热流两个基本方程。需求解θ、θI、Τ三个未知函数,因此,还须补充一个联系方程,即土壤中未冻水含水率θ和温度T的关系(式(9)),水分运动方程、热流方程、联系方程分别为:∂θ∂t=∂∂z[D(θ)∂θ∂z]-∂Κ(θ)∂z-ρΙρw∂θΙ∂t(12)cv∂Τ∂t=∂∂z[Κh∂Τ∂z]+LΙρΙ∂θΙ∂t(13)θ≤θm(T)(14)式中,ρI为冰的密度,0.917g/cm3;ρw为水的密度,1g/cm3;LI为冰的溶解潜热,一般取335J/cm3;θm为相应土体负温下可能的最大未冻水含量,cm3/cm3。2.3双边界段土壤水分通量与基本模式的关系冻结土壤的水热运动状况是由土壤的上下边界的水热条件决定的。热流方程上下边界条件为第一类边界条件,上下边界温度取实测值:Τ(z)={Τtt＞0,z=0Τbt＞0,z=150(15)水流差分方程上边界条件为第一类边界条件,地表含水率较低近似于风干含水率θc=0.05cm3/cm3;下边界有水分通量,为第二类边界条件,方程如下:θ(z)=θct>0,z=0(16)-D(θ)∂θ∂z+Κ(θ)=q(t)t＞0,z=150(17)由于地下水埋深约4m,下边界即1.5m处的土壤水分处于不饱和状态,其水分通量为:q(ˉz)=-Κ(ˉφ)(φ2-φ1Δz+1)(18)其中ˉz=(z1+z2)/2式中,q(ˉz)为边界处单位时间土壤水分通量,cm/h,q>0表明通量向下;q<0表明通量向上;z1、z2为土壤深度,cm;φ1、φ2为z1、z2处基质势;Δz为空间步长;Κ(φ¯)为z¯处非饱和土壤导水率与基质势的关系。基质势φ(单位:cm)与土壤水吸力s(单位:cm)的关系为:s=-φ(19)3中心差分格式由于冻融条件下,水热耦合运移模型的非线性及土壤的非均质和初始边界条件的复杂性,目前用解析方法求解很不理想,采用数值计算方法最有效。有限差分法以差商近似替代微商,将微分方程变成差分方程,组成可直接求解的代数方程组。选用的中心差分格式精度较高,但时间和空间步长较大时,会出现严重的数值振荡,因此应结合实例选取时间与空间步长。本文空间步长Δz、时间步长Δt分别取2.5cm、10min。计算步骤如下。步骤1定义初始含水率、温度,用初始值先代替下一时刻的含水率、温度,如θik+1=θik。步骤2求解水分方程组,获得未冻水含水率预测值。步骤3求解热流方程组,获得温度预测值。步骤4由式(14)校正未冻水含水率,若其值大于式(14),则超出部分为本时段增加的含冰率,即θΙik+1=θik+1-θm+θΙik。步骤5重复步骤2～4,直至前后两次计算的含水率与温度之差小于0.01为止。4模拟结果与实测值对比应用上述参数及数学模型仿真分析土壤冻结过程中的未冻水含水率及温度,结果如图1、2所示。由图可看出:①当土壤深度<30cm时,土壤温度变化剧烈,在周期性的环境温度和太阳辐射的作用下,土壤中各点的温度亦呈周期性变化。50、70、110cm处的土壤温度变化不明显,振幅趋于消失。②在冻结阶段,由于表层土壤未冻水含水率较小,所以在0～10cm处未冻水含水率变化不大;10～50cm未冻水含水率变小,主要原因为温度由正温降至负温使一部分未冻水变成冰,该部分冰的形成由冻结曲线确定;50cm以下深度的未冻水含水率变化不大。综上可知,土壤水分剖面与土壤温度剖面的模拟结果与实测值基本吻合,个别时段模拟结果虽与实测值有一定偏差,但两者差异总体不大。因此,本文应用的数学模型、参数和计算方法基本反映了野外开放条件下水热耦合迁移的基本过程,从而可知在冻结条件下土壤水分的迁移主要受温度控制,为研究季节性冻土的利用、水热耦合迁移的数值模拟提供了可靠的分析工具