渗流作用下砂层冻结模型试验研究

苏彦林,岳祖润,李晓康,张 松,周 圆

(1.石家庄铁道大学土木工程学院,石家庄 050043； 2.石家庄铁道大学道路与铁道工程安全保障省部共建教育部重点实验室,石家庄 050043； 3.石家庄铁道大学研究生学院,石家庄 050043)

1 概述

人工冻结技术是针对松散、软弱含水地层，采用人工制冷进行预加固的施工方法，起到封水、护壁、抵御围压等临时支护作用，具有安全性高、环保节能、灵活高效等优点，广泛应用于矿井工程、隧道及城市地下空间建设[1-3]。在地下水源丰富、水流速度较大的地层进行人工冻结时，水流对冻结壁进行冲蚀，致使冻土帷幕难以封闭或无法达到设计要求厚度，从而导致工程事故发生，例如深圳地铁A标段暗挖隧道采用冻结法施工时，由于冻结区地下暗河水流速度过大，导致冻土帷幕的形成出现问题[4]；东欢沱矿二号井、焦作位村矿副井因水文管水压差异造成地下水流动导致冻结壁不能交圈[5]。

近年来，国内外学者对渗流地层常见冻结工况下冻结壁的扩展规律进行了大量研究。在地下水流速对冻结效果的影响方面：M Vitel[6]通过耦合热液压数值模型提出，当地下水流速超过1～2 m/d时，冻结管之间的土体无法冻结。Endo K[7]由现场实测数据得出，地下水流速达到5 m/d时，地层进行冻结施工时就应借助灌浆来降低水流速，帮助冻土帷幕稳定发展。李方政[8]通过分析北京地铁6号线冻结壁不交圈的地质与环境影响因素，提出在冻结设计时，地下水流速大于5 m/d 必须采取针对性的措施。Ullah S H[9]利用Comsol中的Richards方程与多孔介质传热模型构建了二维有限元计算模型，得出了2.0 m/d的极限冻结流速。芮大虎[10]通过模型试验指出随渗流速度增大冻结壁交圈时间成倍延长，且当流速大于4.3 m/d时，冻结壁交圈困难。刘建刚[11]通过建立渗流场和冻结温度场的耦合数学模型，指出横向水平流速的大小明显影响冻结壁形态和交圈时间，流速大于1.5 m/d 后交圈时间明显增加。在地下水渗流对冻结壁发展及温度场分布规律方面，周晓敏[12]通过室内模型试验研究表明，水流速度是影响冻结壁上游厚度的主要因素。吉植强[13]通过模型试验指出冻结壁厚度受渗流影响不均匀分布，且当渗流到达某一速度后冻结壁不能交圈，甚至不能发展。杨平、皮爱如[14]建立了地下水流动时冻结锋面移动的数学模型，分析了温度场以及渗流场的变化规律。

目前专家学者对地下水流存在的特殊地层冻结研究取得了一些成果，但是针对地下水渗流条件下冻结壁厚度的计算方法则鲜有报导。本文通过模型试验，模拟地下水渗流条件下砂层冻结过程，得到温度场分布规律，同时基于试验数据，对静水条件下的巴霍尔金温度场解析解[15-16]进行修正，得到适用于渗流条件下冻土帷幕厚度的计算方法，为冻结法施工效果评估提供理论指导。

2 模型试验设计

2.1 试验参数缩比

根据渗流地层冻结过程模型试验相似准则[17-18]，采用量纲分析法对各参数进行推导，得到试验相似缩比如下。

(1)几何相似

考虑冻结工程中冻结管尺寸和间距、冻结影响范围、冻结温度场，结合室内试验的可操作性，采用Cl=1∶10的几何缩比，模型试验中选取冻结管尺寸为φ10.8 mm×0.9 mm，模拟工程常用的φ108 mm×9 mm型号冻结管。冻结工程中冻结管间距一般为900 mm，由几何缩比确定试验中冻结管间距为90 mm。冻结影响范围为冻结管间距的3～6倍[19-20]，为尽可能减小边界效应对试验结果的影响，本模型箱设计为直径1 200 mm，高度1 000 mm的圆形箱体。

(2)温度相似

根据科索维奇准则，得到CT=1，即模型温度和原型温度一致，因此，冷源温度、土体温度、循环清水温度均与实际工程一致。

(3)时间相似

(4)冻结管冷媒流量相似

根据相似准则，冻结管中盐水流量CQ=Cl，即模型试验中冻结管内盐水流量1 m3/h相当于工程中10 m3/h。

2.2 模型试验装置

本试验针对富水砂层渗流冻结问题进行研究，试验装置包括模型试验箱、冻结制冷系统、恒温供水系统、温度数据采集系统，装置整体示意如图1所示。

图1 试验装置示意

试验箱体内采用底部供水、顶部出水方式模拟土层中地下水渗流过程，以保证土体内均匀渗流。箱体底部中心处设有进水口，顶部侧壁处设置出水孔，箱内土层自下而上依次为：150 mm厚米字形碎石缓冲层，700 mm厚试验砂土层，150 mm厚碎石缓冲层。砂土层与碎石缓冲层间铺设80目不锈钢滤网，试验砂土物理参数如表1所示。箱体外部包裹50 mm厚度的保温棉进行保温隔热，试验箱内部示意如图2所示。

表1 试验土体物理参数

图2 试验箱体内部示意(单位：mm)

冻结制冷系统采用螺杆式制冷机组配合变频泵，为冻结管提供所需盐水温度和流量。恒温供水系统采用AC200冷浴机对供水箱中清水进行恒温处理，通过变频泵和电磁流量计为试验箱体提供各种渗流速度的恒温清水。温度数据采集系统由DAM温度采集模块、PT100温度传感器组成，对土体温度进行实时采集。

试验箱体内冻结管布置平面如图3(a)所示，3根冻结管采用串联方式连接，横向平行布置在试验砂层1/2高度的平面上。如图3(b)所示，测温点布置在垂直于冻结管的测温面上，各点位综合考虑冻结管间距及冻结影响范围进行优选布置，共计31个。

图3 冻结管和测温点平面布置(单位：mm)

2.3 试验方案

按照GB/T50123—2019《土工试验方法标准》对试验砂土进行处理后，采用分层填土方法，将砂土分4层装填到试验箱中，装填高度为700 mm，每层装填完毕后，向试验箱内充水使砂土保持湿润，降低试样松散性，然后夯实达到试验要求干密度，全部装填完毕后打开供水系统进行渗流循环，使箱体内各点温度保持20 ℃。打开制冷系统对盐水进行降温，当盐水温度达到-20 ℃时，准备进行冻结试验。

试验分为静水冻结和渗流冻结两部分。首先，基于静水冻结试验数据结合巴霍尔金解析解验证模型试验可靠性，进而开展渗流冻结试验，根据试验中不同渗流速度作用下冻结壁发展状况，确定4组渗流冻结工况，渗流速度分别为1，1.5，2.3，2.5 m/d，当砂层上下游冻结壁稳定后停止试验。

3 试验结果分析

3.1 试验结果及可靠性验证

(1)模型试验共进行5次，每次历时约10 h，试验结果汇总于表2。从表2可以看出，当渗流速度逐渐增大时，冻结壁交圈时间变长，交圈位置越偏向下游，轴面与界面交点(T8)处温度降到0 ℃用时越长。当冻结壁厚度基本稳定时，冻结体整体呈上游薄，下游厚的状态，随渗流速度的增大，上下游冻结壁厚度差异性愈加明显。渗流速度为2.3 m/d时，上游冻结壁平均厚度仅占下游厚度的33%，渗流速度为2.5 m/d时，冻结2 h后，冻结锋面热交换处于平衡状态，冻结壁不再发生明显变化，上下游温度场基本维持稳定状态，冻结壁难以交圈，可得到不影响冻结交圈的极限渗流速度介于2.3～2.5 m/d。由于试验土层孔隙率为42.7%，计算得到地下水流速介于5.39～5.85 m/d。DG TJ08—902—2006《上海市旁通道冻结法技术规范》规定地下水流速大于5 m/d时应采取针对措施。可见，本模型试验结果与已有研究成果和施工规范相比差异性不大，考虑到模型缩尺效应的影响及试验过程中产生的误差，本模型得到的极限渗流速度与实际冻结工程中可能存在一定的偏差，但无量级上的差距，可为常见施工工况下冻结法设计提供参考。

表2 试验结果

(2)人工冻结工程中冻结壁厚度是评判冻结效果的一个重要指标。理论研究和实际工程应用表明，在单排管冻结交圈后，巴霍尔金温度场解析解在各种温度场理论中计算结果较为准确[15]，其冻结温度场计算模型如图4所示。巴霍尔金认为在冻土柱交圈之后形成的波浪形冻结帷幕很快就会因为水力冲刷而填平，因此冻结壁侧表面可近似地以平面代替，并提出单排管直线冻结壁交圈后温度场解析公式

(1)

由式(1)，推导出冻结壁厚度计算公式如下

(2)

式中，t(x，y)为冻土区域内计算点的温度，℃；tCT为冻结管外表面的温度，℃；r0为冻结管的外半径；ξ为单管冻土柱半径；x，y均为单管坐标系中计算点坐标；L为冻结管间距。其中

图4 巴霍尔金单排管冻结温度场计算模型

(3)为评估本模型试验的可靠性，以巴霍尔金温度场计算模型为依据，选取静水条件下冻结壁交圈(冻结1.8 h)时、冻结5 h及冻结壁稳定后(冻结9 h)的温度场进行验证。选取主面上T23(0，45)测温点为计算点，此点降温曲线如图5所示。

图5 T23测点降温曲线

通过式(2)，得到巴霍尔金理论计算的冻结壁厚度与模型试验中冻结壁厚度的对比结果，如图6所示，可看出模型试验中冻结壁厚度与巴霍尔金理论计算得到的冻结壁厚度高度吻合且平均误差不超过4%。由此可得，本模型试验模拟现实工况准确性较高。

图6 冻结壁厚度对比

3.2 特征面温度分布规律

为描述流水作用下砂层冻结温度场演变规律，绘制冻结壁交圈时主面和界面上、下游温度分布特征曲线，如图7所示。由图7可得：(1)静水冻结时，上、下游温度曲线在主面和界面上均具有明显对称性，温度梯度基本一致；(2)渗流冻结时，主面、界面上游温度明显高于下游，随渗流速度增大，水流和冻结管之间的热交换逐渐加快，上游温度随之增大，下游温度随之减小；(3)由图7(b)可看出，冻结壁交圈位置处于下游区域，且随渗流速度的增大，交圈位置距轴面与界面交点处距离越远。

图7 主面、界面温度分布特征

3.3 冻结温度场分布

为分析砂层在冻结过程中温度场分布特征，通过各测温点实测数据，绘制冻结壁交圈及稳定时温度场分布云图。

(1)静水条件下冻结温度场云图如图8所示。

图8 静水条件下冻结温度场云图

由图8可得，冻结1.8 h时，冻结壁于轴面与界面交点处交圈，砂层中温度于冻结管列两侧对称分布，此时冻结壁呈现波浪形，受冻结管热交换相互叠加影响，中间管形成的冻土帷幕厚度大于两侧。冻结9 h时，冻结壁基本稳定，冻结管列两侧温度对称性更加明显，冻结体形状近似椭圆形。对比冻结过程中的温度云图，发现随着冻结时间的延长，冻土体波浪形侧表面逐渐变得平滑，主面和界面冻结壁厚度差值逐渐缩小，这种现象与巴霍尔金基础理论基本吻合。

(2)渗流冻结时，以1.5 m/d渗流条件下的冻结温度场为例，如图9所示，冻结2.4 h时冻结壁交圈，冻结管列两侧温度不再对称分布，管列上游温度梯度明显大于下游。冻结过程中水流对上游冻结壁不间断冲刷，冻结管产生的部分冷量被携带到下游，导致下游降温区平均温度低于上游，冻结壁于轴面与界面交点下游位置处交圈，且整体呈现凹槽形。冻结8.6 h时冻结壁基本稳定，管列主面和界面冻结壁厚度基本一致，由于冻结体成型后受水流的冲刷作用和冻结管间温度相互叠加影响，中间冻结管下游冻结壁厚度明显大于两侧冻结管下游冻结壁厚度，冻结体整体近似呈扁桃形。

图9 渗流条件下冻结温度云图

4 渗流条件下冻结壁厚度计算

4.1 冻结壁厚度对比分析

鉴于巴霍尔金温度场计算模型仅适用于静水冻结条件，而关于渗流条件下冻结壁厚度计算的研究则鲜有报道，而渗流条件下冻结壁厚度对冻土帷幕状态评估至关重要。为此，针对渗流冻结时冻结壁厚度计算问题，通过分析渗流作用对管列上、下游冻结壁厚度影响，得到交圈后冻结壁厚度与静水冻结时冻结壁厚度的比值随时间变化曲线。为方便分析与计算，本文采用k值表示渗流条件下冻结壁厚度与静水条件下冻结壁厚度的比值。

(1)如图10(a)所示，渗流冻结过程中，上游冻结壁厚度比值k随冻结时长增长呈现降低趋势，随着渗流速度增大，比值k整体减小。下游冻结壁厚度比值k随冻结时长增长，先增大后减小。在冻结初期，渗流速度越大，单位时间内水流携带的冷量越多，下游冻结壁厚度越大，比值k整体先增大，当冻结壁逐渐趋于稳定时，水流(20 ℃)对冻结壁的冲蚀作用逐渐占主导地位，且随着渗流速度的增大，冲刷强度越大，单位时间内冻结管和水流之间的热量交换加快，比值k整体逐渐减小。

(2)如图10(b)所示，冻结壁发展稳定后，上、下游冻结壁厚度比值k随渗流速度增大，均呈逐渐降低趋势，渗流速度对上游冻结壁厚度影响明显大于下游。

图10 渗流作用对上、下游冻结壁厚度比值影响

4.2 渗流条件下冻结壁厚度计算

(1)基于渗流冻结时上、下游冻结壁厚度与静水冻结时冻结壁厚度随时间和流速变化的比值曲线，拟合得到函数k(τ)和k(v)，对巴霍尔金冻结温度场解析解进行修正，得到渗流条件下上、下游冻结壁厚度的计算公式。为保证冻结壁厚度计算结果更加准确，选取中间冻结管主面测温点T9点和T23点为计算点t(x，y)和t′(x，y)进行分析，得到：①冻结过程中，冻结壁厚度随时间发展的计算公式；②冻结壁稳定时，冻结壁厚度随流速变化的计算公式。

图11(a)为计算点(T9，T23)在静水冻结和渗流冻结时温度差值Δt随时间τ变化的曲线，经过拟合分析，得到温度差值关于冻结时间的函数Δt(τ)。以上游冻结壁厚度计算为例，t(x，y)+Δt(τ)即为某时刻静水冻结时计算点的温度，代入式(2)可得静水冻结时管列两侧冻结壁厚度，其与函数k(τ)的乘积即为渗流条件下上游冻结壁厚度的计算解。因此渗流冻结时，冻结壁交圈后上、下游厚度随时间τ发展的计算公式为

(3)

(4)

式中，ξ、ξ′为上、下游冻结壁厚度；t(x，y)为渗流作用下冻土区域内计算点的温度，℃；(x，y)为坐标系中冻结管周围计算点坐标。

图11(b)为冻结壁发展稳定后，计算点T9、T23在静水冻结时和渗流冻结时温度差值Δt随渗流速度v变化的曲线，经过拟合分析，得到温度差值关于渗流速度的函数Δt(v)，该过程与冻结壁厚度随时间τ发展的计算公式的推导方法一致，可得渗流冻结时，温度场稳定后，上下游冻结壁厚度随流速变化的计算公式为

(5)

(6)

图11 渗流作用对计算点温度影响

(2)为验证上述计算公式的准确性，进行以下拟合分析。

①以1 m/d渗流速度为例，对图10(a)、图11(a)中1 m/d比值曲线及温度差值曲线拟合分析，得到的函数如下

冻结2.5 h时，测温点T9为13.66 ℃，T23为1.69 ℃，将上述函数代入式(3)、式(4)，计算得到上、下游冻结壁厚度分别为27.1，37.6 mm，与试验中上、下游冻结壁厚度28，37.4 mm相比，误差均不超过4%。

②以v=1.5 m/d渗流速度为例，对图10(b)、图11(b)中比值曲线及温度差值曲线拟合分析，得到的函数如下

k(v)=1-0.45v+0.096v2

k′(v)=1-0.24v+0.063v2

Δt(v)=-0.02-9.4v+1.65v2

Δt′(v)=-0.046-2.85v+0.72v2

冻结壁稳定后，测温点T9为7.81 ℃，T23为-2.98 ℃，将上述函数代入式(5)、式(6)，计算得到上、下游冻结壁厚度分别为36.2,62.6 mm，与试验中上、下游冻结壁厚度37，63.1 mm相比，误差均不超过3%。

结合上述算例分析，可知本文提出的渗流条件下冻结壁厚度公式计算结果与试验实测值吻合程度较高。在渗流地层冻结工程中，由于施工工况不同，k值函数及Δt函数可能不同，而本模型试验是模拟冻结法常见施工工况，因此得出的计算结果具有较高的适用性。针对其余渗流地层冻结工况，只需进行模型试验，通过本文提出的计算方法即可得到实际冻结工程施工过程中冻结壁的厚度。

5 结论

(1)静水冻结时，冻结管两侧温度场对称性明显，冻结壁交圈后的温度场与巴霍尔金单排管冻结温度场计算模型吻合程度较好，表明本模型试验研究具有较高的可靠性。

(2)渗流冻结时，上、下游温度场呈现较大差异性，冻结体内温度最低点位于冻结壁交圈位置处，随渗流速度的增大，冻结壁交圈时间越长，交圈位置更加偏向下游；温度场稳定后，冻结壁呈上游薄，下游厚的状态，随渗流速度的增大，这种差异性越加明显；渗流速度为2.3 m/d时，上游冻结壁厚度仅占下游厚度的33%；渗流速度为2.5 m/d时，冻结壁无法交圈，考虑到模型缩尺效应影响，试验结果与实际工程相比可能会存在一定偏差，但无量级差异，可为冻结法设计及施工提供合理性参考。

(3)基于巴霍尔金温度场解析解，提出了适用于渗流条件下冻结壁厚度的计算方法，同时结合模型试验结果，得到了常见冻结工况下冻结壁厚度关于渗流速度及冻结时间的计算公式，该公式可为实际渗流冻结工程中冻土帷幕扩展状态及冻结效果评估提供依据。