基于COMSOL Multiphysics的不同水头压力下水热耦合分析

Hydrothermal Coupling Analysis under Different Head Pressures Based on COMSOL Multiphysics

摘要: 冻土的水热耦合问题一直是冻胀融沉的主要原因之一，涉及到土壤中水分的迁移、热量的传递以及相变过程。在多年冻土和季节冻土区，由于环境因素的影响，冻土路基容易发生不均匀沉降、冻胀融沉等灾害，这些灾害都与水分迁移、相变以及温度变化息息相关。为了解决冻土水热耦合问题，本文通过COMSOL Multiphysic软件建立水热耦合的二维模型，研究不同水头压力作用下同一时间段温度和压力的影响、对流出边界总热通量的影响、对总液态水体积的影响，以及对最低温度的影响进行分析，模拟了3种不同水头压力的工况。结果表明：同一时间段的温度和压力随着水头压力的增加而变化；流出边界离开系统的总热通量也随着水头压力的增大而增大；不同工况下水头梯度越大，冰块融化速度越快，达到总液态水体积最大值的时间越短。

Abstract: The hydrothermal coupling problem of frozen soil has always been one of the main reasons for frost heave and thaw settlement, involving the migration of water in the soil, heat transfer, and phase change process. In permafrost and seasonal permafrost areas, due to environmental factors, permafrost roadbeds are prone to disasters such as uneven settlement, frost heave and thaw settlement, which are closely related to water migration, phase change, and temperature changes. In order to solve the problem of hydrothermal coupling in permafrost, this paper establishes a two-dimensional model of hydrothermal coupling using COMSOL Multiphysic software. The influence of temperature and pressure at the same time period under different head pressures, the influence on the total heat flux at the outflow boundary, the influence on the total liquid water volume, and the influence on the minimum temperature are studied. Three different head pressure working conditions are simulated. The results show that the temperature and pressure during the same time period change with the increase of head pressure; The total heat flux leaving the system at the outflow boundary also increases with the increase of head pressure; The larger the head gradient under different working conditions, the faster the melting rate of ice, and the shorter the time to reach the maximum total liquid water volume.

1. 引言

冻土的冻胀融沉一直是影响冻土区生态环境变化的重要因素之一，在寒冷地区，如果不加以适当的设计和处理，冻胀融沉可能会导致土壤松散、地基沉降、管道破裂等问题，给工程建设和基础设施造成困难，而多年冻土内部的水热过程与冻胀融沉密切相关。路建国等[1]对冻土水热力耦合研究现状及进展进行了阐述。

目前，张明礼等[2]通过COMSOL Multiphysic软件分析了冻土水分对流与温度变化的关系，并考虑了水的相变和对流的温度场与水分场的耦合。尚松浩等[3]根据冻结条件下土壤水分运动和热量迁移的基本方程，推导了冻土水热耦方程并构造了全隐式有限差分计算格式，使模型计算速度得到提高。白青波等[4]基于非饱和土渗流和热传导理论，通过对COMSOL Multiphysic软件的二次开发，实现冻土温度场和水分场全耦合的数值模拟。刘思源等[5]建立考虑水分迁移和相变及初始冰、水含量的路基水–热–力耦合模型，对多年冻土区不同结构路基温度、水分、变形特性以及长期融沉变形预测进行研究。韩昀希等[6]对多场耦合作用下富水砂层人工冻结状态进行了研究，并对冻土的未冻含水率进行实验分析。王蕴嘉等[7]提出了路基覆盖范围冻结锋面下凹扭曲的计算方法，可计算得到填筑完成后不同时刻不同工况下的冻结锋面下凹扭曲曲线，从而对路基填筑短期内冻土温度场的变化规律进行研究。曾照军等[8]研究热传导、相变、温度梯度、水分梯度和重力势的耦合作用，使用FreeFem++有限元软件进行水热耦合的数值模拟，分析了非饱和土水热场的演化特征。申林方等[9]研究表明土体在冻结过程中，未冻区的水分向冻结区移动，并在锋面位置处发生水分积聚。吴东军等[10]进行了人工冻土的冻结试验，研究冷端温度、初始含水率、上覆荷载对其的影响。王青林等[11]通过将渗透系数理论模型用于水热耦合计算，并分析其适用性，反映冻土渗透系数变化规律。

综上所述，国内学者包括蒋瑶钰[12]、郭常凯[13]等人采用COMSOL Multiphysic、FreeFem++等有限元软件分析冻土温度场和水分场耦合问题方面开展了大量的科研工作，但通过COMSOL Multiphysic分析不同水头压力下水热耦合的研究较少。本文研究旨在讨论冻土热液耦合系统的动力学机制，建立二维模型，分析不同水头压力下同一时间段温度和压力的影响、对流出边界总热通量的影响、对总液态水体积的影响，以及对最低温度的影响进行阐述，为多年冻土区工程设计和冻土环境演化数值模拟研究提供借鉴。

2. 数学模型

2.1. 传热方程

张东[14]等研究了非均匀多孔介质等效渗透率的普适表达式，而在局部热平衡假设下，需要描述整个多孔介质的平均温度，我们通常会使用非均质介质热传导方程来考虑多孔介质内部的复杂热传导过程，基于能量守恒和质量守恒以及应用混合规则，传热方程的偏微分形式如下：

${\left(\rho C\right)}_{eq}\frac{\partial T}{\partial t}+\rho C_{w}u\cdot \nabla T+\nabla \cdot \left(-k{}_{eq}\nabla T\right)=Q$ (1)

式中，密度$\rho$ (kg/m³)和恒压热容C的等效值是${\left(\rho C\right)}_{eq}$ ；Q表示热源；k_ep表示有效导热系数；T表示温度；C_w有效流体的恒压热容。

在COMSOL Multiphysic中，通过模拟多孔介质中冰融化的过程来反映冻土随温度变化的热液耦合问题，其中，流体和多孔介质的热特性被称为有效特性，即为有效导热率和有效体积热容：

${\left(\rho C\right)}_{eq}={\theta }_s{\rho }_s{C}_s+{\theta }_f{\rho }_f{C}_f$ (2)

式中：$\rho$ 为密度；C为恒压热容；${\theta }_s$ 为固体体积分数；${\theta }_f$ 为流体的体积分数：${\theta }_f=1-{\theta }_s$ ，其中，多孔介质的结构以及固体和流体的热导率决定了有效导热率。注：将动力黏度和水密度视为恒定且传热方程中不包含热弥散。

2.2. 水迁移方程

冻胀融沉问题一直是冻土区域的核心问题，而冻土的水热迁移能够很好的解释其产生的原因，对研究冻土形成过程中水分迁移规律具有十分重要的意义。非饱和正冻土指在冻土区域，土壤温度低于零度且处于非饱和状态，在这种情况下，土壤中的水分会由于温度过低而结冰，使土壤中只有少量液态水，这种状态称为非饱和正冻土状态。在土势作用下多孔介质的水分迁移过程符合达西定律[14]：

$v=-K\nabla h$ (3)

$K=k\frac{\gamma }{\mu }$ (4)

式中，K为渗透系数；$\nabla h$ 为水力坡降；k为渗透率；v为水流流速；$\mu$ 为流体动力粘滞系数。

利用COMSOL Multiphysics软件模拟水热耦合的过程，其中，根据“达西定律”接口中的“多孔介质”特征包含用于定义存储量S_p；在软件中定义一个质量源，质量源通常指代一种能够引入或移除流体或热量的装置或过程，它可以是液体或气体的注入或抽出，也可以是热量的释放或吸收。而在COMSOL中则定义为由于冰夹杂物融化而产生的额外液态水。

$S_p=S_w\beta {\epsilon }_p$ (5)

$Q_m=S_w{\epsilon }_p\left({\rho }_i-{\rho }_w\right)\frac{\partial S_w}{\partial t}$ (6)

$S_w=S_{w,res}+\left(1-S_{w,res}\right){\theta }_2$ (7)

式中，${\rho }_i$ 和${\rho }_w$ 是冰和液态水的密度。液态水饱和度S_w取决于相变，$S_{w,res}$ 是液态水残余饱和度，${\theta }_2$ 是相变材料节点中定义的平滑阶跃函数。

2.3. 联系方程

传热方程和水分迁移方程中含有多个未知数，若要求解多个方程组则需要建立液态水饱和度与相变之间的关系：

${\theta }_2=\{\begin{array}{c}0;\\ 1;\end{array}\begin{array}{c}大于T_{pc}\\ 大于T_{pc}\end{array}$ (8)

式中，T_pc为相变温度，${\theta }_2$ 为平滑阶跃函数。

3. 材料各属性及工况

3.1. 材料参数

为了在COMSOL Multiphysics软件中模拟冻土融化的过程，以此为背景建立多孔介质中冰夹杂物融化的二维有限元模型。建立一个长3 m高1 m的多孔介质通道，内含30 cm的方形冰块，外界温度为5℃，冰块的温度为−5℃，水流通过多孔介质的孔隙从一端流入另一端，并经过冰块将其热量带走。这涉及温度场与压力场耦合的问题，流体通过对流传热进行传热，它是依靠流体质点的移动进行热量的传递，与流体的流动密切相关。压力场则由通道长度上的水头梯度驱动，通过改变水头梯度来控制压力的变化。各流体、固体的材料属性见表1，其他物理参数见表2。

Table 1. Material Properties

表1. 材料属性

材料属性	冰	水	固体基质
密度	920 kg/m^3	1000 kg/m^3	2650 kg/m^3
恒压热容	2060 J/(kg·K)	4182 J (kg·K)	835 J (kg·K)
导热系数	2.14 W/(m·K)	0.6 W/(m·K)	9 W/(m·K)

Table 2. Other Parameters

表2. 其它参量

参数	值
水的动力黏度	1.793·10^−^3 Pa·s
孔隙率	0.37
有效压缩系数	10^−^8Pa^−^1
残余水饱和度	0.05
阻抗因子	50
固有渗透率	1.3·10^−^10 m^2

3.2. 计算工况

考虑到不同的水头压力对温度场和压力场的影响，利用COMSOL Multiphysics软件对通道长度上不同的水头梯度进行计算分析，分别得到3种状况下同一时间段时，温度和压力的变化，模型域中的最低温度随时间变化，通过流出边界离开系统的总热通量随时间变化以及总液态水随时间的变化。计算数值及工况见表3。

Table 3. Calculation conditions under different head pressures

表3. 不同水头压力下的计算工况

计算工况	值
工况一	通道长度3%的水头压力
工况二	通道长度4%的水头压力
工况三	通道长度6%的水头压力

4. 不同水头压力对热液耦合的影响

4.1. 不同水头压力对同一时间段温度的影响

(a) 工况一

(b) 工况二

(c) 工况三

Figure 1. Temperature distribution surface map

图1. 温度分布表面图

为研究不同水头压力下多孔介质中冰块融化的温度变化，分别绘制了各工况融化9小时后的温度云图，见图1。冰块所处的温度最低，因为流体具有一定的流速对温度场有所影响，涉及到流体动力学、热传递和物质传递等多个物理场过程，流体的流动可以增强冰与周围环境之间的热交换，从而加快冰的融化。流向表面的流体流动比平行于表面的流动产生更高的熔融度，这是因为流体的对流传热效应，所以在通道的下游产生低温。从图1中可以清晰的看出不同工况下同一时间段时，通道下游的温度分布不均匀，沿着流体流动的方向低温区域逐渐扩大，表明不同的水头压力会导致对流传热效应的增加，从而加快冰块的融化。

4.2. 不同水头压力对同一时间段压力的影响

模型模拟了冰块从开始到完全融化的全部过程，较低的温度从通道中对流出来。冰融化会引起液体体积的增加，而液体体积的增加会影响其密度和压缩性，所以压力分布会受到影响。冰融化还会导致流体的扩散，在一定条件下，这种流体扩散可能会改变流体的流动和压力分布。见图2，为不同工况下同一时间的压力分布。工况一最大压力为883 × 10³ Pa，分布在通道的上游，最低压力为−1.15 × 10³ Pa，分布在冰块所在位置；工况二最大压力为1.18 × 10³ Pa，最低压力为−920 Pa；工况三最大压力为1.77 × 10³ Pa，最低压力为−522 Pa。随着水头压力的增加，压力较高的流体会向周围区域扩散，从而使压力分布趋于均匀。

(a) 工况一

(b) 工况二

(c) 工况三

Figure 2. Pressure distribution diagram with ice

图2. 含冰时的压力分布图

4.3. 不同水头压力对总热通量的影响

总热通量是指在一定时间内通过某一特定区域的热量总和，在COMSOL Multiphysics中指的是通过流出边界离开系统的总热通量随时间的变化，见图3。总热通量可以表示热量通过物体或系统的速率，由图可知，工况一离开系统的总热通量在10 h左右急剧增加，当总热通量达到220 W时达到峰值，时间是开始融化后25 h，随后，总热通量的值逐渐减小，直到降为0为止；工况二离开系统的总热通量在9 h左右急剧增加，当总热通量达到270 W时达到峰值，时间是开始融化后18 h，之后总热通量逐渐减小至0；工况三离开系统的总热通量在8 h左右急剧增加，当总热通量达到350 W时达到峰值，时间是开始融化后13 h，之后总热通量逐渐减小至0。由于水头压力的不断增强，流体的流速加快，导致单位时间内通过流出边界的热量较多，而当冰完全融化后，物体处于平衡状态，不再发生热量的传递，所以总热通量都会趋向于0。结果表明不同水头压力加剧了流出边界离开系统的总热通量，且其峰值大小与水头梯度有关，水头梯度越大峰值越高，当冰完全融化后，总热通量趋向与0。

Figure 3. Schematic diagram of the variation of total heat flux over time

图3. 总热通量随时间的变化示意图

4.4. 不同水头压力对总液态水体积的影响

冰夹杂物的融化势必会增加液态水的体积，所以在冰夹杂物融化之前，总液态水的体积是随时间不断增加的，工况一冰块融化后20 h内，总液态水的体积不断增加直到1.11 m³，当完全融化后，总液态水的体积则不再增加，温度的变化不再影响总液态水体积；工况二冰块融化后16 h内，总液态水的体积不断增加直到1.11 m³；工况三冰块融化后14 h内，总液态水的体积不断增加直到1.11 m³，见图4。

Figure 4. Schematic diagram of the change in total liquid water volume over time

图4. 总液态水体积随时间的变化示意图

4.5. 不同水头压力对最低温度的影响

冰块融化的总时间为56 h，温度为−5℃~5℃，在该模型域中，冰块吸收热量温度上升，到达熔点−1℃后温度保持恒定，直到冰块完全融化后继续吸收热量，则液体的温度会继续上升，直到与外界的温度相一致时停止变化，见图5。工况一冰块在18 h后完全融化，工况二和工况三冰块分别在16 h和14 h后完全融化，随后继续吸收热量，达到与外界温度相一致，表明压力水头越大，液体的流速越快，冰块融化的时间也越短。

Figure 5. Schematic diagram of the variation of minimum temperature over time

图5. 最低温度随时间变化的示意图

5. 讨论分析

通过上述理论分析基础以及模拟过程，我们深入探讨了不同水头压力对冰块融化的影响。通过对传热方程、水迁移方程的分析以及相关文献[8]的对照，发现模型对参数水头压力的依赖性最强，各工况在其他条件不变的情况下，只改变水头压力的大小，其结果的变化趋势相同，与吴东军[10]等研究的冻土水热迁移试验及水热耦合数值模拟规律相似，这说明其模拟的结果符合实际，对未来多年冻土水热耦合的研究具有一定的价值。

6. 结论

本文以COMSOL Multiphysics软件建立水热耦合的二维模型，为多年冻土区随气候变化的热液耦合过程提供依据。基于达西定律与包含相变的“多孔介质传热”接口耦合，研究不同水头压力作用下同一时间段温度和压力的影响、对流出边界总热通量的影响、对总液态水体积的影响，以及对最低温度的影响进行分析，具体结论如下：

1) 同一时间段的温度和压力随着水头压力的增加而变化。水头压力越大，对流传热效率越高，沿着流体流动方向的低温区域越大，同时，压力较高的流体会向周围区域扩散，从而使压力分布趋于均匀。

2) 同一时间段流出边界离开系统的总热通量随着水头压力的增大而增大，且不同工况下的水头压力加剧了流出边界离开系统的总热通量。

3) 不同工况下总液态水体积和最低温度也受到影响，水头梯度越大，冰块融化速度越快，达到总液态水体积最大值的时间越短。

NOTES

^*通讯作者。