多圈管冻结壁水热时空分布数值模拟研究

罗阳， 王浩， 倪斌

（1.安徽理工大学土木建筑学院，安徽 淮南 232001；2.安徽寿县农村商业银行股份有限公司，安徽 淮南 232281）

0 引言

人工冻结法是一种使用人工制冷技术对地层进行人工制冷，使得天然岩土变成冻土，形成人工冻结壁，由于冻土的高抗压抗剪性能以及良好的阻水性能，在人工冻结壁内侧形成了良好的施工环境，有利于地下工程的施工［1］。但人工冻结壁［2］的形成过程是一个复杂的多物理场过程，其中主要包括温度场、水分场、应力场的相互作用过程，这个复杂的多物理场过程存在于同一圈不同冻结孔之间，也存在于不同冻结圈之间以及冻结壁与天然岩土层之间。

目前国内人工冻结法的应用主要在一些大型矿井的开挖以及隧道掘进等工程上，人工冻结法根据主冻结圈+辅助冻结圈的冻结形式，从最初的单圈管冻结发展到双圈管、三圈管、四圈管。国内学者［3］对多圈管冻结的研究主要为室内模型试验，通过模型试验研究多圈管冻结壁的形成过程，还有学者通过数值仿真软件研究了多圈管冻结壁温度场发展规律。汪仁和等［4］开展了大型的人工多圈管冻结模型试验，研究了多圈管冻结壁的形成和融化过程中的温度场、水分场、应力场，提出了导水系数是温度梯度的函数的结论。李栋伟等［5］利用ADINA 软件进行了多圈管冻结瞬态温度场的有限元数值仿真，获得了冻结壁的有效厚度、冻结壁的平均温度，与现场数据较为吻合。认为有限元温度场数值仿真可以较好地预测多圈管冻结的温度场。陈军浩等［6］进行了多圈管冻结壁的形成与融化研究，进行了多圈管模型优化设计的室内模型试验。

文中通过对COMSOL Mutiphysics 软件进行了二次开发，实现了多圈管冻结壁的水热耦合数值仿真，通过与现场数据对比，表明正冻非饱和土水热耦合模型可以用于描述多圈管冻结过程中的水热分布。通过对数值模拟结果进行分析，进行了多圈管冻结壁各特征点温度、冰水质量、冻结时间关系分析，以及主界面水热时空分布分析。

1 正冻非饱和土水热耦合模型

1.1 热量迁移方程

基于能量守恒定律和热传导理论，将冰水相变作为内热源处理，建立如下热量迁移方程［7］：

式中，ρ 和ρi为土的密度和冰的密度，kg/m3；C（θ）为体积热容，J/（kg℃）；λ（θ）为导热系数，W/（m℃）；T为土体的瞬态温度，℃；t为时间，s；△为微分算子；θ为体积含水率；θi为孔隙冰体积含量；L为相变潜热，通常取值为334.5kJ/kg。

1.2 水分迁移方程

根据Richard方程［8］，基于质量守恒定律和达西渗流理论，非饱和冻土中的未冻水的水分迁移方程：

式中，D，K 分别为土扩散系数和导水率，均为未冻水含量θu的函数。

1.3 联系方程

通过对徐学组建立的未冻水预测模型变形，获得固液比与温度的关系［9］。

冰的体积含量θi三者的联系方程：

以上，热量迁移方程、水分迁移方程、联系方程3个方程构成的方程组建立了正冻非饱和土的水热耦合模型。

2 多圈管冻结水热耦合模拟

2.1 几何建立

将水热耦合数值模拟理论应用于皖北地区某矿副井的多圈管冻结方案，几何模型为副井平面模型的1/4，冻结孔、冻结孔圈尺寸、和冻结孔位置均按实际工程取值。

2.2 边界条件及物理参数

土体所有外边界假定为绝缘边界，内部冻结孔温度场边界取实际现场循环液去、回路温度的平均值，内部冻结孔温度边界曲线如图1所示。水分场边界假定为绝缘边界，冻结壁周围初始含水率为28%。深厚黏土热量特征参数、土水分特征参数见表1、表2。

图1 冻结孔温度边界

表1 深厚黏土热量特征参数

表2 深厚黏土土水分特征参数

2.3 数值模拟结果及分析

2.3.1 冻结壁厚（180d）。

当达到180d，内边界与辅助冻结圈间距大于冻结壁外边界与主冻结圈的间距，这是因为冻结壁内侧土体在主孔包裹下始终封闭状态，外侧则与外圈的土体发生热交换。

图2 不同时刻冻结壁发展

2.3.2 温度场

图3演示了冻结孔附近土体温度在不同时间内的发展情况：冻结孔降温到指定温度（30d）→主冻结孔圈降温交圈（60d）→辅助孔圈降温交圈（90d）→主、辅助孔温度分布呈规则圆环（120d）→温度分布均匀且有所回升（180d）。

图3 不同冻结时刻温度场分布

在整个过程中可以看出，整个过程中中心土体的温度始终降低直到负温。在冻结的第180d时，主冻结圈和辅助冻结圈间的土体温度有所上升，温度在-15℃左右，这与循环液温度上升，冻结进入养护阶段有关。

2.3.3 水分场

（1） 未冻水含量。图4 演示了冻结孔周围土体含水量随着时间的变化过程，在整个过程可以很明显的看到：①主、辅助孔中间土体含水量最先降低，从开始的14%（30d）降低至8%（60d）直到降至最低含量6%（90d、锯齿状），在冻结的第180d，主、辅助孔之间未冻水含量低于6%土体范围逐渐减小，呈鱼刺状；②主、辅助孔边界土体未冻水含量低于10%的边界土体呈规则环形向内和向外发展。

图4 不同时刻未冻水分布

（2） 含冰量。根据图5 不同时刻冰水分布图分析可知：①在冻结的第30d，冻结孔附近土体含冰量达到24%左右，仅在冻结孔附近很小的范围内有土体冰水含量达到30%以上，在主冻结圈外圈有一环带含冰量小于18%，主冻结圈和辅助冻结圈间的土体含冰量低于10%，呈鱼刺状分布；②在冻结的第60d，主冻结圈和辅助冻结圈间低于18%的土体鱼刺状分布消失，在主冻结圈外和辅助冻结圈内土体呈冰水含量低于18%的环形带，随着冻结的持续，低于18%的环形带向内和向外移动，在主冻结圈外侧形成一条冰水含量为25%的均值环形带；③整个冻结过程中，含冰量大于30%的土体仅存在在冻结孔附近的较小范围内。

图5 不同时刻冰水分布

3 数值结果与现场数据对比

3.1 温度变化

从图6可以看出：①1#、3#测点为主冻结圈外1m处测点，数值计算结果为二者的平均值，变化规律一致，2#测点为主、辅助冻结圈中间测点，数值模拟结果与实测结果在冻结的40～60d 完全重合，曲线变化规律一致；②在试验后期温度均呈上升趋势，数值计算结果低于实测结果，这和计算模型中忽略土体自身温度随季节性变化相关。

图6 1#、3#、2#温度变化

基于数值模拟的计算结果与现场实际情况较吻合，数值模拟的计算过程能够清楚描述多圈管冻结过程中的热时空分布。

3.2 水-热-时间关系

根据图7分析可得：①主、辅助冻结圈中间土体与辅助冻结孔中间土体温度下降过程中出现明显的相变潜热释放过程；②冻结180d 后，土体温度由低到高排序：主冻结圈与辅助冻结圈位置、主冻结孔中间位置、辅助冻结孔中间位置、井帮位置、主冻结圈外1m处。前三的温度值相差不大，均为-15℃左右，井帮位置温度为-10℃左右，主冻结圈外1m 处温度为-6℃左右；③各特征点未冻水含量随时间变化关系为：主、辅助冻结圈位置、主冻结孔中间位置、辅助冻结孔中间位置三处未冻水较早时间内达到平衡状态主冻结圈外1m 处未冻水全过程处于缓慢下降阶段，井帮位置处未冻水变化经历慢-快-慢三阶段；④所有特征点的冰水含量变化均为先减小后增大。开始增大时该处土体温度为冻结温度，即土体降至冻结温度时，土体冰水含量开始增大。主冻结圈外1m处土体冰水质量含量最大，达到23.5%左右，主冻结中间位置土体冰水质量含量最小，为19.5%左右，其他三处为20%左右。

图7 不同特征点水-热-时间关系

3.3 主、界面温度、冰水含量分布

从图8可以看出：①冻结初期，从井筒中心到冻结壁外侧，主面温度：高、低、高、低、高，随着冻结的持续，在主冻结圈和辅助冻结圈之间的土体温度逐渐变得平缓，最终达到和冻结孔一致的温度；②界面温度分布与主面分布变化一致，但界面温度较主面温度分布平缓。

图8 主、界面温度分布

从图9可以看出：①主面和界面的水分分布为高、低、高、低、高、低、高。随着冻结的持续，在冻结壁的外侧出现冰水含量二次高点，这是由于冻结壁向外发展较缓慢，在冻结壁边缘出现水分集聚现象；②同半径距离，主面的冰水含量高于界面，冰水含量大于28%的土体范围仅存在各冻结孔附近直径约0.3m 范围内；③冻结壁两侧的低含水带随冻结壁的增大而向两侧移动；④主、辅冻结圈中间的土体水分在冻结初期呈减小状态，后随着冻结的持续，水分有所增大，在冻结60d后无明显变化。

图9 主、界面冰水分布

4 结语

（1） 基于正冻非饱和土水热耦合模型，进行COMSOL Multiphysics 软件二次开发用于水热耦合计算，能够较好地描述多圈管冻结过程中的水热时空分布，对实际工程有一定指导意义。

（2） 主冻结圈+辅助冻结圈形式的双圈管冻结，主冻结圈首先发生交圈，主冻结圈完全交圈后，辅助冻结孔先于主冻结圈交圈再与临近辅助冻结孔交圈，在交圈的过程中，主、辅冻结圈间会出现局部未冻仓。

（3） 未冻水含量较低的土体首先出现在主冻结孔之间，随着冻结的持续，在主、辅冻结圈间呈锯齿状分布，后由于冻结从积极期进入维护期，锯齿状分布变为鱼刺状分布，范围减小。

（4） 冻结壁范围内各特征点处的冰水含量均经历先减小后增大两个阶段，转折点出现在该处土体温度处于冻结温度，即土体冻结后，土体水分含量开始增大。

（5） 主面温度分布较界面温度分布平缓，在冻结维护期，无论主面还是界面，主、辅冻结圈间的土体温度均与冻结孔处温度一致。

（6） 同半径范围，主面的冰水含量要高于界面，随着冻结的持续，在冻结壁外侧出现冰水含量次高峰，冻结壁两侧的低含水带随冻结壁的增大而向两侧移动，主、辅冻结圈中间的土体水分在冻结初期减小，后随着冻结的持续，水分有所增大，在冻结60d后无明显变化。