REGRESSION ANALYSIS MODEL OF RELATIVE GAS EMISSION RATE OF COALBEDS WITH SAME LEVEL
-
摘要: 利用数理统计学原理,建立了淮南谢二矿-450m水平不同煤层的相对瓦斯涌出量与煤层厚度、煤层结构破碎程度的回归分析模型,表明在同一地质单元同水平条件下,煤层厚度越大,结构越破碎,其相对瓦斯涌出量越大。Abstract: Based on the theory of mathematical statistics,the regression analysis model between relative gas emission rate(q) and the thickness(H) and breakage degree(B) of coalbeds with same level of -450 m in Xiejiaji No.2 coal mine,Huinan coalfield,is estabished.The work reveals that under the conditions of same geological background and with the same level,the thicker the thickness and the higher the breakage degree of coalbeds,the more the relative gas emission rate.It is believed that the regression analysis is a quantitative method to study gas geology.
-
Keywords:
- coalbed /
- gas emission rate /
- regression analysis
-
随着我国基础设施建设的快速发展，富水地层的工程建设项目也日益增多。在该类地区进行工程建设时，如何确保施工安全和高效成为了一个亟待解决的问题[1-4]。人工地层冻结法作为一种有效的地层加固、封水方法，已经在富水地层地下工程施工中得到了广泛应用[5-7]。但如何确定合理的冻结参数、制冷剂，以及准确分析冻结壁的力学特性是影响该工法施工效果的关键因素。
在以往的冻结法施工过程中，多采用单一冷媒进行地层冻结，大量学者针对该类冻结壁的力学特性展开了系统的研究。胡向东[8-9]、荣传新[10]、王彬[11-12]和曹雪叶[13-15]等将冻结壁视作一种随半径变化的功能梯度材料，分别基于M-C(Mohr-Coulomb)强度准则[8-9,12]、D-P(Druker-Prager)强度准则[10-11]和双剪统一强度准则[13-15]，推导得出了单排[8,10,13-15]、双排[9,11]以及三排管[12]所形成的冻结壁的应力解析解，并运用得出的公式对上述冻结壁的应力分布规律以及承载力进行了计算分析。冻结壁并非孤立存在于土体之中，它与邻近的未冻土(岩)体存在相互作用。杨维好等[16-18]在考虑冻结壁与周围岩层相互作用的基础上，推导出了在不同状态下(弹性、弹塑性及塑性)均质冻结壁的应力计算公式。王彬等[19]通过引入应力传递系数来反映开挖卸载过程中的卸荷效应，构建了一个描述冻结壁与周边土体相互作用的力学模型，推导得出了无地下水以及有地下水[20-21]两种状态下非均质冻结壁的弹塑性应力计算公式。
然而，在面临复杂的地质条件和施工要求时，单一的冷媒往往无法满足需求。因此，多冷媒联合应用逐渐成为了研究的热点。多冷媒联合冻结法，指的是在同一个冻结工程中，根据地层的具体特点和施工的实际需求，选择两种或两种以上的冷媒进行组合使用。该技术综合了多种冷媒的优势，不仅提高了冻结效率，还实现了节能环保目标。多冷媒联合冻结技术的应用，为土建工程、矿山工程、环境工程、水利工程、城市管线工程等领域提供了有力的技术支持，确保了工程安全、高效进行。目前，部分学者已经对多冷媒联合冻结技术进行了研究。袁杰[22]结合了液氮冻结和盐水冻结各自的优势，实施了双排孔布局的液氮冻结管，实现了土层的迅速冻结加固。当液氮冻结达到预定的冻结壁厚度时，便中断液氮的供应，并切换至传统的盐水冻结循环，以维持土壤中冻结壁的稳定。这样的过程不仅确保了冻结效果，同时也显著降低了持续使用液氮所带来的高昂成本。在斜井开发领域，叶庆[23]首先推导了单管联合冻结(先盐水后液氮、先液氮后盐水)和双管同时冻结的温度场解析公式；其次，通过数值计算建立了平面导热模型；最后对榆林袁大滩煤矿副斜井第三、四段进行了现场实测研究，进而分析了联合冻结温度场的演变规律。权晓甜[24]在宁波顶管D段管廊改造项目中，应用了液氮与盐水的联合冻结技术，使用COMSOL软件，建立了3种不同冻结顺序的三维有限元模型，以此来分析冻结帷幕的发展过程。董涛等[25]提出了一种创新的多冷媒联合制冷技术，基于珠江三角洲水资源配置项目中的某盾构刀盘更换工程为案例，对多冷媒联合冻结技术的可行性进行了深入分析。
近年来，CO2作为天然制冷剂开始被行业所重视。CO2可以同时作为制冷剂和冷媒，从而减少中间热交换次数，冷媒循环系统简洁高效。并且其黏度低，流动性好，适合长距离恒温输送。与液氮相比，可以循环制冷，具有更高的经济性。因此，在富水地层中可采用盐水−二氧化碳多冷媒的冻结方法。其中，盐水冷媒的优势在于其良好的流动性和热导率，能够有效地在地层中分布和传递热量，形成较为均匀的冻土帷幕；而二氧化碳冷媒的优势在于其低冻结温度和高效的吸热性能，能够快速降低地层温度，提高冷冻速率。通过结合两者的优势，盐水−二氧化碳多冷媒联合冻结方法能够实现更快速、更均匀的冻结效果，提高工程效率，降低工程成本[26-27]。
与单一冷媒形成的冻结壁相比，多冷媒联合冻结帷幕的力学特性显得更为复杂，这主要是由于不同冷媒之间的相互作用和温度梯度导致的力学性能差异。目前，针对这一领域的研究还相对较少，缺乏系统性的深入探讨。笔者针对盐水−二氧化碳多冷媒联合冻结帷幕，深入分析了温度对冻结壁力学参数非均匀性的影响，并结合冻土强度屈服准则，对多冷媒联合双排管非均质冻结壁的力学特征进行了详细的弹塑性分析。这项研究不仅揭示了多冷媒联合冻结帷幕在不同屈服准则下的力学响应，还为进一步优化冻结壁设计、提高工程稳定性和安全性提供了科学依据，对于富水地层多冷媒联合双排管冻结壁的工程实践具有重要的指导意义。
1 冻结壁温度场分布规律
1.1 温度场模型
针对盐水−二氧化碳多冷媒联合双排管冻结壁，选取某段连续的区间作为研究对象，冻结壁局部如图1所示[12]，内侧冻结管用盐水作为低温冷媒，而外侧冻结管用二氧化碳作为低温冷媒。在冻结法施工过程中，为了确保井筒挖掘的安全性和效率，同时为井壁施工创造一个稳定的施工环境，施工方会根据井筒的宽度和直径，对被冻土体进行一定程度的挖掘。因此，冻结壁的实际厚度应当是其原始厚度减去计划挖掘部分的厚度[12]。
已有研究成果表明，距离1/4管距处的特征截面的平均温度能准确反映整个冻结壁的平均温度[28]，为此选取该截面的温度分布曲线代替整个冻结壁的温度曲线，进而开展冻结壁力学特性分析。在非均质冻结壁计算理论中，通常会将冻结壁视为随温度呈一次函数变化或者二次函数变化的非均质材料[29-31]，这样的近似方式不仅可以简化计算，也能够提高冻结壁平均温度计算的准确性和实用性。冻结温度场近似曲线形式需要根据冻结壁实际情况进行选择。
盐水−二氧化碳复合冷媒冻结壁等效温度场分布如图2所示[9]。其中，温度坐标轴的正向表示温度梯度降低的方向。Ⅰ、Ⅱ、Ⅲ区分别对应开挖面至内排冻结管、内排至外排冻结管以及外排冻结管至冻结壁外缘的同心圆区域，全文统一。在“盐水−二氧化碳”复合冻结壁中，两种冷媒的温差较大，且冻结管与土体之间的温度梯度较大，为此选用三段一次函数来近似等效冻结壁截面处温度场的变化规律。
![]()
图 2 盐水−二氧化碳多冷媒冻结壁等效温度场Figure 2. Equivalent temperature field of brine-CO2 frozen wall通过图中4个温度特征点可以确定冻结壁的温度场函数t(R)为：
$$ t(R) = \left\{ \begin{aligned} & {\frac{{{t_1} - {t_{\text{n}}}}}{{{R_1} - {R_{\text{B}}}}}(R - {R_{\text{B}}}) + {t_{\text{n}}}}\qquad {R \in \left[ {{R_{\text{B}}},{R_1}} \right)} \\ & {\frac{{{t_2} - {t_1}}}{{{R_2} - {R_1}}}(R - {R_1}) + {t_1}}\qquad\; {R \in \left[ {{R_1},{R_2}} \right)} \\ & {\frac{{{t_2} - {t_0}}}{{{R_2} - {R_{\text{H}}}}}(R - {R_{\text{H}}}) + {t_0}}\qquad{R \in \left[ {{R_2},{R_{\text{H}}}} \right]} \end{aligned} \right. $$ (1) 1.2 基本假定
为了简化计算得到易于工程应用的弹性、弹塑性和塑性解析解，现对冻结壁力学模型提出以下基本假设：
(1)冻结壁视为材料性质(弹性模量E、黏聚力c)随相对半径呈三段一次函数变化的非均质材料。
(2)当冻结壁由弹性状态进入塑性状态后，认为冻结壁体积不可压缩。
(3)冻结壁材料性质与温度大致呈一次函数的关系。由于当温度变化时，内摩擦角γ和泊松比μ影响很小，因此，本文不作考虑。
(4)考虑到本文中的研究对象具有几何对称性以及荷载和边界条件是轴对称的特性。为了简化计算，可以将冻结壁视为平面轴对称应变问题。
2 冻结壁力学模型
2.1 冻土力学参数
已有相关研究证明冻土的材料性质与温度大致呈一次函数的关系[8-12,32]。因此，当冻结温度场随相对半径呈线性变化时，冻土的弹性模量E和黏聚力c随相对半径r呈线性变化。冻土的弹性模量E(R)和黏聚力c(R)都服从一次函数关系[32]，可设：
$$ \left\{ \begin{aligned} & {E(R) = at(R) + {E_0}} \\ & {c(R) = lt(R) + {c_0}} \end{aligned} \right. $$ (2) 联立式(1)、式(2)，求得弹性模量E(r)和黏聚力c(r)的表达式分别为：
$$ E(r) = \left\{ \begin{aligned} & {{a_1}r + {b_1}}\quad{r \in \left[ {1,{r_1}} \right)} \\ & {{a_2}r + {b_2}}\quad{r \in \left[ {{r_1},{r_2}} \right)} \\ & {{a_3}r + {b_3}}\quad{r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{aligned} \right. $$ (3) $$ c(r) = \left\{ \begin{aligned} & {{l_1}r + {m_1}}\quad{r \in \left[ {1,{r_1}} \right)} \\ & {{l_2}r + {m_2}}\quad{r \in \left[ {{r_1},{r_2}} \right)} \\ & {{l_3}r + {m_3}}\quad{r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{aligned} \right. $$ (4) 式(3)、式(4)中：
$$ {a_1} = a\frac{{{t_1} - {t_{\text{n}}}}}{{{r_1} - 1}}\text{；}{b_1} = a\left( { - \frac{{{t_1} - {t_{\text{n}}}}}{{{r_1} - 1}} + {t_{\text{n}}}} \right) + {E_0}\text{；} $$ $$ {a_2} = a\frac{{{t_2} - {t_1}}}{{{r_2} - {r_1}}}\text{；}{b_2} = a\left( { - \frac{{{t_2} - {t_1}}}{{{r_2} - {r_1}}}{r_1} + {t_1}} \right) + {E_0}\text{；} $$ $$ {a_3} = a\frac{{{t_2} - {t_0}}}{{{r_2} - {r_{\text{H}}}}}\text{；}{b_3} = a\left( { - \frac{{{t_2} - {t_0}}}{{{r_2} - {r_{\text{H}}}}}{r_{\text{H}}} + {t_0}} \right) + {E_0}。 $$ $$ {l_1} = l\frac{{{t_1} - {t_{\text{n}}}}}{{{r_1} - 1}}\text{；}{m_1} = l\left( { - \frac{{{t_1} - {t_{\text{n}}}}}{{{r_1} - 1}} + {t_{\text{n}}}} \right) + {c_0}\text{；} $$ $$ {l_2} = l\frac{{{t_2} - {t_1}}}{{{r_2} - {r_1}}}\text{；}{m_2} = l\left( { - \frac{{{t_2} - {t_1}}}{{{r_2} - {r_1}}}{r_1} + {t_1}} \right) + {c_0}\text{；} $$ $$ {l_3} = l\frac{{{t_2} - {t_0}}}{{{r_2} - {r_{\text{H}}}}}\text{；} {m_3} = l\left( { - \frac{{{t_2} - {t_0}}}{{{r_2} - {r_{\text{H}}}}}{r_{\text{H}}} + {t_0}} \right) + {c_0} 。 $$ 2.2 冻土强度准则
统一强度准则是用于评估和预测冻土在复杂应力状态下的破坏行为的一种屈服准则[21,28-31]。该准则通过引入一个参数
$\xi $ ，将不同的强度准则联系起来，从而描述不同材料的强度特性。以压应力为正，轴对称平面应变状态下的统一强度屈服准则表达式如下[16-18]：
$$ \sigma _\theta ^i = \lambda \cdot \sigma _r^i + \omega \cdot c(r) $$ (5) 式中：上标i=1,2,3为冻结壁受力分析区的3个区域。
其中：λ和ω表达式如下：
$$ \lambda = \frac{{(2 + \xi ) + (2 + 3\xi )\sin \gamma }}{{(2 + \xi )(1 - \sin \gamma )}} $$ (6) $$ \omega = \frac{{4(1 + \xi )\cos \gamma }}{{(2 + \xi )(1 - \sin \gamma )}} $$ (7) 式(6)、式(7)中：
$\xi $ 为统一强度准则的参数，它决定了不同强度准则的权重，取值$0 \leqslant \xi \leqslant 1$ 。当$ \xi ＝0 $ 时，退化为M-C强度准则；当$ \xi ＝1 $ 时，退化为双剪统一强度准则。通过调整$\xi $ 的值，可以描述不同材料的强度特性。在实际工程中，由于被冻土体的性质差别较大，土体冻结之后适用的屈服准则也不相同，为此本研究选取冻土常用的以下4种屈服准则对多冷媒复合冻结壁的力学特性进行计算分析。不同屈服准则中λ和ω的具体表达式见表1[21]。
表 1 λ和ω在不同屈服准则下的表达式Table 1. Expression of λ and ω under different yield criteria屈服准则 λ ω M-C强度准则 $\dfrac{{1 + \sin \gamma }}{{1 - \sin \gamma }}$ $ \dfrac{{2\cos \gamma }}{{1 - \sin \gamma }} $ D-P强度准则 $\dfrac{{\sqrt {3 + {{\sin }^2}\gamma } + \sqrt 3 \sin \gamma }}{{\sqrt {3 + {{\sin }^2}\gamma } - \sqrt 3 \sin \gamma }}$ $\dfrac{{2\sqrt 3 \cos \gamma }}{{\sqrt {3 + {{\sin }^2}\gamma } - \sqrt 3 \sin \gamma }}$ 广义Tresca
准则$\dfrac{{\sqrt {3 + {{\sin }^2}\gamma } + 2\sin \gamma }}{{\sqrt {3 + {{\sin }^2}\gamma } - 2\sin \gamma }}$ $\dfrac{{4\cos \gamma }}{{\sqrt {3 + {{\sin }^2}\gamma } - 2\sin \gamma }}$ 双剪统一强度
准则$\dfrac{{3 + 5\sin \gamma }}{{3(1 - \sin \gamma )}}$ $\dfrac{{8\cos \gamma }}{{3(1 - \sin \gamma )}}$ 3 冻结壁弹塑性区应力计算公式
3.1 冻结壁弹性区应力场求解
当冻结壁的外荷载小于弹性极限承载力时，冻结壁处于弹性状态，需满足下列基本方程[33]：
平衡方程：
$$ \frac{{{\text{d}}\sigma _r^{\text{e}}}}{{{\text{d}}r}} + \frac{{\sigma _r^{\text{e}} - \sigma _\theta ^{\text{e}}}}{r} = 0 $$ (8) 几何方程：
$$ \left\{ \begin{aligned} & {\varepsilon _r^{\text{e}} = {\text{d}}u_r^{\mathrm{e}}/{\text{d}}r} \\ & \varepsilon _\theta ^{\text{e}} = u_r^{\text{e}}/r\end{aligned} \right. $$ (9) 物理方程：
$$ \left\{ \begin{aligned} & {\varepsilon _r^{\text{e}} = \frac{{1 - {\mu ^2}}}{{E(r)}}\sigma _r^{\text{e}} - \frac{{\mu (1 + \mu )}}{{E(r)}}\sigma _\theta ^{\text{e}}} \\ & {\varepsilon _\theta ^{\text{e}} = - \frac{{\mu (1 + \mu )}}{{E(r)}}\sigma _r^{\text{e}} + \frac{{1 - {\mu ^2}}}{{E(r)}}\sigma _\theta ^{\text{e}}} \end{aligned} \right. $$ (10) 式中：上标e为冻结壁弹性区。
$$ \left\{ \begin{aligned} & {\sigma _r^{\text{e}} = \psi /r} \\ & \sigma _\theta ^{\text{e}} = \psi '\end{aligned} \right. $$ (11) 通过式(9)—式(11)，得：
$$ {r^2}\psi '' + r\left( {1 - \frac{{E'(r)}}{{E(r)}}r} \right)\psi ' + \left( {\frac{\mu }{{1 - \mu }}\frac{{E'(r)}}{{E(r)}}r - 1} \right)\psi = 0 $$ (12) 将式(3)中
$E(r)$ 代入式(12)，进一步化简得：$$ {r^2}\psi '' + r\frac{{{b_i}}}{{{a_i} + {b_i}}}\psi ' - \frac{{{b_i}}}{{{a_i} + {b_i}}}\psi = 0 $$ (13) 通过求解式(13)，可得应力函数
$\psi $ 的通解为：$$ \psi = r{M_i}\left( { - {a_i}\frac{1}{r} - \frac{{{b_i}}}{2}\frac{1}{{{r^2}}}} \right) + {N_i}r $$ (14) 式(14)的边界条件为：
$$ \left\{ \begin{aligned} & {\psi = 0}\qquad\;\;\, {\left( {r = 1} \right)} \\ & {\psi = {F_1}{r_1}}\quad\; {\left( {r = {r_1}} \right)} \\ & {\psi = {F_2}{r_2}}\quad\; {\left( {r = {r_2}} \right)} \\ & {\psi = {F_{\text{H}}}{r_{\text{H}}}}\quad{\left( {r = {r_{\text{H}}}} \right)} \end{aligned} \right. $$ (15) 代入应力分量，可得各区域应力分量表达式为：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅰ}})}}} = \frac{{{F_1}{A_2}}}{{{A_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅰ}})}}} = \frac{{{F_1}{A_3}}}{{{A_1}}}} \end{aligned} \right.\quad r \in \left[ {1,{r_1}} \right) $$ (16) $$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅱ}})}}} = \frac{{\left( {{F_2} - {F_1}} \right){B_3} + {B_2}}}{{{B_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}} = \frac{{\left( {{F_2} - {F_1}} \right){B_4} + {B_2}}}{{{B_1}}}} \end{aligned}\right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_1},{r_2}} \right)} \end{array} $$ (17) $$ \left\{\begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_2}){C_3} + {C_2}}}{{{C_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}} = \frac{{\left( {{F_{\text{H}}} - {F_2}} \right){C_4} + {C_2}}}{{{C_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{array} $$ (18) 式(16)—式(18)中：
$$ {A_1} = {a_1}\left( {1 - \frac{1}{{{r_1}}}} \right) + \frac{{{b_1}}}{2}\left( {1 - \frac{1}{{r_1^2}}} \right) \text{；} $$ $$ {A_2} = {a_1}\left( {1 - \frac{1}{r}} \right) + \frac{{{b_1}}}{2}\left( {1 - \frac{1}{{{r^2}}}} \right) \text{；} $$ $$ {A_3} = {a_1} + \frac{{{b_1}}}{2}\left( {1 + \frac{1}{{{r^2}}}} \right) \text{；} $$ $$ {B_1} = {a_2}\left( {\frac{1}{{{r_1}}} - \frac{1}{{{r_2}}}} \right) + \frac{{{b_2}}}{2}\left( {\frac{1}{{r_1^2}} - \frac{1}{{r_2^2}}} \right) \text{；} $$ $$ {B_2} = {a_2}\left( {\frac{{{F_2}}}{{{r_1}}} - \frac{{{F_1}}}{{{r_2}}}} \right) + \frac{{{b_2}}}{2}\left( {\frac{{{F_2}}}{{{r_1}^2}} - \frac{{{F_1}}}{{{r_2}^2}}} \right) \text{；} $$ $$ {B_3} = - \frac{{{a_2}}}{r} - \frac{{{b_2}}}{{2{r^2}}} \text{；} {B_4} = \frac{{{b_2}}}{{2{r^2}}} \text{；} $$ $$ {C_1} = {a_3}\left( {\frac{1}{{{r_2}}} - \frac{1}{{{r_{\text{H}}}}}} \right) + \frac{{{b_3}}}{2}\left( {\frac{1}{{{r_2}^2}} - \frac{1}{{{r_{\text{H}}}^2}}} \right) \text{；} $$ $$ {C_2} = {a_3}\left( {\frac{{{F_{\text{H}}}}}{{{r_2}}} - \frac{{{F_2}}}{{{r_{\text{H}}}}}} \right) + \frac{{{b_3}}}{2}\left( {\frac{{{F_{\text{H}}}}}{{{r_2}^2}} - \frac{{{F_2}}}{{{r_{\text{H}}}^2}}} \right) \text{；} $$ $$ {C_3} = - \frac{{{a_3}}}{r} - \frac{{{b_3}}}{{2{r^2}}} \text{；} {C_4} = \frac{{{b_3}}}{{2{r^2}}} 。 $$ 对于平面轴对称问题有：
$$ u_r^{\text{e}} = \frac{{(1 + \mu )r}}{{E(r)}}\left[ { - \mu \sigma _r^{\text{e}} + (1 - \mu )\sigma _\theta ^{\text{e}}} \right] $$ (19) 根据r1、r2处的位移连续条件可以得到r1、r2处的应变连续条件如下：
$$ \left\{ \begin{aligned} & {u_1^ - = u_1^ + \to {\varepsilon _1^ -} = {\varepsilon _1^ +} } \\ & {u_2^ - = u_2^ + \to {\varepsilon _2^ -} = {\varepsilon _2^ +} } \end{aligned} \right.\begin{array}{*{20}{c}} {}&\begin{gathered} (r = {r_1}) \\ (r = {r_2}) \\ \end{gathered} \end{array} $$ (20) 为了简化计算，由冻结壁的均匀连续性假定，可假设r1和r2处泊松比分别相等。由此可以得出：
$$ \left\{ {\begin{array}{*{20}{l}} { - \mu \sigma _r^{{\text{e({\text{Ⅰ}})}}} + (1 - \mu )\sigma _\theta ^{{\text{e({\text{Ⅰ}})}}} = - \mu \sigma _r^{{\text{e({\text{Ⅱ}})}}} + (1 - \mu )\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}}}&{(r = {r_1})} \\ { - \mu \sigma _r^{{\text{e({\text{Ⅱ}})}}} + (1 - \mu )\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}} = - \mu \sigma _r ^{{\text{e({\text{Ⅲ}})}}} + (1 - \mu )\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}}}&{(r = {r_2})} \end{array}} \right. $$ (21) 将式(16)、式(17)和式(18)代入式(21)，化简后可得F1和F2。
3.2 冻结壁塑性区应力场求解
塑性区的应力平衡方程[34]：
$$ \frac{{{\text{d}}\sigma _r^{\text{p}}}}{{{\text{d}}r}} + \frac{{\sigma _r^{\text{p}} - \sigma _\theta ^{\text{p}}}}{r} = 0 $$ (22) 式中：上标p为冻结壁塑性区。
将强度准则式(5)代入式(22)，可得：
$$ \frac{{{\text{d}}\sigma _r^{\text{p}}}}{{{\text{d}}r}} + \frac{{1 - \lambda }}{r}\sigma _r^{\text{p}} - \frac{\omega }{r}c(r) = 0 $$ (23) 通过解析微分方程式(23)可以得出冻结壁塑性区域内径向应力的通解：
$$ \sigma _r^{\text{p}} = \left\{ \begin{aligned} & {D{l_1}r + P{m_1} + {K_1}{r^{\lambda - 1}}}\quad{r \in \left[ {1,{r_1}} \right)}&{} \\ & {D{l_2}r + P{m_2} + {K_2}{r^{\lambda - 1}}}\quad{r \in \left[ {{r_1},{r_2}} \right)}&{} \\ & {D{l_3}r + P{m_3} + {K_3}{r^{\lambda - 1}}}\quad{r \in \left[ {{r_2},{r_{\text{H}}}} \right]}&{} \end{aligned} \right. $$ (24) 式(24)中：D=ω/(2−λ)，P=ω/(1−λ)。
由r＝1，r＝r1，r＝r2处连续性的应力边界条件可求得：
$$ {K_1} = - (D{l_1} + P{m_1}) \text{；} $$ $$ {K_2} = D({l_1} - {l_2})r_1^{2 - \lambda } + P({m_1} - {m_2})r_1^{1 - \lambda } - (D{l_1} + P{m_1}) \text{；} $$ $$ \begin{split} &{K}_{3}=D({l}_{2}-{l}_{3}){r}_{2}^{2-\lambda }+P({m}_{2}-{m}_{3}){r}_{2}^{1-\lambda }+D({l}_{1}-{l}_{2}){r}_{1}^{2-\lambda }+\\ &\qquad P({m}_{1}-{m}_{2}){r}_{1}^{1-\lambda }-(D{l}_{1}+P{m}_{1})。\end{split} $$ 由此，冻结壁塑性区的环向应力通过下式求得：
$$ \sigma _\theta ^{\text{p}} = \left\{ \begin{aligned} & {\lambda \sigma _r^{{\text{p({\text{Ⅰ}})}}} + \omega {c_1}(r)} \\ & {\lambda \sigma _r^{{\text{p({\text{Ⅱ}})}}} + \omega {c_2}(r)} \\ & {\lambda \sigma _r^{{\text{p({\text{Ⅲ}})}}} + \omega {c_3}(r)} \end{aligned} \right.\begin{array}{*{20}{c}} {} \\ {} \\ {} \end{array}\begin{aligned} & {r \in \left[ {1,{r_1}} \right)} \\ & {r \in \left[ {{r_1},{r_2}} \right)} \\ & {r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{aligned} $$ (25) 3.3 冻结壁弹塑性应力场求解
当冻结壁承受的外荷载大于弹性极限承载力时，冻结壁开始出现塑性区，进入弹塑性状态。此时冻结壁弹塑性力学计算模型如图4所示。
![]()
图 4 不同状态下冻结壁力学计算模型Figure 4. Mechanical calculation model of frozen wall in different states假设冻结壁的塑性半径为rc，在弹塑性分界面(r＝rc)处，冻结壁的径向应力和环向应力满足应力连续条件和屈服准则。冻结壁应力计算公式主要分为以下3种情况：
当满足rc∈[1,r1)时，在[1,rc)区间内，冻结壁处于塑性状态，而在[rc,rH]区间内，冻结壁处于弹性状态。计算模型如图4a所示[12]。
冻结壁在Ⅰ区塑性应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p(I)}}} = D{l_1}r + P{m_1} + {K_1}{r^{\lambda {{ - 1}}}}} \\ & {\sigma _\theta ^{{\text{p(I)}}} = \lambda \sigma _r^{{\text{p(I)}}} + \omega {c_1}(r)}\end{aligned} \right.\begin{array}{*{20}{c}} {r \in \left[ {1,{r_{\text{c}}}} \right)} \end{array} $$ (26) 冻结壁在Ⅰ区弹性应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e(I)}}} = \frac{{({F_1} - {F_{\text{c}}}){G_3} + {G_2}}}{{{G_1}}}} \\ & {\sigma _\theta ^{{\text{e(I)}}} = \frac{{({F_1} - {F_{\text{c}}}){G_4} + {G_2}}}{{{G_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_{\text{c}}},{r_1}} \right)} \end{array} $$ (27) 式(27)中：
$$ {G_1} = {a_1}\left( {\frac{1}{{{r_{\text{c}}}}} - \frac{1}{{{r_1}}}} \right) + \frac{{{b_1}}}{2}\left( {\frac{1}{{r_{\text{c}}^{\text{2}}}} - \frac{1}{{r_1^2}}} \right) ; $$ $$ {G_2} = {a_1}\left( {\frac{{{F_1}}}{{{r_{\text{c}}}}} - \frac{{{F_{\text{c}}}}}{{{r_1}}}} \right) + \frac{{{b_1}}}{2}\left( {\frac{{{F_1}}}{{r_{\text{c}}^2}} - \frac{{{F_{\text{c}}}}}{{r_1^2}}} \right) ; $$ $$ {G_3} = - \frac{{{a_1}}}{r} - \frac{{{b_1}}}{{2{r^2}}} ;\quad {G_4} = \frac{{{b_1}}}{{2{r^2}}}。 $$ 冻结壁在Ⅱ区弹性应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅱ}})}}} = \frac{{\left( {{F_2} - {F_1}} \right){B_3} + {B_2}}}{{{B_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}} = \frac{{\left( {{F_2} - {F_1}} \right){B_4} + {B_2}}}{{{B_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_1},{r_2}} \right)} \end{array} $$ (28) 冻结壁在Ⅲ区弹性应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_2}){C_3} + {C_2}}}{{{C_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_2}){C_4} + {C_2}}}{{{C_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{array} $$ (29) 当满足rc∈[r1,r2)时，在[1,rc)区间内，冻结壁处于塑性状态，而在[rc,rH]区间内，冻结壁处于弹性状态。计算模型如图4b所示[12]。
冻结壁在Ⅰ区塑性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅰ}})}}} = D{l_1}r + P{m_1} + {K_1}{r^{\lambda - 1}}} \\ & {\sigma _\theta ^{{\text{p({\text{Ⅰ}})}}} = \lambda \sigma _r^{{\text{p(I)}}} + \omega {c_1}(r)}\end{aligned} \right.\begin{array}{*{20}{c}} {r \in \left[ {1,{r_1}} \right)} \end{array} $$ (30) 冻结壁在Ⅱ区塑性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅱ}})}}} = D{l_2}r + P{m_2} + {K_2}{r^{{\lambda -1}}}} \\ &{\sigma _\theta ^{{\text{p({\text{Ⅱ}})}}} = \lambda \sigma _r^{{\text{p({\text{Ⅱ}})}}} + \omega {c_2}(r)}\end{aligned}\right.\begin{array}{*{20}{c}} {r \in \left[ {{r_1},{r_{\text{c}}}} \right)} \end{array} $$ (31) 冻结壁在Ⅱ区弹性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅱ}})}}} = \frac{{({F_2} - {F_{\text{c}}}){H_3} + {H_2}}}{{{H_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}} = \frac{{({F_2} - {F_{\text{c}}}){H_4} + {H_2}}}{{{H_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_{\text{c}}},{r_2}} \right)} \end{array} $$ (32) 式(32)中：
$$ {H_1} = {a_2}\left( {\frac{1}{{{r_{\text{c}}}}} - \frac{1}{{{r_2}}}} \right) + \frac{{{b_2}}}{2}\left( {\frac{1}{{r_{\text{c}}^{\text{2}}}} - \frac{1}{{r_2^2}}} \right)\text{；} $$ $$ {H_2} = {a_2}\left( {\frac{{{F_2}}}{{{r_{\text{c}}}}} - \frac{{{F_{\text{c}}}}}{{{r_2}}}} \right) + \frac{{{b_2}}}{2}\left( {\frac{{{F_2}}}{{r_{\text{c}}^2}} - \frac{{{F_{\text{c}}}}}{{r_2^2}}} \right)\text{；} $$ $$ {H_3} = - \frac{{{a_2}}}{r} - \frac{{{b_2}}}{{2{r^2}}}\text{；}{H_4} = \frac{{{b_2}}}{{2{r^2}}}。 $$ 冻结壁在Ⅲ区弹性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_2}){C_3} + {C_2}}}{{{C_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_2}){C_4} + {C_2}}}{{{C_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_2},{r_{\text{H}}}} \right]} \end{array} $$ (33) 当满足rc∈[r2,rH]时，在[1,rc)区间内，冻结壁处于塑性状态，而在[rc,rH]区间内，冻结壁处于弹性状态。计算模型如图4c所示[12]。
冻结壁在Ⅰ区塑性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅰ}})}}} = D{l_1}r + P{m_1} + {K_1}{r^{\lambda - 1}}} \\ & {\sigma _\theta ^{{\text{p({\text{Ⅰ}})}}} = \lambda \sigma _r^{{\text{p(I)}}} + \omega {c_1}(r)}\end{aligned} \right.\begin{array}{*{20}{c}} {r \in \left[ {1,{r_1}} \right)} \end{array} $$ (34) 冻结壁在Ⅱ区塑性区应力场的应力计算公式：
$$ \left\{\begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅱ}})}}} = D{l_2}r + P{m_2} + {K_2}{r^{\lambda - 1}}} \\ & {\sigma _\theta ^{{\text{p({\text{Ⅱ}})}}} = \lambda \sigma _r^{{p} ({\text{Ⅱ}})} + \omega {c_2}(r)} \end{aligned} \right.\begin{array}{*{20}{c}} {r \in \left[ {{r_1},{r_2}} \right)} \end{array} $$ (35) 冻结壁在Ⅲ区塑性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅲ}})}}} = D{l_3}r + P{m_3} + {K_3}{r^{\lambda - 1}}} \\ & {\sigma _\theta ^{{\text{p({\text{Ⅲ}})}}} = \lambda \sigma _r^{{\text{p({\text{Ⅲ}})}}} + \omega {c_3}(r)} \end{aligned} \right.\begin{array}{*{20}{c}} {r \in \left[ {{r_2},{r_{\text{c}}}} \right)} \end{array} $$ (36) 冻结壁在Ⅲ区弹性区应力场的应力计算公式：
$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_{\text{c}}}){I_3} + {I_2}}}{{{I_1}}}} \\ & {\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}} = \frac{{({F_{\text{H}}} - {F_{\text{c}}}){I_4} + {I_2}}}{{{I_1}}}} \end{aligned} \right.\begin{array}{*{20}{c}} {}&{r \in \left[ {{r_{\text{c}}},{r_{\text{H}}}} \right]} \end{array} $$ (37) 式(37)中：
$$ {I_1} = {a_3}\left( {\frac{1}{{{r_{\text{c}}}}} - \frac{1}{{{r_{\text{H}}}}}} \right) + \frac{{{b_3}}}{2}\left( {\frac{1}{{r_{\text{c}}^2}} - \frac{1}{{r_{\text{H}}^2}}} \right); $$ $$ {I_2} = {a_3}\left( {\frac{{{F_{\text{H}}}}}{{{r_{\text{c}}}}} - \frac{{{F_{\text{c}}}}}{{{r_{\text{H}}}}}} \right) + \frac{{{b_3}}}{2}\left( {\frac{{{F_{\text{H}}}}}{{r_{\text{c}}^{\text{2}}}} - \frac{{{F_{\text{c}}}}}{{r_{\text{H}}^2}}} \right) ; $$ $$ {I_3} = - \frac{{{a_3}}}{r} - \frac{{{b_3}}}{{2{r^2}}};{I_4} = \frac{{{b_3}}}{{2{r^2}}}。 $$ 3.4 冻结壁塑性相对半径rc的求解
在r＝rc处，冻结壁的径向应力与环向应力应满足统一强度屈服准则以及应力连续条件，即：
$$ \sigma _\theta ^i = \lambda \cdot \sigma _r^i + \omega \cdot c(r) $$ (38) $$ \sigma _r^{\text{p}} = \sigma _r^{\text{e}} $$ (39) 当
${r_{\text{c}}} \in \left[ {1,{r_1}} \right)$ 时：$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅰ}})}}} = D{l_1}{r_{\text{c}}} + P{m_1} + {K_1}r_{\text{c}}^{\lambda - {\text{1}}}} \\ & \sigma _r^{{\text{e({\text{Ⅰ}})}}} = \\ &\dfrac{{\left({F_1} - {F_{\text{c}}}\right)\left( - \dfrac{{{a_1}}}{{{r_{\text{c}}}}} - \dfrac{{{b_1}}}{{2r_{\text{c}}^{\text{2}}}}\right) + {a_1}\left(\dfrac{{{F_1}}}{{{r_{\text{c}}}}} - \dfrac{{{F_{\text{c}}}}}{{{r_1}}}\right) + \dfrac{{{b_1}}}{2}\left(\dfrac{{{F_1}}}{{r_{\text{c}}^{\text{2}}}} - \dfrac{{{F_{\text{c}}}}}{{r_1^2}}\right)}}{{{a_1}\left(\dfrac{1}{{{r_{\text{c}}}}} - \dfrac{1}{{{r_1}}}\right) + \dfrac{{{b_1}}}{2}\left(\dfrac{1}{{r_{\text{c}}^{\text{2}}}} - \dfrac{1}{{r_1^2}}\right)}} \end{aligned} \right. $$ (40) 当
${r_{\text{c}}} \in \left[ {{r_1},{r_2}} \right)$ 时：$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅱ}})}} }= D{l_2}{r_{\text{c}}} + P{m_2} + {K_2}r_{\text{c}}^{\lambda {{ - 1}}}} \\ & \sigma _r^{{\text{e({\text{Ⅱ}})}}} = \\ &\dfrac{{\left({F_2} - {F_{\text{c}}}\right)\left( - \dfrac{{{a_2}}}{{{r_{\text{c}}}}} - \dfrac{{{b_2}}}{{2r_{\text{c}}^{\text{2}}}}\right) + {a_2}\left(\dfrac{{{F_2}}}{{{r_{\text{c}}}}} - \dfrac{{{F_{\text{c}}}}}{{{r_2}}}\right) + \dfrac{{{b_2}}}{2}\left(\dfrac{{{F_2}}}{{r_{\text{c}}^{\text{2}}}} - \dfrac{{{F_{\text{c}}}}}{{r_2^2}}\right)}}{{{a_2}\left(\dfrac{1}{{{r_{\text{c}}}}} - \dfrac{1}{{{r_2}}}\right) + \dfrac{{{b_2}}}{2}\left(\dfrac{1}{{r_{\text{c}}^{\text{2}}}} - \dfrac{1}{{r_2^2}}\right)}} \end{aligned} \right. $$ (41) 当
$ {r_{\text{c}}} \in \left[ {{r_2},{r_{\text{H}}}} \right] $ 时：$$ \left\{ \begin{aligned} & {\sigma _r^{{\text{p({\text{Ⅲ}})}}} = D{l_3}{r_{\text{c}}} + P{m_3} + {K_3}r_{\text{c}}^{\lambda - {\text{1}}}} \\ & \sigma _r^{{\text{e({\text{Ⅲ}})}}} = \\ &\dfrac{{\left({F_{\text{H}}} - {F_{\text{c}}}\right)\left( - \dfrac{{{a_3}}}{{{r_{\text{c}}}}} - \dfrac{{{b_3}}}{{2r_{\text{c}}^{\text{2}}}}\right) + {a_3}\left(\dfrac{{{F_{\text{H}}}}}{{{r_{\text{c}}}}} - \dfrac{{{F_{\text{c}}}}}{{{r_{\text{H}}}}}\right) + \dfrac{{{b_3}}}{2}\left(\dfrac{{{F_{\text{H}}}}}{{r_{\text{c}}^2}} - \dfrac{{{F_{\text{c}}}}}{{r_{\text{H}}^2}}\right)}}{{{a_3}\left(\dfrac{1}{{{r_{\text{c}}}}} - \dfrac{1}{{{r_{\text{H}}}}}\right) + \dfrac{{{b_3}}}{2}\left(\dfrac{1}{{r_{\text{c}}^{\text{2}}}} - \dfrac{1}{{r_{\text{H}}^2}}\right)}} \end{aligned} \right. $$ (42) 分别联立式(40)、式(41)和式(42)，可以求得冻结壁的塑性区相对半径rc。
4 工程算例
4.1 计算参数
本算例中内外排冻结管尺寸(直径×厚度)均为159 mm×6 mm。温度参数：t0＝−3 ℃，t1＝−25 ℃，t2＝−45 ℃，tn＝−8 ℃。冻结孔设计参数和冻土力学参数分别见表2、表3[21]。
表 2 盐水−二氧化碳复合冷媒冻结孔设计参数Table 2. Design parameters of brine-CO2 freezing hole类别 参数 数值 井筒 井筒内径/m 3 开挖荒径/m 5 内侧冻结管
(低温盐水)布置圈直径/m 14 数量/根 38 管间距/m 1.157 管壁温度/ ℃ −32 外侧冻结管
(二氧化碳)布置圈直径/m 20 数量/根 38 管间距/m 1.653 管壁温度/ ℃ −50 表 3 冻土力学参数Table 3. Mechanical parameters of frozen soil弹性模量
E*/MPa黏聚力
c*/MPa内摩擦角
γ/(°)泊松比
μ−8.074t+35.41 −0.174t+0.803 3.5 0.27 注：*表示根据文献[32]的研究，弹性模量E和黏聚力c与温度之间存在线性关系，这一关系是通过冻土的实验数据拟合得出的，只适用于被冻土体。 将上面的数据分别代入式(3)和式(4)，计算得到非均质冻结壁的材料参数的表达式。
$$ E(r) = \left\{ {\begin{array}{*{20}{l}} {343.145r - 243.143}&{r \in \left[ {1,1.4} \right)}&{} \\ {269.13\dot 3r - 139.527}&{r \in \left[ {1.4,2} \right)}&{} \\ { - 423.885r + 1246.51}&{r \in \left[ {2,2.8} \right]}&{} \end{array}} \right. $$ (43) $$ c(r) = \left\{ {\begin{array}{*{20}{l}} {7.395r - 5.200}&{}&{r \in \left[ {1,1.4} \right)} \\ {5.800r - 2.967}&{}&{r \in \left[ {1.4,2} \right)} \\ { - 9.135r + 26.903}&{}&{r \in \left[ {2,2.8} \right]} \end{array}} \right. $$ (44) 将求得的E(r)、c(r)等参数代入本文公式中，得到该状态下的应力解。在相同条件下，根据表3的数据，计算均质冻结壁的材料参数。将平均温度[30]t＝−23 ℃代入表3，得到弹性模量和黏聚力分别为221.112、4.805 MPa。根据文献[35]的应力公式，求得均质冻结壁的应力计算结果见表4[21]。
表 4 冻结壁应力计算结果Table 4. Stress calculation results of frozen wallMPa 计算准则 弹性极限状态(rc=1) 弹塑性状态(rc=1.2) 弹塑性状态(rc=1.7) 弹塑性状态(rc=2.4) 塑性极限状态(rc=2.8) 均质 非均质 均质 非均质 均质 非均质 均质 非均质 均质 非均质 M-C强度准则 4.456 4.376 6.155 6.153 9.068 9.761 10.991 11.937 11.255 12.170 D-P强度准则 4.453 4.373 6.150 6.149 9.061 9.754 10.983 11.929 11.247 12.161 广义Tresca准则 5.195 5.101 7.198 7.193 10.654 11.447 12.952 14.042 13.269 14.336 双剪统一强度准则 5.942 5.835 8.261 8.251 12.285 13.174 14.980 16.212 15.354 16.575 4.2 冻结壁极限承载力分析
根据4种不同的强度准则，计算得到了2种冻结壁的弹性、塑性极限承载力的计算结果，如图5所示。分析表明，在不同的屈服准则下，冻结壁承受的外荷载存在一定差异。计算结果由小到大依次是D-P准则、M-C准则、广义Tresca准则以及双剪统一强度准则。其中，基于M-C准则和D-P准则的计算结果较为接近，而基于广义Tresca准则以及双剪统一强度准则得出的计算结果依次增大。
![]()
图 5 冻结壁弹性、塑性极限承载力计算结果Figure 5. Calculation results of elastic and plastic ultimate bearing capacity of frozen wall由图5可知，基于不同的屈服准则，在弹性极限状态下(rc=1)，均质冻结壁的弹性极限承载力略高于非均质冻结壁的相应承载力；而在塑性极限状态下(rc=2.8)，均质冻结壁的塑性极限承载力都明显小于非均质冻结壁的塑性极限承载力。通过计算可以得出，考虑非均质特性后，冻结壁的弹性极限承载力降低1.8%，而塑性极限承载力提高8.1%。
考虑冻结壁的非均匀特性后，冻结壁的承载力发生变化的原因如下：弹性极限状态是冻结壁即将进入塑性的状态，在计算过程中认为冻结壁的内缘刚刚进入塑性，对应的荷载就是冻结壁的弹性极限承载力，在考虑冻结壁的非均质特性之后，冻结壁的弹性模量以及黏聚力在冻结壁内缘处小于均质冻结壁，因此对应的弹性极限承载力要小于均质冻结壁。而塑性极限状态表示冻结壁全部区域进入了塑性，此时材料的应力与应变关系变为非线性，材料开始发生不可逆的塑性变形。均质冻结壁在整体达到屈服点时，承载力会迅速下降。相比之下，非均质冻结壁由于不同区域的屈服极限不同，冻结壁中黏聚力的最高值大于均质冻结壁，抵御破坏的能力更强，因此非均质冻结壁对应的塑性极限承载力更大。
4.3 冻结壁应力分布规律分析
选取5种状态下非均质冻结壁与均质冻结壁的应力分布比较：弹性极限状态、弹塑性状态I(塑性区半径为1.2)、弹塑性状态Ⅱ(塑性区半径为1.7)、弹塑性状态Ⅲ(塑性区半径为2.4)、塑性极限状态。
4.3.1 冻结壁弹性极限状态下应力分布
根据4种不同的强度准则，得到2种冻结壁在弹性极限状态下应力分布曲线如图6所示。当相对半径r＝1时，均质冻结壁的径向应力与非均质冻结壁的径向应力之间相差不大，之后随着相对半径的增大，两者之间的差值逐渐增大，而当r>1.4之后，两者之间的差值逐渐减小，并逐渐趋于稳定，但在冻结壁的外边缘，非均质冻结壁的径向应力仍然略小于均质冻结壁的径向应力。而对于环向应力，均质冻结壁的环向应力随着相对半径r的增大而逐渐减小，在冻结壁内壁面(r＝1)处，径向应力达到峰值。对于非均质冻结壁，当1≤r<1.4时，随着相对半径r的增加，环向应力也随之增加，其应力曲线呈现近似抛物线增长，当1.4≤r<2时，曲线随相对半径r的增大呈线性缓慢上升，当r≥2时，环向应力逐渐下降，此时，非均质冻结壁的环向应力最大值出现在冻结壁分区线r＝2的位置。
![]()
图 6 弹性极限状态下冻结壁应力分布曲线Figure 6. Stress distribution curve of frozen wall in elastic limit state4.3.2 冻结壁弹塑性状态下应力分布
2种冻结壁在弹塑性状态Ⅰ应力分布曲线如图7所示。分析结果表明：在两种状态下，冻结壁的径向应力随着相对半径r的增大而增大，之后两者的数值基本相同并趋于稳定。而对于均质冻结壁而言，其环向应力表现出显著的转折点，先是逐渐增加，随后呈现缓慢下降的趋势，且在弹塑性边界达到应力峰值；对于非均质冻结壁，当1≤r<1.2时，环向应力随着相对半径r的增大呈线性急剧增长趋势，当1.2≤r<2时，曲线随r的增大呈现线性缓慢上升趋势，当r≥2时，环向应力逐渐下降。此时，非均质冻结壁的环向应力最大值出现在冻结壁分区线r＝2的位置。
![]()
图 7 弹塑性状态Ⅰ下冻结壁应力分布曲线Figure 7. Stress distribution curve of frozen wall under plastic-elastic state Ⅰ2种冻结壁在弹塑性状态Ⅱ应力分布曲线如图8所示。分析表明，在均质和非均质两种状态下，冻结壁的径向应力均随着相对半径r的增加而上升，但在到达冻结厚度中部(即r＝1.9处)之前，非均质冻结壁的径向应力略小于均质冻结壁的径向应力，而在该点之后，略大于均质冻结壁。对于均质冻结壁，其环向应力特征在于存在一个显著的转折点，表现出先上升，随后逐渐降低的趋势，且在弹塑性界面达到峰值应力；对于非均质冻结壁，环向应力整体呈近似梯形分布。当1≤r<1.7时，环向应力随着相对半径r的增大呈线性急剧增长，当1.7≤r<2时，曲线随r的增大呈现线性缓慢上升，当r≥2时，环向应力急剧下降。此时，非均质冻结壁的环向应力最大值出现在冻结壁分区线r＝2的位置。
![]()
图 8 弹塑性状态Ⅱ下冻结壁应力分布曲线Figure 8. Stress distribution curve of frozen wall under plastic-elastic state Ⅱ2种冻结壁在弹塑性状态Ⅲ应力分布曲线如图9所示。通过分析可以发现，在两种状态下，随着相对半径r的增加，冻结壁的径向应力普遍增大。然而，在到达冻结壁中部(r＝1.9处)之前，非均质冻结壁与均质冻结壁的径向应力差值逐渐减少，该点之后两者之间的差值逐渐增大。对于均质冻结壁，环向应力曲线展现出显著的拐点，先是随着相对半径的增加而增加，达到峰值后逐渐呈现出缓慢下降的趋势，在弹塑性交界区域达到应力峰值；对于非均质冻结壁，环向应力呈现先上升后下降再下降的趋势，且第二段(2≤r<2.4)下降速度比第一段(2.4≤r<2.8)下降得更快。此时，非均质冻结壁的环向应力最大值出现在冻结壁分区线r＝2的位置。
![]()
图 9 弹塑性状态Ⅲ下冻结壁应力分布曲线Figure 9. Stress distribution curve of frozen wall under plastic-elastic state Ⅲ4.3.3 冻结壁塑性极限状态下应力分布
2种冻结壁在塑性极限状态下应力分布曲线如图10所示。在两种状态下，径向应力总体上呈现出随着相对半径增加而增加的趋势。然而，与均质状态相比，非均质状态下的径向应力增长更为迅速。均质冻结壁的环向应力随着相对半径r的增大而逐渐增大，在这种情况下，均质冻结壁的环向应力最大值出现在冻结壁外壁r＝2.8处；而对于非均质冻结壁，环向应力整体呈近似三角形分布。当1≤r<2时，环向应力随着相对半径r的增加而增加，当2≤r<2.8时，环向应力曲线呈现出近似线性的下降趋势。非均质冻结壁的环向应力的峰值出现在冻结壁的分区线，即r＝2的位置。
![]()
图 10 塑性极限状态下冻结壁应力分布曲线Figure 10. Stress distribution curve of frozen wall in plastic limit state4.4 冻结壁塑性区相对半径与承受外荷载之间的关系
随着外荷载的增大，冻结壁逐渐由弹性转变为弹塑性。图11—图12分别展示了冻结壁的塑性区半径随外荷载变化趋势以及均质与非均质冻结壁的计算误差。
![]()
图 11 冻结壁塑性区相对半径和承受外荷载之间的关系Figure 11. Relationship between the relative radius of the plastic zone of the frozen wall and the external load![]()
图 12 均质冻结壁与非均质冻结壁计算结果误差变化Figure 12. Error variation law of calculation results between homogeneous and heterogeneous frozen walls通过对比可以发现：基于M-C准则和D-P准则计算结果较为接近，而基于双剪强度准则计算出的结果最大。冻结壁的承载力随着塑性区半径的增加而增加，但增长速率逐渐降低。当塑性半径1≤rc<1.2时，非均质冻结壁的承载力略小于均质冻结壁的承载力；当塑性半径1.2<rc≤2.8时，非均质冻结壁的承载力明显大于均质冻结壁，此时，冻结壁外荷载误差变化曲线整体呈现先快速上升后轻微下降的趋势，在rc＝2.4处2种计算结果的差值最大。
在均质冻结壁的设计理论中，由于未能全面考虑冻结壁在各区域物理力学性质的差异，会导致对弹塑性状态下冻结壁承载力的计算值偏小，采用该设计理论会让冻结壁设计过于保守，从而会导致大量冷量的浪费，进而会增加后期开挖的难度。相对而言，非均匀冻结壁计算理论充分考虑了冻土强度随温度及区域的变化，采用该计算理论得出的计算结果更加接近实际情况。
在不同状态的冻结壁应力分布规律中，非均质与均质冻结壁的环向应力沿相对半径的分布规律存在显著差异，而径向应力的差别较小，产生该现象的原因如下：本研究中已假设外荷载沿着径向均匀分布在冻结壁的外缘，外荷载使冻结壁整体产生向中心收缩的趋势，径向应力抵御冻结壁的径向收缩，环向应力抵御冻结壁的环向收缩。均质冻结壁的力学参数(弹性模量、黏聚力)处处相等，而非均质冻结壁的力学参数沿着半径方向成三段一次函数变化。在冻结壁内部半径不同的圆环上，力学参数存在较大差异，导致非均质冻结壁的环向应力与均质冻结壁相比存在较大差异。而非均质冻结壁与均质冻结壁的平均温度相等，冻结壁径向整体强度差别较小，外荷载是作用在冻结壁外缘，沿着径向传递到冻结壁内部的，因此，两种计算理论得出的径向应力的变化趋势差别较小。
进一步分析可以发现，对冻结壁的环向应力产生直接影响的是弹性模量和黏聚力，上述2个力学参数与冻结温度成近似线性关系，本研究中将冻结壁特征截面处的温度场等效为三段一次函数形式，在不同区域温度曲线交界处存在转折点，对应的力学参数也会发生变化，因此在区域分界线处，冻结壁的环向应力曲线出现转折点。
5 结 论
(1)针对“盐水−二氧化碳”多冷媒双排管冻结壁，选取距离1/4管距处冻结壁作为特征截面，将该特征截面上的温度分布曲线等效成三段一次函数形式，并将冻结壁视为随温度函数变化的非均质材料，分别基于M-C强度准则、D-P强度准则、广义Tresca准则以及双剪统一强度准则，推导得到该类冻结壁的应力计算公式。
(2)基于均质冻结壁计算理论，弹性极限状态的冻结壁环向应力最大值出现在冻结壁内侧，弹塑性状态的冻结壁的环向应力最大值出现在弹塑性分界面处，塑性极限状态的冻结壁的环向应力最大值出现在冻结壁最外侧；而基于非均质冻结壁计算理论，冻结壁环向应力最大值始终出现在冻结壁分区界线(r＝2)处。
(3)以M-C准则为例，非均质冻结壁与均质冻结壁的弹性极限荷载计算结果分别为4.38、4.46 MPa，塑性极限承载力的计算结果分别为12.17、11.26 MPa。说明考虑冻结壁的非均质特性后，弹性极限承载力降低，而塑性极限承载力提高。
(4)在弹塑性状态下，对应相同塑性区相对半径，非均质冻结壁具有更高的承载力，且这种现象随着塑性区相对半径的增大愈发明显。盐水−二氧化碳多冷媒双排管冻结壁设计时，必须要考虑冻结壁的非均匀特性，否则会让设计过于保守，进而会导致大量冷量的浪费。
(5)本文推导的多冷媒非均质冻结壁弹塑性应力解析解可适用于两种冷媒联合冻结形成的冻结壁的力学特性分析，不仅可用于“盐水−二氧化碳”复合冻结壁，也可用于“盐水−液氮”等复合冻结壁。本公式应用于其他类型的复合冻结壁时，需调整对应特征截面的温度数值以及冻土的基本物理力学参数。
符号注释：
a为弹性模量随温度变化的一次项系数，MPa；ai为弹性模量在Ⅰ、Ⅱ、Ⅲ区随温度变化的统一一次项系数，MPa；a1、a2、a3分别为弹性模量在冻结壁Ⅰ、Ⅱ、Ⅲ区随温度变化的一次项系数，MPa；A1、A2、A3为冻结壁处于弹性状态下在Ⅰ区应力表达式中的待定系数，MPa；b为弹性模量随温度变化的常数项系数，MPa；bi为弹性模量在Ⅰ、Ⅱ、Ⅲ区随温度变化的统一常数项系数，MPa；b1、b2、b3分别为弹性模量在冻结壁Ⅰ、Ⅱ、Ⅲ区随温度变化的常数项系数，MPa；B1、B2、B3、B4为冻结壁在弹性区Ⅱ区应力表达式中的待定系数，MPa；c为黏聚力，MPa；c(R)为整个冻结壁区域Ⅰ、Ⅱ、Ⅲ随半径变化的黏聚力，MPa；c(r)为整个冻结壁区域Ⅰ、Ⅱ、Ⅲ随相对半径变化的黏聚力，MPa；c0为初始黏聚力，MPa；c1(r)、c2(r)、c3(r)分别为冻结壁在I、Ⅱ、Ⅲ区随相对半径变化的黏聚力，MPa；C1、C2、C3、C4为冻结壁在弹性区Ⅲ区应力表达式中的待定系数，MPa；D为冻结壁在塑性区应力表达式中的一次项待定系数；E为弹性模量，MPa；E(R)为整个冻结壁区域Ⅰ、Ⅱ、Ⅲ随半径变化的弹性模量，MPa；E(r)为整个冻结壁区域Ⅰ、Ⅱ、Ⅲ随相对半径变化的弹性模量，MPa；E0为初始弹性模量，MPa；F1为冻结壁Ⅰ区和Ⅱ区分界处的相互作用力，MPa；F2为冻结壁Ⅱ区和Ⅲ区分界处的相互作用力，MPa；FH为外部荷载，MPa；Fc为弹性区与塑性区分界处的相互作用力，MPa；G1、G2、G3、G4为冻结壁处于弹塑性状态下在Ⅰ区应力表达式中的待定系数，MPa；h1为内排冻结管之间的管间距，m；h2为外排冻结管之间的管间距，m；H1、H2、H3、H4为冻结壁处于弹塑性状态下在Ⅱ区应力表达式中的待定系数，MPa；I1、I2、I3、I4为冻结壁处于弹塑性状态下在Ⅲ区应力表达式中的待定系数，MPa；K1、K2、K3分别为冻结壁在Ⅰ、Ⅱ、Ⅲ塑性区应力表达式中的(λ−1)项待定系数，MPa；l为黏聚力随温度变化的一次项系数，MPa；l1、l2、l3分别为黏聚力在冻结壁Ⅰ、Ⅱ、Ⅲ区随温度变化的一次项系数，MPa；L为内、外侧冻结管之间的排距，m；m为黏聚力随温度变化的常数项系数，MPa；m1、m2、m3分别为黏聚力在冻结壁Ⅰ、Ⅱ、Ⅲ区随温度变化的常数项系数，MPa；Mi、Ni为应力函数通解的待定系数，MPa；P为冻结壁在塑性区应力表达式中的常数项待定系数；相对半径r＝R/RB为冻结壁截面内的任意一点，其中，r1＝R1/RB，r2＝R2/RB，rH＝RH/RB，rc＝Rc/RB；R为冻结壁截面上任意一点的半径，m；R0为冻结壁内径，m；RB为开挖荒径，m；R1为内侧冻结管中心至井壁中心的距离，m；R2为外侧冻结管中心至井壁中心的距离，m；RH为冻结壁外径，m；Rc为冻结壁塑性区半径，m；t为冻结壁的温度， ℃；t0为冻土的结冰温度， ℃；t1为距离盐水冻结温度场主面1/4管距处冻结壁的最低温度， ℃；t2为距离二氧化碳冻结温度场主面1/4管距处冻结壁的最低温度， ℃；tn为冻土开挖后内壁温度， ℃；
$u_r^{\mathrm{e}} $ 为冻结壁弹性区径向位移；$u_1^{-} $ 为冻结壁r1处的内侧位移；$u_1^{+} $ 为冻结壁r1处的外侧位移；$u_2^{-} $ 为冻结壁r2处的内侧位移；$u_2^{+} $ 为冻结壁r2处的外侧位移；γ为内摩擦角，(°)；δ为内侧冻结管到内侧冻结壁边界的距离，m；$\varepsilon _r^{\mathrm{e}} $ 为冻结壁弹性区径向应变；$\varepsilon _\theta ^{\mathrm{e}} $ 为冻结壁弹性区环向应变；$\varepsilon _1^{-} $ 为冻结壁r1处的内侧应变；$\varepsilon _1 ^{+} $ 为冻结壁r1处的外侧应变；$\varepsilon _2^{-} $ 为冻结壁r2处的内侧应变；$\varepsilon _2 ^{+} $ 为冻结壁r2处的外侧应变；ζ为外侧冻结管到外侧冻结壁边界的距离，m；λ为小主应力项的屈服强度系数；μ为泊松比；$\xi $ 为统一强度准则参数；σr、σθ分别为冻结壁的径向、环向应力，MPa；$\sigma _r^{{\text{e}}} $ 、$\sigma _\theta^{{\text{e}}} $ 分别为冻结壁弹性区的径向、环向应力，MPa；$\sigma _r^{{\text{p}}} $ 、$\sigma _\theta^{{\text{p}}} $ 分别为冻结壁塑性区的径向、环向应力，MPa；$\sigma _r^{{\text{e({\text{Ⅰ}})}}} $ 、$\sigma _r^{{\text{e({\text{Ⅱ}})}}} $ 、$\sigma _r^{{\text{e({\text{Ⅲ}})}}} $ 分别为冻结壁在Ⅰ、Ⅱ、Ⅲ弹性区的径向应力，MPa；$\sigma _\theta ^{{\text{e({\text{Ⅰ}})}}} $ 、$\sigma _\theta ^{{\text{e({\text{Ⅱ}})}}} $ 、$\sigma _\theta ^{{\text{e({\text{Ⅲ}})}}} $ 分别为冻结壁在Ⅰ、Ⅱ、Ⅲ弹性区的环向应力，MPa；$\sigma _\theta ^{i} $ 为大主应力，MPa；$\sigma _r ^{i} $ 为小主应力，MPa；$\sigma _r^{{\text{p({\text{Ⅰ}})}}} $ 、$\sigma _r^{{\text{p({\text{Ⅱ}})}}} $ 、$\sigma _r^{{\text{p({\text{Ⅲ}})}}} $ 分别为冻结壁在Ⅰ、Ⅱ、Ⅲ塑性区的径向应力，MPa；$\sigma _\theta ^{{\text{p({\text{Ⅰ}})}}} $ 、$\sigma _\theta ^{{\text{p({\text{Ⅱ}})}}} $ 、$\sigma _\theta ^{{\text{p({\text{Ⅲ}})}}} $ 分别为冻结壁在Ⅰ、Ⅱ、Ⅲ塑性区的环向应力，MPa；ψ为应力函数，MPa；ω为黏聚力项的屈服强度系数。
下载:
计量
- 文章访问数: 25
- HTML全文浏览量: 11
- PDF下载量: 1
