冻结法施工中围岩的冻胀应力、冻胀应变和冻胀潜势

Freezing method construction of surrounding rock freezing expansion stress freezing expansion strain and freezing expansion potential

摘要

为预测和控制冻结法施工中围岩的冻胀行为，依据弹塑性体在具有膨胀趋势时可蓄积变形能的理论，分析具有膨胀趋势受约束的冻土具有该类变形能，并在约束不足或解除时表现为冻胀应力和冻胀应变。为此，在弹性变形能理论基础上，提出冻胀潜势的概念及其表述。该概念以冻胀应力和冻胀应变函数为参数，表征围岩因冻结而蓄积在土体中的变形能。并在此基础上，对某地铁沿线的粉质黏土进行系列冻胀试验。经计算得到了不同温度时的冻胀应力、冻胀应变，当冻结温度为-10℃时，冻胀应力、冻胀应变和冻胀率分别为0.338 MPa、0.446%、3.13%，计算结果与工程实际需求基本吻合。冻胀潜势可用来描述岩土体的整体冻胀行为，也对丰富和发展冻土理论和冻结法技术提供新的思路和途径。

关键词

冻结法施工 / 冻胀应力 / 冻胀应变 / 冻胀潜势 / 冻胀率

Key words

Freezing method construction / frost heave stress / frost heave strain / frost heave potential / frost heave rate

0 引言

人工冻结法施工是岩土体加固的一种有效方式，多用于隧道联络通道施工、软弱土层开挖和高含水率地下土层渗透等工程中^[1-4]。但在冻结和解冻时，土体会不可避免地产生冻胀和融沉，进而会影响支护结构的稳定性、安全性和耐久性。

为此众多学者对预测和控制冻结法施工中围岩的冻胀变形问题进行了广泛研究。蔡海兵等^[5-6]通过理论与数值模拟先后提出了不同冻结阶段隧道冻胀位移预测方法与隧道水平冻结期地层位移的热力耦合数值分析方法。Yang等^[7]提出了一种针对地表竖向冻胀位移的水热力耦合预测模型。吴东军等^[8]通过建立土体冻融试验装置研究了温度、含水率对水分迁移的影响。Yue等^[9]通过采用有限元分析方法得到了人工地面冻结期间冻胀压力和隧道变形随冻结温度、时间的规律。崔宏环等^[10]研究了季节区粉质黏土路基的冻胀特性，开发了冻胀模型。董建华等^[11]研究了不同冻胀模式下隧道施工引起的地表竖向位移。

土体的冻胀力和冻胀变形与其含水率、外部荷载和土体密实度等因素密切相关^[12-14]。增大外部荷载会减小冻胀变形，但同时会增大冻胀力。同样，减小外部荷载能够降低冻胀力，但会诱发较大的冻胀变形，进而引发解冻后的过大融沉。可见，冻胀力与冻胀变形的研究均不能直接反映土体的冻融力学性能。因此，建立冻胀力与冻胀变形之间的关系，从能量角度评价土体的力学状态，进而构建工程冷源与土体冻结能量之间的关系，对冻结工程具有理论指导意义。

为此，本研究在弹塑性变形能理论和能量守恒理论基础上，建立了基于冻胀应力和冻胀变形的冻胀潜势概念。根据不同冻结温度条件下的冻胀力和变形量确定其冻胀潜势，可据此评价冻土的综合力学性能。在此基础上，给出了轴对称冷源状态下冻胀潜势与冻胀应力和冻胀应变之间的关系。

1 基于弹性理论的围岩冻胀分析

1.1 计算假定

将联络通道、周围冻结区和未冻结区视为2个轴对称弹性受力体，如图1所示。其中，隧道及联络通道分别视为位于无限大土体中的圆形孔和矩形孔。为建立理论化冻胀潜势模型，作如下假设。

1）隧道及周围土体为均匀各向同性材料。

2）隧道围岩为平面应变问题。

3）冻结区水分分布均匀。

4）只考虑冻胀效应对应力和变形的影响。

5）不考虑冻结管之间的相互作用。

若联络通道周围土体冻结时向衬砌和未冻结区传递的冻胀位移分别为Δh ₁和Δh ₂，σ ₁和σ ₂分别为Δh ₁和Δh ₂对应的冻胀力，如图2所示。图2中a、b和c分别为中心点到衬砌内壁、外壁和围岩内壁的距离。联络通道衬砌为仅受上覆土体重力和冻结区土体冻胀力作用的厚壁圆筒，属轴对称问题。

1.2 围岩冻胀应力计算分析

在冻胀力（σ ₁）作用下，可基于弹性理论得到联络通道衬砌径向位移（ u_{\mathrm{1}}），即

u_{\mathrm{1}} ( r ) = - \frac{b^{\mathrm{2}}{\sigma }_{\mathrm{1}}}{E_{\mathrm{1}}(b^{\mathrm{2}}-a^{\mathrm{2}})} \left[(\mathrm{1}-\mathrm{2}{v}_{\mathrm{1}})(\mathrm{1}+v_{\mathrm{1}})r+\frac{(\mathrm{1}+v_{\mathrm{1}})a^{\mathrm{2}}}{r}\right]。（1）

式中：E ₁和v ₁分别为联络通道衬砌的弹性模量和泊松比；r为任意点到隧道中心的距离。

在衬砌外径处，即当r=b时，其位移（ {\delta }_{\mathrm{1}}）为

{\delta }_{\mathrm{1}} = - \frac{b{\sigma }_{\mathrm{1}}}{E_{\mathrm{1}}(b^{\mathrm{2}}-a^{\mathrm{2}})} ( \mathrm{1} + v ) \left[(\mathrm{1}-\mathrm{2}v)b^{\mathrm{2}}+a^{\mathrm{2}}\right]。（2）

冻结土体可等效为σ ₁与σ ₂作用下的轴对称问题，其联络通道衬砌径向位移u ₂（r）为

u_{\mathrm{2}} ( r ) = \frac{(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})(\mathrm{1}+v_{\mathrm{2}})(b^{\mathrm{2}}{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}}{c}^{\mathrm{2}})r}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})} + \frac{(\mathrm{1}+v_{\mathrm{2}})({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}})b^{\mathrm{2}}{c}^{\mathrm{2}}}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})r}。（3）

式中，E ₂和v ₂分别为冻结土体弹性模量和泊松比。在冻结内壁处，即当r=b时其位移（ {\delta }_{\mathrm{11}}）为

{\delta }_{\mathrm{11}} = \frac{(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})(\mathrm{1}+v_{\mathrm{2}})(b^{\mathrm{2}}{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}}{c}^{\mathrm{2}})b}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})} + \frac{(\mathrm{1}+v_{\mathrm{2}})({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}})b{c}^{\mathrm{2}}}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})}。（4）

在外壁处，即当r=c时其位移（ {\delta }_{\mathrm{12}}）为

{\delta }_{\mathrm{12}} = \frac{(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})(\mathrm{1}+v_{\mathrm{2}})(b^{\mathrm{2}}{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}}{c}^{\mathrm{2}})c}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})} + \frac{(\mathrm{1}+v_{\mathrm{2}})({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}})b^{\mathrm{2}}c}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})}。（5）

未冻结区土体可等效为仅受内壁σ ₂作用的轴对称厚壁圆筒，得其内壁位移（ {\delta }_{\mathrm{2}}）为

{\delta }_{\mathrm{2}} = \frac{c}{E_{\mathrm{3}}} ( \mathrm{1} + v_{\mathrm{3}} ) {\sigma }_{\mathrm{2}} （6）

式中，E ₃和v ₃分别为未冻结土体弹性模量和泊松比。

在衬砌和冻结土体接触面及土体冻结区和未冻结区接触面处，应满足以下各式，即

- {\delta }_{\mathrm{1}} + {\delta }_{\mathrm{11}} = \mathrm{\Delta } h_{\mathrm{1}}；

{\delta }_{\mathrm{2}} - {\delta }_{\mathrm{12}} = \mathrm{\Delta } h_{\mathrm{2}}；

\mathrm{\Delta } h_{\mathrm{1}} = \frac{\mathrm{\Delta }{V}_{\mathrm{1}}}{\mathrm{2}\mathrm{\pi }b} = \frac{\eta \mathrm{\pi }\left[{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}-b^{\mathrm{2}}\right]}{\mathrm{2}\mathrm{\pi }b} = \frac{\eta \left[{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}-b^{\mathrm{2}}\right]}{\mathrm{2}b}；

\mathrm{\Delta } h_{\mathrm{2}} = \frac{\mathrm{\Delta }{V}_{\mathrm{2}}}{\mathrm{2}\mathrm{\pi }c} = \frac{\eta \mathrm{\pi }\left[c^{\mathrm{2}}-{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}\right]}{\mathrm{2}\mathrm{\pi }c} = \frac{\eta \left[c^{\mathrm{2}}-{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}\right]}{\mathrm{2}c}。（7）

式中，η为冻结土体冻胀率；ΔV ₁和ΔV ₂分别为冻结土体与衬砌和非冻结土体接触面的膨胀量。

因此，冻结土体的冻胀应力（ \sigma）为

\sigma = \frac{A}{B}。（8）

式中：A为不同因素导致的冻胀变形总和；B为冻胀过程中的弹性约束。

A = A_{\mathrm{1}} + A_{\mathrm{2}} + A_{\mathrm{3}}， A_{\mathrm{1}}和 A_{\mathrm{2}}分别为不同材料参数导致的变形量； A_{\mathrm{3}}为某种体积变形的影响，其计算公式分别为

A_{\mathrm{1}} = \frac{\eta bc(\mathrm{1}+v_{\mathrm{2}})(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})}{\mathrm{2}{E}_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})k} \left[c^{\mathrm{2}}-{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}\right]；

A_{\mathrm{2}} = \frac{\eta bc(\mathrm{1}+v_{\mathrm{2}})}{\mathrm{2}{E}_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})k} \left[c^{\mathrm{2}}-{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}\right]；

A_{\mathrm{3}} = \frac{\eta }{\mathrm{2}b} \left[{(\frac{b+c}{\mathrm{2}})}^{\mathrm{2}}-b^{\mathrm{2}}\right]。（9）

B = B_{\mathrm{1}} + B_{\mathrm{2}} + B_{\mathrm{3}} + B_{\mathrm{4}}， B_{\mathrm{1}}为冻胀土体自身的弹性特性引起的约束； B_{\mathrm{2}}为不同部位的变形差异而产生的额外约束； B_{\mathrm{3}}为不同泊松效应引起的横向约束； B_{\mathrm{4}}为土体与外界的接触应力差异产生的约束，其计算公式分别为

B_{\mathrm{1}} = \frac{(\mathrm{1}+v_{\mathrm{2}})b}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})} \cdot \left[(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})b^{\mathrm{2}}+c^{\mathrm{2}}\right]；

B_{\mathrm{2}} = - \frac{(\mathrm{1}+v_{\mathrm{2}}{)}^{\mathrm{2}}{b}^{\mathrm{3}}{c}^{\mathrm{3}}}{{E_{\mathrm{2}}}^{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}}{)}^{\mathrm{2}}k} \cdot \left[(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}}{)}^{\mathrm{2}}+\mathrm{1}\right]；

B_{\mathrm{3}} = - \frac{\mathrm{2}(\mathrm{1}+v_{\mathrm{2}}{)}^{\mathrm{2}}(\mathrm{1}-\mathrm{2}{v}_{\mathrm{2}})b^{\mathrm{3}}{c}^{\mathrm{3}}}{{E_{\mathrm{2}}}^{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}}{)}^{\mathrm{2}}k}；

B_{\mathrm{4}} = \frac{(\mathrm{1}+v_{\mathrm{1}})b\left[(\mathrm{1}-\mathrm{2}{v}_{\mathrm{1}})b^{\mathrm{2}}+a^{\mathrm{2}}\right]}{E_{\mathrm{1}}(b^{\mathrm{2}}-a^{\mathrm{2}})}。（10）

式中， k = \frac{c}{E_{\mathrm{3}}} ( \mathrm{1} + v_{\mathrm{3}} ) + \frac{(\mathrm{1}+v_{\mathrm{2}})(\mathrm{1}-v_{\mathrm{2}})c^{\mathrm{3}}}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})} + \frac{(\mathrm{1}+v_{\mathrm{2}})b^{\mathrm{2}}c}{E_{\mathrm{2}}(c^{\mathrm{2}}-b^{\mathrm{2}})}

根据式（8）—式（10）可知，冻胀应力与衬砌内径（a）、外径（b）和冻结区外径（c）有关。在材料参数方面，冻胀应力与不同材料的弹性模量（E）和泊松比（v）也有关。同时由于黏性土颗粒间咬合能力较小，在冻结过程主要由胶结力和冻胀力控制变形^[15]。因此，黏性土不会像砂土那样容易发生脆性破裂。

1.3 围岩冻胀应变计算分析

在冻结过程中围岩受温度和周围约束共同作用^[16]。因此，冻胀总应变增量由2部分组成。即

1）主应力方向上冻胀应力引起的弹性应变增量。

2）冻结过程中孔隙水含量与体积变化引起的应变增量。

根据弹性理论，3个方向上的主应变分别为

{\epsilon }_{\mathrm{1}} = \frac{{\sigma }_{\mathrm{1}}-v({\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{2}})}{E} + \mathrm{\Delta } w \cdot \eta。（11）

{\epsilon }_{\mathrm{2}} = \frac{{\sigma }_{\mathrm{2}}-v({\sigma }_{\mathrm{3}}+{\sigma }_{\mathrm{1}})}{E} + \mathrm{\Delta } w \cdot \eta。（12）

{\epsilon }_{\mathrm{3}} = \frac{{\sigma }_{\mathrm{3}}-v({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}})}{E} + \mathrm{\Delta } w \cdot \eta。（13）

式中：σ ₁、σ ₂、σ ₃分别为3个主应力；ε ₁、ε ₂、ε ₃分别为3个主应力方向上的应变；v为冻土的泊松比；Δw为冻土含水量变化量；η为冻土冻胀率。

冻胀率指土体在冻结前后体积之差与冻结前的体积之比，即

{\eta }_{\mathrm{1}} = \frac{\mathrm{\Delta }V}{V} \times \mathrm{100} \%；

{\eta }_{\mathrm{2}} = \frac{\mathrm{\Delta }H}{H} \times \mathrm{100} \%。（14）

式中：η为冻胀率，%；ΔV和H分别为试样冻结前后的体积差和试样高度。

基于土体各向同性假设，3个方向上的应变可表示为

\epsilon = \frac{\sigma -\mathrm{2}{v}_{\mathrm{2}}\sigma }{E_{\mathrm{2}}} + \mathrm{\Delta } w \cdot \eta。（15）

式中：ε为每个主应力方向上的应变；σ为主应力；v ₂为冻结土体的泊松比；Δw为冻结土体含水量增量；E ₂为冻结土体弹性模量。

2 冻胀潜势的表现与计算分析

冻胀潜势的宏观表现是当减小或解除约束时，冻结土体在冻胀应力的方向上发生膨胀做功现象。其中冻胀应力不是常数，其大小随冻胀变形的变化而改变。根据能量守恒定律可知，贮存在物体中的冻胀潜势（U）应等于冻结土体在3个主应力方向上所做的功（W）。根据弹性理论，三维应力状态下单位土体的应变能为

W = \frac{{\sigma }_{\mathrm{1}}{\epsilon }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}{\epsilon }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}{\epsilon }_{\mathrm{3}}}{\mathrm{2}}。（16）

对于各向同性土体，若只考虑冻胀作用，则3个主方向上的应力和应变是一致的。因此式（16）化简为

W = \frac{\mathrm{3}}{\mathrm{2}} \sigma \cdot \epsilon。（17）

式中，σ和ε为各向同性冻结任一方向上的应力和应变。

因此，体积为V的土体其冻胀潜势为

U = W \cdot V = \frac{\mathrm{3}}{\mathrm{2}} \sigma \cdot \epsilon \cdot V。（18）

由于σ和ε与冻结温度密切相关，因此冻胀潜势和温度也是密切相关的。另外，σ与ε之间的力学关系可通过弹性模量表达，不同温度下的冻结土体与弹性模量呈线性关系。因此可以建立冻胀潜势与温度的关系。

3 案例分析

3.1 工程概况

某市轨道交通1号线一期工程（试验段）TJSG-2标联络通道里程为K33+000位于粉质黏土层。取样地点位于1^#联络通道与右线盾构隧道交汇处，深度为13.71 m，其物理力学参数见表1。该联络通道各参数如图2所示，各参数物理意义见表2。

3.2 冻胀力和冻胀变形测试

试样冻胀变形测试试验在高精度低温试验箱内进行，为了能够获取试样的实际温度变化，采用K型热电偶温度线布置在试样的不同高度进行监测。在试样的上方安装YWC-位移传感器监测冻胀位移。整个试验监测过程由电脑操作控制和记录，如图3所示。

试样制作为直径100 mm、高150 mm的圆柱体，压实度为90%。试验步骤如下。

1）将有机玻璃筒（直径100 mm、高210 mm）内壁均匀涂抹凡士林后放入试样，沿有机玻璃筒轴线均匀布置5个热敏电阻预埋孔。

2）试样放入有机玻璃筒后将热敏电阻通过预留孔插入试样内部。

3）用保温材料对有机玻璃筒进行保温，并在其上固定位移计。

4）将装置调整完毕之后放入冻结箱内。冻结温度由常温逐渐降至-3、-5、-8、-10、-15 ℃，每个冻结温度维持冻结72 h。

5）记录不同温度不同冻结历时的冻胀变化量。

3.3 温度对冻胀潜势的影响

将冻结后的试样在压力机上进行无侧限抗压试验^[18]，加载速度控制为0.01 mm/min，从而得到其应力应变关系，如图4所示。

通过应力-应变曲线图得出：冻结黏土试样在无侧限抗压试验中呈现应变硬化型。土体在刚加载阶段处于弹性变形阶段，随应力增加，变形表现为塑性屈服。

弹性模量与物质构成相关，不同温度时冻土内部水、冰的比例不同，弹性模量也将随之变化。冻结后孔隙水的结冰增强了土颗粒间的连接作用，从而改变了冻土的弹性模量^[19]。因此，温度降低后的微观表现为孔隙水结冰，而宏观表现则为弹性模量的增大。

冻结黏土的弹性模量指单轴抗压情况下最大抗压强度一半与其应变的比值^[20]。根据图4应力应变测试结果计算出不同温度下冻结土体的弹性模量，进一步线性拟合得到土体的弹性模量与不同温度基本呈线性关系，相关系数为0.989，温度越低，弹性模量越大，如图5所示。因此，根据土体在某温度下弹性模量和冻胀潜势与应力应变之间的关系，进行推导换算后代入图5中的拟合公式可知，冻胀潜势与温度呈负相关关系。即土体的温度越低则产生的冻胀潜势越大，则土体发生的冻胀程度就越明显。

3.4 结果与讨论

按照施工方案要求，当冻结壁温度达到-10 ℃时进行土体开挖。由上述可知，在-10 ℃下土体弹性模量E ₂取值为45.24 MPa，经计算参数k等于0.527× 10⁶ m³/N。当冻结温度为-10 ℃时，土体冻胀率为3.13%。将表2数据代入到式（8）可以计算出冻胀应力为0.338 MPa。根据冻胀防治要求，在冻结壁附近未冻土两侧各布设了2个卸压孔，每个卸压孔的最大压力控制在0.4 MPa以下。按模型计算出的冻胀应力0.338 MPa，小于最大控制压力，满足冻胀防治要求。

当土体冻胀率为3.13%时含水率为21.14%，由表1中初始含水量16.44%可得Δw为4.7%。将表2数据代入式（15）计算出冻胀应变值为0.446%。根据融沉控制要求，联络通道沉降累计大于1.0 mm时应进行融沉补偿注浆。根据模型计算0.446%应变会造成8.028 mm冻胀，略大于结构沉降控制量。因此，在对围岩冻结和解冻后，要注意冻结壁的冻胀变形以及沉降变化，及时采取防护措施，保证施工与支护结构的安全。

将应力、应变及体积代入到式（18）可以得出冻胀潜势为1.37×10⁶ J。该冻胀潜势表征了受约束冻结土体在-10 ℃温度下的做功现象，土体内部发生能量变化，并且以能量大小可体现土体的冻胀程度。受温度梯度影响下，土体的冻胀来源于原位水分冻结或是迁移水冻结。当冻结锋面移动速度超过水分迁移速度，土层中水分在迁移前原位冻结，冻胀的水分来源于初始含水量，当冻结锋面发展速度小于水分迁移速度，水分迁移会加剧冻结锋面区域冻胀量。因此，通过区分冻结区域与非冻结开挖区域，给出开挖过程由边界条件改变引起的应力、应变及冻胀潜势释放规律。可据此对土体冻胀力与冻胀变形，冷源功率和冻胀力等问题进行反向验证。

4 结论

本研究通过对某市联络通道冻结施工过程中冻胀应力与冻胀应变的理论计算，结合冻胀变形试验测试，主要得出以下结论。

1）通过计算冻胀力与冻胀应变，给出了基于冻胀应力和冻胀变形函数的冻胀潜势表示方法，即 U = W \cdot V = \frac{\mathrm{3}}{\mathrm{2}} \sigma \cdot \epsilon \cdot V，可用以表征围岩因冻结蓄积的能量。在-10 ℃温度下，经计算得到冻胀潜势为1.37×10⁶ J，进而从能量角度可构建工程中冷源输出功率与土体冻结能量之间的关系，为冻结法施工工程提供理论参考。

2）不同温度冻结条件下土体弹性模量随温度降低而逐步增大，但冻土温度与弹性模量间呈现线性的关系，相关系数为0.989。冻土应力的增加源于水相变成冰后的体积膨胀所造成的密实度增加（孔隙含量减小）。

参考文献

[1] 陈湘生.地层冻结法[M].北京：人民交通出版社，2013.

[2] CHEN X S.Artificial ground freezing[M].Beijing：China Communications Press，2013.

[3] 董新平，于少辉，张毅豪，地铁联络通道冻结法施工中涌水成因及防治[J].地下空间与工程学报，2022，18（1）：322-329，340.

[4] DONG X P， YU S H， ZHANG Y H，et al.The causes and prevention measures of water inflow at cross passage located in sandy gravel layers and constructed by artificial ground freezing[J].Chinese Journal of Underground Space and Engineering，2022，18（1）：322-329，340.

[5] 鲁先龙，陈湘生，陈曦.人工地层冻结法风险预控[J].岩土工程学报，2021，43（12）：2308-2314.

[6] LU X L， CHEN X S， CHEN X.Risk prevention and control of artificial ground freezing （AGF）[J].Chinese Journal of Geotechnical Engineering，2021，43（12）：2308-2314.

[7] 胡俊，刘勇，曾晖.新型管幕冻结法不同管幕填充形式的温度场数值对比分析[J].森林工程，2015，31（6）：135-141，145.

[8] HU J， LIU Y， ZENG H.Numerical analysis on the temperature field of a new type of pipe-roof freezing method with different fill-in procedures[J].Forest Engineering，2015，31（6）：135-141，145.

[9] 蔡海兵，彭立敏，郑腾龙.隧道水平冻结施工引起地表冻胀的历时预测模型[J].岩土力学，2012，33（6）：1761-1768.

[10] CAI H B， PENG L M， ZHENG T L.A duration prediction model of surface frost heave induced by tunnelling with horizontal freezing method[J].Rock and Soil Mechanics，2012，33（6）：1761-1768.

[11] 蔡海兵，程桦，姚直书，基于冻土正交各向异性冻胀变形的隧道冻结期地层位移数值分析[J].岩石力学与工程学报，2015，34（8）：1667-1676.

[12] CAI H B， CHENG H， YAO Z S，et al.Numerical analysis of ground displacement due to orthotropic frost heave of frozen soil in freezing period of tunnel[J].Chinese Journal of Rock Mechanics and Engineering，2015，34（8）：1667-1676.

[13] YANG P， KE J M， WANG J G，et al.Numerical simulation of frost heave with coupled water freezing，temperature and stress fields in tunnel excavation[J].Computers and Geotechnics，2006，33（6/7）：330-340.

[14] 吴东军，朱杰，李锐阳.人工冻结深厚表土冻胀特性及水热迁移的试验研究[J].煤炭技术，2018，37（5）：49-52.

[15] WU D J， ZHU J， LI R Y.Experimental study on freezing expansion and water-heat transference of artificially frozen deep overburden soil[J].Coal Technology，2018，37（5）：49-52.

[16] YUE F T， LYU S G， SHI R J，et al.Numerical analysis and spot-survey study of horizontal artificial ground freezing in tunnel connecting passage construction[J].Applied Mechanics and Materials，2015，744：969-977.

[17] 崔宏环，闫子麟，何静云，季冻区粉质黏土路基的冻胀特性的模型试验研究[J].科学技术与工程，2020，20（2）：814-820.

[18] CUI H H， YAN Z L， HE J Y，et al.Model experimental study on frost heaving characteristics of silty clay subgrade in seasonal frost region[J].Science Technology and Engineering，2020，20（2）：814-820.

[19] 董建华，吴晓磊，师利君.不同冻胀模式下冻结隧道施工引起的地表竖向位移时空解[J].中国公路学报，2020，10（17）：1-12.

[20] DONG J H， WU X L， SHI L J.Spatial and temporal solutions of vertical displacement of ground surface caused by freezing tunnel construction under different frost heave modes[J].China Journal of Highwan and Transport，2020，10（17）：1-12.

[21] 王思文，姜国辉，李玉清，不同系统单向冻结条件下土体冻胀率与冻结速率分析[J].水电能源科学，2021，39（10）：164-167.

[22] WANG S W， JIANG G H， LI Y Q，et al.Analysis of frost heave rate and freezing rate of soil under unidirectional freezing conditions of different systems[J].Water Resources and Power，2021，39（10）：164-167.

[23] 侯兆领.压实度及含水率对盐渍土冻胀规律影响试验研究[J].公路，2023，68（3）：316-320.

[24] HOU Z L.Research on the influence of compactness and moisture content on frost heaving law of saline soil[J].Highway，2023，68（3）：316-320.

[25] 杨平，王岩梓，刁鹏程.上覆荷载对重塑粉质黏土补水冻胀特性影响[J].郑州大学学报（工学版），2022，43（6）：83-89.

[26] YANG P， WANG Y Z， DIAO P C.Influence of load on frost heaving characteristics of remolded silty clay[J].Journal of Zhengzhou University （Engineering Science），2022，43（6）：83-89.

[27] 刘振亚，刘建坤，李旭，毛细黏聚与冰胶结作用对非饱和粉质黏土冻结强度及变形特性的影响[J].岩石力学与工程学报，2018，37（6）：1551-1559.

[28] LIU Z Y， LIU J K， LI X，et al.Effect of capillary cohesion and ice cementation on strength and deformation ofunsaturated frozen frozen soil[J].Chinese Journal of Rock Mechanics and Engineering，2018，37（6）：1551–1559.

[29] 曹善鹏，夏才初，周舒威，考虑寒区隧道围岩冻结温度渐变的冻胀力解析解[J].同济大学学报（自然科学版），2024，52（3）：360-369.

[30] CAO S P， XIA C C， ZHOU S W，et al.Analytical solution forfrost heave of cold region tunnel caused by non-uniform frost heave of surrounding rock considering freezing temperature gradient[J].Journal of Tongji University （Natural Science Edition），2024，52（3）：360-369.

[31] 张帅伟，李顺群，程学磊，基于荷载结构法的冻结壁内力与变形研究[J].建井技术，2024，45（2）：90-94.

[32] ZHANG S W， LI S Q， CHENG X L，et al.Research on internal force and deformation of frozen wall based on load structure method[J].Mine Construction Technology，2024，45（2）：90-94.

[33] 中华人民共和国住房和城乡建设部. GB/T 50123—2019 土工试验方法标准 [S].北京：中国建筑工业出版社，2019.

[34] Ministry of Housingand Urban-Rural Development of the People's Republic of China. GB/T 50123—2019，standard for geotechnical testing method [S].Beijing：China Architecture& Building Press，2019.

[35] 李寒.滨海人工冻土压折力学特性及影响因素分析[D].福州：福建工程学院，2024.

[36] LI H.Analysis of mechanical properties and influencing factors of artificial frozen soil in coastal areas[D].Fuzhou：Fujian University of Technology，2024.

[37] 张雅琴，杨平，江汪洋，粉质黏土冻土三轴强度及本构模型研究[J].土木工程学报，2019，52（S1）：8-15.

[38] ZHANG Y Q， YANG P， JIANG W Y，et al.Study on triaxial strength and constitutive model of frozen silty clay[J].ChinaCivil Engineering Journal，2019，52（S1）：8-15.

基金资助

国家自然科学基金项目(52478354)

甘肃省重点研发计划项目(23YFGA0038)

天津市研究生科研创新项目资助(2022SKYZ329)