非饱和正冻土一维水热耦合模型与试验

明锋，李东庆

(中国科学院 寒区旱区环境与工程研究所 冻土工程国家重点实验室，甘肃 兰州，730000)

摘要: 利用Clapeyron方程来描述冻土内温度、冰压力和水压力的关系，根据质量守恒原理和能量守恒原理，建立一维非饱和冻土水热耦合模型。采用数值模拟和室内试验相结合的方法，对建立的水热模型进行检验。研究结果表明：随着冻结时间的增加，冻土的温度变化由剧烈转成平稳，最后沿高度线性分布；冻土段含水量增大，未冻土段含水量减小；在冰透镜体的存在的地方，孔隙比明显增大，含水量明显增大；模拟结果与实测结果较吻合，验证了该模型的正确性。

关键词：非饱和正冻土；水热耦合；孔隙比；数学模型

中图分类号：TU445            文献标志码：A         文章编号：1672-7207(2014)03-0889-06

Modeling and experimental investigation of one dimension coupled moisture and heat in unsaturated freezing soil

MING Feng, LI Dongqing

(State Key Laboratory of Frozen Soil Engineering, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, China)

Abstract: Clapeyron equation was applied in freezing soil to describe the relationship among temperature, water pressure and ice pressure when ice and water coexist in phase equilibrium. According to the mass and energy conservation principle, a mathematical model of coupled moisture and heat in unsaturated freezing soil was established. The methods of numerical simulation and indoor experiments were applied to verify the correctness of this model. Results show that as freezing time increases, temperature fluctuation turns from being severe into smooth and the final temperature shows obvious linear distribution by height. Water content increases in frozen soil area, and decreases in unfrozen area. Wherever the ice lens body exists, void ratio tends to increase and water content increases substantially. The simulation results and measured results can match each other well and the effectiveness of this model is fully verified.

Key words: unsaturated freezing soil; coupling of moisture and heat field; void ratio; mathematical model

中国的多年冻土面积占国土面积的22.3%，分布较广。冻土区的冻融作用使得构筑物出现裂缝、倾斜、沉陷、结构断裂等现象，给道路和构筑物造成很大的危害。冻土冻胀的机理研究一直是冻土界注重的问题之一^[1-3]。Harlan^[4]根据未冻水动力学原理，提出考虑水分迁移和冰水相变的一维非线性水热耦合方程。该方程可用来求解温度和水分的变化。Konrad等^[5]通过大量实验，提出分凝势的概念，定义为水分迁移速率与温度梯度的比值。Shen等^[6]提出简化的水热力三场耦合模型，并进行数值模拟。Fremond等^[7]提出热力学水、热、力三场耦合模型，该模型考虑由冻胀引起的孔隙变化，但是不能很好地解释冻土的成冰机制。李洪升等^[8]将冻土视为空间弹性体，建立冻土水热力耦合作用的数学模型，并给出离散方程与解法。周家作等^[9]通过引入Heaviside 阶梯函数和Z函数, 提出开放系统下饱和土在冻结过程中热质迁移的数学模型，并对一个侧面绝热的圆柱土样的热质迁移过程进行数值模拟。李东庆等^[10]考虑基质势和压力势与温度的关系，给出季节性冻土水-热-三场耦合的数学模型，并将耦合问题归结为求解一个非线性、非稳态温度场微分方程。赵刚等^[11]利用室内水分迁移试验装置，研究初始含水量及温度模式对水分迁移的影响。国内外学者已进行大量饱和土的温度场、水分场、应力场及其耦合研究，但在实际中，土体多处于非饱和状态，在浅层地表更是如此。本文作者综合考虑水、热、力的相互关系，提出非饱和土的水热耦合数学模型，并通过数值模拟和室内试验对模型进行验证。

1  数学模型

1.1  基本假定

为充分考虑冻土中水、热耦合过程的主要因素，合理简化模型，作如下基本假定：

(1) 土体均匀连续，且为各向同性体。

(2) 水分迁移均以液态形式进行，不考虑冰晶的迁移。

(3) 土颗粒和冰晶不可压缩，忽略冰颗粒的压融。

(4) 水分迁移、空气流动均符合达西定律。

(5) 冻土体内的空气与外界相通。

(6) 计算过程中温度无损失。

图1所示为非饱和冻土组相图，其中，e为孔隙比，S_sat为饱和度，S_i为含冰率。

图1  非饱和冻土组相图

Fig. 1  Phase diagram of unsaturated freezing soil

冻土中未冻水含量与温度保持动态平衡关系，含冰率是温度的函数^[12]：

        (1)

式中：为温度；为系数；为冻结温度。

1.2  水分迁移方程

1.2.1  冻土渗透系数

根据文献[8]，非饱和土体的渗透系数k₀与干容重和饱和度之间有如下关系：

       (2)

式中：k₀为土体的渗透系数，m/s；土体的干容重，g/cm³。

冻土渗透系数k与温度θ有关^[13]。当冰透镜体形成之后，在土中就形成一个隔水层，水分无法穿过冰透镜体进入到冻结区。因此，冻土渗透系数k可改写成θ的函数：

      (3)

x轴方向从未冻区指向冻结区，x_s是冰透镜体形成的位置。采用有限单元法计算时，先依次判别求解域内每个节点上的孔隙比e是否达到临界孔隙比e_s。如果达到e_s，则标记该点为x_s，然后，在x≥x_s范围内，将渗透系数强制改变为0。

1.2.2  水分迁移驱动力

在考虑外部荷载的未饱水开放系统中，土水势可用下式表示^[10]：

                (4)

式中：为转化系数，P_w为水压力，G为重力势。

Clapeyron方程可来描述相平衡状态下，压力与温度的关系：

            (5)

式中：L为单位质量的水(或冰)冻结(或融化)所释放(或吸收)的热量；和分别为孔隙水和孔隙冰的绝对压力；为绝对温度；为水冻结的绝对温度。

一点处的总应力P由孔隙冰压力P_i和孔隙水压力P_w组成，它们之间的关系可用下式表示：

                      (6)

式中：为孔隙水压力系数；P为总应力。

联立式(5)和(6)，并用相对压力和摄氏温度表示P_w，则有

          (7)

对于任意截面，其总应力等于外荷载。

1.2.3  水分迁移方程

任意时刻，微小单元dx内冰和水的质量分别为m_i和m_w，由图1可得

，    (8)

根据质量守恒有

          (9)

将式(8)带入式(9)，可得

                (10)

考虑土颗粒不可压缩性，所以在上式两边同时除以dx_s。通过整理可得水分迁移方程：

               (11)

1.3  热扩散方程

微元dx内冰体积含量为，单位时间微元体内由于冰体积含量增加所释放的潜热为。根据能量守恒，建立考虑相变的热扩散方程：

    (12)

式中：C为土体的体积热容量；为土体的导热系数。

土体的体积热容量和导热系数分别为

   (13)

          (14)

式中：C_s，C _i，C_w和C_g分别为土颗粒、冰、水和气体的体积热容量。λ_s，λ_i，λ_w和λ_g分别为土颗粒、冰、水和气体的导热系数。

将各项参数带入式(13)中，整理可得：

    (15)

由水分迁移方程(式(11))及热流方程(式(15))构成本问题完整的数学描述，由此建立考虑荷载作用的非饱和冻土的水热耦合数学模型。

1.4  冰分凝条件

对于与外界连通的非饱和土，当空气体积减小到0以后，冻土内才会有冻胀的可能。O’Neill^[14]认为当孔隙压力等于总应力时，冰透镜体产生。这一方法得到Thomas等^[15]的进一步发展。当冻胀力大于上覆压力时，土颗粒将被推开，留出的孔隙将由分凝冰填充。本文采用O’Neill的计算理论，来判定冰透镜体的产生。

由水热耦合方程和冰分凝条件，通过初始条件、边界条件及其他所需参数，就可以算出温度、含水量、冻胀变形和冰透镜体的分布及其随时间的变化规律。

2  结果分析

将一个直径10 cm、高度10 cm的圆柱形土样，放入1 ℃环境中的恒温箱中恒温24 h，让土样温度均匀分布。然后，将土样顶端温度降为-3 ℃，恒温；底端温度固定为1 ℃，侧面绝热，冻结100 h。初始空隙比为0.6，初始饱和度为0.72。土样顶端瞬时冻结，孔隙水只发生原位冻结。土样底端采用无压补水。对土样顶部施加轴向应力100 kPa，土样完全侧限。由于水分和热量只在试件轴向传递，变形只在轴向发生，本问题可简化成一维水热耦合问题。

为能使数值模拟与室内试验具有良好的可比性，模拟条件与试验条件一致。模拟中，初始条件和边界条件如下：，，= 0.6， ℃， ℃，，℃。

数值模拟采用有限元分析软件进行求解，沿土样高度划分120个单元，时间步长为3 min，总计算时间为100 h。计算参数及取值如表1所示。

为直观地表达冰透镜体的生长过程及分布情况，图2(a)所示为冻结35，60和100 h时的冰透镜体柱状分布图。从图2可见：冰透镜体总是呈不连续层状分布，这是符合实际情况的。在35 h冰透镜体规模较小，随着时间延长，冻结锋面下移，下层冰透镜体生长。由于底部有水分补给，在孔隙压力和温度梯度的作用下，水分不断被吸收上来，使土样产生冻胀，土样高度逐渐增大。从图2(b)可见：土样顶端冻结后产生网状裂隙，裂隙中充填冰，裂隙越密集表明含水量越大，在2.5~4.5 cm高度位置范围内，裂隙最密集；在4 cm处有冰透镜体产生，可知该位置处含水量最大。冰透镜体位置往上，主要是原位冻结，因此，裂隙密度逐渐减小，而在融土内没有裂隙产生。

表1  参数表

Table 1  Parameters table

图2  不同时刻冰透镜体分布图(白色为冰透镜体)

Fig. 2  Histogram of ice lenses distribution at different times (white parts are ice lenses)

图3  各点温度随时间变化曲线

Fig. 3  Temperature curves with time at different positions

为分析试样各点温度的变化规律，沿土样高度方向均匀布置10个温度探头。图3所示为2个不同测点的温度随着时间的变化曲线。从图3可见：各点温度在冻结初始阶段，温度迅速降低；随着冻结时间增长，降温速率有所减缓，最后达到稳定的变化速率；靠近冷端的点温度下降迅速，远离冷端的点温度下降较为缓慢。因此，各点温度变化达到稳定所需要的时间不同，表现为随着该点与冷端距离增大而耗时更长。温度达到稳定状态，也只是相对稳定，存在一定的波动。

图4  不同时刻土中的温度剖面

Fig. 4  Temperature profiles at different times

从图4可见：在开始的1~10 h内温度变化幅度较大；随着冻结时间的增加，土样的温度分布逐渐趋于稳定，大约在25 h后温度分布进入似稳定状态，最后温度随着土样高度近似呈线性分布。由此可见：随着冻结时间的增加，冻结深度不会一直增加，而是将达到一个稳定值。随着冻结时间的增长，冻结锋面下移，冻土段长度逐渐增加而未冻土段长度缩短。由于土样制作的不均匀性以及温度探头存在标定误差，因此，实测温度并没有严格按照高度呈线性分布。数值模拟是在理想情况下计算的，因此，温度是随高度线性分布的。

图5  冻胀量随时间变化

Fig. 5  Frost heave with time

图5所示为土样冻胀量随时间的变化。在冻结初期，由于是非饱和土体，在较长一段时间里，气体被慢慢排出，因此，不存在冻胀现象。在降温的同时进行加压，土样有一个被压缩的过程，因此，在冻结前期土体不但没有冻胀，反而处于被压缩状态。此外，由于土样具有热胀冷缩的性质，因此，在初始一段时间内，冻胀量为负值。在冻结13 h后，土样开始产生冻胀；在冻结后期，冻胀量随时间几乎线性增长。从图5可见：2条冻胀量曲线变化趋势基本一致，数值上也差别较小，验证数学模型的正确性。

图6  冻结100 h时的含水质量分数剖面

Fig. 6  Water content profile after freezing 100 h

图6所示为冻结100 h后，含水量沿土样高度的分布情况。水分迁移是引起冻土体冻胀、融沉的主要影响因素。冻土中水分场的分布是研究冻土稳定性和预防冻胀、融沉的主要参考因素。试验结束后，将土样沿高度方向每1 cm划分1份(最后不足1 cm按1份考虑)，测定其含水量。从图6可看到：在冻结锋面附近，含水量较其他地方大。这与图2所示的冰透镜体的分布规律是一致的。实测含水量在底部较模拟值增大，是因为该位置处于补水端，浸泡在水中，导致实测含水量偏大。

不同时刻孔隙比沿高度变化见图7。由图7可见：在未冻土段的孔隙比较初值有所减小，冻土段孔隙比较初始孔隙比明显增大，尤其是在2.5~4.5 cm处。由于未冻土段的孔隙压力为负，有效应力增大，受到冻胀力挤压而产生了压密过程，导致孔隙比减小。而冻土段，由于孔隙水原位冻结和迁移水形成了冰透镜体，冰透镜体要生长，必然要推开土颗粒，从而增大土颗粒的间隙，导致孔隙比增大。从图7还可见：随着冻结时间的增长，孔隙比逐渐增大，最大孔隙比的位置也逐渐往暖端移动，最后在冻结锋面处达到最大。在靠近暖端位置，土样被压密，孔隙减小，这与图6中该位置的含水量较低是吻合的。

图7  不同时刻孔隙比沿高度变化

Fig. 7  Porosity ratio profiles at different time

3  结论

(1) 根据质量守恒原理和能量守恒原理，利用Claypeyron方程，建立以孔隙比为变量的非饱和冻土水热耦合模型。

(2) 冻土段含水量增大，未冻土段含水量减小。约在35 h形成冰透镜体。在冰透镜体的存在的地方，孔隙比明显增大。

(3) 随着冻结时间的增加，冻结锋面逐渐下移，最后温度逐渐趋于稳定，温度沿土样高度近似呈线性分布。

(4) 从温度分布、含水量分布以及冻胀方面看，数值模拟结果与试验结果在趋势上具有良好的一致性，验证了模型的正确性。

参考文献：

[1] Grenier S, Konrad J M, LeBouf D. Dynamic simulation of falling weight deflectometer tests on flexible pavements using the spectral element method: Forward calculations[J]. Canadian Journal of Civil Engineering, 2009, 36(6): 944-956.

[2] 武建军, 韩天一. 饱和正冻土水-热-力耦合作用的数值研究[J]. 工程力学, 2009, 26(4): 246-251.

WU Jianjun, HAN Tianyi. Numerical research on the coupled process of moisture-heat-stress fields in saturated soil during freezing[J]. Engineering Mechanics, 2009, 26(4): 246-251.

[3] 马巍, 王大雁. 中国冻土力学研究50年回顾与展望[J]. 岩土工程学报, 2012, 34(4): 625-640.

MA Wei, WANG Dayan. Studies on frozen soil mechanics in China in past 50 years and their prospect[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(4): 625-640.

[4] Harlan R L. Analysis of coupled heat-fluid transport in partially frozen soil[J]. Water Resources Research, 1973, 9(5): 1314-1323.

[5] Konrad J M, Morgenstern N R. Prediction of frost heave in the laboratory during transient freezing[J]. Canadian Geotechnical Journal, 1982, 19(3): 250-259.

[6] SHEN Mu, Ladanyi B. Modeling of coupled heat, moisture and stress filed in freezing soil[J]. Cold Regions Science and Technology, 1987, 14(3): 237-246.

[7] Fremond M, Mikkola M. Thermo mechanical modeling of freezing soil[C]// Proceedings of the Sixth International Symposium on Ground Freezing. Yu X, Wang C, ed. Rotterdam: Kalke Ma, 1991: 17-24.

[8] 李洪升, 刘增利, 梁承姬. 冻土水热力耦合作用的数学模型及数值模拟[J]. 力学学报, 2001, 33(5): 621-629.

LI Hongsheng, LIU Zengli, LIANG Chengji. Mathematical model for coupled moisture、heat and stress field and numerical simulation of frozen soil[J]. Acta Mechanica Sinica, 2001, 33(5): 621-629.

[9] 周家作, 李东庆, 房建宏, 等. 开放系统下饱和正冻土热质迁移的数值分析[J]. 冰川冻土, 2011, 33(4): 791-795.

ZHOU Jiazuo, LI Dongqing, FANG Jianhong, el at. Numerical analysis of heat and mass transfers in saturated freezing soil in an open system[J]. Journal of Glaciology and Geocryology, 2011, 33(4): 791-795.

[10] 李东庆, 周家作, 张坤. 等. 季节性冻土的水-热-力建模与数值分析[J]. 中国公路学报, 2012, 25(1): 1-7.

LI Dongqing, ZHOU Jiazuo, ZHANG Kun, et al. Modelling and numerical analysis of moisture, heat and stress in seasonal frozen soil[J]. China Journal of Highway and Transport, 2012, 25(1): 1-7.

[11] 赵刚, 陶夏新, 刘兵. 重塑土冻融过程中水分迁移试验研究[J]. 中南大学学报(自然科学版), 2009, 40(2): 519-525.

ZHAO Gang, TAO Xiaxin, LIU Bing. Experimental research on water migration in remoulded soil during freezing and thawing process[J]. Journal of Central South University (Science and Technology), 2009, 40(2): 519-525

[12] Tice A R, Anderson D M, Banin A. The prediction of unfrozen water contents in frozen soils from liquid limit determinations[R]. U.S. Army Corps of Engineers, Cold Regions Research & Engineering Laboratory, 1976: 76-83.

[13] Gilpin R R. A model for the prediction of ice lensing and frost heave in soils[J]. Water Resources Research, 1980, 16(5): 918-930.

[14] O’Neill K, Miller R D. Exploration of a rigid ice model of frost heave[J]. Water Resources Research, 1985, 21(3): 281-296.

[15] Thomas H R, Cleall P, Li Y C, et al. Modelling of cryogenic processes in permafrost and seasonally frozen soils[J]. Geotechnique, 2009, 59(3): 173-184.

(编辑  邓履翔)

收稿日期：2013-02-20；修回日期：2013-04-15

基金项目：国家重点基础研究发展计划(“973”计划)项目(2012CB026102)；国家自然科学基金资助项目(41271080)；中国科学院西部计划项目(KZCX2-KB3-19)

通信作者：明锋(1986-)，男，重庆人，博士研究生，从事寒区岩土工程研究；电话：0931-4967118；E-mail: mf0329@163.com