在弹塑性力学的理论中,对围岩弹塑性的分析计算通常采用圆形断面进行简化; 但是在实际工程中,许多隧道的形状都是非圆形,因此对于非标准圆形断面形状的隧道的弹塑性就变得十分必要.朱大勇等人基于复变函数提出了非圆形巷道的粘弹性解析解^[1]; 蔡晓鸿等人根据复变函数保角变换提出了二向不等围压和内压作用下椭圆形洞室的围岩应力计算方法^[2]; 非圆形隧道的弹塑性解答大多是都是基于复变函数进行计算的,计算规则复杂,不易应用.同时在分析隧道周边围岩的稳定性时,应用最多的是Mohr-Coulomb强度准则,虽然其在土体计算时吻合较好,但不能反映岩体及岩石的非线性破坏特性,而Hoek-Brown则可以描述岩体的非线性变形,充分考虑岩体裂隙及结构面信息,压力释放等因素,在描述岩体的变形方面具有明显的优势.虽然有部分学者采用Hoek-Brown强度准则^[3-6],但都是基于圆形隧道断面来进行分析计算的,几乎没有将Hoek-Brown准则应用于非圆形断面隧道的弹塑性解答的; 同时绝大部分专家学者在采取岩石的峰后应力 - 应变的本构关系时也都是基于理想弹塑性模型很少采用应变软化模型来进行分析计算,而应变软化模型更能代表实际宏观观测的结果^[7-8].针对以上问题,本文提出一种基于Hoek-Brown准则采用应变软化模型来对非圆形隧道断面进行弹塑性分析的简便计算方法,通过将非圆形隧道断面等效成圆形隧道断面,在MATLAB中编程来实现Hoek-Brown强度准则的应用; 对比分析在采用外接圆圆心法,大小半径之和法以及当量半径法进行等代时计算得到的弹塑性解答与现场工程实际的弹塑性解答的吻合程度,得出在工程实际中应用Hoek-Brown强度准则计算围岩稳定性时的最佳等代方法.

Hoek-Brown屈服准则由于其可以很好的反应岩石和岩体的非线性破坏特征以及结构面、应力状态对强度的影响,并且能够解释最大主应力和最小主应力对强度的影响,可延用到破碎岩体和各向异性岩体等情况,是为数不多的非线性准则之一,从而得到了广泛应用,但没有考虑应变软化区域的影响,与实际情况有偏差.本文基于Hoek-Brown准则的应变软化模型,推导圆形隧道塑性区半径的表达式,并应用于非圆形隧道.首先将模型简化为圆形断面,隧道初始应力场均匀分布.围岩发生屈服状态时满足Hoek-Brown准则:

σ₁=σ₃+σ_ci(m_b(σ₁)/(σ_ci)+s)^a(1)

式中:σ₁、σ₃为第一和第三主应力,σ_ci为岩块的单轴抗压强度; m_b、s、a均为材料参数,可以通过下面公式计算得出

{m_b=m_iexp((GSI-100)/(28-14D))

s=exp((GSI-100)/(9-3D))

a=1/2+1/6(e^-(GSI)/(15)-e^-(20)/3)(2)

其中:GSI为地质强度指标; D为扰动参数.应变软化模型的一个主要特点是失效准则和塑性势函数不仅取决于应力张量σ_ij,而且取决于软化参数η; 即失效准则方程里的各个参数不是一个定值,而是随着材料的变形阶段不同而变化的,是关于软化参数η的函数.

因此本文将失效准则定义为

f(σ_r,σ_θ,η)=0(3)

在弹性阶段,软化参数η=0,对应于弹性阶段失效准则f(σ_r,σ_θ,0)=0; 当进入软化阶段以后(0<η<η^*),对应于软化阶段失效准则f(σ_r,σ_θ,η)=0; 进入残余阶段后(η>η^*),则对应残余阶段失效准则f(σ_r,σ_θ,η^*)=0.软化参数η^*通过实验测定.采用极坐标形式的应变定义软化参数η为

η=ε^p_θ-ε^p_r(4)

应变软化材料的本构方程通过塑性增量理论得到,而塑性增量理论可以通过塑性势函数得到,定义塑性势函数为

g(σ_r,σ_θ,η)=0(5)

根据塑性增量理论

式中:λ^.是一个未知的塑性乘数,式(6)～(7)是塑性阶段的本构方程,称为流动法则.目前在进行围岩弹塑性分析时认为岩体体积保持不变,而岩体在变形时体积会增大,产生扩容,对于岩体的扩容的考虑主要通过采用相关联流动法则和将扩容通过剪胀角来体现的^[9],因此本文采用相关联流动法则建立采用Hoek-Brown准则应变软化模型的塑性势函数:

g(σ_r,σ_θ,η)=f(σ_r,σ_θ,η)

=σ_θ-σ_r-σ_ci(η)(m(η)(σ_r)/(σ_ci(η))+s(η))^a(η)

=0(8)

其中:σ_θ为环向应力; σ_r为径向应力; σ_ci为岩石单轴抗压强度; m为量纲为一的材料常数,用来反应岩石材料的软硬程度; s反应岩体的破碎程度,参数a用来调整失效包络曲线.在本文的所将建立的模型中,需要获取的参数为弹性阶段的弹性模量E和泊松比ν; 软化阶段参数,包括损伤模量M和剪胀角ψ和残余阶段参数η^*.定义峰值参数σ_cp,s_p,a_p,m_p,残余参数σ_cr,s_r,a^{_</sup>r,m_r,,其中损伤模量M根据岩体质量和侧限应力水平求得,计算得到:}

M=-ω·E(9)

其中,ω取决于峰值地质强度指标GSI^peak和侧限水平应力,可有下式计算得到

当(σ₃)/((s_p)^1/2σ_ci)≥0.1时,

残余阶段参数即临界软化参数η^*,根据公式(4)可得

η^*=ε^p₁-ε^p₃

其中,ε^p₁,ε^p₃分别表示最大最小塑性应变,经过推导可得

ε^p₁=(σ^p₁-σ^r₁)/E-(σ^p₁-σ^r₁)/M(12)

ε^p₃=-(11+sinψ)/(21-sinψ)ε^p₁(13)

其中,σ^p₁和σ^r₁分别为在某一侧压下σ₁的峰值和残余值,ψ为剪胀角,本文采用HOEK和Brown提出的方法确定剪胀角即:

ψ=(5GSI-125)/(1 000)φ(14)

式中,φ为摩擦角.

假设地基材料为各向同性材料,且符合Hoek-Brown失效准则; 此时地基可视为半无限的弹塑性体.设原岩应力为σ₀,在这样的地基材料中开挖一个硐室半径为 b,长度较长的地下圆形隧道结构; 并且基于这种圆形结构对围岩进行弹塑性分析,将其看作平面应变问题来求解.围岩的弹塑性分析,则就归结为弹塑性力学中的平衡微分方程、几何方程和物理方程这三大微分方程的求解,由于塑性区的解答对于围岩稳定性的影响最显著所以本文仅给出塑性区的应力场应变场及位移场的解答,计算示意图如图1所示.

图1 围岩弹塑性区域划分示意图
Fig.1 Schematic diagram of elasto-plastic zone division of surrounding rock

图2 塑性区域微元圆环划分示意图
Fig.2 Schematic diagram of ring division of plastic zone elements

对于轴对称问题,围岩在极坐标下的平衡方程可表示为

(dσ_r)/(dρ)+(σ_r-σ_θ)/ρ=0(15)

将上式变换应用到第i圈,列出平衡方程,根据差分法,将微分方程用微元的代数方程代替,求出近似的数值解答如下

dσ_r=σ_r(i)-σ_r(i-1)(16)

dρ=ρ_i-ρ_i-1(17)

在塑性区,岩石材料已经屈服,其应力状态始终满足失效准则方程(8),则第i圈的应力状态满足

σ_r(i)-σ_θ(i)=

-σ_ci(η_i-1)(m(η_i-1)(σ_r(i)^-)/(σ_ci(η_i-1))+s(η_i-1))^a(η)(18)

σ_r(i)^-=(σ_r(i)+σ_r(i-1))/2)(19)

将式(15)中的半径ρ用圆环半径表示,即

ρ=(ρ_(i)+ρ_(i-1))/2(20)

将公式(16)～(20)带入到公式(15)中,可得

ρ_(i)=

(2σ_ci(η_i-1)(m(η_i-1)(σ_r(i)^-)/(σ_ci(η_i-1))+s(η_i-1))^a(η)+Δσ_r)/(2σ_ci(η_i-1)(m(η_i-1)(σ_r(i)^-)/(σ_ci(η_i-1))+s(η_i-1))^a(η)-Δσ_r)ρ_(i-1)(21)

当i=1时,则有ρ₍₀₎=R_p

当i=n时,则有ρ_(n)=b

将初始值和已知条件带入式(21),经过n次迭代和变换可求得塑性区半径R_p的表达式为

R_p=

b/(∏ⁿ_i-1(2σ_ci(η_i-1)((σ_r(i)^-)/(σ_ci(η_i-1))+s(η_i-1))+Δσ^a(η)_r)/(2σ_ci(η_i-1)(m(η_i-1)(σ_r(i)^-)/(σ_ci(η_i-1))+s(η_i-1))-Δσ^a(η)_r))(22)

对于平面应变问题,对几何方程进行迭代消去位移变量得到

(dε_θ)/(dρ)+(ε_θ-ε_r)/ρ=0(23)

应变由弹性应变和塑性应变两部分组成,将其表示为

ε_θ=ε^p_θ+ε^e_θ,ε^{_</sup>r=εp_r+εe_r(24)}

将其24带入到23中得到

(dε^p_θ)/(dρ)+(ε^p_θ-ε^p_r)/ρ=-(dε^e_θ)/(dρ)-(ε^e_θ)/(dρ)-(ε^e_θ-ε^e_r)/ρ(25)

经过推导计算可得第i个圆环与第i-1个微元之间的应力表达和应变表达,表达式如下:

应力场的表达式为

{σ_r(i)=σ_r(i-1)+Δσ_r

σ_θ(i)=σ_θ(i-1)+H(σ_r(i)^-,η_(i-1))(26)

Δσ_r=(p_i-σ_R)/n=σ_r(i)-σ_r(i-1)(27)

H(σ_ci(η)^-,η_(i-1))=σ_ci(η)(m(η)(σ_r)/(σ_ci(η))+s(η))^a(η)(28)

应变场的表达式为

{ε_r(i)=ε_r(i)+Δε^e_r(i)+Δε^p_r(i)

ε_θ(i)=ε_θ(i)+Δε^e_θ(i)+Δε^p_θ(i)(29)

{Δε^p_θ(i)=

(-(Δε^e_θ(i))/(Δρ_(i))-(1+ν)/E(H(σ_r(i)^-,η_(i-1)))/(ρ_(i)^-)-1/(ρ_(i)^-)(ε^p_θ(i-1)-ε^p_r(i-1)))/(1/(Δρ_(i))+(1+k_(i-1))/(ρ_(i)^-))

Δε^p_r(i)=-k_(i-1)Δε^p_θ(i)(30)

{Δε^e_θ(i)=(1+ν)/E[-νΔσ_r(i)+(1-ν)Δσ_θ(i)

Δε^e_r(i)=(1+ν)/E[(1-ν)Δσ_r(i)-νΔσ_θ(i)(31)

{Δσ_r(i)=σ_r(i)-σ_r(i-1)

Δσ_θ(i)=σ_θ(i)-σ_θ(i-1)(32)

k(σ_r,η)=

1+m(η)·a(η)·[m(η)(σ_r)/(σ_ci(η))+s(η)^a(η)-1](33)

Δε^p_θ(i)、Δε^p_r(i),Δε^e_θ(i)、Δε^e_r(i)分别为塑性应变增量和弹性应变增量,E为弹性模量,ν为泊松比,Δρ为相邻两个微元的半径差.

经过n次迭代循环,即可求出弹塑性区域内的应力场和应变场,同时位移场也可由下列方程迭代求得:

u=ε_θ·ρ(34)

将上述表达式在,MATLAB环境下编程求解即可得到具体的应力应变及位移.

(1)本文采用考虑应变软化特性的Hoek-Brown强度准则,在相关联流动法则下推导得出了圆形隧道应力场、应变场及塑性区半径的表达式,实现Hoek-Brown在实际工程中的简便应用,充分考虑了因塑性而引起的岩体的扩容、岩体的结构面信息及扰动情况,使计算结果更加接近围岩的实际变形.

(2)以重庆兴隆隧道为依托结合现场声波测试结果得到围岩松动圈范围,并通过公式等效计算围岩塑性区范围,通过MATLAB自编程分析在选取不同等代圆计算方法下计算得到的塑性区半径及塑性区范围发现,当隧道形状为多心圆时,应用Hoek-Brown准则分析其塑性区范围时,选取当量半径法分析得到的结果与现场实测结果相差最为接近,误差在5.3%左右.

(3)通过计算在不同等代方法下隧道围岩的位移及软化区域发现,采用外接圆法计算的得到的位移最大,采用高度和跨度的四分之一的等效方法计算得到的位移最小,而采取当量半径法计算得到的位移居中,同时发现计算半径选取的越大计算得到的软化区域范围越小,半径越小,岩体的软化特性越明显.