(安徽理工大学 土木建筑学院，安徽 淮南 232001)

1 研究背景

冻结法凿井是采用人工制冷的方法，在井筒周围的含水岩土层中形成封闭的冻结壁，以此来抵抗水压力，隔绝地下水与井筒的联系，确保井筒掘砌安全的一种特殊工法[1-2]。随着冻结深度的增加，冻结壁的外荷载也逐渐增加，目前我国冻结法施工穿过冲积层的最大厚度达750 m，这对冻结壁的强度和稳定性都提出了更高的要求[3-4]，单排管形成的冻结壁已经无法满足实际的施工要求，因此必须要采用多排管冻结的施工方案。在多排管冻结研究领域，汪仁和和曹荣斌[5]运用ANSYS有限元计算程序对双排管冻结壁的形成及其表征变化进行了详细的计算分析，得出了双排管冻结下冻结时间缩短、冻结效率提高、冻结壁平均温度下降等特性，并探讨了双排管冻结下冻结壁温度的简化计算方法。肖朝昀等[6]从冻结深度、厚度方向上分析多排管局部冻结排内和排外的温度发展特性，并分析计算出了积极冻结期排内冻结壁交圈时间以及发展速度。胡向东等[7-11]对多排管冻结温度场进行了系统研究，并提出了直线双排管冻结壁平均温度的等效梯形、等效抛物弓模型以及三排管冻结的梯形-抛物弓等效温度场模型，并验证了其精度。胡向东等[12]以及王彬等[13-14]将双排管冻结形成的冻结壁等效成由材料性质沿着径向呈梯形分布的功能梯度材料(FGM)构成的无限长厚壁圆筒，分别结合摩尔-库伦屈服准则以及Druker-Prager屈服准则对该类冻结壁进行弹塑性分析，并对冻结壁设计理论提出了具有参考价值的建议。

近些年在深厚冲击层冻结法凿井过程中，为了应对较大流速的地下水对井壁掘砌施工带来的危害，三排管冻结方案得到了越来越广泛的运用。三排管冻结壁的力学特性分析是决定该类冻结方案设计以及进一步推广的主要因素，传统的三排管冻结壁计算理论将冻结壁视为均质材料，该计算方法忽略了冻结壁在不同冻结区间的冻土的物理力学参数的差异性，因此得出的结果与冻结壁实际的受力情况存在一定的差异。

本文将基于三排管冻结的梯形-抛物弓等效温度场模型，根据冻结壁的弹性模量和黏聚力与温度之间的一次函数关系，并结合摩尔-库伦准则，对三排管非均质冻结壁力学特性进行弹塑性分析。

2 温度场模型

从三排管冻结壁中选取一段连续的区间作为研究对象，并将其简化为长方形区间,如图1所示。其中ζ为冻结管到冻结壁边界的距离，h为管间距，L为排距，冻结壁内缘一小部分被冻土体以及冻结壁包围的未冻土体在冻结完成后会被一起一次性移除，从而为井壁的施工提供一个稳定的空间环境，因此冻结壁的实际范围为被冻土体的厚度减去冻结壁拟开挖控制厚度，如图1所示。

图1 三排管冻结壁局部示意图Fig.1 Local sketch of triple-row-piped frozen soil wall

三排管冻结温度场特征截面处的温度分布如图2所示。其中：T0=-3 ℃为冻土的结冰温度；T1=-8 ℃为冻土开挖后内壁的温度；Tk为冻结管外壁温度；Ta为距离冻结温度场主面1/4管距处的温度场的最低温度；Tb为冻结管在该特征截面对应位置处的温度；RB为冻结壁内壁至开挖部分中轴线的距离(简称内半径)；R1—R3分别为第1排—第3排冻结管至开挖部分中轴线的距离；RH为冻结壁外壁至开挖部分中轴线的距离(简称外半径)。

图2 三排管冻结温度场特征截面温度分布示意图Fig.2 Schematic diagram of temperature distributionof characteristic section of temperature field oftriple-row piped freezing

已有的研究成果表明，三排管冻结温度场具有明显的周期性，并且距离主面1/4管距处的截面可以比较全面而准确地反映整个冻结壁的温度分布规律[7-11]，因此选取该截面为特征截面来对冻结壁进行力学特性分析。

通过梯形-抛物弓的温度场模型来代替距离主面1/4管距处截面的温度场[10]，如图3所示。通过图中5个温度特征点可以确定温度场模型中2段一次函数曲线以及1段二次函数曲线的表达式。为了便于计算，将冻结壁分为区间Ⅰ，Ⅱ，Ⅲ。

图3 等效温度场分布示意图Fig.3 Diagram of equivalent temperaturefield distribution

3 力学模型

通过分析课题组多年来总结得到的不同矿区400 m深度处的冻土的力学参数试验结果，可以发现，冻土的弹性模量E、黏聚力c与温度之间存在近似的一次函数关系，如图4所示。

图4 冻土的弹性模量、黏聚力随温度变化曲线Fig.4 Curves of elastic modulus versus temperatureand cohesive force versus temperature

图5 冻结壁受均布荷载示意图Fig.5 Schematic diagram offrozen wall under uniformload

因此在计算中假设冻结壁的弹性模量以及黏聚力与温度呈线性关系；而泊松比与内摩擦角受温度影响较小，故忽略两者随温度的变化[15-16]。

冻结壁在均布荷载作用下的力学计算模型如图5所示。

图5中：p为外荷载；σp为弹塑性区分界处塑性区与弹性区间的相互作用力；R为冻结壁的任一点半径；Rp为塑性区半径。冻结壁内任一点相对半径r=R/RB，则半径R1和R2对应的相对半径分别为r1=R1/RB，r2=R2/RB。塑性区相对半径ρ=Rp/RB，外相对半径rH=RH/RB。

根据弹性模量E、黏聚力c与温度之间的一次函数关系，基于梯形-抛物弓的温度场模型的三排管非均质冻结壁的E(r)与c(r)的表达式为：

(1)

(2)

式中a1,a2,a3,b1,b2,b3,c2,l2,m1,m2,m3，n1,n2,n3为待定常数，由冻土力学参数与温度的关系决定。

4 应力求解

4.1 弹性状态下应力求解

假设冻结壁处于平面应变状态，将冻结壁等效为周围受均布荷载的厚壁圆筒[17]。

设应力与应力函数φ之间的关系为

(3)

式中：σr为径向应力；σθ为环向应力。

结合弹性力学中平衡方程、几何方程以及物理方程[17]有

(4)

将式(1)中E(r)代入式(4)，进一步化简得

(5)

其中：

解变系数微分方程得应力函数φ的通解为

(6)

式中i1，i2，i3为待定常数。

式(6)的边界条件为：

代入边界条件，得各分区(r∈[1,r1)，r∈[r1,r2)，r∈[r2,rH])的应力表达式分别为：

(7)

(8)

(9)

其中：

B3=σr1a2lnr2-σr2a2lnr1，

在平面轴对称应变问题[15]中，有：

(10)

(11)

(12)

式中：ur为径向位移；M,N为待定参量。

联立式(10)—式(12)可得

(13)

(14)

(15)

式中上标(1)—上标(3)分别表示冻结壁所对应的区域Ⅰ，Ⅱ，Ⅲ。

联立式(14)和式(15)可求得σr1与σr2。

4.2 弹塑性状态下应力求解

当冻结壁的外加荷载p超过冻结壁的弹性极限荷载时，冻结壁进入弹塑性状态，沿着冻结壁的径向，由内到外依次为塑性区与弹性区。

4.2.1 塑性区应力场求解

假定冻结壁的塑性区体积不可压缩，则将摩尔-库伦准则对应的公式进行化简[18-21]，即

σθ=Aσr+Bc(r) 。

(16)

平衡方程为

(17)

将式(16)代入式(17)得

(18)

通过求解式(18)得冻结壁塑性区径向应力的通解为

(19)

式中C1,C2,C3为待定常数。

其中：

分别利用r=1,r=r1以及r=r2处的连续条件，计算得出冻结壁塑性区径向应力的特解为

(20)

其中，

依据摩尔-库伦准则可得冻结壁塑性区环向应力解为

(21)

4.2.2 弹性区应力场求解

(22)

其中：

区域Ⅱ和区域Ⅲ的弹性区应力依据式(8)和式(9)进行计算。

当塑性区半径ρ∈[r1,r2)时，区域r∈[1,ρ)上的冻结壁全部处于塑性状态，依据式(20)计算其径向应力，并结合式(21)计算对应的环向应力；在区域r∈[ρ,r2)上弹性区应力表达式为

(23)

其中：

区域Ⅲ依据式(9)进行应力计算。

当塑性区的半径ρ∈[r2,rH]时，区域r∈[1,ρ)上的冻结壁处于塑性状态，依据式(20)与式(21)计算其应力；在区域r∈[ρ,rH]上冻结壁的应力计算公式为

(24)

其中：

4.2.3 塑性区半径的求解

在r=ρ处，径向应力和环向应力满足摩尔-库伦准则以及应力的连续条件，即

(25)

区域Ⅰ(ρ∈[1,r1))、区域Ⅱ(ρ∈[r1,r2))、区域Ⅲ(ρ∈[r2,rH])的ρ可分别由式(26)—式(28)求得，即：

(28)

其中，

H1=2cosφ(m1ρ+n1)·

5 工程算例

淮南某矿区进风井采用人工地层冻结法进行施工，在深度为400 m处为砂质黏土层，冻结壁的内径RB=5 m，对应的井帮温度为T1=-8 ℃，外径RH=13.5 m，对应的温度T0=-3 ℃。特征截面处冻结管G1与G3对应位置处的冻结壁温度Tb=-28 ℃，最低温度Ta=-32 ℃，冻结壁平均温度为-23.2 ℃，泊松比μ=0.3，内摩擦角为3°。通过数值计算可以得到温度场的曲线方程为

(29)

将不同负温条件下测得的该矿区的砂质黏土的力学参数进行线性拟合，可得冻土的弹性模量E和黏聚力c分别为：

E=-8.092T+6.4 ,

(30)

c=-0.215T-0.308 。

(31)

式中弹性模量E和黏聚力c的单位均为MPa。

将温度场的曲线方程代入式(30)和式(31)中可得：

(33)

冻结壁外载p一般按照重液公式进行计算，即

p=0.013h。

(34)

式中h为计算深度(m)。

由推导的三排管非均质冻结壁的应力计算公式，可得到非均质冻结壁的弹性极限承载力(即外加荷载p=3.29 MPa时)与塑性极限承载力(即外加荷载p=10.32 MPa时)作用下的应力分布。根据冻结壁平均温度相等的原则计算出对应均质冻结壁的弹性模量与黏聚力，根据参考文献[22]的推导结果，通过计算可得均质冻结壁的弹性极限荷载p=4.26 MPa，而塑性极限荷载为p=10.34 MPa，可以发现非均质冻结壁的弹性极限承载力比均质冻结壁承载力小29.5%，并分别求得对应的应力分布。另外，也求出当冻结壁承受外荷载分别为5.05,9.10,10.00 MPa时的三排管非均质冻结壁与均质冻结壁的应力分布，如图6所示。

图6 三排管非均质冻结壁与均质冻结壁应力分布Fig.6 Stress distributions of heterogeneous frozen soilwall and homogeneous frozen soil wall in triple-rowpiped freezing

由图6 (a)可发现，当冻结壁处于弹性极限状态时，径向应力随着相对半径r的增大不断增大;由于均质冻结壁的弹性极限承载力大于非均质冻结壁的弹性极限承载力，所以该状态下，均质冻结壁的径向应力与环向应力均大于三排管非均质冻结壁的计算结果;均质冻结壁的环向应力随着相对半径r的增大而减小;非均质冻结壁环向应力在区域Ⅰ随着相对半径r增大而增大，在区域Ⅱ、区域Ⅲ则随着相对半径r的增大而减小。

如图6(b)所示，当处于塑性极限状态时，冻结壁的径向应力随r的增大而增大;当冻结壁相对半径r<1.8时，均质冻结壁的径向应力大于三排管冻结壁的径向应力，当r>1.8时均质冻结壁的径向应力小于三排管冻结壁的径向应力，但是数值相差不大；均质冻结壁的环向应力随r的增大而线性增大；三排管冻结壁的环向应力在区域Ⅰ、区域Ⅱ随r增大而增大，在区域Ⅲ随r的增大而减小；三排管非均质冻结壁中部的环向应力大于均质冻结壁的环向应力，而靠近内缘与外缘部分的环向应力则小于均质冻结壁的环向应力。

当冻结壁承受的外荷载分别为5.05,9.10,10.00 MPa时，三排管非均质冻结壁与均质冻结壁的应力分布分别如图6(c)—图6(e)所示。对应的非均质冻结壁的塑性区相对半径分别为1.2,1.8,2.5，均质冻结壁的塑性区相对半径则分别为1.12,1.85,2.23。冻结壁的径向应力随着r的增大而增大，且三排管非均质冻结壁与均质冻结壁之间计算结果基本相同；均质冻结壁的环向应力最大值出现在弹塑性分界线上，三排管冻结壁的环向应力的最大值则需分情况讨论：当塑性区相对半径ρ≤1.4 时，环向应力的最大值出现在r=1.4处(如图6(c))；当塑性区相对半径1.4<ρ<2.1时，环向应力的最大值出现在弹塑性分界线处(如图6(d))；当塑性区相对半径ρ≥2.1时，环向应力的最大值出现在r=2.1处(如图6(e))。

冻结壁的承载力与塑性区相对半径的关系如图7所示。由图7可见，三排管非均质冻结壁与均质冻结壁的承载力都随着塑性区相对半径的增大而增大，并且增大的速度逐渐变慢；对应相同的塑性区相对半径，在区域Ⅰ(1

图7 冻结壁承受外荷载与塑性区相对半径之间的关系Fig.7 Relationship between external load of frozenwall and relative radius of plastic zone

为了保证冻结壁的整体稳定性以及降低其内部冻结管的断管率，冻结壁的塑性区应当被控制在一个较小的范围内。由于均质冻结壁的设计方法并没有考虑冻结壁在不同区间物理力学参数的差异性，因此在塑性区控制范围内的承载力计算结果偏大，基于该理论的冻结壁设计方法存在一定的风险性。为了解决该问题，通常会在设计中引入安全系数，而安全系数的选取存在较大的不确定性，安全系数过大会造成建井成本较高，而安全系数过小又无法保证冻结壁的安全性,因此对冻结壁的受力情况进行进一步精确的计算是解决冻结壁设计难题的重要手段。本文基于三排管冻结温度场的梯形-抛物弓的温度场模型，并根据冻结壁的弹性模量、黏聚力与温度之间的一次函数关系对冻结壁的受力机理进行了比较准确的分析，其计算结果对三排管冻结壁的设计具有重要参考意义。

6 结 论

基于三排管冻结温度场的抛物弓-等腰梯形的温度场模型，推导出了更加符合实际工况的三排管非均质冻结壁的应力计算公式，结合淮南某矿区冻土的力学参数，对该矿区进风井冻结壁的受力机理进行了计算分析可得：

(1)三排管非均质冻结壁的径向应力随着相对半径r的增大而增加，其环向应力在区间Ⅰ，Ⅱ，Ⅲ呈现不同的变化规律。当塑性区相对半径ρ≤1.4 时，环向应力的最大值出现在r=1.4处；当塑性区相对半径1.4<ρ<2.1时，环向应力的最大值出现在弹塑性分界线处；当塑性区对半径ρ≥2.1时，环向应力的最大值出现在r=2.1处。

(2)依据均质冻结壁计算方法，该冻结壁的弹性极限承载力计算结果为4.26 MPa，而考虑冻结壁的非均匀性后，该冻结壁的弹性极限承载力计算结果为3.29 MPa。由此可见，均质冻结壁计算理论关于弹性极限承载力的计算值偏大，基于该计算结果的冻结壁设计方法存在一定的风险性，而本文提出的三排管非均质冻结壁的计算理论更加符合冻结壁实际的受力情况，其计算结果将为深厚冲积层多排管冻结壁的设计提供一定的理论参考。