基于对数应变考虑渗流影响的软弱围岩圆形隧洞弹塑性解

于旭光

(唐山工业职业技术学院，河北 唐山 063299)

0 引言

随着南水北调工程的建设和运营，出现了大量深埋水工圆形隧洞。深埋水工圆形隧洞在开挖过程中会出现围岩卸荷、应力重分布，进而围岩会出现一定的塑性区。同时，由于深埋水工圆形隧洞在运行期存在一定的内水压力，当内水压力过大时，围岩主应力顺序将发生改变，进而会重新产生塑性区，著名的Fenner公式[1]和Kastner公式[2]将不再适用。针对上述问题，大量学者对此进行了研究。李宗利等[3]采用理想弹塑性模型并基于Mohr-Coulomb准则推导了考虑渗流场影响的深埋圆形隧洞弹塑性解；刘成学等[4]采用理想弹塑性模型并基于Mohr-Coulomb准则以及考虑应力重分布推导了渗流场下深埋圆形隧洞弹塑性解；王睢等[5]采用理想弹塑性模型并基于DP准则推导了考虑渗流场影响下的深埋圆形隧洞弹塑性解；张常光等[6]采用弹脆塑性模型并基于统一强度理论推导了考虑渗流场影响下的深埋圆形隧洞弹塑性解；张黎明等[7]采用应变非线性软化模型推导了考虑渗流场影响下的深埋圆形隧洞弹塑性解。

从文献[3-7]来看，其研究成果都是基于小应变假设，未采用大变形理论分析方法。当(p0-pi)>2G(p0、pi、G分别为初始地应力、支护力、剪切模量)时，根据小应变假设求解的弹性位移可能大于开挖半径，因此，上述研究成果求得的位移仅适用于硬岩工程，对软岩工程不适用。由于在大变形理论中，对数应变(即真实应变)能真实反映变形的积累过程,故在固体力学中应用广泛[8-11]。为此，国内外学者展开了大量研究。蒋斌松等[12]基于对数应变推导了圆形巷道围岩的弹性解，但未考虑渗流影响和围岩处于弹塑性状态时的弹塑性解；Papanastasiou等[13]基于Mohr-Coulomb准则和DP准则推导了理想弹塑性模型的圆形巷道扩张和收缩的大变形弹塑性解；Park[14]采用弹性区小应变假设，推导了基于Mohr-Coulomb准则的理想弹塑性模型球形或圆形岩土体的大应变相似解；Guan等[15]在围岩弹性区基于小应变假设、塑性区采用对数应变推导了基于Mohr-Coulomb准则的应变软化模型圆形隧道的大应变数值解；Zhang等[16]在围岩弹性区、塑性区均采用对数应变推导了基于Mohr-Coulomb准则、Hoek-Brown准则的应变软化模型圆形隧道的大应变数值解。

可以看出，文献[12-16]推导了软弱围岩的弹性解或弹塑性解，仍有3点不足之处：1)除了文献[13]，其他文献均采用未考虑中间主应力的Mohr-Coulomb准则或Hoek-Brown准则而使计算结果偏于保守；2)均是考虑施工期(径向应力为第一主应力)时的解，而未求解深埋水工圆形隧洞在运行期由于内水压力过大，切向应力为第一主应力时的解，也未求解隧洞运行一段时期后放水检修时检修期的解；3)未考虑深埋、富水环境下渗流场对软弱围岩圆形隧洞的影响。

综上所述，本文首先采用能考虑中间主应力和最小主剪应力影响的三剪应力统一强度理论[17-18]，基于拉格朗日坐标下的对数应变推导施工期、运行期和检修期3种工况下理想弹塑性模型软弱围岩的弹塑性解；然后，通过实例验证本文求解方法的正确性；最后，分析弹性模量、泊松比、孔隙水压力和强度准则4个参数对塑性区厚度和洞壁处径向位移的影响。

1 三剪应力统一强度理论

三剪应力统一强度理论是高江平等[17-18]根据菱形十二面单元体所有应力分量所建立的，于旭光等[19-20]经过推导得到平面应变状态下的表达式为：

(1)

式中：σ1为第一主应力；σ3为第三主应力；φ0为内摩擦角；c0为黏聚力；b为中间主剪应力τ12及其作用面上正应力σ12综合影响的作用系数，0≤b≤1；c为最小主剪应力τ23及其作用面上正应力σ23综合影响的作用系数，0≤c≤1。

需特别指出：当b=0、c=0时，三剪应力统一强度理论对应为Mohr-Coulomb准则；当b=1、c=0时，对应为双剪应力准则；当b=1、c=1时，对应为三剪应力准则；当b、c分别取0～1的其他值时，可得到一系列新的强度准则。

2 渗流场求解

(a)施工期

施工期的轴对称恒定渗流连续微分方程为：

(2)

将边界条件H|r=ri=hi，H|r=r0=h0代入式(2)可得到任意r处水头

(3)

式中：hi为r=ri处水头；h0为r=r0处水头；β=r0/ri。

3 圆形隧洞围岩弹塑性解

3.1 基本方程

平衡微分方程为：

(4)

式中：ξ为岩石等效孔隙水压力系数；γw为水的重度；σr和σθ分别为径向和切向有效应力，以压应力为正，拉应力为负(本文中所有应力均为有效应力)。

径向位移u(施工期)、u*(运行期)和应变(径向应变εr，切向应变εθ)之间的关系可以用对数函数的形式表示如下[16,21]。

1)对于施工期或检修期

(5)

(6)

2)对于运行期

(7)

(8)

特别说明：式(7)—(8)与式(5)—(6)不同，是因为本文假定指向洞内位移为正，指向围岩位移为负，因此有正负号之分。

忽略隧洞围岩开挖前初始地应力p0所引起的变形，各向同性材料的平面应变本构方程有如下2种表述方式[22]。

对于应变-应力公式(Hooke)：

(9)

(10)

式(9)—(10)中：E为弹性模量；ν为泊松比。

对于应力-应变公式(Lame)：

σr-p0=2Gεr+λ(εr+εθ)；

(11)

σθ-p0=2Gεθ+λ(εr+εθ)。

(12)

3.2 施工期围岩弹塑性解

在施工期，隧洞围岩由于开挖卸载而出现一定范围的塑性区，此时切向应力σθ、径向应力σr、轴向应力σz均为压应力，径向应力σr为第一主应力，有-σr>-σθ，即|σr|<|σθ|，因此式(1)变为：

σθ=Aσr+B。

(13)

因此，施工期求解的是第一主应力为径向应力时的隧洞围岩弹塑性解。

3.2.1 弹性区应力和位移

施工期围岩弹性区应力和位移可仿照文献[16,21]的求解方法。将式(5)—(6)代入式(11)—(12),同时联立式(4)，得到：

(14)

由于

(15)

(16)

将式(5)—(6)代入式(11)，并将式(11)两端进行求导得到：

(17)

将式(15)—(16)代入式(17)得到：

(18)

将式(18)代入式(14)得到：

(19)

式中：

联立式(14)和式(19)，得到：

dσr=f(X)dX。

(20)

对于计算区域r=r0内的圆形隧洞，径向位移随半径的增大逐渐减小，且在r=r0(可看成无限远)处不发生位移。因此，当r=r0时，X=1。并且根据X表达式，可以看出X取值区间为(0,1]。根据f(X)表达式，由洛必达法则可求得f(1)的值为：

(21)

对式(20)进行积分，并以r=r0，σr=p0为边界条件，可得到弹性区径向应力表达式为：

(22)

将式(11)减去式(12)，联立式(5)—(6)，可得到弹性区切向应力表达式为：

σθ=σr-2GlnX。

(23)

在弹塑性交界面处，将式(22)—(23)代入式(13)，可得到关于与施工期围岩塑性区半径rp相对应的Xrp的函数表达式为：

(24)

将式(19)进行积分，可以确定弹性区任意位置r处的表达式为：

(25)

将式(5)减去式(6)，得到：

εr=εθ+lnX。

(26)

将式(26)代入式(11)，并联立式(6)，得到弹性区径向位移u的表达式：

(27)

3.2.2 塑性区应力和位移

在围岩塑性区，将式(13)代入式(4)，并结合边界条件r=ri，σr=pi，得到塑性区径向和切向应力表达式为：

(28)

(29)

在围岩塑性区，其体积具有明显剪胀特性，即塑性体积应变不等于0，根据线性非关联流动法则，有

(30)

式中：χ为剪胀特性参数，χ=(1+sinψ)/(1-sinψ)，ψ为剪胀角。

由式(30)可得

(31)

将式(5)—(6)代入式(31)，可得

(32)

式(32)也可改写成如下表达式：

(33)

首先将式(28)和(29)分别代入式(9)—(10)，然后代入式(33)，最后积分式(33)并以弹塑性交界面r=rp处的径向位移urp为边界条件，整理得到塑性区径向位移u的表达式为：

(34)

式中：

为了方便计算，将式(34)中指数函数展开成幂级数，其表达式为：

(35)

将式(35)代入式(34)，即可得到塑性区径向位移u的表达式为：

(36)

在弹塑性交界面r=rp处，径向应力相等(即式(22)与式(28)相等)，可得到rp/ri计算表达式为：

(37)

由于洞壁处径向位移可定义为初始位置ra与当前半径ri之差，即

ra-ri=u|r=ri。

(38)

当r=ri时，将式(27)、式(37)、式(38)代入式(36)，得到ri计算表达式为：

{1-[T(rp)]k+D5}。

(39)

将ri值代入式(28)、式(29)、式(36)可得到施工期围岩塑性区应力和径向位移。

3.3 运行期围岩弹塑性解

σr=Aσθ+B。

(40)

因此，运行期求解的是第一主应力为切向应力时的隧洞围岩弹塑性解。

3.3.1 弹性区应力和位移

仿照3.2.1节的求解方法，运行期围岩弹性区径向应力和切向应力的表达式与式(22)、式(23)相同。但需注意：式(22)和式(23)中的X采用X*来代替。

(41)

3.3.2 塑性区应力和位移

(42)

(43)

根据线性非关联流动法则，有

(44)

由式(44)可得

(45)

将式(7)—(8)代入式(45)，可得

(46)

(47)

(48)

式中：

同样，将式(48)中指数函数展开成幂级数，其表达式为：

(49)

将式(49)代入式(48)，即可得到塑性区径向位移u*的表达式为：

(50)

(51)

(52)

(53)

3.4 检修期围岩弹塑性解

由于运行期围岩出现塑性变形，应力和应变关系呈现非线性，叠加原理不再适用，为求解隧洞围岩因正常检修放空之后的应力和位移，可采用有限环变位协调法进行求解[23]。

(54)

(55)

式中：Ec、νc、rc分别为衬砌的弹性模量、泊松比、内半径；Δpc、Δph、Δpi分别为衬砌的内半径处内水压力变化量、外半径表面力的变化量、围岩与衬砌相互作用力变化量。

(56)

4 实例验证

为了验证本文方法的正确性，选取文献[16]中的软弱围岩，具体参数如下：初始开挖半径ra=1 m，初始地应力p0=1 MPa，弹性模量E=26 MPa，泊松比ν=0.3，黏聚力c0=0.233 MPa，内摩擦角φ0=40°，剪胀角ψ=20°。除以上工程参数外，另附加β=50，ξ=1，h0-hi=0 m，b=0，c=0。采用MATLAB自编程序(程序中将区间[0,1]分为105等分，在每个子区间利用辛普森公式计算以保证计算精度)进行计算，并与文献[16]计算结果进行对比。当内压力pi取不同值时，图2给出了塑性区厚度Qp与pi/p0关系曲线以及pi/p0与ua/ra关系曲线。其中：Qp=rp-ri，pi/p0为施工期内压力与初始地应力的比值，ua/ra为施工期洞壁处径向位移与初始开挖半径比值。经过与文献[16]的计算结果吻合，验证了本文方法正确性。

(a)Qp与pi/p0关系曲线

5 参数分析

从第3节推导中可以看出，本文解与小应变解不同，弹塑性交界面处的应力不仅依赖于强度准则，还依赖于弹性模量和泊松比，进而将影响软弱围岩塑性区径向位移和半径。由于检修期围岩应力状态与施工期类似，应力、位移和塑性区半径变化规律也类似。下面仅取弹性模量、泊松比、孔隙水压力和强度准则4个参数对施工期和运行期2种工况下塑性区厚度和洞壁处径向位移进行分析。在本节分析中，除另附加β*=30外，未特殊说明时其他参数与第4节相同。

5.1 弹性模量的影响

(a)施工期

(a)施工期

从图4可以看出：在施工期围岩塑性区厚度随弹性模量的增大呈现先增大后逐渐趋于一个稳定值(即小应变解0.148 3 m)的趋势，并且当E≥80p0时本文解和小应变解近似相等；而运行期围岩塑性区厚度随弹性模量的增大呈现逐渐减小后逐渐趋于一个稳定值(即小应变解0.110 2 m)的趋势，并且当E≥200p0时本文解和小应变解近似相等。因此，对于弹性模量较小的软岩来说，在施工期时采用本文解求解的塑性区半径比小应变解求解的塑性区半径小，可在初期支护时减小锚杆的长度；而在运行期时求解的塑性区半径大，可采取增加衬砌刚度方法。

(a)施工期

从图5可以看出：在施工期和运行期2种工况下，洞壁处径向位移随弹性模量的增大逐渐减小；并且当E≥50p0时，2种工况下本文解和小应变解得到的洞壁处径向位移几乎相等。

5.2 泊松比的影响

图6 施工期Qp、ua/ra与υ关系曲线

图7 运行期与υ关系曲线

从图6和图7可以看出：塑性区厚度随泊松比增大几乎无变化，其值大约为0.143 m(施工期)、0.119 4 m(运行期)；洞壁处径向位移随泊松比增大呈线性增大，泊松比从0.1增大到0.4时，洞壁处径向位移分别增大了18.63%(施工期)、28.00%(运行期)。

5.3 孔隙水压力的影响

图8 施工期塑性区厚度Qp与pi/p0的关系曲线

从图8和图9可以看出：1)在施工期围岩塑性区厚度随孔隙水压力(h0-hi=0、100、200 m)的增大而增大，孔隙水压力对塑性区的发展起增强作用；而在运行期围岩塑性区厚度随孔隙水压力的增大而减小，孔隙水压力对塑性区的发展起抑制作用。2)塑性区厚度随内压力变化分为3个阶段，这与5.1节所述相同，但是施工阶段由塑性转变为弹性的临界内压力变化显著，与h0-hi=0、100、200 m相对应的临界内压力分别为0.177、0.233、0.277 MPa；随着内压力增大，由弹性再次转变为塑性并且主应力发生变化的临界内压力分别为1.80、2.35、2.74 MPa。

图9 运行期塑性区厚度与关系曲线

5.4 强度准则的影响

图10和图11分别示出了施工期、运行期在不同强度准则下对塑性区厚度的影响。

从图10和图11可以看出：不同强度准则下的塑性区厚度随内压力的增大呈近似线性减小(施工期)、线性增大(运行期)的趋势，且斜率近似相等。不同强度准则下对应的临界内压力差异很大，在施工期，当b=1、c=1时，对应的临界内压力比b-0、c=0时施工期减小97.18%、运行期增大19.44%。实际中应考虑强度准则效应，可以提高围岩的承载能力。

图11 运行期不同强度准则对塑性区厚度的影响

不同强度准则下施工期、运行期洞壁处径向位移与内压力的关系曲线和塑性区厚度与内压力的关系曲线变化规律相似，在此不再赘述。

6 结论与讨论

1)本文根据三剪应力统一强度理论，基于对数应变并考虑渗流影响推导了施工期、运行期和检修期3种工况下软弱围岩的弹塑性解，该解可以充分考虑大变形的影响，因而本文解具有更广泛的适用性。

2)围岩塑性区厚度随弹性模量的增大呈现施工期逐渐增大、运行期逐渐减小后逐渐趋于一个稳定值(即小应变解)的趋势，洞壁处径向位移在施工期和运行期2种工况下均随弹性模量的增大逐渐减小；施工期和运行期围岩塑性区厚度随泊松比增大几乎无变化，而洞壁处径向位移随泊松比增大呈线性增大；围岩随内压力的增大经历塑性、弹性、塑性3个阶段，对应的塑性区厚度和洞壁处径向位移随孔隙水压力的增大而增大(施工期)、减小(运行期)；不同强度准则下的塑性区厚度和洞壁处径向位移均随内压力的增大呈近似线性减小(施工期)、线性增大(运行期)的趋势，且斜率近似相等。

3)检修期围岩应力状态与施工期类似，因此塑性区厚度和洞壁处径向位移变化规律相类似。

4)本文研究的围岩为理想弹塑性模型，未考虑围岩流变、材料的不均匀性及不同力学模型(弹脆塑性模型、应变软化模型)的影响，并且本文解给出的是未在衬砌外半径处设置环向排水槽、纵向排水管的弹塑性解，这些内容将在后续做进一步研究。