DOI: 10.11817/j.issn.1672-7207.2016.12.033
考虑水-力耦合效应时软化围岩解析
邹飞1, 2, 3，邹金锋2
(1. 贵州省交通规划勘察设计研究院股份有限公司，贵州 贵阳，550081；
2. 中南大学 土木工程学院，湖南 长沙， 410075；
3. 贵州省交通建设工程质量监督局，贵州 贵阳， 550000)
摘要：基于线性Mohr-Coulomb(M-C)强度准则，改进考虑水-力耦合作用下应变软化围岩深埋圆形隧道的应力与位移求解的逐步位移法。该方法全面考虑水-力耦合作用下围岩强度和变形参数的劣化、剪胀角和塑性区内弹性应变的变化。基于平面应变的假设和改进的逐步位移法，将整个塑性区分为n个同心圆环，以弹塑性交界面处的应力应变作为塑性区的初始值，获得塑性区内的应力及位移解。同时，对水-力耦合作用下的强度及变形参数进行分析。研究结果表明：在考虑水-力耦合作用下，应力减小，收敛和塑性半径均变大，位移与塑性半径提高幅值分别为27.00%和3.17%。
关键词：水-力耦合；M-C强度准则；应变软化；逐步位移法；深埋圆形隧道
中图分类号：TU457 文献标志码：A 文章编号：1672-7207(2016)12-4216-08
Solution for strain-softening surrounding rock incorporating hydraulic-mechanical coupling
ZOU Fei1, 2, 3, ZOU Jinfeng2
(1. Guizhou Transportation Planning Survey & Design Academe, Guiyang 550081, China;
2. School of Civil Engineering, Central South University, Changsha 410075, China;
3. Construction Engineering Quality Supervision Bureau of Guizhou Province, Guiyang 550000, China)
Abstract: An improved numerical stepwise procedure for the stress and displacement of a deep-buried circular opening excavated in a strain-softening rock mass under hydraulic-mechanical coupling was proposed. It followed the generalized Mohr-Coulomb (M-C) failure criterion. This method considered the deterioration of the strength, deformation, dilation angle and the variation of elastic strain in the plastic region. The total plastic region was presumably divided into n concentric annuli based on the assumption of plane strain and improved numerical stepwise procedure. The stress and displacement on each ring were expressed by the recursive relation from the initial value of the outermost annulus. Meanwhile, parametric studies were also conducted to highlight the influence of hydraulic-mechanical coupling on stress and displacement. The results show that the stress confinement is lower, and displacement and plastic radius are higher than those obtained when hydraulic-mechanical coupling is not considered. Displacement and plastic radius increase by 27.00% and 3.17%, respectively.
Key words: hydraulic-mechanical coupling; Mohr-Coulomb failure criterion; strain softening; numerical stepwise procedure; deep-buried circular opening
在工程实践中，在富水软化围岩中开挖隧道不可避免，了解富水软化围岩的应力-应变状态对工程设计与建设具有重要的理论意义和工程应用价值。对水-力耦合作用下的软化围岩受力特性及其稳定性分析也是富水软化围岩工程设计的难点和重点问题之一。虽然渗透力以及水-力耦合作用对软化围岩的应力和位移的影响至今没有明确界定，但工程实践表明，渗透力对水下隧道围岩及其支护系统有显著影响。虽然许多学者对隧道围岩的应力和位移进行了大量研究，但考虑水-力耦合作用隧道围岩尤其是软化围岩的理论解并不多见。BROWN等[1]提出了基于HOEK-BROWN屈服准则的软化隧道围岩的理论解。在BROWN解的基础上，PARK等[2-3]考虑了塑性区内弹性应变和剪胀角为变量的因素，改进了BROWN的解，获得了更严格的理论解。然而，这些研究并未考虑渗透力和水-力耦合作用。LEE等[4]提出了基于逐步应力法的软化围岩应力位移解析方法。该方法将塑性区分为许多同心圆环，每环的应力和应变增量可以通过应力平衡方程和相容方程解得。WANG等[5-7]提出了分析圆形隧道应变软化围岩应力应变的新方法，该方法针对线性Mohr-Coulomb准则的弹-脆-塑性岩体，获得弹—脆—塑性应变软化围岩的理论解。ZHANG等[8]提出利用多步弹—脆—塑性近似模拟软化围岩应力应变特性的分析模型。然而，他们均没有考虑渗透力的影响。FAHIMIFAR等[9]提出了获得地下水下开挖的圆形隧道简化闭合解的方法。KOLYMBAS等[10]推导出定常流动时地下水渗入隧道渗流量及渗透力解析式。ALONSO 等[11]建立了软化围岩应力应变的通用求解方法。然而，上述研究也没有考虑渗透力和水-力耦合的影响。现有的研究主要集中在圆形隧道的理论或数值解[12]，且都只是考虑了渗透力的影响[13-16]，很少考虑渗透力作用下软化围岩的应力应变变化特性，而且未综合考虑强度和变形劣化以及在塑性区剪胀角和弹性应变的变化。为此，本文作者综合考虑强度和变形劣化以及在塑性区剪胀角和弹性应变的变化，同时考虑水-力耦合对软化围岩应力应变的影响规律，进而改进现有的软化围岩求解过程中的逐步位移法，利用改进的逐步位移法，针对线性Mohr-Coulomb软化围岩，实现富水软化围岩应力应变的求解。
1 求解方法
1.1 问题定义
隧道围岩受力状态及计算分析模型如图1所示，图1中：σθ和σr分别为大、小主应力；σc为完整岩石的单轴抗压强度；R为塑性区半径；RS为软化区半径。在连续、均匀、各向同性的弹性介质中开挖半径为r0圆形深埋隧道。隧道围岩初始地应力为σ0，孔隙水压力为pw(r，θ)，开挖后沿半径方向均匀作用在隧道壁表面内部支护力为pin。在极坐标系下，由于考虑到力学分析的轴对称条件，可不考虑重力场，围岩的应力和位移仅是半径r的函数。在水力分析中，围岩周围的孔隙水压力pw(r，θ)及渗透力Fr的分布是非轴对称的，即孔隙水压力pw(r，θ)和渗透力Fr均是r和θ的函数。因此，由力学分析得到的是任意方向的应力应变，而由水力分析得到的是某一特定方向上的应力应变。
图1 隧道围岩受力状态及计算分析模型
Fig. 1 Axisymmetric calculation model and stress distribution of rocks around tunnel
1.2 基本假设
实际工程中工况复杂，为了便于力学分析，根据取重舍轻的原则进行如下假设：隧道为深埋隧道，隧道围岩初始地应力各个方向的差异性较小，同时考虑到力学分析轴对称条件，隧道围岩受到的初始地应力可简化为各个方向相同的压力σ0；隧道围岩视为均质、各向同性的、连续的透水介质。在水-力分析(渗透力和孔隙压力分析)中，孔隙水压力是径向(r)和方向角(θ)的函数，即pw= pw(r，θ)。
1.3 孔隙水压力和渗透力
孔洞周围不同区域的渗透网如图2所示，在地下水位线以下开挖具有恒定水头的圆形隧道，孔隙水压力用伯努利方程确定。基于稳定流体状态，并以隧道埋深作为基准水平线。当塑性区半径与地下水位到隧道深度之比(R/h1)足够小时，可假设流网在塑性区沿着径向方向变化。弹性区与塑性区圆形边界面上的水头处处相等，假定向内渗流速率为正。
图2 孔洞周围不同区域的渗透网
Fig. 2 Seepage flow net in different regions around opening
基于KOLYMBAS等[10]提出的围岩中的孔隙水压力公式，可推导出隧道弹性空间任何1点的孔隙水压力公式：
(1)
式中：
；
。当θ=0°时，有
(2)
式中：R为隧道塑性半径；h1为地下水位到隧道的深度；(r, θ)为极坐标；
为在隧道壁上最终孔隙压力；
为围绕隧道壁在弹性-塑性交界面的最终水头；
为水容重。利用式(2)可得径向渗透力表达式为
(3)
式中：ζ为Biot–Willis多孔弹性耦合常数，在本文中假定为1。
1.4 水-力耦合方程
在考虑水-力耦合的软化围岩隧道稳定性分析中，根据BROWN等[1]的研究成果，在弹性阶段围岩的渗透系数为常量k0r，在围岩的塑性范围内，其渗透系数可假设同围岩的变形有关，即
(4)
式中：kr为裂隙岩体的次生渗透系数；k0r为岩体的初始渗透系数；
或
，为体积应变；η为比例常数(耦合常数)。
由FAHIMIFAR等[9]的研究成果可知，孔隙水压力[pw(r，θ)]和在水-力耦合作用下塑性区渗流量(q)可分别表示为：
(5)
(6)
1.5 强度和变形参数的演变方程
根据ALEJANO[11]的研究成果，可采用塑性剪应变描述应变软化岩体的强度和变形参数的软化规律：
(7)
式中：γp为塑性偏应变；
和
分别为大、小主塑性应变。
利用塑性剪切应变描述围岩物理参数的双线性函数为
(8)
式中：ωp为强度参数，本文考虑黏聚力C，内摩擦角φ，剪胀角ξ随塑性偏应变γp的变化；
为临界塑性偏应变；下标p和r分别表示峰值和残余值；上标p和e分别表示塑性区和弹性区。
1.6 Mohr-Coulomb屈服准则
Mohr-Coulomb屈服准则为
(9)
式中：
和
分别为最大主应力与最小应力；
；
。
2 弹性区内的应力与应变求解
2.1 平衡方程与应力边界条件
考虑渗透力时，极坐标中的轴对称平衡方程为
(10)
式中：
和
分别为径向和环向应力。边界条件可表示为：
(11)
(12)
2.2 应力和应变求解
随着隧道开挖，径向应力在弹性区外边界处从初始值(
)开始减小。此外，初始岩石变形是弹性的。在径向对称的条件下，弹性应力-应变和位移-应变关系可以表示为：
(13)
(14)
式中：E为弹性模量；v为泊松比；εr和εθ分别为径向应变和环向应变。
当考虑平面应变条件时，应力-应变方程可以用胡克定律得到：
(15)
(16)
式中： 
和
分别为径向和环向弹性应变。将式(3)，(13)和(16)代入式(10)得
(17)
式中：u为径向位移；pw为在径向距离r的孔隙水压力。
式(17)是线性差分方程，因此，可以应用叠加原理进行求解。总位移可以分解成uw和ub共2部分，分别表示由内部渗流体的力和边界的压力导致的位移，即
以及
导致的位移。
求解式(17)可得r处的径向位移和应变为：
(18)
(19)
(20)
式中：
；
；
；
；C1和C2为积分常数，可由边界条件获得。
(21)
(22)
将式(19)和(20)代入式(15)和(16)，可得应力
和
。
由边界的压力引起的位移和应力可由下式计算 获得：
(23)
(24)
(25)
总的位移和应变则可由下式表示：
(26)
(27)
(28)
3 塑性区内的应力和位移求解
在应变软化岩体中，应力和位移很难获得解析解，尤其是考虑水-力耦合作用时。但可以构造考虑水-力耦合作用的逐步位移迭代法(该耦合作用引入了强度破坏、变形、剪胀角以及塑性区弹性应变的变化)，实现富水软化围岩的应力位移求解。具体过程如下。
将总塑性区分为n个同心环，如图3所示。
图3 塑性区的应变软化模型
Fig. 3 Strain-softening rock model in plastic zone
第j个环的外径为
，内径为
。其中，每个圆环厚度一般不相等，因为每个圆环的厚度是由平衡方程、边界条件等决定的。在围岩处于弹塑性临界状态时，最外部环的外径为
，相应的应力可以写为
(29)
因此，最外层环的应力和位移由式(26)~(28)确定。在轴向对称情况下，环可以足够薄，应变-位移关系可以由BROWN等[1]提出的求解方法获得：
(30)
方程(30)可转化为
(31)
半径r(j)处的位移可表示为
(32)
半径r(j)处的位移以及塑性半径(R)的位移可利用逐步迭代法及MATLAB编程获得。假设任意足够小的环向应变增量为
，则第j环的应变表示为：
(33)
(34)
(35)
式中：塑性区第j环的弹性应变
和
可通过胡克定律确定；β为剪胀系数。
对于第j环，引入水-力耦合作用的平衡微分方程可由方程(9)和(10)求解。
对于第j环，考虑水-力耦合作用的平衡差分方程可由式(3)，(5)，(9)和(10)获得，具体如下：
(36)
式中：
；
；
；
；
；
。
环向应力可以通过M-C屈服准则公式得到：
(37)
每个环上的应力通过利用递归关系从最外侧环的初始值逐步获得，而最外环的初始值可以通过方程式(27)和(28)得到。当每个环厚度足够小时，就可以得到应变软化区域内的应力应变曲线。
4 验证
为了验证本文计算方法正确性，将本方法数值计算结果与ALONSO等[11]的计算结果进行对比分析，取参数如下：
隧道开半径r0=3.0 m，弹性模量E=10 GPa，泊松比v=0.25，初始地应力σ0=20 MPa，黏聚力峰值和残余值Cp=1.0，Cr=0.7，内摩擦角峰值和残余值φp=30°，φr=22°，剪胀角峰值和残余值ξp=ξr=3.75°，临界塑性偏应变
=0.004，水力耦合常数η=1×105。计算结果如图4所示。
从图4可见：本方法不考虑水力耦合作用时所得结果与ALONSO等[11]的计算结果基本一致，例如不考虑水力耦合作用时计算的临界支护力pic=9.134 MPa，该结果与ALONSO等[11]计算结果一样，由此可以证明本方法的正确性。在考虑水力耦合作用时，位移增大。例如考虑水力耦合pi/σ0为28.9，6.7，2.8时，(u/b)(2G)(σ0-pic)分别为28.9，6.7，2.8；当在不考虑水力耦合作用时，(u/b)(2G)(σ0-pic)分别为28.9，6.7，2.8，均小于考虑水-力耦合时相应值，这也验证了在考虑水力耦合下圆形隧道收敛性增大的结论。
图4 软化围岩的地面响应曲线
Fig. 4 Ground reaction curves of softening rock
5 算例分析
为了验证本文方法的正确性，同时分析相关参数对理论解的影响规律，取如下参数进行计算(参考FAHIMIFAR等[4]的计算结果)：r0=3.0 m，σ0=20 MPa，pin=1.0 MPa，cp=1.0，cr=0.7，φp=30°，φr=22°，ξp=5°，ξr=3.75°，rw=0.01，h1=400 m，
，kc=10-6 m/s，γ=0.004，η=1×105。本文主要分析考虑水-力耦合和渗透力对本文理论解的影响，同时分析影响软化特性的相关参数对理论解的影响规律。
5.1 渗透力的影响
为了分析渗透力的影响，考虑渗透力作用下与简化条件下应力与位移的分布如图5所示。
从图5可见：考虑渗透力的应力和位移比不考虑渗透力时的大；考虑渗透力后，塑性半径与位移都相应增大。其原因可能是渗透力使围岩的强度减弱，导致塑性区扩大。例如，在隧道内壁处，考虑渗透力时隧道内壁处位移u=0.043 m，此时，塑性半径R= 8.770 m；不考虑渗透力时，位移u=0.034 m，相对应的塑性半径R=8.500 m。而位移与塑性半径提高的幅值分别为3.17%及27.00%。因此，在隧道设计中，渗透力和水-力耦合作用不能被忽视。
5.2 应变软化的影响
为了分析考虑水-力耦合作用时应变软化参数对软化围岩应力与位移的影响，分别取γp为0.004，0.008以及0.012进行对比，结果如图6所示。
从图6可见：软化参数取值变化时，位移与塑性半径的变化很大，这说明软化参数是决定软化围岩塑性区范围的关键；随着软化参数增大，位移与塑性半径均减小。例如，当r/r0=1时，若γp=0.004，则位移u=0.043 m，塑性半径R=8.770 m；若γp=0.008，则位移u=0.033 m，塑性半径R=7.990 m；若γp=0.012，则位移u=0.026 m，塑性半径R=7.390 m。
5.3 强度参数的影响
为了分析考虑水-力耦合作用时强度参数对位移与应力的影响，对于C分以下情况讨论：1) Cp=Cr=0.7；2) Cp=1，Cr=0.7；3) Cp=Cr=1。对于φ分以下3种情况讨论：1) φp=φr=30°；2) φp=30°，φr=22°；3) φp=φr=22°。结果分别如图7与图8所示。
从图7与图8可见：当强度参数C增大时，塑性区径向应力与环向应力增大，位移与塑性半径随着强度参数C的增大不断减小；当强度参数φ增大时，塑性区径向应力与环向应力均增大，位移与塑性半径随着强度参数φ的增大不断减小。因此，考虑水-力耦合作用时，强度参数C与φ对应力与位移的影响趋势相同，且对应力与位移均具有显著影响。
图5 不同渗透力作用下的应力与位移分布
Fig. 5 Distribution of displacement and stress under different seepage forces
图6 不同软化参数下应力与位移分布
Fig. 6 Distribution of stress and displacement under different critical values of strain-softening parameters
图7 C不同时应力与位移分布
Fig. 7 Distributions of stress and displacement under different strength parameters C
图8 φ不同时的应力与位移分布
Fig. 8 Distributions of stress and displacement under different strength parameters φ
6 结论
1) 重构了考虑强度和变形参数劣化、剪胀角和塑性区弹性应变为变量的应变软化围岩的演化方程，同时在该方程中引入水-力耦合作用，提出了基于重构的应变演化方程的应变软化围岩应力位移求解的逐步位移法。该方法基于线性M-C岩体，可有效求解圆形隧道应力位移，尤其是塑性区、软化区及残余区的应力与位移。
2) 在考虑水-力耦合作用下，应力约束减小，隧道围岩的收敛、位移和塑性半径均增大。渗透力、应变软化参数、强度参数对围岩的收敛有显著影响。
3) 在实际隧道工程中，该方法也可用于估算隧道围岩的加固程度，也可用于现场监测结果对比。
参考文献：
[1] BROWN E T, BRAY J W, LADANYI B, et al. Ground response curves for rock tunnels[J]. Journal of Geotechnical Engineering, 1983, 109(1): 15-39.
[2] PARK K H, TONTAVANICH B, LEE J G. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses[J]. Tunnelling and Underground Space Technology, 2008, 23(2): 151-159.
[3] PARK K H. Similarity solution for a spherical or circular opening in elastic-strain softening rock mass[J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 71: 151-159.
[4] LEE Y K, PIETRUSZCZAK S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass[J]. Tunnelling and Underground Space Technology, 2008, 23(5): 588-599.
[5] WANG Shuilin, YIN Xiaotao, TANG Hua, et al. A new approach for analyzing circular tunnel in strain-softening rock masses[J]. International Journal of Rock Mechanics and Mining Sciences, 2010, 47(1): 170-178.
[6] WANG Shuilin, ZHENG Hong, LI Chunguang, et al. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(1): 67-76.
[7] WANG Shuilin, YIN Shunde, WU Zhenjun. Strain-softening analysis of a spherical cavity[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(2): 182-202.
[8] ZHANG Qiang, JIANG Binsong, WANG Shuilin, et al. Elasto-plastic analysis of a circular opening in strain-softening rock mass[J]. International Journal of Rock Mechanics and Mining Sciences, 2012, 50: 38-46.
[9] FAHIMIFAR A, ZAREIFARD M R. A theoretical solution for analysis of tunnels below groundwater considering the hydraulic-mechanical coupling[J]. Tunnelling and Underground Space Technology, 2009, 24(6): 634-646.
[10] KOLYMBAS D, WAGNER P. Groundwater ingress to tunnels–the exact analytical solution[J]. Tunnelling and Underground Space Technology, 2007, 22(1): 23-27.
[11] ALONSO E, ALEJANO L R, VARAS F, et al. Ground response curves for rock masses exhibiting strain-softening behaviour[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27(13): 1153-1185.
[12] ZOU Jinfeng, SU Yu. Theoretical solutions of a circular tunnel with the influence of the out-of-plane stress based on the generalized Hoek–Brown failure criterion[J]. International Journal of Geomechanics, 2015, 16(3): 06015006.
[13] 王星华, 章敏, 王随新. 考虑渗流及软化的海底隧道围岩弹塑性分析[J]. 岩土力学, 2009, 11(5): 3267-3272.
WANG Xinghua, ZHANG Min, WANG Suixin. Theoretical solutions and analysis of submarine Tunnel in strain-softening surrounding rock considering the seepage force[J]. Rock and Soil Mechanics, 2009, 11(5): 3267-3272.
[14] 黄阜, 杨小礼. 考虑渗透力和原始 Hoek-Brown 屈服准则时圆形洞室解析解[J]. 岩土力学, 2010, 31(5): 1627-1632.
HUANG Fu, YANG Xiaoli. Analytical solution for circular cavern considering the seepage force and original Hoek-Brown yield criterion[J]. Rock and Soil Mechanics, 2010, 31(5): 1627-1632.
[15] ZOU Jinfeng, LI Shuaishuai, XU Yuan, et al. Theoretical solutions for a circular opening in an elastic-brittle-plastic rock mass incorporating the out-of-plane stress and seepage force[J]. KSCE Journal of Civil Engineering, 2016, 20(2): 687-701.
[16] 邹金锋, 李帅帅, 张勇, 等. 考虑轴向力和渗透力时软化围岩隧道解析[J]. 力学学报, 2014, 46 (5): 747-755.
ZOU Jinfeng, LI Shuaishuai, ZHANG Yong, et al. Theoretical solution and analysis of circular tunnel for the strain-softening rock mass considering the axial in situ stress and seepage force[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(5): 747-755.
(编辑 陈灿华)
收稿日期：2016-03-14；修回日期：2016-05-10
基金项目(Foundation item)：国家重点基础研究发展计划(973计划)项目(2013CB036004)；国家自然科学基金资助项目(51208523) (Project(2013CB036004) supported by the National Basic Research Program(973 Program) of China; Project(51208523) supported by the National Natural Science Foundation of China)
通信作者：邹飞，博士，高级工程师，从事岩石动力学和边坡极限分析研究；E-mail：68870417@qq.com
