人工冻结是在地铁与隧道等地下联络通道施工中常见的施工方法^[1]。人工冻结法主要是通过人为制冷装置暂时性地改变岩土体的化学物理性质,形成有足够强度与密闭性的冻结帷幕,待施工结束后,撤出冷冻装置,恢复冻结岩土体^[2-3]。冻结过程中,地下水的流动性会对冻结壁形成一定的干扰甚至阻止其造成稳定区间,对冻结帷幕安全性造成一定的威胁,因此针对地下水的流动性对冻结区域的影响规律进行研究,以有效防止冻结帷幕受到破坏,具有重要的意义^[4]。通过对地下水流动条件下地下联络通道的冻土帷幕的温度场与渗流场的数值模拟,可以分析冻结区域内稳定性的变化规律^[5]。

针对冻结法施工中出现的复杂问题,国内外研究者开展了大量的研究并获得丰硕的成果。王朝晖等^[6]通过渗流状态下的土层液氮冻结室内试验,发现水流对液氮冻结有重要影响。刘泉声等^[7]通过公式推导出冻结过程中岩土体的温度场、渗流场与应力场数值耦合应用控制方程,并与有限元数值模拟的冻结效果做比较分析,发现冻结过程对岩体渗透系数有较为明显的影响,并获得在低温情况下冻胀效应的温度场、渗流场及应力场的分布规律。刘子昕^[8]通过运用有限元软件探讨氮冻结方式下的二维温度场变化,发现冻土的单轴抗压强度和冻结温度之间存在明显的线性相关,进而以此探究土层中液氮冻结形成的温度场变化规律。Bronfenbrenner等^[9]通过建立数值计算模型和室内试验的方法,研究了冻土在特定温度下的冻结壁厚度及强度的变化规律,从而验证了试验模型Stefan准则的一般形式。

综观现有研究,鲜见涉及地下水渗流情况下耦合温度场变化的冻结规律,故本研究在前人研究的基础上,以福建省福州市地铁5号线某区间联络通道在地下水渗流作用下的冻结过程为例,建立温度场与渗流场数值模拟耦合模型,结合实际工程情况,分析不同渗流状态下的冻结帷幕演变规律,归纳出地下水渗流对地下联络通道冻结区域的影响规律,以期为地下联络通道建设过程中的安全与经济分析提供参考。

地下隧道联络通道人工冻结施工存在水冰相变问题、内部热源问题、斯特凡问题,是涉及移动边界条件的复杂工程。在冻结区域中,温度的变化将导致渗流场中水流的物理性质参数变化,尤其是渗透系数的变化; 而渗流场内水流的热传导效应、流体的压缩性能与黏性损失又会反过来影响温度场的物理性质参数^[10-12]。冻结温度场通过热传导把能量传递到土体中,而地下渗流场又会改变能量传递的方向及位置,进而改变其温度场的分布性质,从而得出渗流场与温度场之间的数值变化耦合关系^[13-15]。温度场与渗流场耦合关系如图1所示。

图1 温度场与渗流场耦合关系
Fig.1 Coupling relation diagram of temperature field and seepage field

基于能量守恒原理与孔隙介质对流换热原理,运用显热容法研究相变过程,推导出温度场控制方程^[16-17]

式(1)中:Q_t为热源,W/m³; C_eq为等效体积显热容,J/(kg·K); T为温度,K; t为时间,s; ρ _l为流动水的密度,kg/m³; L为单位质量的冰水相变释放的潜热值,取值为334 J/kg; θ _l为液态水的含量,%; C_l为水的比热容,J/(kg·K); u为地下水渗流速度矢量,m/s; λ_eq为等效的导热系数,W/(m·K); Δ为哈密顿算子。

结合达西定律^[18],可推导出渗流场控制方程

Δ(ρ u)=Q_m (2)

u=-K/ηΔp (3)

式(2)～(3)中:ρ为水的密度,kg/m³; Q_m为源项,W/m³; K为冻结过程中土体的渗透率,m/d; η为水的黏度,kg/(m·s); p为孔隙水的压强,Pa。

考虑冻结中存在水冰相变因素,可修正土层的渗透系数k(T)^[18]为

k(T)=k_u×H(T)+k_f×[1-H(T)] (4)

式(4)中:k_u为未冻结时土层的渗透率,m/d; H(T)为二维阶跃平滑函数H关于温度T的海维赛德函数; k_f为冻结后土层的渗透率,取值为10^-25 m/d,即表示当土层下降到冻结温度时不透水。

以福州市轨道交通5号线某区间地下联络通道为例,依据详细的地质勘查报告可知,联络通道的地层主要为淤泥质土层(夹砂层)和中砂层。上部主要为淤泥质土层(夹砂层),下部主要为中砂层。孔隙承压水层主要分布于隧道洞身范围和隧道下方,含水层岩性主要为细砂、中砂,属中等强透水层。本工程使用冻结法加固土层,并依据实际情况采用冻结管87根。

为了更好地研究地下水渗流对冻结过程中地下联络通道的影响,对冻结过程中的影响因素进行一定的简化和假设,以达到降低收敛难度的目的。现做出以下假设:

1)假设冻结区域的岩土层为淤泥质土,土体在冻结过程中均质且连续;

2)地下水渗流过程中仅考虑渗流场影响冻结温度,且渗流场流动规律符合达西定律;

3)冻结管周围的影响区域至恒温边界,忽略冻结管内部盐水的流动损失,冻结管管道外侧温度等同于盐水温度,并直接赋予冻结管壁温度参数。

本研究以地下水渗流作用下的福州地铁5号线某区间联络通道冻结过程为例,运用COMSOL Multiphysics软件建立基于温度和渗流控制方程的数值模拟计算模型,结合实际岩土层的各项物理参数,分析地下水渗流作用对冻结区域产生的温度渗流影响。数值模拟布置87根冻结管,设计冻结联络通道三维数值模型规格为40 m×20 m×30 m。数值模拟计算模型及其网格分布如图2所示。

图2 数值模拟计算模型及其网格分布
Fig.2 Numerical simulation calculation model and grid distribution

冻结区域土层所在层位起始温度约为293.15 K,将在联络通道两侧的管壁设置为热传导边界,其他边界为热绝缘边界; 将联络通道两侧的管壁设置为透水边界,其他边界设置为不透水边界。

本模型网络划分采用自由四面体和扫掠结合的网格方法,划分为161 165个四面体单元、21 216个边界单元、8 760个边单元和682个顶点单元,平均单元质量为0.256 8。

在模型的模拟过程中,主要的参数设置来源于前期现场勘察和试验,模型基本参数见表1。

表1 模型基本参数
Table 1 Model basic parameters

根据实际地下联络通道的工程施工情况,对冻结区域进行数值耦合模拟,此时地下水渗流方向顺着冻结管方向前进,模拟参数取值见表2。

盐水温度下降采用内插函数intl(t)表示,可得冻结管盐水降温图(图3)。

表2 模拟参数取值
Table 2 Trial parameter value

图3 冻结管盐水降温图
Fig.3 Salt water cooling diagram of freezing pipe

实际工程中地下水先到的一端为渗流上游,则另一端为渗流下游,将A-A截面(简称A截面)、B -B截面(简称B截面)及实际测温点位置作为分析对象,可得获得A、B截面测温点及其工况结构示意图(图4)。

图4 A、B截面测温点及其工况结构示意图
Fig.4 Schematic diagram of temperature measurement points and working conditions of section A and section B

将工程实测数据与数值模拟结果做对比分析,得到各位置实测温度和模拟温度的对比图(图5)。

图5 各位置实测温度和模拟温度的对比
Fig.5 Comparison chart of measured temperature and simulated temperature at various locations

由图5可知,距隧道中心线越近,温度变化越大; 距中心线越远,受冻结施工的影响越小,则温度变化越小。由于受环境、施工和人为因素等影响,数值模拟结果和现场实测数据结果存在偏差,但总的趋势还是比较接近,从而证明了本数值模型的可靠性。

根据工程数值模拟结果,得出冻结时间为5、15、27、40 d时受地下水影响下的冻结壁变化情况,由此得到实际工程参数下冻结壁随冻结时间的变化情况(图6)。

图6 实际工程参数下冻结壁随冻结时间的变化情况(单位:m)
Fig.6 Freezing wall variation with freezing time under actual engineering parameters(unit: m)

由图6(b)可知,当冻结时间小于15 d时,地下联络通道上下部分的冻结帷幕均未全部冻结交圈,随后从左右两端开始冻结壁逐渐交圈,左右侧壁顶部也慢慢交圈。依据图6整体变化情况,冻结壁随着冻结时间的增大而逐渐交圈,其厚度也逐渐增加,侧面隧道顶部位置冻结区域整体发展良好。

提取截面A、B中冻结时间为5、15、27、40 d的冻结壁温度分布图,获得截面A和截面B冻结壁温度分布(图7)。由图7可知,冻结时间大致分为2个阶段,第一阶段冻结时间为t≤27 d的时候,温度由正温度向负温度过渡阶段,冻结壁中大多数的自由水在土体中结冰并释放出大量潜热,进一步产生冻胀,其变化表现为急剧增长; 第二阶段冻结时间为27 d<t≤45 d,冻结壁初具形态,自由水大量减少,土体冻胀变化趋缓。在地下水流动过程中,由于渗流作用把大部分热量带入下游,位于截面B的冻结厚度比截面A的冻结厚度大,冻结影响范围较大,冻结情况符合实际。在地下水渗流冲刷作用下冻结过程中冻结壁的发育不对称,随着时间的增加,冻结壁厚度缓慢增加,在渗流水域上游部分(即截面A)比渗流水域下游部分(即截面B)冻结壁发育速度迟缓,如图7(b)与(f)所示。冻结过程中,由于截面A顶部和底部存在双排冻结管,因此提供的冷源比截面B中的多,又因为冻结过程中地下水的渗流作用,导致整体的冻结效果是渗流上游部分比渗流下游部分冻结速度慢。

选取实际工程中冻结时间为27 d的时间点,取沿着冻结管方向侧壁剖面的纵截面,得出静水状态下和实际地下水渗流状态下的模拟结果,进一步分析地下水渗流对联络通道的冻结区域影响的规律,由此获得第27 d时在静止和渗流状态下纵截面温度场对比图(图8)。

图7 截面A和截面B冻结壁温度分布
Fig.7 Freezing wall temperature distribution of section A and section B

图8 第27 d时在静止和渗流状态下纵截面温度场对比
Fig.8 Comparison of longitudinal section temperature field under static and seepage conditions on the 27th day

由图8可知,离冻结管距离越远,土体冻结的温度下降越慢。图8(a)中最外圈表示在区域内冻结温度相等,且最外圈的等温线到冻结管表面的距离也比较均匀,即表示在静水状态下冻结温度场在冻结管附近分布较为均匀。由图8(b)可知,最外圈的等温线范围蔓延到地下隧道渗流下游附近,说明地下水流动导致冻结区域沿着冻结管打孔方向延伸,且冻结壁发展较为缓慢,这是由于地下水渗流引导温度场冷源偏向下游转移。由图8可知,该位置截面第27 d时已经基本上形成较为良好的冻结帷幕,其基本冻结变化规律和截面A与截面B相吻合。

提取模型在第0、15、27、40 d实际工程参数下各个时期渗流场分布图(图9),图中箭头表示渗流的速度和方向,颜色的不同表示渗流的速度大小不同。

图9 实际工程参数下各个时期渗流场分布
Fig.9 Distribution map of seepage field in different periods under actual engineering parameters

由图9(a)与(e)可知,由于地下联络通道的渗流过程为三维渗流场且渗流不对称,导致地下水渗流速度也不对称。处于隧道顶部和底部的渗流速度最大为1.92 m/d(渗流速度等于地下水流速度与孔隙率的乘积),由此可以推测存在隧道初始渗流场的隧道渗流现象,尤其处于隧道顶部和底部位置渗流速度有变快迹象,此刻的流速约为初始流速的1.92倍; 其他位置渗流速度均较小,并且隧道中心线中点位置流速更低,甚至无流动。综上可得,若地下水渗流的初始速度为1 m/d,则在隧道影响附近的冻结区域地下水流速范围是0～1.92 m/d。另外,地下水渗流过程中由于存在绕流的现象,导致在隧道口处的渗流变化最大,因此2个隧道口和联络通道交接处的冻结厚度最为薄弱,需要在设计施工过程中采取加固措施。

由图9(b)与(f)可知,冻结时间为15 d时,就冻结整体而言渗流场尚未发生较大变化,顶部和底部也还未封闭交圈。地下水基本上顺着隧道顶部外部流过,从上游的外部穿过联络通道的内部,并从隧道内部渗流至下游外壁流出。从图9(b)与(f)分布图可以看出,此时渗流场还未发生明显的变化。

由图9(c)与(g)可知,冻结时间为27 d时,冻结帷幕基本上完成交圈,联络通道的内部渗流场发生较大程度的变化,冻结帷幕中地下水无法渗入,底部存在较薄的冻结壁,侧壁的地下水仅仅顺着冻结壁流动。

由图9(d)与(h)可知,冻结时间为40 d时,冻结帷幕已经形成封闭式完整体,地下水渗流已无法进入联络通道内部,其仅仅顺着外侧冻结帷幕流动。

综上所述,地下水渗流状态下温度场变化会引起渗流场变化,并且在冻结帷幕的变化过程中,渗流场也会形成自身的内部自由边界面,不断随冻结帷幕的发育情况而变化。地下水渗流场在联络通道冻结区域内会不断随冻结帷幕发生变化,且离冻结区域的距离越近受影响程度越大。

本研究从模拟地下水渗流对地下联络通道冻结影响入手,建立以能量守恒和达西定律为基础的温度场与渗流场的数值方程,借助COMSOL数值模拟软件完成温度场与渗流场之间相互作用的耦合模型,并结合实际工程数据,分析在地下水渗流作用下对联络通道冻结区域影响的变化规律,得出如下结论:

1)通过实际工程实测数据和数值模拟结果的数据对比,得出二者较为吻合,验证了本模型的可靠性;

2)在地下水渗流的影响下,分析地下联络通道冻结壁的演变过程,得出靠近渗流上游越近,冻结壁冻结越迟缓,靠近渗流下游越近,冻结壁冻结越快;

3)通过对渗流场演变分析得出,随着在地下水渗流状态中的冻结时间变化,可将冻结壁发育分为前后期; 冻结前期,地下水对冻结区域有较大的影响; 冻结后期,地下水渗流对其影响较小甚至无影响。

本研究通过探讨冻结区域地下水冻结壁厚度变化情况,可为类似工程提供参考,从而保证其建设过程中的安全性和经济性。