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

基于中心差分格式的冻结过程模拟 季节性冷冻土壤是一种特殊的土壤水系统。土壤冻结过程非常复杂,伴随物理、化学、力学现象和过程,包括水分、热量迁移、水分的相变和盐分的积聚。国内对冻土水热耦合迁移的数值模拟问题研究较晚,起步于20世纪80年代末期,但发展较快。雷志栋等采用Harlan模型模拟了水平、垂直土柱的冻结过程,并定性分析了土壤初始含水率对土壤冻胀量的影响;尚松浩等模拟了室内土柱冻结试验和田间的水热状况,改进了求解方法使计算速度加快;郑秀清对山西省汾河灌区土壤冻融过程水热动态进行仿真分析;李扬基于非饱和土多孔介质理论的水热迁移耦合模型,采用有限元法求解,模拟了长春市季节性冻土水分迁移,经试验验证,效果较好。鉴此,本文采用Matlab7.1软件编制计算程序,应用中心差分格式并用线性化迭代法求解冻土一维垂直水热耦合迁移数学模型,对沈阳农业大学水利学院综合试验基地土壤在冻结过程中的水、热动态进行了模拟,为揭示季节性冻土水热耦合迁移状况提供依据。 1 数据采集和参数测量 1.1 温度测区温度和含水率测量 2007年11月26日～12月8日在沈阳农业大学水利学院综合试验基地采集野外数据。测区为一个由不透水保温板围成的3 m×3 m×26 m(长×宽×高)试验区。测区埋设TRM-ZS1气象及生态环境监测系统探头,可收集气象数据及地表至地下1.5 m处土壤各剖面的温度数据;1.5 m深的冻土器可测量冻深;在固定位置安放中子仪测管和TDR测管测量总含水率和未冻水含水率;测区土壤为潮棕壤土,0～20 cm的土样干容重为1.39 g/cm3,饱和含水率0.48 cm3/cm3;20～150 cm的土样干容重为1.47 g/cm3,饱和含水率0.45 cm3/cm3。 1.2 土壤水吸力与d的关系 土壤非饱和扩散率D(θ)(单位:cm2/min)通过室内水平土柱吸渗法确定。 当土层深度为0～20 cm时: {D(θ)=0.0984e7.482θθ≤0.38D(θ)=0.0005e22.144θθ＞0.38(1) 当土层深度为20～150 cm时: D(θ)=0.000 6e18.422θ(2) 土壤非饱和导水率K(θ)(单位:cm/min)与D(θ)土壤水吸力s(θ)的关系为: Κ(θ)=-D(θ)dθds(3) s(θ)(单位:cm)由压力仪法测定。当土层深度为0～20 cm时: s(θ)={5×108e-41.379θθ＞0.336166374θ2-17403θ+4554θ≤0.336(4) 当土层深度为20～150 cm时: s(θ)={2×107e-35.119θθ＞02-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.5 w+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.244 5、0.562 0、2.180 0、0.025 0。 土壤体积比热容具有可加性,是土壤中各组成物热容量之和: cv=cvsxs+cvwxw+cvlxvl+cvaxva(11) 式中,xs、xw、xvI、xva分别为单位体积土壤矿物、水、冰、空气所占体积;cvs、cvw、cvI、cva分别为相应组成物的体积比热容,J/(cm3·K),数值分别为2.223、4.187、1.730、0.0