0 引言

近些年来，隧道建设在交通、水利水电等领域发挥着日益重要的作用^[1-2]，然而隧道开挖往往会破坏地层的原始平衡状态，造成围岩应力重分布，进而产生一定范围的塑性区。在塑性区内，岩石往往会出现剪胀效应及应变软化现象，因此需要考虑其对隧道围岩塑性区及变形等力学行为的影响^[3]。剪胀是指围岩在剪应力作用下，岩石矿物颗粒排列发生变化从而产生塑性变形，导致体积膨胀的现象。大量试验及工程实践表明，岩石破坏与其剪胀机制密切相关^[4]，这又直接影响到隧道围岩的稳定性。而应变软化则是指在围岩应力应变关系曲线峰后阶段，岩石强度等参数随塑性应变的增加而逐渐降低，直至达到残余值的现象。

首先，剪胀角是评价岩石峰后剪胀行为的重要参数，因此许多学者研究了其不同取值方法对隧道围岩力学特性的影响。潘继良^[5]取剪胀角为小于岩石内摩擦角的一个常量，认为剪胀角不会影响围岩塑性区半径，但会影响塑性区的位移分布，若不考虑剪胀将会低估围岩变形。范文^[6]认为随着岩石剪胀角的增加，隧道围岩塑性残余区半径增大，但增幅较小，而洞壁位移增幅则较大。王凤云^[7]认为剪胀系数随围岩变形衰减，通过采用半数值方法，基于统一强度理论建立了同时考虑岩石剪胀和软化效应的隧道围岩力学模型。王凤云^[8]、高永涛^[9]按照Hoek和Brown^[10]工程实际经验分类方法，将岩体质量分为较好、中等和较差，并取其剪胀角分别为内摩擦角的1/4，1/8和0。以上都是基于常量剪胀角模型建立的力学模型，而高启程^[11]通过对饱水砂岩的卸围压试验发现岩石剪胀角随其塑性剪切应变的增加先降低后增加，而非定值。曹朔^[12]研究了不同岩石剪胀角对隧道围岩洞壁位移的影响，发现其影响程度较大。由此可以看出应考虑隧道开挖过程中岩石剪胀角的变化，否则得到的隧道围岩应力、位移等将与实际情况相差较大^[13]。大量试验研究也表明，岩石剪胀角并非定值，而是与临界软化系数和围压等密切相关，呈现非线性变化规律。Pourhosseini^[14]基于砂岩、粉砂岩和泥岩数据归纳总结出了初始剪胀角及其演化公式，并采用数值模拟与三轴试验结果对其进行了验证。Cai^[15]提出了考虑围压和临界软化系数的岩石峰后剪胀角模型，该模型形式简洁，参数少且使用方便。

其次，关于岩石应变软化对隧道围岩力学特性的影响，江权^[16]认为高应力洞室开挖过程中，围岩峰后弹性模量、黏聚力等会随着应变的增加而降低，而内摩擦角则相反，并由此建立了硬岩参数劣化的弹塑性本构模型，并将其应用于某二级水电站厂房开挖数值模拟。王俊杰^[17]通过考虑深埋软岩隧道塑性区围岩的强度及刚度参数的劣化效应、流变特性及支护结构施作时机等建立了深埋软岩隧道复合式衬砌力学模型，计算结果与实测吻合很好。郑可跃^[18]基于三维H-B强度准则建立了考虑围压和中主应力的深埋软岩隧道弹塑性解，并认为二者将分别加剧和降低围岩的软化程度，并进而影响隧道的挤压变形。余俊^[19]考虑黏聚力弱化，建立了用内摩擦角和黏聚力表征的岩石峰后应力应变方程，并基于大理岩试验结果验证了峰后内摩擦角和黏聚力呈二次函数跌落的假设。而余敏^[20]则认为岩石强度衰减的原因主要在于黏聚力线性跌落，而内摩擦角不变，并考虑岩石剪胀性质，推导出围岩应力、塑性区半径和洞壁位移的解析式。Lee^[21]假设强度参数为临界软化系数的线性函数，采用应力释放的观点将隧道围岩塑性区分为若干圆环，建立了一种半解析半数值的方法，得到了隧道围岩应力及位移分布规律。Zareifard^[22]考虑了隧道开挖扰动对围岩的损伤，基于Hoek-Brown准则建立了隧道围岩线性及非线性软化模型，得到了围岩损伤区的半理论半数值解。王凤云^[8]根据统一强度理论及塑性增量的非关联流动法则，提出了考虑中主应力及岩石应变软化的围岩特征曲线有限差分计算方法。高永涛^[9]归纳总结了屈服准则统一方程，基于线性软化及常量剪胀模型，得到了4种屈服准则下的围岩应力、位移、塑性区分布规律，探讨了各准则的差异性，并优选出了适用于某实际工程的屈服准则。王永红^[23]认为围岩峰后强度参数服从线性软化规律，基于考虑中主应力的广义张-朱强度准则建立了相应的隧道力学模型，并对围岩应力及变形进行了深入研究。

总之，可以看出目前常用的岩石剪胀及软化模型存在着参数多、适用性差等问题，同时以多适用于硬岩的理想弹塑性和弹脆性模型为主，然而在实际工程中，很多隧道围岩为中等强度岩体或软岩，因此必须同时考虑岩石应变软化及剪胀特性对隧道围岩力学行为的影响。为此，拟首先根据Lee^[21]提出的围岩应力释放法，并基于目前常用的M-C屈服准则，建立考虑围岩峰后阶段剪胀角非线性衰减和强度参数软化的力学模型，并验证其合理性；而后采用Alejano^[24]的围岩参数对本研究提出的软化模型的两种特例，即理想弹塑性和弹脆性模型进行对比验证。最后讨论目前常用的5种剪胀模型和3种本构模型对围岩变形及塑性区分布的影响规律，分析其差异性。

1 隧道围岩力学模型 1.1 基本假设

建立如图 1所示隧道围岩力学模型，并做如下假设：(1)隧道开挖断面为圆形；(2)围岩初始状态为弹性状态，且围岩为均质、连续的各向同性体，围岩屈服或破坏满足M-C屈服准则；(3)隧道长度远远大于其截面尺寸，可视为平面应变问题；(4)隧道为深埋隧道，处于静水应力场，且不考虑重力作用；(5)隧道内壁受均匀恒定的支护力。

如图 1所示，隧道开挖后，围岩自隧道中心向外依次为塑性区(这里进一步细分为软化区和残余区)和弹性区。分别用e、p和s代表弹性区、塑性软化区及塑性残余区的物理量。p₀为初始地应力；p_i为隧道内壁施加的支护力；r_p为塑性区半径，即弹塑交界面处的半径；r_s为塑性区残余半径，即塑性软化区与残余区交界面处的半径；r₀为隧道半径。

1.2 平衡方程与几何方程

根据上述假设，以隧道中心为极点建立极坐标，其平衡方程为：

(1)

式中σ_r，σ_θ分别为围岩径向和环向应力。

假设隧道围岩变形为小变形，其几何方程为：

(2)

式中，ε_r, ε_θ分别为围岩径向和环向应变；u为围岩径向位移。

围岩应力边界条件为：

(3)

1.3 围岩软化本构模型

围岩在峰值强度后将发生劣化，其强度参数如黏聚力c和内摩擦角φ会随着软化系数η的增加而降低，目前软化系数η定义通常为^[25]：

(4)

式中，ε_r^p, ε_θ^p分别为围岩塑性径向和环向应变。

目前研究多认为围岩强度在软化阶段会随软化系数呈现线性衰减，即：

(5)

式中，ω_peak为围岩峰值强度参数，如c_peak，φ_peak等；ω_res为围岩残余强度参数，如c_res，φ_res等；η^*为围岩由软化阶段进入残余阶段的临界软化系数，上述参数均可通过如图 2所示的单轴或三轴压缩试验曲线获得。

围岩本构模型对其变形计算结果影响很大，η^*的不同取值代表不同的本构模型。当η^*=0时，代表弹脆性模型(EBP)；当η^*→∞时，代表理想弹塑性模型(EPP)；当0 < η^* < ∞时，代表不同软化程度的应变软化模型(ESS)。上述3种本构模型对应的强度参数变化如图 2所示。

1.4 塑性流动法则

随着隧道的开挖，围岩应力不断释放，部分围岩逐步由弹性区进入塑性区。进入塑性区后围岩塑性应变采用塑性流动法则进行计算，即：

(6)

式中，dλ为与塑性变形过程相关的系数, 在塑性区内，dλ>0；在弹性区内dλ=0；Q为塑性势函数，当其与屈服函数相等时，称为关联流动法则；当其与屈服函数不相等时，称为非关联流动法则。岩土材料一般采用非关联流动法则，以此来修正因岩土材料破裂膨胀而产生的实际结果与采用关联流动法则计算结果的偏差。Q通常表示为：

(7)

式中，K_Ψ为剪胀系数，K_Ψ=(1+sinΨ)/(1-sinΨ)；Ψ为岩土材料剪胀角，其表示岩石在剪切过程中提交变化与应力变化的关系，可通过三轴压缩试验或直剪试验获得。

将式(7)代入式(6)，得到满足非关联流动法则的径向与环向塑性应变增量之间的关系为：

(8)

1.5 剪胀模型

由于岩土材料的剪胀特性导致传统塑性位势理论并不能很好地描述其变形机制^[26]，因此剪胀系数应由塑性位势理论和非关联流动法则得到，而目前研究多采用与M-C准则形式相同的塑性势函数下的剪胀系数K_Ψ，显然不甚合理^[5]。剪胀角通常被视为评估岩石体积应变和膨胀特性的重要参数，Cai^[2]认为剪胀角与围压σ₃和软化系数η相关，提出了考虑剪胀角衰减的剪胀模型，并基于11类岩石的剪胀数据证明其剪胀模型更符合实际情况，且计算参数少, 数学形式简单，更能反映深部岩体的剪胀特性，因此更适用于本研究。其表达式为：

(9)

式中，k为围岩影响系数，硬岩和软岩时分别取0.017和0.004 1；σ_ci为完整岩石单轴抗压强度。

2 围岩应力场与位移场 2.1 弹性区应力场与位移场

如图 1所示，隧道围岩弹性区的应力σ_r^e、σ_θ^e，应变ε_r^e、ε_θ^e和位移u_r^e分别为^[27]：

(10)

式中，σ_rp为弹塑性交界面处的径向应力；r为围岩任一点到隧道中心的距离；E, v分别为围岩弹性模量和泊松比。

2.2 塑性区应力场与位移场 2.2.1 塑性区应力场

假设隧道围岩应力满足σ_θ>σ_z>σ_r，且弹塑性交界面处的应力满足连续性假定，因此其同时满足式(10)及M-C准则, 即，由此可得弹塑性交界面处的径向应力σ_rp为：

(11)

式中，A=(1+sinφ_peak)/(1－sinφ_peak)；B=(2c_peakcosφ_peak)/(1－sinφ_peak)；c, φ分别为岩石黏聚力和内摩擦角。

相应的弹塑性交界面处的环向应力σ_θp可由M-C准则求得。

当围岩进入塑性区后，围岩强度参数随着临界软化系数而减小，利用解析法无法反映其变化过程。本节拟基于Lee^[21]提出的应力释放模型，求出考虑岩石应变软化和剪胀效应的隧道围岩应力与位移的半数值半解析解。已知弹塑性交界面处只有弹性应变，而无塑性应变。因此以弹塑性交界面为外边界以一定步长向内逐步递推，依次计算塑性区每一步长的应力应变。具体方法如图 3所示。首先将围岩塑性区划分为n个同心圆环，其中，弹塑性交界面处的径向和环向应力、半径、塑性径向及环向应变分别为σ_r(0)=σ_rp，σ_θ(0)=σ_θp，r₍₀₎=r_p，ε_r(0)^p=ε_θ(0)^p=0；隧道内壁径向应力和半径分别为σ_r(n)=p_i和r_(n)=r₀。第i环内外边界的径向和环向应力为σ_r(i)、σ_θ(i)和σ_r(i-1)、σ_θ(i-1)。

第i环由两个归一化半径的圆界定，即ρ_(i)和ρ_(i-1)，其表达式为：

(12)

可知第一环外边界ρ₍₀₎=1。

已知隧道内壁处支护力为p_i，弹塑性交界面径向应力σ_r(0)，设相邻圆环同边界应力差Δσ_r为固定值，则应力差表达式为：

(13)

因此，第i环径向应力为已知量。以圆环内边界为例，第i环的径向应力为：

(14)

当求出每一环的径向应力σ_r(i)后，代入M-C准则即可求出每一环的环向应力σ_θ(i)。每一环的环向应力受制于强度参数ω_(i)，而其又与软化系数η_(i)有关，故需首先求出软化系数的值。但软化系数求解的前提是环向应力已知，因此存在应力和软化系数迭代求解的矛盾。而当圆环数量足够多时，将使用上一环的临界软化系数η_(i-1)代替本环临界软化系数η_(i)而产生的误差足够小，由此便可求出每一环的环向应力。其表达式为：

(15)

2.2.2 塑性区位移场

在弹塑性交界面处，径向总应变ε_r(0)和环向总应变ε_θ(0)为：

(16)

塑性区隧道围岩应变由弹性应变分量和塑性应变分量组成，因此第i环的径向应变ε_r(i)和环向应变ε_θ(i)可表示为：

(17)

式中，ε_r(i)，ε_θ(i)分别为第i环的径向及环向应变；ε_r(i-1)，ε_θ(i-1)分别为第i-1环的径向及环向应变; Δε_r(i)^e，Δε_θ(i)^e分别为第i环的径向及环向弹性应变增量；Δε_r(i)^p，Δε_θ(i)^p分别为第i环的径向及环向塑性应变增量。

第i环径向及环向弹性应变增量Δε_r(i)^e和Δε_θ(i)^e分别为：

(18)

式中，Δσ_r(i)为相邻两圆环的径向应力差；Δσ_θ(i)为相邻两圆环的环向应力差，Δσ_θ(i)=σ_θ(i)－σ_θ(i-1)

若圆环数量n足够大，可采用差分形式将式(1)近似表示为：

(19)

将式(19)整理为相邻归一化半径相互关系的表示形式：

(20)

将式(2)几何方程整理，并将式(12)归一化半径代入，可得到变形协调方程：

(21)

同样将式(21)表示为差分形式，同时将总应变分为弹性应变和塑性应变：

(22)

将式(8)亦用差分形式表示为：

(23)

式中，。

由弹性本构方程可得第i环弹性应变差值为：

(24)

将式(23)、(24)分别代入式(22)，并化简可求出每一环塑性环向应变增量Δε_θ(i)^p：

(25)

将式(25)代入式(23)可求出每一环塑性径向应变增量Δε_r(i)^p，将式(18)，(23)和(25)代入式(18)可求得每一环总应变。

由式(20)，可从第一环逐步向内推可以求出每一环对应的归一化半径，以此可求出围岩塑性区峰值半径为：

(26)

式中ρ(n)为隧道内壁处归一化半径。

围岩洞壁处位移可表示为：

(27)

围岩软化系数η的定义见式(4)，由此可得第i环的临界软化系数η_(i)为：

(28)

上述求解流程如图 4所示。

3 模型验证

根据2.3节内容可知η^*的不同取值代表不同的本构模型，因此该方法可求解3种岩石本构模型，即理想弹塑性模型(EPP)、应变软化模型(ESS)和弹脆性模型(EBP)，下面分别采用上述3种模型进行对比分析，以验证本方法的正确性。计算过程中令圆环个数n=5 000，p_i=0 MPa，即无支护。

(1) ESS模型验证

选取Lee等^[21]的围岩参数见表 1，采用M-C屈服准则进行计算，且认为剪胀角为定值。由图 5所示计算结果可知由本研究求得的隧道内壁位移数值解、塑性区峰值和残余半径的数值解与Lee等结果十分一致，证明了软化模型的合理性。

(2) EBP模型验证

选取Alejano^[24]的围岩参数如表 2所示，采用M-C屈服准则进行计算，且认为剪胀角为定值。由上文可知，若令临界软化系数η^*=0，模型可转化为弹脆性模型。由图 6所示计算结果可知由本研究求得的围岩应力和洞壁位移数值解与Alejano^[24]结果十分一致，说明了弹脆性模型的合理性。

(3) EPP模型验证

选取Alejano^[24]的围岩参数见表 3，采用M-C屈服准则进行计算，且认为剪胀角为定值。为便于对比验证，令临界软化系数η^*=10⁵，模型可转化为理想弹塑性模型。由图 7所示计算结果可知由本研究求得的围岩应力和洞壁位移曲线数值解与Alejano等^[24]结果十分一致，证明了理想弹塑性模型的合理性。

4 不同剪胀模型计算结果对比

为了探究不同剪胀模型对围岩塑性区范围和位移的影响，分别考虑如下5种剪胀模型: (1)Hoek^[10]提出的将剪胀角ψ分别取为岩体内摩擦角φ的1/4，1/8和0倍; (2)Alejano^[28]提出的峰值剪胀角及剪胀角衰减模型；(3)Cai^[15]提出的剪胀模型。在M-C屈服准则基础上，考虑了围岩峰后强度参数衰减的软化效应，并将以上剪胀模型分别代入计算，比较各剪胀模型对围岩影响。本节采用Alejano^[28]的围岩参数，如表 4所示。

由图 8所示计算结果可知：(1)从围岩塑性区位移u分布曲线来看，随着剪胀角的增加，洞壁处的围岩位移u₀呈增长趋势。当Ψ分别为0，φ/8，φ/4，Alejano剪胀角及本研究剪胀角时，洞壁位移分别为0.250，0.299，0.364，0.270 m和0.297 m，其中Ψ=0时位移最小，而Ψ=φ/4时位移最大。与Ψ=0相比，当Ψ分别为φ/8，φ/4，Alejano剪胀角及本研究剪胀角时，洞壁位移分别增加了19.6%，45.7%，6.9%和19.1%。同时与Ψ=φ/8相比，本研究剪胀角计算位移随着距隧道内壁距离的减小而呈先略微增大而后减小的趋势，但整体上二者计算结果非常接近。(2)从剪胀角变化曲线来看，相对于Hoek^[10]的常剪胀角模型，Alejano剪胀角模型^[28]和本研究剪胀角模型可以定量反映围岩变形过程中剪胀角变化情况。Alejano剪胀角模型中，剪胀角在塑性区内持续减小，至洞壁处，剪胀角减小为0°。本研究剪胀角模型中剪胀角会表现出先略微增加后快速下降的现象，至洞壁处，剪胀角减小为0°，这与周辉^[29]试验结果相符。在弹塑性交界面处，本研究与Ψ=φ/8计算的峰值剪胀角非常接近，而Alejano^[28]计算的剪胀角偏小，约为φ的8.9%。

表 5为5种剪胀角模型对应的围岩塑性区半径。整体上，随着Ψ的增加，塑性区峰值半径r_p及残余半径r_s均呈现增长趋势，即剪胀角的增加导致围岩塑性区及残余区范围增加。由Ψ=0剪胀角模型计算得到的塑性区半径最小，而Ψ=φ/4剪胀角模型计算得到的塑性区半径最大。与Ψ=0相比，当Ψ分别为φ/8，φ/4，Alejano剪胀角及本研究剪胀角时，r_p分别增加1.2%，2.5%，0.7%和1.5%；而r_s分别增加5.9%，12.0%，3.1%和7.4%。可见本研究与Ψ=φ/8计算的r_p，r_s最为接近。

5 不同本构模型计算结果对比

为了探究不同本构模型对围岩应力、位移及塑性区范围影响，分别考虑以下3种本构模型: (1)EPP本构模型；(2)ESS本构模型；(3)EBP本构模型。在M-C屈服准则基础上，考虑围岩峰后剪胀效应，比较上述模型对围岩影响的差异性。本节采用Alejano^[28]的围岩参数，如表 4所示。

由图 9所示计算结果可知：(1)由图 9(a)可知在3种本构模型下，围岩位移曲线差别较大，整体来看EBP模型对应的位移最大，而EPP模型对应的位移最小。在洞壁处，由EPP，ESS和EBP等3种模型计算得到的位移分别为0.140，0.297 m和0.370 m。与EPP本构模型相比，ESS和EBP模型洞壁位移分别增加112.1%和164.3%，可见仅考虑围岩理想弹塑性变化状态，会低估围岩变形量，易出现隧道塌方等安全隐患。而仅考虑围岩弹脆性的EBP模型则会高估围岩变形量，导致支护结构工程量增加，造成浪费。(2)由图 9(b)可知在3种本构模型下，围岩径向应力σ_r和环向应力σ_θ峰后强度跌落状态不同。EPP模型下围岩峰后强度参数不变，峰后强度跌落缓慢，不存在残余区；ESS模型下围岩峰后强度参数会逐渐减小至残余值，峰后强度跌落稍快，会逐步进入残余区，在残余区内应力与EBP模型应力值相同；EBP模型下，围岩峰后强度参数立即跌落到残余值，峰后强度跌落极快。(3)由EPP，ESS和EBP等3种模型计算得到的围岩塑性区半径r_p分别为8.60，13.49 m和16.09 m，可以看出与EPP本构模型相比，ESS和EBP模型对应的r_p分别增加56.8%和87.1%，可见EPP模型会极大地低估围岩塑性区范围。总之，选择合适的本构模型可以更精确预测围岩应力、位移及塑性区范围。

6 结论

(1) 针对目前的岩石剪胀或软化模型无法很好地描述中、低强度岩石力学行为对隧道围岩力学特性影响的不足，基于围岩应力释放法及M-C屈服准则，建立了同时考虑岩石峰后阶段剪胀角非线性衰减和强度软化的隧道围岩力学模型。同时基于前人研究成果，对由本模型的两个特例即理想弹塑性和弹脆性模型计算得到的隧道围岩特征曲线和塑性区半径等与前人研究成果进行了对比，吻合很好，说明本模型的合理性。

(2) 对比常用的5种剪胀角模型对隧道围岩位移和塑性区半径的影响，同时给出了围岩峰后剪胀角变化规律。

(3) 对比了3种本构模型对隧道围岩应力、位移和塑性区半径的影响差异。其中，EPP本构模型会低估围岩塑性区范围及变形能力，易出现隧道塌陷的安全隐患。EBP本构模型会高估围岩塑性区范围及变形量，导致支护结构工程量增加，造成浪费。

(4) 本研究很好地考虑了岩石剪胀及软化特性对隧道围岩力学行为的影响，可应用于中、低强度隧道围岩力学行为分析。但是限于篇幅与实测数据的原因，对该模型合理性的进一步验证仍是今后的一个主要工作。