非饱和土体水平冻结过程中的水热耦合作用机制研究

申林方，王志良，李腾风，孟鸿斌

(昆明理工大学 建筑工程学院，云南 昆明 650500)

季节性冻土地区土体的冻胀融沉病害严重影响铁路路基结构的稳定性以及列车行驶的安全性[1-3],是目前工程界亟须解决的问题。这些病害与土体中的水分迁移、温度变化及冻结相变密切相关,是典型的水热耦合问题。一方面土体冻结时未冻区的水分向冻结锋面迁移,并伴随着相变潜热的释放,从而改变土体中的温度分布及冻结锋面的移动;另一方面土体温度场的演化又会引起水分冻结相变、特征参数的变化及迁移运动[4]。因此,有必要考虑冻结相变作用研究非饱和土体的水热耦合作用机制。

针对冻结相变作用下的土体水热耦合问题,相关学者展开了大量的研究工作。Harlan[5]根据土体冻结相变过程中的能量变化以及水分运移规律,最早提出了描述土体冻结过程的水热耦合模型。在此基础上,Taylor等[6]通过引入阻抗系数来表征冰对土体导水能力的阻碍作用,研究了水热耦合作用下冰阻抗对水分迁移规律的影响。由于土体含相变水热耦合模型的非线性、初始条件及边界条件的复杂性,采用解析法或半解析法求解非常困难,需借助数值计算方法进行计算[7]。为此,文献[8-9]对求解区间进行网格划分和时间离散,采用有限差分法求解了Harlan模型;文献[10-11]采用消除源项的方法降低温度场与水分场演化方程间的耦合性,基于全隐式有限差分计算格式求解了冻土的水热耦合方程,并分析了阻抗系数对土体水热耦合演化规律的影响。然而,上述研究未考虑冻结锋面位置附近,未冻区和已冻区未冻含水率与温度间导数不连续的问题。为此,文献[12-13]对冻结锋面位置处的特殊节点进行优化处理,基于有限容积法建立一种求解水热耦合模型的新数值计算方式。为简化计算,该方法将非饱和土体在负温条件下的最大未冻含水率作为冻结区土体的未冻含水率,这显然与实际情况存在一定的误差,且会影响冰组分的计算。

为此,在进行演化方程组解耦、冻结锋面位置精细化处理的同时,采用土体冻结特征曲线改进冻结区未冻含水量计算,通过引入总含水量优化冰组分生成,提出一种新的水热耦合作用模型。然后,采用有限差分法对该模型进行数值求解,并结合试验数据验证本模型的准确性和有效性。最后,以张掖壤土为研究对象,分析讨论冻结时间、初始含水量、温度梯度以及冰阻抗等因素,对水平冻结条件下非饱和土体水热耦合作用机制的影响。

1 水热耦合数学模型

1.1 基本方程

忽略水气运移、热量对流、冰透镜体形成以及外荷载的作用,得到在水平冻结作用下(不考虑重力方向的渗流)非饱和土体的水热耦合控制方程[6]。

温度场:

(1)

水分场:

(2)

式中:T为土体的温度,℃;θw、θI分别为非饱和土体的体积未冻含水量和含冰量,%;Dw为水分扩散系数,m2/s;Cv为体积比热容,kJ/(kg·℃);λ为体积热传导系数,W/(m·℃);ρI、ρw分别为冰、水密度,kg/m3;L为相变潜热,kJ/kg;t为时间,h。

土体冻结后冰的生成阻断了部分水流通道,从而影响水分迁移。为了表征冰组分对冻土水分扩散系数的影响,Taylor等[6]基于试验和数值模拟结果提出了阻抗因子的概念。考虑冰阻抗的影响,水分扩散系数Dw可表示为

(3)

式中:D为未冻区的水分扩散系数,m2/s;I0为阻抗因子,I0=1010θI。

考虑土体各组分影响,其体积比热容Cv及体积热传导系数λ可表示为[14-15]

Cv=ρsCs(1-φ)+ρwCwθw+ρICIθI+ρaCa·

(φ-θw-θI)

(4)

(5)

式中:φ为土体的孔隙率,下标s、w、I、a分别表示土颗粒、未冻水、冰、空气。

上述方程有温度、未冻含水量及含冰量三个待求变量,而仅有两个方程,故需补充一个联系方程——土体冻结特征曲线,并由此可建立非饱和土体冻结区的未冻含水量与温度的曲线方程为

θw≤θmax(T)

(6)

式中:θmax(T)为非饱和土体在负温条件下对应的最大未冻含水量。

1.2 方程的改进

由温度场演化方程式(1)和水分场演化方程式(2)可知,相关参数及源项间均存在耦合作用,而联系方程式(6)又取决于试验结果,在求解过程中需考虑数值计算和试验数据之间的整合问题。为此,本文在对演化方程组进行解耦的同时,对含冰量计算进行了优化处理。

1.2.1 演化方程组解耦

为降低温度场和水分场间的耦合作用,将式(1)与式(2)联立,消除温度场演化方程中的源项,即:含冰量变化率∂θI/∂t。同时,在求解过程中,将未冻含水量θw作如下变换

(7)

则温度场的演化方程可变换为[13]

(8)

1.2.2 含冰量优化计算

联系方程式(6)所对应的未冻含水量主要为土颗粒周围未冻结的水分,其数值只能用来校核冻结区的未冻含水量,而不能直接求解含冰量的变化。为较好地处理冰组分的生成过程,本文引入总含水量θaw,其计算过程为

(9)

含冰量的计算过程为:根据式(2)计算土体的未冻含水量,然后结合式(6)对其进行校核,并根据式(9)进行总含水量计算,最后求得的含冰量为总含水量与未冻含水量之差,即

(10)

因此,将式(2)、式(6)、式(8)～式(10)联合,建立新的水平冻结作用下非饱和土体水热耦合模型。

2 计算模型的离散化处理

2.1 控制方程的离散过程

为了求解所建立的非饱和土体水热耦合模型,采用隐式中心差分格式,对其进行数学离散。以水分场演化方程式(2)为例,在节点i位置的离散形式可表示为

(11)

式中:Δt为时间步长;Δx为网格间距;上标n为离散的时步;下标p、s分别表示向前、向后半个节点位置处对应参数的平均值。

将式(11)进行整理,可得

(12)

式中:αi、βi、γi、hi分别为节点i位置处的参数,满足

(13)

同理,可得到温度场演化方程式(8)的离散形式为

(14)

(15)

总含水量θaw求解的离散形式可表示为

(16)

2.2 冻结锋面的精细化处理

求解离散方程时,在冻结锋面处非饱和土体的温度、总含水量、未冻含水量等参数需借助于联系方程式(6)求得。同时,未冻水转换成冰使得总含水量和未冻含水量开始发生分化,而水冻结成冰的瞬间释放相变潜热,因此在冻结锋面位置处其水热演化机制极其复杂,相应的离散过程需进行单独讨论。

为此,假设If位置处于未冻区,而If+1位置则为冻结区。在冻结锋面处,由于If节点的未冻含水量与温度无函数关系,而If+1位置处则需引入特征曲线,因此对于dθw/dT在(If)p位置处或者(If+1)s处没有明确定义。故dθw/dT在(If)p位置处或者(If+1)s处,进行离散展开得

(17)

将式(17)代入式(8)可得

(18)

对式(18)进行简化,可得

(19)

在冻结锋面位置,无法通过式(19)直接求得未冻含水量和温度两个变量,故需要补充相关方程。为此,在迭代过程中,对dθw/dT进行展开,并根据牛顿切线法得[16]

(21)

式中:*表示本时步上次迭代的计算结果。

联立式(19)和式(21),可以求得n+1时刻冻结锋面位置处的未冻含水量及温度。

2.3 求解流程

本文冻结相变作用下非饱和土体水热耦合模型的数值计算步骤如下:

Step1根据初始条件及边界条件,结合式(14)、式(19)及式(21)求解非饱和土体本时刻的温度及冻结锋面处未冻含水量。

Step2根据Step1求得的温度,结合联系方程式(6)的特征曲线,对非饱和土体冻结锋面位置处的未冻含水量进行校核:当其数值不大于本时刻温度所对应的θmax时,未冻含水量的数值不变;而当超过θmax时,未冻含水量则由θmax代替。然后将校核后的未冻含水量代入式(12),求解整个计算域内的未冻含水量。

Step3根据本时刻的温度,结合联系方程式(6),对Step2求得的未冻含水量进行再次校核,确保数值的准确性。然后将校核后的未冻含水量代入式(16),求得本时刻总的含水量。同时根据式(10),求得本时刻的含冰量。

Step4由于上述求解过程中的参数信息并非由本时刻的温度、未冻含水量、含冰量求得,因此,需要重复步骤Step1～Step3更新本时刻的计算结果,并得到本时刻第一次迭代的温度、未冻含水量、总含水量、含冰量。然后再次重复步骤Step1～Step3,并根据相邻两次迭代计算结果的相对误差评判其收敛性,如相对误差小于控制误差(本文取为10-4)时,本时刻迭代计算结束;否则,继续迭代计算。

Step5重复步骤Step1～Step4,直至达到预定冻结时间。

3 算法验证

结合文献[8]中的试验算例,验证本文计算模型的有效性。在封闭条件下对石英粉进行水平冻结试验,试样直径为0.1 m,长为0.3 m。在计算过程中,将水平方向0.3 m均匀划分为300个网格,计算时步设为2 s。初始时刻石英粉温度为20 ℃,含水率为15.59%。试验过程中,模型左端为恒温不透水边界,右端为冷源边界,其温度为

(22)

冻结区未冻含水量与温度间的关系为

(23)

水分扩散系数为

(24)

体积比热容为[8]

Cv=ρd(Cs+4.184mw+2.10mI)

(25)

式中:ρd为石英粉的干密度,ρd=1 330 kg/m3;Cs为石英粉的骨架质量比热容,Cs=837 J/(kg·℃);mw、mI分别为质量含水量及含冰量。

石英粉的导热系数为[17]:冻结区λf=6.875 kJ/(m·h·℃),未冻结区λu=12.396 kJ/(m·h·℃),冻结锋面位置处取为二者的平均值。

图1及图2分别为石英粉冻结过程中温度场和水分场的本文数值解与试验结果对比。从图1和图2可以看出,在不同时刻土体温度场及水分场的分布规律均与石英粉的试验结果吻合较好,这充分证明了本数值方法的有效性,而两者间的误差主要是由于数值计算的参数取值所致。

图1 温度场的本文数值解与试验结果对比

图2 水分场的本文数值解与试验结果对比

4 分析讨论

为研究非饱和土体在水平冻结过程中,冻结时间、初始体积含水量、温度梯度及冰阻抗等因素对温度场及水分场演化规律的影响,基于本文计算模型,针对张掖壤土的冻结演化过程进行了数值计算。计算模型长为0.15 m,相应的划分成150个网格。土体的初始温度为5 ℃,初始体积含水量为25%。右端为恒温不透水边界,左端为冷源边界,冷源温度从初始温度开始每小时降温1 ℃,直至达到设定T0,然后保持恒温T0不变。根据文献[10],土体的孔隙率φ为0.444,张掖壤土各组分的热物理参数见表1。

表1 非饱和土体各组分的热物理参数

根据文献[10]的试验结果,最大未冻含水量和温度间的关系为

(26)

水分的扩散系数取为

Dw=2.03×G7.35/I0

(27)

式中:G为土体饱和度。

4.1 不同冻结时刻的水热耦合演化规律研究

为研究温度场及水分场随时间的发展演化规律,设土体初始温度为5 ℃,冷端达到恒温为T0=-5 ℃时,计算得到了冻结12、24、72、144 h后的温度场及水分场分布见图3,冻结锋面随时间的演化规律见图4。从图3可以看出,土体温度随时间的推移逐渐降低,且在冻结前期发展较快,随后逐渐趋于平缓。由于冻结相变作用,未冻区的水分会向冻结锋面迁移,并在锋面位置处明显发生水分积聚的现象,从而导致冻结区的总含水量明显大于未冻区。同时,由图4可知,冻结锋面在初始阶段发展较快,随着时间的推移,其推进速度逐渐减缓,从而导致锋面在相应位置的持续时间较长,使得未冻区水分有充足的时间向锋面位置迁移,水分积聚的现象愈加明显,这也是大多数冻土工程产生破坏的主要原因。

图3 不同时刻的温度场及水分场分布

图4 冻结锋面位置随时间的演化规律

4.2 初始含水量对水热耦合演化规律影响

为研究土体初始含水量对水热耦合作用机制的影响,在保持土体孔隙率不变的情况下,研究了冷端达到恒温为T0=-5 ℃、初始体积含水量θw0分别为22%、25%、28%时,不同时刻土体温度及体积含水量的变化规律。

图5为不同初始含水量情况下不同时刻土体的温度及体积含水量分布趋势。由图5可知,在相同时间内初始体积含水量越高,其温度下降速度越快,且随着时间的推移,该差别愈发明显。这主要是由于土体冻结过程中热扩散系数发生变化,在初始阶段体积含水量为22%、25%、28%的土体的热扩散系数分别为1.92×10-7、2.03×10-7、2.15×10-7m2/s,由此可见,虽然热扩散系数随初始含水量的增加而增大,但三者间的差距较小。随着时间的推移,当t=144 h时,三者在冻结区的热扩散系数分别为3.21×10-7、3.84×10-7、4.51×10-7m2/s,相比于初始阶段其差异显著增大,从而使得温度场的发展趋势也产生明显不同。同时,在三种不同初始含水量的情况下,冻结区的体积含水量均明显大于未冻区,且随着土体初始含水量的增加,冻结区的体积含水量增加更为显著。当t=144 h时,冻结区内的平均体积含水量比初始阶段分别增加了4.65%、6.20%、7.38%;而在未冻区则下降了3.68%、6.40%、9.47%。这主要是由于在外界条件以及土体孔隙率不变的情况下,冻结区未冻水和温度间的关系特征曲线是唯一的,故土体初始体积含水量越大,冻结区内其体积含水量偏离特征曲线的程度也越剧烈,水分场受到的影响也越显著。

图5 不同初始含水量情况下土体温度场及水分场的分布

4.3 温度梯度对水热耦合演化规律影响

在保持土体孔隙率不变的情况下,初始体积含水量取为θw0=25%,冷源温度T0分别取为-3、-5、-7 ℃时,计算土体在冻结过程中温度场和水分场的演化情况,见图6。由图6可知,不同冷源温度作用下,土体中温度场的分布规律较为相似,且冷源温度越低,土体温度下降速度越快,冻结锋面的发展速度也越快。同时,冷源温度越高,土体冻结速度越慢,水分迁移过程也就更为充分,导致相同时间内冻结区的平均体积含水量越高,其水分积聚现象越显著,这与文献[18]的试验结果相一致。当冷源温度T0分别为-3、-5、-7 ℃时,冻结144 h后冻结区的平均体积含水量比初始时刻分别增加了8.12%、6.22%及5.27%,同时,冻结锋面位置距离冷源分别为0.075、0.095、0.106 m。

4.4 冰阻抗对水热耦合演化规律影响

在冷源作用下,土体内部水分逐渐冻结成冰,进而阻碍该区域的水分迁移。为此,针对考虑冰阻抗[I0=10(10θI)]和不考虑冰阻抗(I0=1)两种情况,研究冷端达到恒温T0=-5 ℃时、初始体积含水量分别为22%、25%、28%时,土体中温度场及水分场的演化规律。

图7为考虑冰阻抗与不考虑冰阻抗两种情况下冻结144 h后土体温度场及水分场的分布规律。由图7可知,冰阻抗通过影响冻结锋面处的温度,进而改变土体中的整体温度分布,但温度总体变化数值不大。这主要是由于冰阻抗的作用影响冻结锋面附近的水分迁移,从而改变了土体的热扩散系数,影响其温度场的演化。在冻结144 h后,考虑阻抗作用和不考虑阻抗作用两种情况下土体的热扩散系数分别为3.84×10-7、4.02×10-7m2/s,两者的差别并不显著,因此冰阻抗作用对温度场的影响相对较小。对于水分场而言,无论是否考虑冰阻抗的影响,由于冻结区的未冻含水量较小,水分均由未冻区向冻结锋面处迁移,从而导致该位置处的含水量显著增加。如果不考虑冰阻抗影响,其冻结锋面附近的水分迁移能力会增强,从而导致体积含水量有所增大,并对冻结锋面位置产生一定影响。这主要是由于土体在冻结过程中其含冰量的积累,会使冰阻抗因子呈指数形式的增大,而水分扩散系数则以相应的形式减小,导致锋面位置处的水分迁移严重受阻。

图7 冰阻抗对土体温度场及水分场的影响

5 结论

(1)针对计算模型进行演化方程组解耦、冻结锋面位置精细化处理的基础上,采用土体冻结特征曲线改进冻结区的未冻含水量,并引入总含水量优化含冰量的计算,提出了一种新的水平冻结水热耦合计算模型。由验证实例可以看出,本文计算结果与试验结果吻合较好,证明了该计算方法的准确性和有效性。

(2)在土体孔隙率不变的情况下,土体初始含水量越大,其冻结区的体积含水量增加越显著,同样未冻区的含水量减少量也越大。当初始体积含水量分别为22%、25%、28%时,土体冻结144 h后,冻结区内的平均体积含水量比初始阶段分别增加了4.65%、6.20%、7.38%,而在未冻区则下降了3.68%、6.40%、9.47%。

(3)在不同冷源温度作用下,土体中温度场的变化规律较为一致,但冷源温度越低,土体温度下降速度越快,冻结锋面的发展速度也越快。同时,冷源温度越高,土体冻结速度越慢,水分迁移过程也就越充分,相同时间内冻结区的平均体积含水量越高,其水分积聚现象越显著。当冷源温度分别为-3、-5、-7 ℃时,冻结144 h后冻结区的平均体积含水量比初始时刻分别增加了8.12%、6.22%、5.27%。

(4)冰阻抗通过影响冻结锋面附近的水分迁移作用,而导致其体积含水量及热扩散系数发生变化,并对冻结锋面位置产生一定影响。同时,冰阻抗对冻结锋面附近的温度影响最大,而对冻结区及未冻区的影响较小,但相比于含水量,冰阻抗对温度场的影响并不显著。