基于改进西原模型的圆形隧洞位移解析解

2022-04-02刘咏梅

成都理工大学学报(自然科学版) 2022年2期

刘咏梅

(中铁二十三局集团第三工程有限公司，成都 611137)

深埋洞室的位移受围岩流变效应的影响，呈现出显著的时间相关性。洞室位移的时间特性对洞室的长期稳定性(如深部盐岩地下储气库[1]、核废料地下储存室[2]、地下工程支护结构[3])有显著的影响。因此，对地下洞室黏弹塑性位移的分析具有重要的理论与工程应用价值。

针对该问题，国内外学者进行了大量的研究。魏昂[4]采用摩尔库伦脆塑性非关联流动模型和宾海姆黏弹塑性模型，得到了任意截面隧道围岩位移的解析解；经纬等[5]考虑围岩流变、塑性扩容及中间主应力的影响，得到了基于蠕变终止轨迹的巷道围岩黏弹塑性解析解；夏晓等[6]以Levy-Mises本构关系为基础，选用基于D-P屈服准则的西原模型，推导了围岩黏弹塑性变形计算公式；A.R.Kargar[7]基于Burgers模型，推导了隧道在有支护及无支护条件下的黏弹塑性位移计算公式；P.Nomikos等[8]基于Burgers本构，得到了围岩具有流变效应的围岩支护特性曲线解析分析模型；卞跃威等[9]推导了考虑围岩塑性软化及体积膨胀效应的圆形隧洞黏弹塑性解。

近年来，基于分数阶微积分理论的元件模型被较多地应用于对岩石流变特性的描述。如：徐国文等[10]采用分数阶黏壶及NVPB模型分别取代传统西原模型黏弹性体中的牛顿黏壶及黏塑性体，提出了改进的西原模型；丁靖洋等[11]根据Weibull损伤理论，提出了一个考虑元件损伤的三元件蠕变模型，对岩盐的蠕变特性进行描述。分数阶流变模型具有元件少、参数意义明确等特性；同时，还可以完整地反映岩石流变的三个阶段。

目前，基于分数阶理论流变模型对隧洞围岩黏弹塑性分析的研究成果较少。鉴于此，本文首先根据Riemann-Liouville分数阶微积分理论，采用损伤分数阶Abel黏壶取代传统西原模型黏弹性体及黏塑性体中的牛顿黏弹，提出改进的西原模型；然后基于改进西原模型，推导出圆形隧道位移的黏弹塑性解；并以金瓶岩隧道为例，对解析解的有效性进行验证。

1 流变力学模型

1.1 损伤Abel黏壶

分数阶微积分是指运算阶次为任意实数或复数的微分或积分运算[8]。其中函数f(t)的Riemann-Liouvilleγ阶积分为

(1)

其微分形式为

(2)

Abel黏壶元件基于Riemann-Liouville微积分定义，其能够较好地描述材料的黏性特征

(3)

式中：ηγ为Abel黏壶的黏滞系数，ε为应变，t为时间。值得注意的是，当γ=0或1时，Abel黏壶元件就退化成理想Hooke弹性元件与牛顿体黏性元件。

岩石的流变过程伴随着微裂纹的产生与扩展。裂纹的出现会导致岩石内部出现损伤。采用如下的表达式对Abel体的损伤行为进行表征

ηγ=ηγ(1-S)

(4)

式中S为损伤变量, 0

S=1-e-ζ t

(5)

式中ζ岩石的材料参数。

结合公式(3)～(5)，得到损伤Abel黏壶的本构关系

(6)

通过对式(6)进行Laplace变换，得到其蠕变方程

(7)

1.2 损伤蠕变模型的建立

采用损伤Abel黏壶取代传统西原模型黏弹性体及黏塑性体中的黏性单元，提出改进的西原模型(图1)。

图1 流变模型Fig.1 Creep constitutive model

黏弹性状态方程为

(8)

黏弹塑性状态方程为

(9)

式中：σ、ε分别为总应力和总应变；σ1、σ2、σ3分别为1、2、3部分元件(图1)的应力；ε1、ε2、ε3分别为1、2、3部分元件的应变；E1、E2为弹性体和黏弹性体的弹性系数；γ、β为Abel损伤体的分数阶阶次；λ、α为Abel损伤体的损伤系数。

当σ0<σs：

a.弹性体的应力应变关系由式(8)中的第1式描述。

b.根据式(8)中的第2式可得[11]

(10)

c.当σ0≥σs，根据式(9)中的第3式可得[12]

(11)

综上所述，模型的本构方程为

(12)

(13)

1.3 蠕变模型验证

徐飞[13]通过室内分级加载蠕变试验，得到了千枚岩单轴压缩三阶段蠕变特性及长期强度(12 MPa)。本文以σ=13.5 MPa荷载下的试验结果为基础，对蠕变模型进行验证。

采用最小二乘法[12]对蠕变参数进行模拟，结果如表1所示。从图2可以看出，试验曲线与计算曲线吻合较好，均表现出明显的蠕变三阶段特性，即衰减蠕变(A点之前)、常速蠕变(A—B)及加速蠕变(B点之后)。可见，与同类流变元件组合模型相比，该模型可以采用较少的元件反映岩石三阶段流变过程。同时，可以定量地表征岩石流变全过程的损伤特性以及损伤对流变行为的影响。

表1 千枚岩流变参数Table 1 Creep parameters of phyllite

图2 流变模型验证Fig.2 Validation of the creep constitutive model

2 隧洞位移解析解

2.1 基本假定

本文解析解的推导前提为：①隧道断面为圆形且处于平面应变状态，同时原岩应力场为静水应力场(σ0)；②围岩的蠕变在隧道瞬时弹塑性变形完成后发生；③岩体应力以受压为正。

2.2 隧洞围岩瞬时弹塑性应力场

深埋隧道围岩特性具有显著的非线性特性，因此，采用Hoek-Brown屈服准则对围岩蠕变变形之前的弹塑性应力场进行描述

(14)

式中：σ1为最大主应力；σ3为最小主应力；σc为岩石的单轴抗压强度；m、s为与岩体性质相关的参数，即

(15)

(16)

式中：IGSI为地质强度因子；M2为开挖扰动系数。

根据已有的隧洞围岩瞬时弹塑性应力场推导结果，隧洞开挖后，围岩的弹性与塑性应力为

(17)

(18)

(19)

(20)

(21)

(22)

2.3 蠕变位移解

2.3.1 黏弹性区

在围岩的黏弹性区，塑性体流变元件未被触发。岩体的流变特性由元件1与元件2的组合描述。位移与应变的关系为

(23)

式中：εθ、εr分别为围岩的切向与径向应变；u为围岩的径向位移。

求解式(23)得到

(24)

2.3.2 黏塑性区

黏塑性区位移与应变的关系为

(25)

(26)

式中φ为岩体内摩擦角。

将式(21)和(22)代入式(25)可得

(27)

式中:

(28)

(29)

弹塑性截面位移边界条件为

(30)

[kψ(3-sinψ)-(3+sinψ)]

(31)

式中:

(32)

3 模型验证

金瓶岩隧道位于成兰铁路松潘段，隧道单洞长度为 12 773 m，隧道最大埋深为791 m，开挖断面的最大跨度和高度分别为13.7 m与11.5 m(图3)。隧道全长范围内广泛分布着软弱千枚岩，开挖过程中，频繁遭遇喷层开裂、初支钢架扭曲破坏、二衬开裂等大变形灾害。岩石单轴抗压强度为15 MPa，IGSI=50，mi=9。岩石的流变参数见表1。原岩地应力场简化为静水压力场，即σ0=13.97 MPa。

图3 隧道开挖断面Fig.3 Tunnel excavation section

隧道采用三台阶法进行开挖。开挖过程中，对拱顶沉降与拱腰收敛进行了监测。计算曲线与测试曲线如图4所示。可以看出，计算曲线的变化趋势及量值与拱顶、拱腰位移曲线的平均值较为吻合，验证了本文理论推导结果的可靠性。需要说明的是，对于监测曲线而言，围岩净空变形随时间增长迅速，且未呈现收敛趋势；随着掌子面推进，拱顶沉降和洞周收敛呈现明显的三台阶波动增长。而理论推导未考虑围岩的开挖过程，因此曲线未呈现阶跃式上升的特征。

图4 位移曲线Fig.4 Displacement curves

4 模型参数敏感性分析

4.1 分数阶阶次的影响

不同分数阶阶次γ、β时，围岩的位移曲线如图5所示。可以看出，随着阶次数值的增大，围岩的蠕变变形量值增加。这是因为对于Abel损伤体而言，分数阶阶次越高，元件的黏性特性就越明显。随着时间的增加，同一时间不同分数阶阶次对应的围岩位移之间的差值也增加。同时，对于任一时间点，不同分数阶阶次对应的围岩位移之间的差值随着阶次的增加而减小。可见，分数阶阶次与围岩的变形之间具有较为复杂的非线性关系。

图5 分数阶阶次的影响Fig.5 Influence of fractional order

4.2 损伤参数的影响

不同损伤系数(α及λ)时，围岩位移的变化如图6所示。可以看出，围岩的位移随着损伤系数量值的增加而增大。这是因为对于Abel损伤体而言，损伤系数量值越大，元件的黏性特性就越明显。随着时间的增加，同一时间不同损伤次数值对应的围岩位移之间的差值也增加。同时，对于任一时间点，不同损伤次数对应的围岩位移之间的差值随着阶次的增加而减小(特别是黏弹性体的损伤系数λ)。可见，黏塑性体的损伤对围岩变形的影响更为显著。