CN116680832B - 一种考虑结构面蠕变特性的隧道围岩稳定性评价方法 - Google Patents

Abstract

本发明提出一种考虑结构面蠕变特性的隧道围岩稳定性评价方法，包括：S1.对隧道围岩区域离散化，构建隧道‑围岩结构的数值模型；S2.初始化所述数值模型；S3.预设蠕变时间，基于初始化后的所述数值模型，获取结构面节点的法向力；S4.基于所述法向力，获取岩体区域单元不平衡力；S5.判断所述岩体区域单元不平衡力，完成隧道围岩稳定性评价。本发明基于西原蠕变模型，采用有限差分方法，实现了含结构面的隧道围岩蠕变应力‑变形模拟，通过观察变形量值，分析隧道围岩稳定性；该方法能够模拟和评价岩质隧道等工程长期变形特征，为工程稳定性评价提供了重要思路。

Description

一种考虑结构面蠕变特性的隧道围岩稳定性评价方法

技术领域

本发明属于隧道围岩稳定性评价技术领域，尤其涉及一种考虑结构面蠕变特性的隧道围岩稳定性评价方法。

背景技术

在高应力作用下，由于岩体自身的蠕变特性，山岭隧道围岩运营期可能会发生较大的变形，造成隧道破坏。岩体由岩石和结构面组成，结构面的蠕变将显著影响岩体的变形和破坏，尤其在裂隙发育岩层中。因此，结构面的蠕变研究尤为重要。

目前大量的研究集中在岩石蠕变特性上，成果也较为丰富。早在1939年，Griggs(Griggs D.Creep of rocks.J Geol[M].1939:225-251)就已经对砂岩、石灰岩、页岩试件以及云母等矿物晶体进行大量蠕变试验，并提出了相应的蠕变本构模型。后来，国内外学者不断丰富了岩石蠕变模型。并将岩石蠕变本构模型写成有限元或者有限差分格式，方便开展隧道等岩石工程的数值评价。例如现有的商业软件，ANSYS，ABAQUS，FLAC3D等软件都有岩石蠕变计算模块。关于结构面蠕变本构模型近些年的研究越来越多，但相应的数值模拟的实现目前仍然是空白。例如，

CN202111150288.7公开了一种岩石非线性蠕变模型建立方法，包括进行室内岩石单轴压缩试验，获取岩石的平均压缩强度；采用分级加载方式进行室内岩石蠕变试验；绘制室内岩石蠕变试验获得的岩石全应变-时间曲线和岩石在不同应力水平下的分级应变-时间曲线分类；进行蠕变模型的元件种类识别，构建初始蠕变模型；拟合岩石加速蠕变阶段蠕变速率与对应的蠕变时间历程之间的非线性关系；最终获得岩石非线性蠕变模型的本构方程和蠕变方程。然而该方法未能实现岩石蠕变的数值仿真。

CN201910628858.5公开了一种改进Burgers岩石剪切蠕变模型的方法，包括以下步骤：S1、针对蠕变状态下岩石抗剪强度的时间差异性，构建一个由一个塑性元件和一个粘性元件并联而成描述岩石加速蠕变阶段力学特性的非线性元件；S2、采用广泛适用于流变领域的Kachanov定律，描述岩石加速蠕变阶段时效损伤函数D(t)；S3、基于时效损伤函数D(t)，确定抗剪强度τd的具体表达式τd(t)；S4、将残余强度τr引入抗剪强度函数τd(t)，构建考虑残余强度的修正抗剪强度函数τd(t)；S5、将修正抗剪强度函数τd(t)代入步骤S1的非线性元件中，确定元件本构方程γ(t)；S6、将元件本构方程γ(t)引入经典Burgers模型本构方程，获得精准反映全蠕变过程力学行为的改进Burgers模型本构方程γB(t)。该方法同样未能实现岩石蠕变的数值仿真。

CN202210179692.5公开了一种模拟巷道围岩蠕变失稳的离散元方法及系统，包括根据应力–应变曲线对压缩劈裂模型进行反演得到合理接触面和合理块体的参数值；根据岩石蠕变曲线对根据参数值建立的第一单轴蠕变数值模型进行反演得到合理块体蠕变参数值；根据岩石蠕变曲线对根据合理接触面和合理块体蠕变的参数值构建的第二单轴蠕变数值模型进行反演得到合理块体蠕变参数退化方程；对基于合理块体蠕变参数退化方程构建的仿真模型进行网格划分；对划分后的仿真模型进行模拟得到各网格的累计张拉裂隙长度和累计剪切裂隙长度；基于累计张拉裂隙长度和累计剪切裂隙长度确定待模拟巷道围岩蠕变失稳的范围。然而这种方法无法精确分析巷道或隧道的应力-变形情况，具有一定的局限性。

发明内容

为解决上述技术问题，本发明提出了一种考虑结构面蠕变特性的隧道围岩稳定性评价方法，基于西原蠕变模型，采用有限差分方法，实现了含结构面的隧道围岩蠕变应力-变形模拟，通过观察变形量值，分析隧道围岩稳定性；该方法能够模拟和评价岩质隧道等工程长期变形特征，为工程稳定性评价提供了重要思路。

为实现上述目的本发明提出了一种考虑结构面蠕变特性的隧道围岩稳定性评价方法，包括：

S1.对隧道围岩区域离散化，构建隧道-围岩结构的数值模型；

S2.初始化所述数值模型；

S3.预设蠕变时间，基于初始化后的所述数值模型，获取结构面节点的法向力；

S4.基于所述法向力，获取岩体区域单元不平衡力；

S5.判断所述岩体区域单元不平衡力，完成隧道围岩稳定性评价。

可选地，对隧道围岩区域离散化包括：

采用四边形单元对岩体区域离散，一维线单元对结构面离散，结构面单元的节点与岩体单元的节点共用。

可选地，初始化所述数值模型包括：

设置参数和边界条件，初始化所述数值模型的位移场和应力场；其中，所述参数包括：结构面弹性模量、和黏性系数、塑性极限、岩体的弹性模量和泊松比；所述边界条件包括：位移边界和力学边界；所述位移场包括：结构面上的位移和岩体的位移；所述应力场包括：结构面上的应力和岩体的应力。

可选地，获取岩体区域单元不平衡力包括：

基于所述法向力，获取结构面节点的切向力；

基于所述法向力和切向力，获取岩体单元节点的节点力；

基于所述节点力，获取所述岩体区域单元不平衡力。

可选地，所述法向力为：

F_n＝-k_n×u_n

其中，F_n为法向力，k_n为法向刚度，u_n为法向位移。

可选地，获取结构面节点的切向力包括：

判断所述法向力是否超过最大塑性应力，如果超过，则按第一预设方程，计算所述切向力，否则，按第二预设方程，计算所述切向力；

所述第一预设方程为：

其中，F′为切向力，X和Y按照下式计算：

u′为蠕变位移，u⁰为初始位移，F⁰为初始剪切力，k为弹簧元件弹性模量，Δt为蠕变时间，η为蠕变变形速率；

所述第二预设方程为：

其中，X′和Y′按照下式计算：

F_s为塑性元件受到的力。

可选地，获取岩体单元节点的节点力包括：将结构面节点的法向力和切向力添加到数值模型中相应的岩体单元节点上，获取新的所述节点力。

可选地，获取岩体区域单元不平衡力包括：

基于所述节点力，获取岩体区域节点速度；

基于所述节点速度，获取岩体区域节点的应变增量；

基于所述应变增量，获取岩体区域节点的应力增量和总应力；

基于所述总应力，获取所述岩体区域单元不平衡力。

可选地，岩体区域节点的所述速度为：

其中，为在t时步l节点在i方向上的速度，F_i ^l(t)为在t时步l节点在i方向的不平衡力分量，m^l为l节点的集中质量；

岩体区域节点的所述应变增量为：

其中，Δε_ij为应变增量，u_i,j为位移u_i在x_j方向上的偏导，u_j,i为位移u_j在x_i方向上的偏导；

岩体区域节点的所述应力增量为：

其中，Δσ_ij为应力增量，G为剪切弹性模量，ε_ij为单元ij的应变，E为弹性模量，ε_kk为单元上平均应力，δ_ij为Kronecker符号；

岩体区域节点的总应力为：

其中，σ_ij为总应力；

所述岩体区域单元不平衡力为：

其中，F^l为岩体区域单元不平衡力，n⁽¹⁾为单元体1面的法向向量，S⁽¹⁾为单元体1面的面积，n⁽²⁾为单元体2面的法向向量，S⁽²⁾为单元体2面的面积。

可选地，判断所述岩体区域单元不平衡力包括：

判断所述岩体区域单元不平衡力是否大于预设的允许容差，如果是，减少预设蠕变时间，返回S3，重新计算，否则，进行下一时间步计算；

判断下一时间步是否大于所述预设蠕变时间，如果不大于，返回S3，进行下一轮计算；如果大于，则结束计算，获得岩体和结构面的变形和应力；

基于计算的岩体和结构面的变形和应力，获取隧道周围围岩区域最大变形值；

判断所述最大变形值是否小于预设的允许值，如果小于所述允许值，则隧道是稳定的，否则，隧道不稳定。

与现有技术相比，本发明具有如下优点和技术效果：

(1)能够实现复杂裂隙岩体中隧道、边坡等工程性开挖长期变形模拟以及稳定性评价；(2)网格划分简单，划分速度快，避免了有限元计算对网格要求高的不足；(3)计算原理简单，编程容易，计算速度快。

附图说明

构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1为本发明实施例的岩体不平衡力计算示意图。

图2为本发明实施例的岩体区域和结构面单元划分示意图；

图3为本发明实施例的隧道围岩结构面蠕变有限差分求解流程示意图；

图4为本发明实施例的隧道几何模型和网格示意图；其中，(a)为隧道几何模型，(b)为网格化的隧道几何模型；

图5为本发明实施例的数值模拟与监测对比示意图；

图6为本发明实施例的竖向变形云图示意图；其中，(a)为蠕变1年，(b)为蠕变2年，(c)为蠕变3年；

图7为本发明实施例的五元件黏弹塑性力学模型示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

在介绍本实施例的具体步骤前，先介绍本实施例的基本假设和西原结构模型；

基本假设

(1)岩体应力-应变关系满足弹性本构关系；

(2)结构面应力-应变关系满足西原蠕变关系；

(3)岩体和结构面变形满足小变形假定。

西原结构模型

西原模型能够较全面地反映结构面的弹-黏弹-黏塑性蠕变变形特性。西原模型由一个弹簧、一个Kelvin体和一个黏塑性体串联组成，其力学模型如图7所示，u表示位移，k₁、k₂分别表示弹簧的弹性模量，η、η₂分别表示Kelvin体、黏塑性体中阻尼器的牛顿粘滞系数，σ₁、σ₂分别表示单独一个弹簧和Kelvin体的应力，σ_s表示塑性元件的屈服极限，ε₁、ε₂、ε₃分别表示单独一个弹簧、Kelvin体和黏塑性体的应变。本构方程分为两种情况：

当σ₀<σ_s时，模型退化为广义Kelvin模型：

当σ₀≥σ_s时：

可见，k₁、k₂的取值共同影响着最终蠕变变形的大小，k₁的取值直接影响着瞬时变形的大小，η决定了该蠕变曲线的斜率，即蠕变变形的速率，σ_s决定结构面是否进入塑性。

西原结构面模型的有限差分格式推导

结构面法向力按照弹性计算，即：

F_n＝-k_n×u_n (3)

结构面切向力按照西原模型进行计算，其位移增量为弹簧体、Kelvin体和黏塑性体的位移增量和，对于其中的弹簧体：

上式中，u₁ ⁰、u₁′分别表示弹簧体初始位移、经过时间步Δt后的位移，F₀、F′分别表示初始剪切力和经过时间步Δt后的剪切力。

对于Kelvin体，整理可得：

上式中，u₂ ⁰、u₂′分别表示Kelvin体初始位移、经过时间步Δt后的位移。

对于黏塑性体，整理可得：

上式中，u₃ ⁰、u₃′分别表示黏塑性体初始位移、经过时间步Δt后的位移。F_S分别表示黏性元件、塑性元件所受的力。

当结构面未进入塑性，即σ₀<σ_s时，结合式(4)、(5)，西原模型退化为广义Kelvin模型，剪切力的表达式为：

其中，

当结构面进入塑性，即σ₀≥σ_s时，结合式(4)、(5)、(6)，可得西原模型中剪切力的表达式：

其中，

岩体变形有限差分格式推导

岩体区域离散成四边形单元，分析时以节点为计算对象，力和质量集中在节点上，通过运动方程在时域内进行求解。节点运动方程可表示如下，

式中，为在t时步l节点在i方向上的速度，F_i ^l(t)为在t时步l节点在i方向的不平衡力分量，m^l为l节点的集中质量，在分析中，采用虚拟质量以保障数值稳定。采用中心差分来近似，得到：

f_i ^l(t)为非黏性阻尼，在求解静态问题时，使的系统的振动衰减至达到平衡状态，表示为，

式中，α为阻尼系数，默认值为0.8。

采用速率来求某一时步的单元应变增量，

由弹性本构方程求得应力增量Δσ_ij，

叠加各时步的应力增量，得到最新时步的总应力，

利用每个单元的应力，计算出每个节点的不平衡力F^l(t)，如图1所示，

式中，n为单元边的法向量，S为单元边面积或者长度，n⁽¹⁾为单元体1面的法向向量，S⁽¹⁾为单元体1面的面积，对于四边形单元来说指的是2面边长，n⁽²⁾为单元体2面的法向向量，S⁽²⁾为单元体2面的面积。如果所有节点的不平衡力小于容差，迭代结束。

岩体和结构面应力和变形耦合求解

岩体区域离散成四边形单元，结构面离散成一维线单元。结构面的单元与岩体区域单元共用节点，如图2所示。在计算中，需要将结构面上的不平衡力添加到对应的岩体节点上，以完成迭代求解，整个求解过程如图3所示。

如图3所示，本发明实施例提供了一种考虑结构面蠕变特性的隧道围岩稳定性评价方法，实现方法包括以下步骤：

S1.对模拟研究区域离散化。采用四边形单元对岩体区域离散，一维线单元对结构面离散，结构面单元的节点与岩体单元的节点共用；

S2.输入参数：输入结构面弹性模量k₁、k₂和黏性系数η、η₂，塑性极限σ_s以及岩体的弹性模量E和泊松比μ；

S3.输入模型的边界条件，包括位移边界和力学边界；

S4.初始化模型的位移场和应力场。包括结构面上的位移和应力，以及岩体的位移、速度和应力；

S5.获取时间步，按照公式(3)计算结构面节点的法向力或法向应力；

S6.判断法向应力σ是否超过最大塑性应力σ_max。如果超过，那么按照公式(7)计算结构面节点的切向力，如果没有超过，那么按照公式(10)计算结构面节点的切向力；

S7.将结构面节点的法向和切向力添加到相应的岩体单元节点上，得到新的节点结点力(不平衡力)；

S8.按照公式(14)计算岩体区域节点的速度；

S9.按照公式(16)计算岩体区域节点的应变增量；

S10.按照公式(17)和(18)计算岩体区域节点的应力增量和总应力；

S11.按照公式(19)计算岩体区域单元不平衡力；

S12.判断单元最大不平衡力是否大于允许容差，如果是，减少时间步长，返回S5，重新计算，如果不是，那么进行下一时间步计算。

S13.判断下一时间步是否大于设置的时间，如果不是，那么返回S5，进行下一轮计算；如果是则结束计算。

S14.观察隧道周围围岩区域最大变形值，是否小于规范要求的允许值，如果小于该值则，表明隧道在给定运行时间内是稳定的，如果大于该值，则表明隧道是不稳定的。

某次具体实施案例：

(1)某隧道蠕变计算

工程背景：某隧道为越岭隧道，位于宜都市红花套渔洋溪村和刘家村之间，进洞口端位于渔洋溪村，出洞口端位于刘家村，最大埋深约238m，水平长度为1011.2m，隧道属于中长隧道。隧道上方有一断层F1，在地表的延伸距离可达1km以上，如图4所示。隧道出现失稳的可能性大。

实施方法：

(1)建立数值模型。根据工程背景介绍，简化隧道几何结构，建立隧道-围岩结构模型，如图4(a)所示。模型左右边界长度近似取30倍跨度，上边界取自地表，下边界自隧道中心向下取25倍隧道高度，最终建立的模型整体尺寸为100m×338m×1m(宽×高×厚)。采用四边形单元对岩体区域离散，一维线单元对断层离散，结构面单元的节点与岩体单元的节点共用，共划分四边形单元8620个，节点32624个，其中，断层单元32个，如图4(b)所示；

(2)输入参数：输入结构面弹性模量k₁、k₂和黏性系数η、η₂，塑性极限σ_s以及岩体的弹性模量E和泊松比μ；

(3)输入模型的边界条件，包括位移边界和力学边界；

(4)施加自重和构造应力，模拟应力场，包括结构面上的位移和应力，以及岩体的位移、速度和应力；

式中，λ＝0.3.σ_z和σ_H分别为自重应力和构造应力，H为埋深。

(5)获取时间蠕变时间0.1年，按照公式(3)计算结构面节点的法向力或法向应力；

(6)判断法向应力σ是否超过最大塑性应力σmax。如果超过，那么按照公式(7)计算结构面节点的切向力，如果没有超过，那么按照公式(10)计算结构面节点的切向力；

(7)将结构面节点的法向和切向力添加到相应的岩体单元节点上，得到新的节点结点力(不平衡力)；

(8)按照公式(14)计算岩体区域节点的速度；

(9)按照公式(16)计算岩体区域节点的应变增量；

(10)按照公式(17)和(18)计算岩体区域节点的应力增量和总应力；

(11)按照公式(19)计算岩体区域单元不平衡力；

(12)判断单元最大不平衡力是否大于允许容差，如果是，减少时间步长，返回(5)，重新计算，如果不是，那么进行下一时间步计算。

(13)判断下一时间步是否大于设置的时间，如果不是，那么返回(5)，进行下一轮计算；如果是则结束计算；获得岩体和结构面的变形和应力。

这种计算会得到所有节围岩节点的变形，找出节点变形最大值，作为最大变形值。

(15)观察隧道周围围岩区域最大变形值，是否小于规范要求的允许值，如果小于该值则，表明隧道在给定运行时间内是稳定的，如果大于该值，则表明隧道是不稳定的。

(16)再重新输入蠕变时间分别为0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2和3年，重复(5)～(15)计算隧道围岩的变形。

(17)选取洞顶点作为监测点，获取开挖后的竖向位移变化，并将本方法计算结果与监测资料进行对比。图5绘制了数值计算与监测结果的对比。可见在监测日期内，该模拟结果在数值上与实测变形非常接近，趋势上与测点的位移实际过程线吻合，初期变形较快，随着时间增长变形量增加，变形速率趋于0mm/d，有稳定的趋势。图6为不同蠕变时间隧道变形结果；其中，图6(a)为蠕变1年，图6(b)为蠕变2年，图6(c)为蠕变3年。

以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (8)

1.一种考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，包括：

S1.对隧道围岩区域离散化，构建隧道-围岩结构的数值模型；

S2.初始化所述数值模型；

S3.预设蠕变时间，基于初始化后的所述数值模型，获取结构面节点的法向力；

S4.基于所述法向力，获取岩体区域单元不平衡力；

获取岩体区域单元不平衡力包括：

基于所述法向力，获取结构面节点的切向力；

基于所述法向力和切向力，获取岩体单元节点的节点力；

基于所述节点力，获取所述岩体区域单元不平衡力；

获取结构面节点的切向力包括：

判断所述法向力是否超过最大塑性应力，如果超过，则按第一预设方程，计算所述切向力，否则，按第二预设方程，计算所述切向力；

所述第一预设方程为：

其中，F′为切向力，X和Y按照下式计算：

u′为蠕变位移，u⁰为初始位移，F⁰为初始剪切力，k为弹簧元件弹性模量，Δt为蠕变时间，η为蠕变变形速率；

所述第二预设方程为：

其中，X′和Y′按照下式计算：

F_s为塑性元件受到的力；

S5.判断所述岩体区域单元不平衡力，完成隧道围岩稳定性评价。

2.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，对隧道围岩区域离散化包括：

采用四边形单元对岩体区域离散，一维线单元对结构面离散，结构面单元的节点与岩体单元的节点共用。

3.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，初始化所述数值模型包括：

设置参数和边界条件，初始化所述数值模型的位移场和应力场；其中，所述参数包括：结构面弹性模量、和黏性系数、塑性极限、岩体的弹性模量和泊松比；所述边界条件包括：位移边界和力学边界；所述位移场包括：结构面上的位移和岩体的位移；所述应力场包括：结构面上的应力和岩体的应力。

4.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，所述法向力为：

F_n＝-k_n×u_n

其中，F_n为法向力，k_n为法向刚度，u_n为法向位移。

5.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，获取岩体单元节点的节点力包括：将结构面节点的法向力和切向力添加到数值模型中相应的岩体单元节点上，获取新的所述节点力。

6.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，获取岩体区域单元不平衡力包括：

基于所述节点力，获取岩体区域的节点速度；

基于所述节点速度，获取岩体区域节点的应变增量；

基于所述应变增量，获取岩体区域节点的应力增量和总应力；

基于所述总应力，获取所述岩体区域单元不平衡力。

7.根据权利要求6所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，岩体区域的所述节点速度为：

其中，为在t时步l节点在i方向上的速度，F_i ^l(t)为在t时步l节点在i方向的不平衡力分量，m^l为l节点的集中质量；

岩体区域节点的所述应变增量为：

其中，Δε_ij为应变增量，u_i,j为位移u_i在x_j方向上的偏导，u_j,i为位移u_j在x_i方向上的偏导；

岩体区域节点的所述应力增量为：

其中，Δσ_ij为应力增量，G为剪切弹性模量，ε_ij为单元ij的应变，E为弹性模量，ε_kk为单元上平均应力，δ_ij为Kronecker符号；

岩体区域节点的总应力为：

其中，σ_ij为总应力；

所述岩体区域单元不平衡力为：

其中，F^l为岩体区域单元不平衡力，n⁽¹⁾为单元体1面的法向向量，S⁽¹⁾为单元体1面的面积，n⁽²⁾为单元体2面的法向向量，S⁽²⁾为单元体2面的面积。

8.根据权利要求1所述的考虑结构面蠕变特性的隧道围岩稳定性评价方法，其特征在于，判断所述岩体区域单元不平衡力包括：

判断所述岩体区域单元不平衡力是否大于预设的允许容差，如果是，减少预设蠕变时间，返回S3，重新计算，否则，进行下一时间步计算；

判断下一时间步是否大于所述预设蠕变时间，如果不大于，返回S3，进行下一轮计算；如果大于，则结束计算，获得岩体和结构面的变形和应力；

基于计算的岩体和结构面的变形和应力，获取隧道周围围岩区域最大变形值；

判断所述最大变形值是否小于预设的允许值，如果小于所述允许值，则隧道是稳定的，否则，隧道不稳定。