绝热边界作用下人工冻结温度场解析解

摘要

隧道地下工程采用冻结法施工时经常遇到地连墙等构筑物影响冻土帷幕发展的问题，给冻结方案设计和封水效果评估带来困难，危及后续开挖作业安全。为探明地连墙等绝热边界作用下冻结温度场的分布状态，掌握冻结壁局部受限时的发展规律，建立了半无限平面内包含直线绝热边界的温度场数学模型。通过热势函数叠加结合镜像法求解了单根、2根以及3根冻结管的温度场解析解，并采用数值模拟验证了其准确性和适用性。结果表明：解析解与稳态数值解吻合较好，与瞬态数值解的误差随冻结时间延长而减小，冻结第50天3种模型的误差分别为0.39℃、0.17℃和0.06℃，均能控制在0.5℃以内；等温线在绝热边界处与边界垂直，只存在平行于边界的热流而没有法向热流，绝热边界对于其与冻结管之间区域的温度降低更加有利；随着绝热边界与冻结管距离d的增加，对冷量传递的阻断作用降低，冻结模型逐渐向无限大平面模型转化，实际工程应根据冻结壁设计厚度，合理匹配管间距l和边界距离d的取值。

关键词

隧道工程; 人工冻结法; 温度场; 解析解; 数值模拟; 绝热边界

随着我国城市化进程进一步加快，城市建设规模不断扩大，有限的地面空间很难满足日益增长的生产生活需要，开发利用地下空间逐渐成为城市可持续发展的必经之路。特别是在北京、上海、深圳等超一线城市，地下空间的利用程度越来越高，大量建成的地下构筑物对新建工程提出了更高的要求。而人工地层冻结法自1862年首次在英国南威尔士的矿井中获得成功以来，因冻土帷幕具有良好的封水效果和力学强度，并且在保护周边环境方面存在优越性，目前已广泛应用于地下工程建设之中，如联络通道、盾构进出洞、地基加固和地下水污染防控等^［1-4］。作为一种特殊的岩土加固手段，冻结法中冻土帷幕的力学性质、冻土厚度、冻结封水路径等重要参数都依赖于温度分布。获得冻土温度场是准确判断冻结壁性状的基础，对于冻结法的科学设计与安全施工尤为重要。天然土体通过人工制冷技术循环低温冷媒形成冻土，必然包含土中水由液态产生相变转变为固态的过程，因此精确求解冻结温度场原则上必须要处理包含相变潜热释放的瞬态导热问题^［5］。长期以来由于数学上的困难，精确解一直局限于具有简单边界的单管冻结问题，目前仍未有突破。考虑到土体降温速度一般都很慢，采用稳态算法近似求解冻结过程中特定时刻的温度场不会造成较大误差^［6-7］。结合冻结管的不同布置形式，国内外学者基于稳态假设在温度场解析方面取得了系列成果：ТРУПАК^［8］和БАХОЛДИН^［9］提出了单管和双管冻结温度场计算方法；Sanger ^［10］求解了单排管解析解；Tobe等^［11］得到了直线上布置多根冻结管的温度场；周晓敏等^［12］和荣传新等^［13］通过模型试验研究了渗流条件下多管冻结的温度分布，文献［14-15］给出了冻结壁形状的解析表达；胡向东对双排冻结管的解答进行了补充和优化^［16］，结合水力学中群井干扰的结论求解了任意少量管^［17］、广义排管^［18］和圈管^［19］冻结温度场解析解。

然而，已有这些研究成果都是基于热传导发生在无限大平面内的假设进行的，没有考虑因导热边界限制而产生的冻土非对称发展问题。随着地下工程建设条件的越发严苛，新建工程采用冻结法施工时不可避免地紧邻地铁、地下车站、地下商场等构筑物，冻结管附近常常存在地下连续墙等影响传热，从而改变冻结温度场的分布^［20］。此外，施工中为了减少低温盐水的冷量损失，冻结孔附近的地连墙表面常会铺设一层或多层保温材料以保证墙外土体的冻结效果。对于这种工况，敷设隔热层的地连墙将抑制低温盐水的冷量耗散，可以近似地看成绝热边界，对冻结壁的发展具有不可忽视的影响，此时温度场计算若继续采用无限大平面假设会产生很大误差。

为了给这一实际工程问题提供理论指导，本文基于绝热边界在热传导中的作用特点建立冻结模型，结合势函数叠加法和镜像法等数学方法，尝试推导了直线绝热边界附近布置少量冻结管的稳态温度场解析表达，并利用数值手段验证了其准确性和适用性，以期为冻结设计与施工提供依据。

1 热传导控制方程及解的叠加原理

人工冻结过程中盐水冷量在土体中的传递方式主要以传导传热进行，根据傅里叶定律，有

$$q_x=-k\frac{\partial T}{\partial x}$$

（1）

式中：q_x为热流密度，x为横坐标；k为导热系数，T为温度。定义热势Φ =kT，则式（1）变为

$$q_x=$$$$-\frac{\partial \Phi }{\partial x}$$

（2）

另外，冻结管为细长构件，轴向尺寸远大于横截面尺寸，因此可将3维空间导热问题简化为二维平面问题。结合能量守恒定律，平面内热传导控制方程为

$$\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}-Q+\rho c\frac{\partial T}{\partial t}=\mathrm{0}$$

（3）

式中：q_y为热流密度；Q为从外界吸收或放出的热量；ρ为密度；c为比热。假设土体为各向同性材料，计算区内不考虑系统与外界的热量交换，在稳态条件下将式（2）代入式（3），可进一步在极坐标系中化简控制方程，为

$$\frac{\mathrm{d}}{\mathrm{d}r}(r\frac{\mathrm{d}\Phi }{\mathrm{d}r})=\mathrm{0}$$

（4）

式（4）实际上是二维平面内的Laplace方程，很容易证明，若Φ₁和Φ₂是该方程的解，那么Φ=Φ₁+Φ₂也一定是方程的解，即稳态热传导控制方程满足解的叠加原理。

2 无限大平面多管冻结的热势分布

在式（4）控制方程前提下，为研究含绝热边界的半无限平面问题，首先研究无限大平面的冻结模型。以最简单的单管冻结模型为例，如图1所示。因其温度分布在圆周上具有对称性，故可简化为沿半径方向的一维冻结问题。直接对式（4）积分，得

图1  无限大平面单管冻结模型

Fig. 1  Temperature model with one freezing pipe

$$\Phi =C_{\mathrm{1}}\mathrm{l}\mathrm{n}r+C_{\mathrm{2}}$$

（5）

再代入冻结管和冻土温度2个边界条件：① 冻结管表面，r = r₀，T=T_f，Φ = Φ_f。② 冻土边界处：r = ξ，T=T₀，Φ = Φ₀。解得

$$\Phi =\frac{{\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}}}{\mathrm{l}\mathrm{n}\frac{\xi }{r_{\mathrm{0}}}}\mathrm{l}\mathrm{n}\frac{\xi }{r}-{\Phi }_{\mathrm{0}}$$

式中：r₀为冻结管半径；ξ为冻土厚度；r为极径；T_f 和Φ_f分别为冻结管表面温度和热势；T₀和Φ₀分别为土体结冰温度和冻土边界的热势。由Φ =kT可得无限大平面单管冻结温度场分布表达式为

$$T=\frac{T_{\mathrm{f}}-T_{\mathrm{0}}}{\mathrm{l}\mathrm{n}\frac{\xi }{r_{\mathrm{0}}}}\mathrm{l}\mathrm{n}\frac{\xi }{r}+T_{\mathrm{0}}$$

（6）

式（6）与Трупак^［8］得到的经典单管解是完全相同的。当平面内布置n根冻结管时，由式（5）和上述解的叠加原理可以写出冻结区域内任一点的热势为

$$\Phi ={\displaystyle \sum _{i=\mathrm{1}}^n}{C}_{\mathrm{1}i}$$ lnri+C2

（7）

联立n根冻结管表面温度和冻土边界处温度条件，即可求解C₁₁、C₁₂、… 、C_1n 及C₂共计n+1个待定系数，从而根据热势定义获得温度场表达式。

3 存在绝热边界的温度场解析解

3.1　绝热边界附近布置单根冻结管

如图2所示，距离无限长绝热边界d处存在一根冻结管，以绝热边界为x轴、过冻结管中心并垂直于绝热边界的直线为y轴，建立平面直角坐标系。记冻结管中心到冻土边界的距离为ξ，稳态情况下冻结管表面温度维持恒定。

图2  单管绝热边界冻结模型示意

Fig. 2  Diagram of adiabatic model with one freezing pipe

绝热边界的存在将无限大平面模型变成包含一根冻结管的半无限平面模型，直接求解该模型比较困难。根据绝热边界无热流通过的特点，采用镜像法将单管模型映射成以绝热边界为对称轴的双管模型（图3）。因镜像管和原始冻结管的管壁温度相同，由对称性知，镜像模型在x轴上沿y方向的热流量恒为零，从而其上半平面内的温度场和原问题一致。

图3  双管镜像冻结模型

Fig. 3  Diagram of mirror model with two freezing pipes

设冻结管半径为r₀、管壁温度为T_f、土体结冰温度为T₀，由式（7）可知，映射后的模型中平面内任意一点P的热势为

$$\Phi =C_{\mathrm{11}}\mathrm{l}\mathrm{n}{r}_{\mathrm{1}}+C_{\mathrm{12}}\mathrm{l}\mathrm{n}{r}_{\mathrm{2}}+C_{\mathrm{2}}$$

（8）

分别代入冻土边界和冻结管边界条件。冻土边界为

$${\Phi }_{\mathrm{0}}=C_{\mathrm{11}}\mathrm{l}\mathrm{n}\xi +C_{\mathrm{12}}\mathrm{l}\mathrm{n}(\xi +\mathrm{2}d)+C_{\mathrm{2}}$$

（9）

冻结管边界为

$${\Phi }_{\mathrm{f}}=C_{\mathrm{11}}\mathrm{l}\mathrm{n}{r}_{\mathrm{0}}+C_{\mathrm{12}}\mathrm{l}\mathrm{n}(\mathrm{2}d-r_{\mathrm{0}})+C_{\mathrm{2}}$$

（10）

$${\Phi }_{\mathrm{f}}=C_{\mathrm{11}}\mathrm{l}\mathrm{n}(\mathrm{2}d+r_{\mathrm{0}})+$$$$C_{\mathrm{12}}\mathrm{l}\mathrm{n}$$$$r_{\mathrm{0}}$$$$+C_{\mathrm{2}}$$

（11）

由于实际工程中冻结管半径r₀与其间距2d相比非常小，对温度场的影响可以忽略，因此可将式（10）和（11）简化为

$${\Phi }_{\mathrm{f}}=C_{\mathrm{11}}\mathrm{l}\mathrm{n}{r}_{\mathrm{0}}+C_{\mathrm{12}}\mathrm{l}\mathrm{n}\mathrm{2}d+C_{\mathrm{2}}$$

（12）

$${\Phi }_{\mathrm{f}}=C_{\mathrm{11}}\mathrm{l}\mathrm{n}\mathrm{2}d+$$$$C_{\mathrm{12}}\mathrm{l}\mathrm{n}$$$$r_{\mathrm{0}}$$$$+C_{\mathrm{2}}$$

（13）

联立式（9）、（12）、（13），可以解得

$$C_{\mathrm{11}}=C_{\mathrm{12}}=\frac{{\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}}}{\mathrm{l}\mathrm{n}\frac{\mathrm{2}d{r}_{\mathrm{0}}}{\xi (\xi +\mathrm{2}d)}}$$

$$C_{\mathrm{2}}={\Phi }_{\mathrm{0}}-\frac{{\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}}}{\mathrm{l}\mathrm{n}\frac{\mathrm{2}d{r}_{\mathrm{0}}}{\xi (\xi +\mathrm{2}d)}}\mathrm{l}\mathrm{n}\xi (\xi +\mathrm{2}d)$$

代入式（8），可得热势分布函数表达式为

$$\Phi ={\Phi }_{\mathrm{0}}+\frac{\mathrm{l}\mathrm{n}\frac{r_{\mathrm{1}}{r}_{\mathrm{2}}}{\xi (\xi +\mathrm{2}d)}}{\mathrm{l}\mathrm{n}\frac{\mathrm{2}d{r}_{\mathrm{0}}}{\xi (\xi +\mathrm{2}d)}}({\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}})$$

式中：$r_{\mathrm{1}}=\sqrt[]{x^{\mathrm{2}}+{(y-d)}^{\mathrm{2}}}$；$r_{\mathrm{2}}=\sqrt[]{x^{\mathrm{2}}+{(y+d)}^{\mathrm{2}}}$。由Φ =kT，可得温度分布函数，为

$$T=T_{\mathrm{0}}+\frac{\mathrm{l}\mathrm{n}\frac{r_{\mathrm{1}}{r}_{\mathrm{2}}}{\xi (\xi +\mathrm{2}d)}}{\mathrm{l}\mathrm{n}\frac{\mathrm{2}d{r}_{\mathrm{0}}}{\xi (\xi +\mathrm{2}d)}}(T_{\mathrm{f}}-T_{\mathrm{0}})$$

（14）

由式（14）很容易验证

$$\left.\frac{\partial T(x,y)}{\partial x}\right|$$ y=0 =0

即满足绝热边界处y方向上无热流的条件。在式（14）中限制y≥0，即得所求绝热边界附近存在单根冻结管的稳态温度场解析解。

另外可以发现，原始冻结管和镜像管在无限大平面温度场中的权重相同，即C_1i=C_1j。实际上，若设冻结管圆周上的总热流量为q_c，则由式（2）可知

$$q_c=\mathrm{2}\mathrm{\pi }rq=-\mathrm{2}\mathrm{\pi }r\frac{\partial \Phi }{\partial r}$$

两边积分可得

$$\Phi =-\frac{q_{\mathrm{c}}}{\mathrm{2}\mathrm{\pi }}\mathrm{l}\mathrm{n}r+C$$

与式（5）对比可以发现，冻结管权重系数C₁实际上反映的是冻结管内的热流量情况-q_c/（2π）。对于映射后的无限大平面多管冻结问题，在管壁温度相同的条件下，若第i根和第j根冻结管坐标也满足对称条件，则其在群管温度场中的权重相同，即有C_1i=C_1j，这一结论将为后续更复杂问题的求解带来便利。

3.2　绝热边界附近布置2根冻结管

以冻结管平行绝热边界作为算例，假设距离d处布置有2根冻结管，间距为2l，如图4建立平面直角坐标系，记冻结管到冻土边界的距离为ξ。

图4  双管绝热边界冻结模型示意

Fig. 4  Diagram of adiabatic model with two freezing pipes

参考单管问题求解方法，采用镜像法将半无限平面内的双管冻结模型映射为无限大平面内的四管冻结模型，如图5所示。

图5  四管镜像模型示意

Fig. 5  Diagram of mirror model with four freezing pipes

因映射后的4根管关于x轴和y轴两两对称，故各冻结管在温度场分布中的权重相同，由3.1节可知，热势通解中C₁₁、C₁₂、C₁₃、C₁₄应相等，记为C₁，则图5中任意一点热势可表示为

$$\Phi =C_{\mathrm{1}}(\mathrm{l}\mathrm{n}{r}_{\mathrm{1}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{2}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{3}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{4}})+C_{\mathrm{2}}$$

（15）

代入冻土边界条件Φ（l， ξ+d）= Φ₀和冻结管边界条件Φ（-l， -r₀+d）= Φ_f，并联立解得

$$C_{\mathrm{1}}=\frac{{\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}}}{\alpha }$$

C₂=Φ_f-$\frac{{\Phi }_{\mathrm{f}}-{\Phi }_{\mathrm{0}}}{\alpha }$ln(2d·r₀·2l·$\sqrt[]{\mathrm{4}{d}^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}$ )

代入（15），同时考虑热势定义可得

T=T0+$$\frac{\beta }{\alpha }$$ （Tf-T0)

(16)

其中

α=ln $\frac{\mathrm{2}d·r_{\mathrm{0}}·\mathrm{2}l·\sqrt[]{\mathrm{4}{d}^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}}{\xi (\xi +\mathrm{2}d)·\sqrt[]{{\xi }^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}·\sqrt[]{{(\xi +\mathrm{2}d)}^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}}$

$\beta$=ln $\frac{r_{\mathrm{1}}·r_{\mathrm{2}}·r_{\mathrm{3}}·r_{\mathrm{4}}}{\xi (\xi +\mathrm{2}d)·\sqrt[]{{\xi }^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}·\sqrt[]{{(\xi +\mathrm{2}d)}^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}}$

式中：r₁、r₂、r₃、r₄分别表示待求温度点（x， y）到4根冻结管边界点的距离。容易证明式（16）也满足x轴上沿y方向无热流的绝热边界条件，限制y≥0即得图4模型的稳态温度场解析解。

3.3　绝热边界附近布置3根冻结管

如图6所示，假设距绝热边界d处分布3根冻结管，冻结管间距为l，记中间冻结管到冻土帷幕边界距离为ξ。同样采用镜像法将包含绝热边界的三管问题映射成无限大平面的六管问题，如图7所示。

图6  三管绝热冻结模型示意

Fig. 6  Diagram of adiabatic model with three freezing pipes

图7  六管镜像冻结模型

Fig. 7  Diagram of mirror model with six freezing pipes

由6根冻结管的几何关系可知，位于中间列的冻结管（管2、5）和位于两侧的冻结管（管1、3、4、6）热流量不同，因此它们在温度场分布中的权重也不同，设两侧冻结管的权重为C₁₁、中间冻结管权重为C₁₂，则图7中任意点热势的表达式为

$$\begin{array}{l}\Phi =C_{\mathrm{11}}(\mathrm{l}\mathrm{n}{r}_{\mathrm{1}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{3}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{4}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{6}})+\\ \hspace{1em}\hspace{1em}{C}_{\mathrm{12}}(\mathrm{l}\mathrm{n}{r}_{\mathrm{2}}+\mathrm{l}\mathrm{n}{r}_{\mathrm{5}})+C_{\mathrm{2}}\end{array}$$

（17）

代入冻土边界条件Φ（0， ξ+d） = Φ₀、冻结管1边界Φ（-l， -r₀+d） = Φ_f和冻结管2边界Φ（0， -r₀+d）= Φ_f，联立并考虑Φ =kT，解得

$$\begin{array}{l}T=T_{\mathrm{0}}+\frac{({\alpha }_{\mathrm{22}}-{\alpha }_{\mathrm{12}})(T_{\mathrm{f}}-T_{\mathrm{0}})}{{\alpha }_{\mathrm{11}}{\alpha }_{\mathrm{22}}-{\alpha }_{\mathrm{12}}{\alpha }_{\mathrm{21}}}{\beta }_{\mathrm{1}}+\\ \frac{({\alpha }_{\mathrm{11}}-{\alpha }_{\mathrm{21}})(T_{\mathrm{f}}-T_{\mathrm{0}})}{{\alpha }_{\mathrm{11}}{\alpha }_{\mathrm{22}}-{\alpha }_{\mathrm{12}}{\alpha }_{\mathrm{21}}}{\beta }_{\mathrm{2}}\end{array}$$

（18）

其中

$${\alpha }_{\mathrm{11}}=\mathrm{l}\mathrm{n}\frac{\mathrm{2}d\cdot r_{\mathrm{0}}\cdot \mathrm{2}l\cdot \sqrt[]{\mathrm{4}{d}^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}}}{({\xi }^{\mathrm{2}}+l^{\mathrm{2}})\cdot \left[{(\xi +\mathrm{2}d)}^{\mathrm{2}}+l^{\mathrm{2}}\right]}$$

$${\alpha }_{\mathrm{12}}=\mathrm{l}\mathrm{n}\frac{l\cdot \sqrt[]{\mathrm{4}{d}^{\mathrm{2}}+l^{\mathrm{2}}}}{\xi (\xi +\mathrm{2}d)}$$

$${\alpha }_{\mathrm{21}}=\mathrm{l}\mathrm{n}\frac{l^{\mathrm{2}}\cdot (\mathrm{4}{d}^{\mathrm{2}}+l^{\mathrm{2}})}{({\xi }^{\mathrm{2}}+l^{\mathrm{2}})\cdot \left[{(\xi +\mathrm{2}d)}^{\mathrm{2}}+l^{\mathrm{2}}\right]}$$

$${\alpha }_{\mathrm{22}}=\mathrm{l}\mathrm{n}\frac{\mathrm{2}d\cdot r_{\mathrm{0}}}{\xi (\xi +\mathrm{2}d)}$$

$${\beta }_{\mathrm{1}}=\mathrm{l}\mathrm{n}\frac{r_{\mathrm{1}}{r}_{\mathrm{3}}{r}_{\mathrm{4}}{r}_{\mathrm{6}}}{({\xi }^{\mathrm{2}}+l^{\mathrm{2}})\cdot \left[{(\xi +\mathrm{2}d)}^{\mathrm{2}}+l^{\mathrm{2}}\right]}$$

$${\beta }_{\mathrm{2}}=\mathrm{l}\mathrm{n}\frac{r_{\mathrm{2}}\cdot r_{\mathrm{5}}}{\xi (\xi +\mathrm{2}d)}$$

式中：r₁、r₂、r₃、r₄、r₅、r₆表示平面内任一点（x，y）到6根冻结管边界点的距离。同样，当y=0时很容易验证式（18）对y的偏导数为零，也即绝热边界上沿y方向无热流通过，满足绝热条件。因此，在式（18）中令y>0即得图6所示3根管模型的稳态温度场解析解。

4 解析解的准确性和适用性

4.1　准确性的稳态数值验证

在解析推导过程中，为方便计算采取了一定的简化处理，如将冻结管壁的圆周温度近似用管壁上一点进行代替，这可能会使温度场产生误差。为检验解析解的准确性以及求解方法的正确性，采用COMSOL Multiphysics数值软件进行热力学稳态模拟对比。针对单管、双管和三管模型分别计算，共计4组，其中单管模型中考虑了d＜ξ和d≥ξ 2种不同冻结情况。参考实际冻结工程，各组参数如表1所示，其中d为边界距离，ξ为冻土厚度，l为管间距。冻结管表面温度T_f=-30℃，土体结冰温度T₀=0℃，冻结管半径r₀=0.054m。

表1  稳态数值计算模型参数

Tab. 1  Calculation parameters adopted in steady-state numerical simulation

冻结模型	d/m	ξ/m	l/m
单管	0.5	1.5
	2.0	1.5
双管	0.5	1.5	0.5
三管	0.5	1.5	0.8

由于绝热边界的存在，冻土帷幕形状是事先未知的，所以在数值模型中不能预设出一条特定的冻土边界。采用的方法是在离冻结管较远处设置一圆形恒温边界，通过调节恒温边界上的温度荷载T_b使冻土边界条件点达到结冰控制温度T₀。以单管数值模型为例（图8），通过调节施加在半径20m的圆形恒温边界上的温度荷载T_b来满足冻土边界条件，使坐标（0， ξ）处的温度为T₀；绝热边界通过设置热绝缘实现（即图中半圆的直径部分）。因数值计算范围为靠近冻结管附近很小的区域，在采用平面自由三角形单元进行网格划分时，对冻结管处的网格进行加密处理。根据稳态导热原理可知，这种情况下的温度场分布与解析解的冻结模型是十分接近的。

图8  数值计算模型示意

Fig. 8  Schematic diagram of numerical calculation model

再将表1冻结参数分别代入单管、双管和三管模型的解析表达式（14）、（16）和（18），计算得到各模型的温度场解析结果。选择经过冻结管中心并平行于y轴的主面、经过相邻冻结管连线中点并平行与y轴的界面以及经过冻结管中心平行于x轴的轴面3个特征面进行解析和稳态数值结果的对比，4组参数下的对比曲线如图9所示，其中离散点代表数值结果，连续线代表解析结果。

图9  解析解和数值模拟计算结果对比

Fig. 9  Comparison of analytical results and numerical simulation

从图9可以看出，无论是单管模型、双管模型还是三管模型，解析解与稳态数值解的温度曲线都吻合得很好，特征面上温度分布的最大误差在0.01℃以内。这说明在稳态条件下本文所获得的解析解具有非常高的精确性，微小的误差可能是由于网格划分的精细程度不足造成。由于冻结管截面尺寸相对较小，采用冻结管表面一点的温度作为管壁边界条件不会造成较大误差。同时，对于冻结管热流量不同的多管非对称问题，采用差异化的温度权重系数进行求解是可行的。

4.2　适用性的瞬态数值验证

在人工地层冻结法中，土的冻结过程是带有土中水相变、包含冻结锋面移动的瞬态传热问题，虽然稳态解析解应用于瞬态冻结过程的计算和预测已有较多文献论证，但大多都是针对无限大平面进行的，对于绝热边界半无限平面采用镜像法求解的结果依然缺乏依据。本文采用瞬态数值模拟的方法对稳态解析解的适用性进行验证，土体的材料参数参照《长江口土层人工冻土物理力学性能试验报告》进行选取，所用到的土体热物理性质参数如表2所示。

表2  瞬态模拟土体热物理性质参数

Tab. 1  Thermal properties of soil adopted in transient simulation

土壤状态	温度/℃	导热系数/（W·（m·K） -1）	密度/（kg·m-3）	比热/（J·（kg·K）-1）	焓值/（kJ·m-3）
冻土	-28	1.74	1 820	1 560	0
冻土	-10	1.74	1 820	1 560	55 714.4
冻土	-4	1.74	1 820	1 560	115 302.2
冰点	-2.1	1.74	1 820	1 560	199 526.7
融土	18	1.42	1 870	1 650	247 954.0

瞬态模拟的冻结管位置和冻结管半径等参数与表1一致，在冻结过程早、中、后期分别选取第10天、第30天和第50天的温度数据与解析解进行对比。限于篇幅，这里仅选取不同时期单管模型（d≥ξ）的主面、双管模型的界面和三管模型的轴面温度作为代表，对比结果如图10所示。

图10  解析解与瞬态模拟结果的对比（第10天、30天和50天）

Fig. 10  Comparison of steady-state analytical solutions and transient simulation results (10^th, 30^th, and, 50^th days)

由图可以看出，在冻结早期第10天，解析结果与瞬态数值结果有一定差异，绝热边界处最大差值为单管模型的2.4℃；冻土边界（冰点）处温差很小，误差约为0.02℃；模型边界为正温区域，误差约为3.5℃。这是因为在冻结初期，低温盐水与高温土体之间温差大，冻结管和土体之间剧烈的热交换使得地层处于快速降温状态，此时的温度场分布尚未达到稳定状态，从而与稳态解析结果之间存在较大差异。随着冻结时间的延长，盐水与周边土体之间的温差逐渐减小，冻结管与地层间的热交换速度降低，冻结温度场逐渐向稳态转化，解析解和瞬态数值解之间的差异也逐渐减小。至冻结中期第30天，绝热边界处最大误差缩小至0.76℃，冻土边界处最大误差0.006℃，模型边界正温处最大误差0.90℃，解析精度足以满足实际工程的应用需求。至冻结第50天，解析和瞬态曲线几乎完全重合，单管、双管和三管冻结模型最大误差分别为0.39℃、0.17℃和0.06℃，均在0.5℃内。说明解析解能够很好适用于瞬态传热过程，特别是在冻结中后期，温度场表达式在整个计算域上都有足够高精度供实际工程采用。

4.3　基于解析解的温度分布云图

在验证解析解的准确性和适用性之后，为更加直观显示绝热边界作用下冻结温度场的分布情况，采用MatLab软件对单管、双管和三管冻结模型的温度场解析函数进行云图绘制，各模型的二维及三维温度场分布如图11所示。

图11  基于解析解的温度场分布云图

Fig. 11  Cloud maps of temperature field based on analytical formulas

由图可见，温度场等值线随着与冻结管距离的增加而逐渐稀疏，说明距离冻结管越近，温度梯度越大；而远离冻结管的温度梯度逐渐减小，这符合冻结温度分布的一般规律。其次，由傅里叶定律可知热流量总是从高温区域流向低温区域，热流方向垂直于等温线。二维云图中绝热边界处等温线总是与边界正交，因此热流量必然平行于绝热边界，而无垂直，解析云图完全符合绝热边界作用下的热流特征。另外可以发现，相较于冻结管外侧浅色区域，冻结管与绝热边界之间被更加低温的深色区域覆盖，说明相比于无限大平面空间，绝热边界的存在使得其与冻结管之间的温度降低更加有利。

4.4　讨论

（1）关于冻结管和绝热边界相对位置。由镜像求解方法可知，随着d的增大，绝热边界的影响逐渐减小，冻土帷幕的形状也发生变化。以单管模型为例，冻土帷幕从d较小时的半椭圆形状，到d较大时冻土轮廓在绝热边界上产生凹陷，再到当d继续增大至一定值时，冻土边界必将与绝热边界分离。冻土帷幕在绝热边界上形成的冻结壁胶结长度与冻结管相对绝热边界距离具有密切的关系，冻土帷幕在边界上的最大分布、开始凹陷和开始脱离3种临界极限状态均可以通过温度场解析表达式进行计算。

（2）关于土体结冰温度。实际工程中因土体存在不同程度的含盐量，结冰温度往往低于0℃，有时甚至达到-2℃以下，如果依然按照0℃的结冰温度来进行冻结壁设计必然会造成工程隐患。也正是考虑到这一点，在本文解析推导部分，针对性地将结冰温度设为变量T₀，而并没有限定其为0℃，因此解析解针对不同结冰温度的地层均是适用的。在稳态数值模拟部分，选择具有代表性的纯水结冰温度0℃作为算例，其目的仅是验证解析解的准确性和精确性，并不表明实际土体的冰点。因为如果解析推导无误，那么对于任何T₀值的选取，解析结果和数值结果都应当是相吻合的。而且，瞬态模拟以上海地区冰点为-2.1℃的土体作为算例也验证了对不同结冰温度的适用性。

（3）关于稳态解析解的应用。在冻结设计中，一般根据冻结壁所受外荷载（包括水土荷载、地面超载等），由力学计算先拟定出冻结壁的厚度值，再根据该设计厚度来布置、优化、调整冻结管的布置。本文解析解显式给出了温度分布与冻土厚度ξ以及其他冻结参数之间的定量关系，对于设计人员在拟定冻结壁厚度后合理布置冻结管具有参考价值。在施工过程中，通过现场监测，由任一测点温度结合解析解就可以精确反算冻土帷幕的厚度，这相比于目前采用平均降温速率进行估算的方法具有更高的准确性。因此，获得温度场解析解对于研究绝热边界作用下冻结管的合理布置参数、判断是否形成可靠的冻土封水路径具有一定的工程指导意义。

5 结论

以城市地下工程冻结法施工附近存在地连墙，特别是采用冻结管进行外侧局部冻结、内侧大面积保温为工程背景，进行冻结温度场的解析研究，获得的主要结论如下：

（1）通过对实际工程问题的抽象，建立了绝热边界作用下单管、双管和三管冻结的热传导数学模型，并在稳态假设下结合镜像法求解了温度场分布。冻结管在温度场中的权重系数由其热流量决定，解析解在绝热边界处满足无垂直热流通过的边界条件。

（2）利用和解析模型相同的冻结参数进行稳态导热数值计算，考虑了冻结管和绝热边界相对位置关系，通过主面、界面和轴面3个特征面的温度曲线对比，验证了稳态解析表达式的精确性和势函数镜像求解方法的正确性。

（3）结合上海长江口土体热物理参数进行瞬态导热数值计算以反映实际冻结过程。根据不同冻结时期的对比结果发现，随着冻结过程的进行，解析解从冻结早期在冰点处精确逐渐发展到整个温度计算域上都较为吻合，冻结50天3种模型的最大误差均能控制在0.5℃以内，解析表达式具有较高的精度。

应指出的是，因实际工程中地连墙保温效果常不能达到理想绝热状态，地连墙和内部高温气流间不可避免存在热量交换，因此，更加精确的温度场计算需将理论模型中的绝热边界条件转换为第3类传热边界条件进行求解。本文基于绝热边界获得的温度场解析解是相对理想情况下的结果，在实际应用中需根据边界条件变化和实测温度数据合理调整。

作者贡献声明

洪泽群：数据处理及数值分析。

付硕任：图表呈现，论文返修。

胡向东：指导理论分析，论文初稿撰写。

王宝生：解析解的适用性验证。