补水条件下路基土体冻结过程中的温湿度分布研究

2022-05-17鲜少华李清鹏罗振先

水利与建筑工程学报 2022年2期

鲜少华，李清鹏，罗振先

(1.武汉市政工程设计研究院有限责任公司，湖北武汉 430023；2.河南大学土木建筑学院，河南开封 475004)

当季节冻土区地下水埋深较浅时，在冬季冻结过程中水分会在温度梯度的作用下迁移至路基土体中，引起冻胀变形，继而产生路面开裂破坏；融化过程中，大量的迁移水不易消散，造成土体强度降低，进而引起路面沉陷等病害。温度梯度的变化导致水分迁移的过程即是一种水热耦合的过程。水热耦合理论一直是国内外学者的研究热点，Harlan[1]通过对比分析非饱和土与部分冻结土中的水分迁移规律，最早建立了描述水热迁移的数学模型。随后，Taylor等[2]、Jame等[3]、Newman等[4]不同程度地改进和发展了水热迁移模型。Mu等[5]引入固体力学理论，在计算体应变时考虑了水分冻胀引起的土体应变，提出了水热力三场耦合模型。毛雪松等[6]利用水热力耦合模型计算了青藏公路唐南段路基温度场、湿度场和应力分布规律，得到的结果能够体现道路冻胀的主要诱因。张玉芝等[7]联合土体冻结过程中的热传导方程和热弹性力学理论建立了高铁路基的变形和内力计算模型，预测了路基服役期内在冻融作用下可能产生的冻胀量和融沉量。邰博文等[8]在冻土水热耦合方程的基础上通过土体冻胀率与含冰量的关系改进了土体冻胀模型，并利用寒区路基变形监测数据验证了模型的准确性。Zheng等[9]将一维冻胀模型扩展到了二维情形并考虑了冻结速率和约束应力的影响。Li等[10-11]在水热力三场耦合模型将冻融循环作用对土体的影响考虑了进去，并利用室内试验和现场监测进行了验证。褚志成等[12]借助温度场-应力场耦合模型研究了多年冻土区边坡在冻融作用下的局部稳定性。吕龙等[13]在Harlan水热耦合模型的基础上考虑水分迁移和冰水相变，认为水分迁移是导致基坑工程冬季冻胀破坏的主要原因之一。

综上，国内外诸多研究者在土体冻胀机理和水热耦合理论方面已经进行了大量研究工作，并且取得了一定的成果，但上述研究均未考虑有外界水分持续补给的情况。鉴于此，基于现有的水热耦合理论，针对季节冻土区路基填土的水热迁移问题，采用室内试验和数值计算的手段，研究有外界持续补水条件下土体内部的温度场和湿度场变化规律。

1 水-热耦合理论

1.1 基本假定

基于现有理论[14]，作以下基本假定：(1) 土体均匀连续，属于各向同性体；(2) 土中无盐分影响；(3) 冻结区和未冻区水分仅以液态形式进行迁移；(4) 土颗粒不可压缩；(5) 已冻土和未冻土为弹性体；(6) 无溶质迁移；(7) 水分迁移符合Darcy定律。

1.2 温度场方程

忽略对流传热和热辐射，考虑土体中冰水相变的热传导方程表示为：

(1)

1.3 水分场基本方程

利用水动力学模型来描述土体冻结过程中水分的迁移过程。对于饱和土和非饱和土，水分迁移方程均可表示如下：

(2)

式中：θw为水分总体积含量，即θw=θu+θiρi/ρw，其中θu为未冻水的体积含量；ρw为水的密度；kx，ky为水平向和垂直向的渗透系数；Ψ为土水势，参考文献[3]将土体冻结和融化过程中的土水势均定义为Ψ=Ψm+Y，其中Ψm为基质势，Y为重力势。

(3)

1.4 水热联系方程

上述公式(1)和公式(3)共有T，θu和θi三个未知量，只有两个方程，因此需要补充一个联系方程才能解方程组。联系方程为未冻水含量与温度变化的关系，如下所示：

θu=θ0|Ttrans|B|T|-B

(4)

式中：θ0为初始含水率Ttrans为冰水相变温度；B为与土质相关的常数。

给定式(1)、式(3)和式(4)边界条件和初始值就组成了路基工程水热耦合的数学模型。

将未冻水含量θu对时间和空间求偏导得：

(5)

将水分迁移方程代入热传导方程并联系上式可消去含冰率θi：

(6)

式(6)实现了方程组的解耦。将上式与水热联系方程(4)联立即可求得温度和水分的分布。

冻结区和未冻区之间有一个固态和液态分界界面，而且这个界面不固定。为了将该界面简化为固定状态，引入Heaviside阶梯函数H(T)[15]利用COMSOL Multiphysics软件内部自带的表达式表示如下：

H(T)=flc2hs(T-Ttrans,dT)

(7)

式中：Ttrans为相变点温度；dT为转变温度间隙取Ttrans=-0.3℃，dT=0.3℃，建立的函数变化曲线如图1所示。

图1 阶梯函数曲线

根据式(7)，在整个求解区域内表示未冻水含量为：

θu=θ0|Ttrans|B|T|-B[1-H(T)]+θ0H(T)

(8)

容积热容为：

C=Cf+(Cu-Cf)H(T)

(9)

式中：Cu为未冻土容积热容量；Cf为已冻土容积热容量。

渗透系数为：

k=kf+(ku-kf)H(T)

(10)

式中：ku为未冻土渗透系数；kf为已冻土渗透系数。

扩散系数为：

D=Df+(Du-Df)H(T)

(11)

式中：Du为未冻土水分扩散系数；Df为已冻土水分扩散系数。

参考以往研究成果[2,10]，定义阻抗系数I:

I=1010θi

(12)

那么渗透系数和水分扩散系数定义如下：

(13)

(14)

(15)

(16)

2 土体开放系统单向冻结试验

为验证土体开放系统中水热耦合理论的有效性，选取路基填土制作圆柱形试样，开展开放系统条件下单向冻结室内试验。试样直径为38 mm，高度为76 mm。试样初始含水率为20.2%，初始干密度为1.57 g/cm3。试样制作完成后，用防水薄膜将试样密封，仅保留底面以方便补水。试样放置在圆柱形筒中，并安装在如图2所示的冻胀装置中。该装置能够从上向下单向冻结试样，最低温度可达-40℃。此外，冻结过程中保持试样底部恒温在1℃左右，方便从外部给试样补水。沿试样高度方向每隔1.5 cm插入一个小型热电阻，实时监测试样内部温度变化。冻结试验共进行8 h，试样取出后，将土柱沿高度方向分层切片取样烘干测含水率。试样不同高度处的温度随冻结时间的变化及温度分布如图3所示。可以看到试样在恒温冻结的情况下，前2 h温度剧烈变化，之后逐渐稳定，冻结5 h后试样全部达到0℃以下。并且随着冻结时间的增加温度分布呈现由非线性分布向线性分布过渡的趋势，表明温度梯度逐渐稳定。

图2 土体冻融循环试验原理图

3 水热耦合理论的验证

3.1 计算模型及参数

本次研究借助于第二节中开放系统土柱单向冻结试验的温度和水分变化结果验证水热耦合理论在开放系统季节冻土区应用的合理性。采用有限元软件COMSOL Multiphysics根据圆柱土样尺寸建立有限元模型如图4所示，计算模型为二维轴对称模型，其中高度为7.6 cm，半径为1.9 cm。

图3 冻结过程中不同时刻试样的温度分布

图4 有限元计算模型

将试验开始时室内温度作为试样的初始温度，即T0=17.5℃。试样初始含水率为20.2%，因此初始体积含水率设定为θ0=0.32。

试样顶部即冷端(AD)温度边界条件按照试样实际的温度变化(冷却温度为-30℃时)进行设置。根据图5可以发现试样顶部(h=7.6 cm)温度变化趋势为先线性下降然后逐渐稳定，为简化计算过程，将试样顶部温度变化拟合为分段函数即式(17)，拟合结果如图5所示。

(17)

图5 试样顶部温度变化

AB为轴对称边界，CD为绝热边界，BC为暖端恒温边界，温度为试验过程中的实际温度，Tb=2℃。对于水分边界，AD和CD均为不透水边界，AB为轴对称边界，BC为补水边界，含水率恒定且为试样的饱和体积含水率，即θs=0.42。参考徐学祖等[16]总结的亚黏土热参数，综合取值如表1所示[17]。那么可以得到土体的未冻水含量随温度的变化曲线如图6所示。

表1 土体热参数取值

图6 未冻水含量随温度的变化曲线

COMSOL Multiphysics软件自带有系数型偏微分方程接口，可以将上述计算条件、参数及水热耦合方程导入软件中对方程重新编译并进行计算。

3.2 结果与分析

根据室内试验结果，试样冻结持续5 h后，内部温度全部达到零下并逐渐稳定。因此以试样冻结时间分别为1 h、2 h、4 h和5 h的温度变化、冻结锋面位置和水分分布对计算结果进行验证。

(1)试样温度变化对比。试样在冻结过程中的温度变化云图如图7所示。可以明显看到随着冻结时间的持续试样冻结区逐步向下发展。数值计算时试样底部保持恒温2℃，因此图7中冻结5 h后试样的底部仍高于0℃。而试验过程中虽然底板温度也保持在2℃左右，但是冻融循环试验系统中酒精循环液的冷却效率明显高于其保持恒温的效率，因此出现5 h后试样被完全冻结。

图7 试样冻结过程中的温度变化云图(单位：℃)

假定试样冻结锋面处土体的温度为0℃，那么试样在不同冻结时间计算冻结锋面高度和实际冻结锋面高度对比如图8所示，其中试样实际的冻结锋面用白线标记。可以看到冻结持续1 h后冻结锋面的计算高度和实际高度均为55 mm左右；冻结2 h后高度均为30 mm左右；冻结4 h后冻结锋面的高度出现一定差异，实际高度为20 mm左右，而计算高度为15 mm左右；冻结5 h后室内试验的试样已经被完全冻结，而模拟结果的冻结锋面的高度为5 mm左右。因此，在试样完全冻结以前冻结锋面的计算高度基本上和实际高度相差不大。

试样不同高度处的计算温度和实测温度随冻结时间的变化对比如图9所示。可以看到，计算温度和实测温度变化趋势一致，前2 h温度迅速降低，冻结5 h后逐渐稳定，且数值也比较接近。上述结果表明利用水热耦合理论得到的开放系统土体冻结过程中的温度计算结果是可信的。

图8 试样冻结过程中冻结锋面的变化

图9 试样不同高度计算温度和实测温度变化对比

(2)试样含水率变化对比。图10为不同冻结时间试样总体积含水率的计算值与实测值的对比。由于室内试验中5 h后试样即被完全冻结，而数值计算中温度边界条件的影响，土柱不会完全冻结。因此仅研究冻结时间分别为1 h、2 h和4 h时试样内部的总体积含水率分布状况。所述总体积含水率包含了未冻水含水率和冰晶融化为液态水后的体积含水率。从图10中可以看到计算含水率与实测含水率仅在试样顶部分布差异较大其余位置数值大小接近，整体的变化趋势基本一致，均呈现试样上部含水率高于下部并在中部发生突变，并且计算得到的含水率突变位置与实际的突变位置基本一致。将图10与图8对比，可以发现含水率突变位置即冻结锋面位置，冻结锋面以上为冻结区，冻结锋面以下为未冻区。由于土体中的水分是由未冻区向冻结区迁移的，因此试样上部含水率要高于下部。上述结果表明，利用水热耦合理论计算土体的含水率变化是符合实际的，可以应用于季节冻土区受地下水迁移影响的岩土体的温度和含水率计算中。