基于改进西原模型巷道围岩流变参数反演分析1)

2020-03-16高子璐荣传新

力学与实践 2020年1期

高子璐荣传新唐彬

(安徽理工大学土木建筑学院，安徽淮南232001)

对岩土工程的安全性和稳定性来说，研究岩石的蠕变特性有着重要的意义。传统本构模型采用线性元件组合，不能很好地描述岩石蠕变的三阶段，尤其是加速蠕变阶段。因此国内外学者[1-2]引进了非线性流变模型。目前建立非线性流变模型的方法有很多。一种是对传统线性元件进行非线性改进。如周家文等[3]构造出一个非线性函数，带入到广义Bingham 模型方程中，形成了一个新的非线性蠕变模型，但模型中的元件均为理想线性元件，无法完整表述岩石的全过程流变，因此不具备普遍适用性。另一种是建立应变与时间的非线性关系来描述加速蠕变阶段的特征。如徐卫亚等[4-5]将非线性元件与线性五元件模型串联，建立一个新的岩石非线性流变模型(河海模型)；张治亮等[6]将非线性流变与广义Kelvin 模型串联，建立一个岩石非线性黏弹塑性流变模型，用于反映岩石亚稳定蠕变和加速蠕变阶段的特征；蒋海飞等[7]采用幂函数与对数函数混合方程对加速阶段进行拟合，通过类比提出一个非线性黏性元件，将其与Burgers 模型串联，得出一个新的六元件非线性黏弹塑性蠕变模型，用于描述加速蠕变阶段蠕变特性；杨逾等[8]以Burgers 模型为基础，基于Lemaitre 原理建立了改进的Burgers 非线性蠕变损伤模型，并以砂岩为研究对象进行试验，验证了模型的合理性；张亮亮等[9]通过引入作用函数表征衰减蠕变阶段，提出一种非线性的改进模型，将其与宾汉姆模型串联，得出一个新的蠕变模型。

图1 改进西原体模型

在众多传统模型当中，西原体模型虽然对岩石稳定蠕变及其变形阶段的模拟取得了良好的效果，但是，它们在模拟加速蠕变时仍存在一定不足，需要通过对西原体模型的改进，建立出能更有效地反映加速蠕变阶段的模型。曹树刚等[10]对西原正夫模型进行改进，使之适用于描述岩石流变3 个阶段，尤其是非线性的加速蠕变变形；佘成学[11]在西原模型的基础上，引进岩石时效强度理论和损伤理论，建立了岩石非线性黏弹塑性蠕变模型。范翔宇等[12]基于煤岩蠕变力学实验对西原模型进行改进与分析。王军宝等[13]基于岩石非线性流变力学理论提出一个能够描述岩石加速蠕变的非线性黏滞体，与基本弹性体和塑性体进行组合建立一个新的四元件非线性黏弹塑性蠕变模型。田小朋等[14]以西原蠕变模型为基础添加含水元件建立与含水量相关的泥岩蠕变模型，同时选取实际油田中一井网为研究对象建立注采数值模型，通过控制注采比来调节地层中的含水量，验证了改进泥岩蠕变模型的正确性与合理性；杨秀荣等[15]考虑时效劣化和含水弱化对岩石蠕变参数的劣化效应，对西原模型进行参数修正，提出了一种新的蠕变损伤模型。本文通过引入一个新的非线性流变元件并将其串联在西原体模型上，完成对西原体模型的修正，建立了新的蠕变本构方程。该模型可以充分反映岩石蠕变的三个阶段，特别是加速蠕变阶段的特征，并且结构简单、参数少，能够更好地描述岩石三阶段的蠕变特性。

1 改进西原体模型的建立

根据岩石的流变力学特性以及对蒋海飞等[7]提出的加速蠕变拟合曲线进行修改，由此引入一个新的非线性流变元件，该流变元件应力应变关系为

式中，σ为应力，η为材料的黏性参数，˙ε为应变速率，n为蠕变指数，t为流变时间，σb为中间应力，σs为极限应力。

将新的非线性流变元件与西原体模型串联，组成六元件非线性黏弹塑性模型，即改进的西原体模型，如图1 所示，从左到右依次为元件1，元件2，元件3，元件4。六元件非线性黏弹塑性模型满足以下条件：

(1)当0＜σ≤σb时，新的非线性流变元件等同于牛顿体，改进的西原体模型等同于广义Kelvin 体。模型的状态方程为

式中，σ与ε分别为模型的总应力与总应变；σ1和σ2分别为元件1 和元件2 对应部分的应力；ε1和ε2分别为元件1 和元件2 对应部分的应变；E1和E2为岩石材料的弹性参数；η1为岩石黏性参数。

整理式(2)，得本构方程为

(2)当σb＜σ≤σs时，新的非线性流变元件等同于牛顿体，改进的西原体模型等同于Burgers 模型。模型的状态方程为

式中，σ与ε分别为模型的总应力与总应变；σ1，σ2和σ4分别为元件1，元件2 和元件4 对应部分的应力；ε1，ε2和ε4分别为元件1，元件2 和元件4 对应部分的应变；E1和E2为岩石材料的弹性参数；η1和η3为岩石黏性参数。

整理式(4)，得本构方程为

(3)当σ ＞σs时，新的非线性流变元件等同于非线性黏弹塑性体，改进的西原体模型状态方程为

整理式(6)，得本构方程为

由上面各式假定t= 0 时施加瞬时应力，则σ=σ0代入边界条件：t= 0，˙ε= 0，¨ε= 0，解微分方程，可求得非线性黏弹塑性模型一维状态下的蠕变方程

将式(8)两边分别对时间进行一次和二次求导，可以得到下列结果

(1) 当0＜σ≤σb时

(2)当σb＜σ≤σs时

(3)当σ ＞σs时

很显然，当σ≤σs时，˙ε恒大于0，¨ε ＜0，当t→∞时，则˙ε趋于某一非零恒定值。这表明，该模型可以很好地反映岩石的衰减和稳定蠕变过程。当σ ＞σs时，˙ε恒大于0，当t→∞时，则˙ε趋于某一非零恒定值。该模型可以很好地反映岩石蠕变全过程曲线中的衰减、稳定和加速3 个蠕变阶段。如果n ＞1，随着时间t的增长，¨ε可以小于、等于或大于0，对应于岩石的衰减、稳定和加速3 个阶段，且蠕变曲线的形态与n取值有关。

2 圆形巷道围岩径向位移

由于内部岩体被开挖的过程中产生了巷道变形，巷道半径为R0，围岩可以视为均质的各向同性体，巷道远处的边界应力为各项等压应力，采用改进的西原体模型描述圆形巷道的变形特性，计算模型如图2所示。

图2 圆形轴对称巷道计算模型

当0＜σ≤σb时，在等围压状态下，由圆形轴对称巷道(图2)径向位移表达式为

式中，P0为原岩应力；E和μ分别为岩体的弹性模量与泊松比。

由式(3)得到岩石的本构关系

对上式进行Laplace变换

对式(12)进行Laplace变换可以得到

将式(17)代入式(16)得到

对式(18)进行Laplace逆变换可以得到

当σb＜σ≤σs时，同理可得

对式(20)进行Laplace逆变换可以得到

当σ ＞σs时，同理可得

式中，M=tn+nlnt，A=E1+E2，B=E1E2，C=η2η3。

式(24)中的参数识别可采用优化反演法[16]。以实际位移和计算位移残差平方和最优值设计目标函数，通过不断优化位移残差的最小值得到蠕变参数值，由此得到基于改进的西原体模型下的圆形巷道围岩径向位移变化情况。

3 工程应用及参数反演

山西某煤矿中8 号采区的在建南回风巷道，采用全断面隧道掘进机(Tunneling Boring Machnie，简称TBM) 法进行巷道掘进。巷道轴向截面为圆形断面，断面直径为4.52 m，断面面积为16.12 m2，竖向地应力为P= 14.3 MPa。巷道地层所在的主要层位是砂岩层，主要组成成分是细粉砂岩。弹性模量范围为24.38～46.71 GPa，泊松比范围为0.121～0.36。假设巷道周边足够范围内都位于巨厚砂岩层中，在巷道掘进的影响范围内都不考虑围岩分层。

巷道开挖完成后，选取施工巷道靠近掘进面的断面进行围岩变形监测，断面1 沿巷道圆周选取5 个点，每点打孔锚固钢片，在钢片上留永久标记作为观测点，巷道收敛监测点布置如图3所示。

图3 巷道断面收敛监测测点布置示意图

取其中较有代表性的曲线作为目标蠕变曲线，进行参数反演。基于最小二乘法原理，对给定的数据点{(Xi,Yi)}(i=1,2,3,···,n)，在取定的函数W中，求拟合函数f(x) ∈W，使误差平方和∑[f(Xi-Yi)2]最小。假定W是式(24)提出的位移与时间的非线性关系式

已给出(W,T)的n对监测值(wi,ti)，要求待定参数(E1,E2,η1,η2,η3,n),通过使

平方和最小，得到模型相应参数值。由于现场条件限制，与5 测点相关的收敛数据均无法测得。通过巷道实际监测数据以第一次测量数据为基准，对后续测量数据取差进行记录，可得1-2，1-3，1-4监测点之间的收敛值，如图4 所示，将收敛值加上基准值即1-3监测点之间的实际距离d13，1-4 监测点之间的实际距离d14等。由监测布置图3 可见，监测点3到圆形巷道中心的距离d3近似由图3 可见，1-5点的距离近似为2d3，2-5 的距离近似等于1-4 测点之间的距离，由此得测点2到圆形巷道直径的竖向距离d26=d12·d14/2d13。根据勾股定理，可以算出监测点2 竖向投影点6 到监测点5 的水平距离d65，因此监测点2 到圆心的距离由此算出监测点2 和3 到圆心的距离值d2和d3。读取第一次测量数据值，计算初始时刻监测点2 和3 到巷道半径的距离，并以该值为初始值。将d2，d3的值分别减去巷道半径值和初始值，即测点2 和3 的径向位移值，如表1 所示。采用Origin 软件将测点2 和测点3 位移值的平均值进行非线性拟合，对蠕变参数进行初始识别并赋值，通过上述算法反复迭代反演，由此获得蠕变参数如表2 所示。将所得的蠕变参数代入到式(24)中，计算得到巷道围岩径向位移，如图5 所示。由图5 可见，实际监测结果与理论计算结果吻合较好，说明该模型可较好地模拟其围岩的蠕变特性。在t≤25 d 时，巷道径向位移速率随着时间的增长而不断增加；当t ＞25 d 时，巷道径向位移变化速率随时间的增长而趋于稳定。