一维非饱和土热一湿耦合的数值解

吴礼舟，黄润秋

(成都理工大学地质灾害防治与地质环境保护重点实验室，成都610059)

地表浅层土体是工程活动的主要对象，由于受到诸多因素的影响，土体中的温度场、水分场等均是变化的，而且相互影响。土体中的温度变化会引起水分迁移使含水量变化，含水量的变化又会引起土的导热系数等发生变化，从而影响传热过程及温度分布。因此，研究浅层土体水(水分场)、热(温度场)问题时，必须考虑其相互影响[1-2]。

非饱和土耦合现象具有较广的工程背景，因而引起众学者的兴趣[2-5]。岩土深处核废料处置过程中的放热现象亦是近年来人们关注的问题之一。路基中的温度场、水分场是动态变化的。路基中热量的差异和改变引起水分的迁移与转化，同时，路基中的水通过改变土的热特性来影响土壤的温度[3]。白冰指出在给出非饱和介质热能-水分传输的耦合质量控制方程和基于Fourier热传导定律的热能平衡方程的基础上，对热能传输及水分迁移的基本特征和机理进行了分析[4]。Wu L Z利用傅立叶积分变换开展了非饱和土变形-渗流耦合的解析研究[5]。张玉军从水-气-热三场耦合的角度上分析了孔隙介质中各组分的渗流与变形问题[6]。

基于非饱和土水-热耦合的理论模型，运用美国的商业软件FlexPDE软件分析一维非饱和土热-湿耦合解，计算非饱和土的温度场和水分场。试图分析非饱和土热-湿耦合作用下非饱和土水分迁移及温度分布的特点与规律。

1 理论模型

为了较好的分析非饱和土热-湿耦合效应，作出以下假设:非饱和土是均匀，介质不发生变形。该理论模型的基本假定:(1)土体均匀连续，且为各向同性体;(2)非饱和土中水分迁移均以液态形式进行，不考虑水的固相;(3)土体为弹性体，且已完全固结;(4)非饱和流体符合Darcy定律。

非饱和土介质中的热-湿耦合的控制方程为[7]:

式中 :θ——体积含水量 ;T ——温度;t——时间;c——比热 ;ρ——密度 ;ρ c——体积热容;λ——热传导率;L——潜热;Dθv——等温下的气扩散系数;DT——热湿扩散率;Dθ——等温条件的扩散系数;K——非饱和水力传导系数。

对于一维问题，耦合控制方程，式(1)和式(2)可写成:

不考虑重力项，式(4)右边第3项为0，即式(4)变为:

若不考虑温度对含水量的影响，一维的水分控制方程为:

不考虑重力项，式(6)写成:

假设含水量和渗透系数是压力水头ψ的函数:

式中 :θs——饱和含水量;Ks——饱和状态时的渗透系数;α——去饱和系数。

其扩散系数为:

初始条件表达为

下边界条件为

上边界条件为

2 算例分析

由于非饱和土热-湿耦合的解析解求解很困难，这里采用了商业软件FlexPDE(PDE-Partial differential equation)。FlexPDE使用有限元方法获得数值解。它能用于分析求解各类偏微分方程(组)，通过建立耦合的数学模型，并能有效的分析耦合各场的分布状态 ，程序编写遵照 coordinates、variables、definitions、equations、boundaries、monitors、time 等流程 。其中variables是定义方程组的变量，equations是描述耦合方程组，time是运行求解。求解耦合问题非常方便，可灵活分析用户建立的耦合模型。

这里对有限长度的一维非饱和土热-湿耦合问题进行分析 。l=50 cm ，θs=0.4，α=0.01/cm ，Ks=0.4×10-4cm/s。初始含水量θi=0.2，初始温度 Ti=10℃;z=0 处 θ0=0.2，T0=10℃;z=l处均为Neumann 边界条件 ，∂θ/∂z=0.5×10-4cm-1，∂T/∂z=0.01℃/cm，其他相关系数见表1[7]。

表1 非饱和土热-湿耦合的物理参数值

图1 非饱和土热-湿耦合中温度的变化

图1显示的是非饱和土热-湿耦合中不同时刻下的温度分布。图2为非饱和土热湿耦合过程中的含水量随时间的变化。图1和图2中，地表边界发生热量交换，同时出现降雨入渗。图1表明温度随时间在增大，从图2可看出，含水量变大，地表的含水量变化更为明显。图1和图2采用的耦合模型来源于式(3)和式(4)。

图3为非饱和土热-湿耦合与非耦合情况下含水量随时间的分布。虚线是没有考虑温度对含水量影响的分布图，而实线表示考虑非饱和土热-湿耦合效应。从图3可以看出，温度对水分迁移是有影响的。在地表，不考虑热-湿耦合时的水分迁移比考虑耦合效应时快，而在土体深处，没有考虑热-湿耦合的水分运移比考虑耦合效应的快。

图2 非饱和土热-湿耦合中含水量的变化

图3 耦合与非耦合情况下含水量的分布

图4 耦合中重力对含水量的影响

图5 等温下的气扩散系数对耦合作用下温度的影响

图4展示的是非饱和热-湿耦合中重力对含水量分布的影响。图中线条表示考虑了重力势，非饱和土热-湿耦合模型来源于式(4);图中点代表没有考虑重力势，来源于式(5)。在入渗初期t=200 s重力对非饱和土热-湿没有任何影响，随着入渗继续，重力对含水量分布显示有非常微弱的影响，如t=1000 s时，考虑重力势时的含水量微微大于没有考虑重力时的含水量。在分析非饱和土湿-热耦合时，入渗较小时，重力对含水量分布的影响基本可忽略不计，印证了文献[7]的结论。

图6 等温下的气扩散系数对耦合作用下含水量的影响

图5和图6分别为等温下气扩散系数对耦合作用下温度和含水量分布的影响。无论是否考虑等温下的气扩散系数，温度随时间的分布几乎没有变化;同样，气扩散系数对耦合作用下含水量分布也没有任何影响。

3 结论

(1)当非饱和土耦合问题较为复杂时，采用解析方法显得困难时，采用FlexPDE软件分析非饱和土热-湿耦合问题是可行的，也是有效的。

(2)非饱和土热-湿耦合数值分析表明温度对水分迁移是有影响的，分析非饱和土水分迁移时应考虑非饱和土热-湿耦合效应。

(3)分析非饱和土湿-热耦合时，降雨入渗较小时，重力对含水量分布的影响甚微。等温下的气扩散系数对非饱和土湿-热耦合几乎没有任何影响，耦合分析时可略去其影响。

[1] Lee J，Kim J T，Chung I M ，et al.Analytical solution for heat and moisture diffusion in layered materials[J].Canadian Geotechnical Journal，2010 ，47(6):595-608.

[2] 王铁行，李宁，谢定义.土体水热力耦合问题研究意义、现状及建议[J].岩土力学，2005，26(3):488-493.

[3] 毛雪松，王秉纲，胡长顺，等.冻土路基水热迁移问题的理论模型及数值模拟[J].中外公路，2006，26(1):23-26.

[4] 白冰，刘大鹏.非饱和介质中热能传输及水分迁移的数值积分解[J].岩土力学，2006，27(12):2085-2089.

[5] Wu L Z，Zhang L M.Analytical solutions to 1D coupled water infiltration and deformation in unsaturated soils[J].International Journal for Numerical and Analytical Methods in Geomechanics，2009，33(6):773-790.

[6] 张玉军.核废料处置概念库近场热-水-应力耦合模型及数值分析[J].岩土力学，2007，28(1):17-22.

[7] Basha H A，Selvadurai A P S.Heat-induced moisture transport in the vicinity of a spherical heat source[J].InternationalJournalfor Numericaland Analytical Methods in Geomechanics，1998，22(12):969-981.