CN113177248A - 隧道围岩破裂碎胀大变形失稳灾变过程2d-fdem数值模拟方法 - Google Patents

Abstract

本发明公开了一种隧道围岩破裂碎胀大变形失稳灾变过程2D‑FDEM数值模拟方法，模拟步骤包括建模及赋参、地应力施加、地应力平衡和核心材料软化卸荷。针对上述模拟步骤，本发明提出了：①一种节点力反向施加技术和坐标转换方法，使得任意倾角的地应力能够施加于具有任意形状边界的模型上；②一种新型迟滞阻尼模型及相应的临界阻尼系数取值方法，使得地应力加载阶段的动能能够被快速耗散掉；③一种核心材料逐步软化法以模拟隧道准静态开挖卸荷过程并能够反映掌子面径向支撑效应逐渐消失的过程。本发明提供的方法实现了隧道围岩破裂碎胀大变形失稳灾变全过程的模拟，能够捕获围岩裂隙场、应力场和位移场的孕育演化全过程。

Description

隧道围岩破裂碎胀大变形失稳灾变过程2D-FDEM数值模拟 方法

技术领域

本发明属于岩石力学与岩石工程领域，具体涉及隧道围岩破裂碎胀大变形失稳灾变过程的有限元-离散元耦合数值模拟方法(the combined Finite-Discrete ElementMethod,FDEM)。

背景技术

在深埋地层中开掘隧道后将造成围岩的切向应力升高而径向应力降低，即隧道开挖从力学角度而言是一个围岩切向加载而径向卸载的过程。当升高的切向应力超过岩体强度且隧道表面无支护结构时，围岩便不可避免发生破裂，隧道表面围岩的破裂使其承载失效造成集中切应力向深部转移进而使得破裂不断向围岩深处扩展，直至在深处达到极限平衡状态。隧道围岩的变形全过程包括弹塑性连续变形、围岩裂纹的萌生-扩展-贯通非连续变形和破碎块体间的剪切滑移及块体翻转大运动。岩体弹性恢复和块体间错位不咬合造成的空隙将使围岩发生体积膨胀现象，上述过程称为围岩破裂碎胀大变形失稳灾变全过程。

由于岩体材料介质的复杂性、赋存环境的复杂性和和围岩破裂碎胀大变形失稳灾变全过程的复杂性，理论解析显得无能为力。此外，室内模型试验和现场测试存在研究成本高昂、研究周期漫长、研究成果只对特定工程有效等不足，更重要的是，模型试验和现场测试难以获取岩石破裂碎胀大变形全过程，包括围岩裂隙场、应力场和位移场孕育演化全过程。因此，数值模拟成为最佳研究手段，包括连续性方法(如有限元法FEM、有限差分法FDM、边界元法BEM和无网格法MM等)、非连续性方法(如颗粒元法PFC、通用离散单元法UDEC和不连续变形法DDA等)和耦合方法(如数值流形法NMM和有限元-离散元耦合方法FDEM 等)。

其中，FDEM方法将材料划分为三角形单元，在两三角形单元边界插入初始无厚度的四边形节理单元，三角形单元只发生弹性连续变形，材料的塑性及断裂失效通过四边形节理单元来体现，当节理单元断裂失效后，其两侧的三角形单元由完整材料的粘结关系转化为块体间的接触关系。因此，FDEM能够模拟岩体从弹塑性连续变形至断裂失效非连续变形再至破碎块体间的接触挤压效应，能够反映真实的含有粗糙度和开度等几何特征的岩石裂纹，且计算效率高、参数标定较为便捷，能够适用于真实案例，如隧道开挖和边坡滑移的模拟。然而，现有的2D-FDEM隧道开挖模拟存在如下不足之处：①采用施加边界荷载的方法施加地应力只能加载水平与垂直地应力，难以适用于倾斜地应力的施加，且对任意形状的模型边界无能为力；②模型边界加载后将产生动能，只有在动能耗散完毕后才能进行隧道开挖的模拟，然而现有技术对动能的耗散效率是极低的；③核心材料的卸荷速率、卸荷曲线和临界动能具有任意性，难以确保准静态开挖卸荷过程。

因此，为解决上述存在的不足，本发明提出了隧道开挖围岩破裂碎胀大变形失稳灾变过程二维FDEM模拟方法，包括地应力施加技术、动能快速耗散技术及准静态开挖卸荷模拟技术。

发明内容

为了解决上述技术问题，本发明的目的是提出一种隧道围岩破裂碎胀大变形失稳灾变过程二维有限元-离散元耦合数值模拟(FDEM)方法，该方法能够对具有任意形状的模型边界施加任意倾角的地应力，动能耗散效率极高且能够实现准静态开挖卸荷过程，进而实现隧道围岩破裂碎胀大变形失稳灾变过程的分析模拟预测。

本发明所采用的技术方案如下：

隧道围岩破裂碎胀大变形失稳灾变过程2D-FDEM数值模拟方法，包括如下步骤：

步骤1.根据模拟需求，建立计算模型并进行网格划分，将模型分为三个区域：核心材料区、网格细化区和远场边界区，模型尺寸>20倍洞径，如图1所示。模型中包括三角形单元和四边形节理单元，分别对三角形单元和四边形节理单元赋予力学参数，对节点和模型边界赋予对应属性；

步骤2.根据所需地应力，采用下式计算出所有节点的节点力，将节点力反向施加至对应节点上，此时模型边界自由；

步骤3.在模型中施加步骤2的节点力后，模型将产生大量动能，需待系统动能耗散完毕后方可进行开挖的模拟，即为地应力加载平衡阶段；采用迟滞阻尼模型和临界阻尼系数将动能快速耗散掉；

步骤4.待步骤3中动能耗散完毕后，固定模型边界以保持住地应力，迟滞阻尼系数β设定为实际阻尼值：β＝(0.005～0.01)×β_c；将核心材料的弹性模量和粘滞阻尼分成至少5万次逐步软化；当上一步软化造成的围岩扰动达到平衡，即不再有新的裂纹产生后方可进行下一步的软化以确保准静态开挖卸荷过程。

进一步，所述步骤1中力学参数通过单轴压缩、直接拉伸和三轴压缩室内试验和相应的模拟试验进行标定得到。

进一步，所述步骤2中，节点力包括三角形三个节点的力(节点按逆时针排序)分别为：

式中：f_x0、f_y0、f_x1、f_y1、f_x2和f_y2分别为三角形三个节点(如图2所示)在x方向和y方向的节点力；x₀、y₀、x₁、y₁、x₂和y₂分别为三个节点的坐标；

σ_xx、σ_yy、τ_xy分别为所需地应力根据坐标换算后的x、y方向正应力与剪应力：

式中：σ_h、σ_v分别为水平和垂直地应力，以压应力为负；l₁、m₁、l₂、m₂分别为x轴与y轴在h轴、v轴的投影，即夹角余弦，如图3所示。

进一步，所述步骤3中，所述迟滞阻尼模型如下：

式中：m为节点质量，c为迟滞阻尼，β为迟滞阻尼系数，E为弹性模量，ρ为材料密度，F_in为节点内力，F_ex为外荷载，x为节点位移。

进一步，所述步骤3中，地应力加载带来的系统动能通过将迟滞阻尼系数β取为临界阻尼系数β_c来快速耗散，临界阻尼系数β_c可通过悬臂梁简谐振动模拟试验得到。

更进一步，采用任意阻尼系数β得到振幅-时间曲线，临界阻尼系数β_c根据下式得到：

式中：δ为在迟滞阻尼β作用下的相邻振幅比，即：

式中：u为第i个振幅，u_i+1为下一个振幅。

进一步，所述步骤4中，所述动能耗散完毕的判断方法为模型动能达到极小值，该极小值的范围为0.001-0.01kJ。

进一步，所述步骤4中，隧道掌子面对其前、后方一定范围内的岩体存在径向支撑效应，因此在二维平面应变条件下，可采用核心材料逐步软化技术模拟掌子面径向支撑效应逐渐消失的过程。将核心材料(如图1所示)的弹性模量和粘滞阻尼

h为网格尺寸) 分成至少5万次逐步软化。当上一步软化造成的围岩扰动达到平衡，即不再有新的裂纹产生后方可进行下一步的软化以确保准静态开挖卸荷过程。采用图4实现核心材料准静态开挖卸荷过程。

本发明的有益效果：

与现有隧道掘进围岩变形数值模拟相比，本发明实现了隧道围岩由弹塑性连续变形至裂纹萌生-扩展-交汇贯通非连续变形，再至破碎块体间的接触挤压、剪切滑移及块体自身的翻转大运动全过程模拟，能够捕获围岩裂隙场、应力场和位移场孕育演化全过程，实现具有任意模型边界和任意倾角地应力的施加，耗能效率高，能够确保准静态开挖卸荷模拟过程。

附图说明

图1为隧道开挖2D-FDEM模拟模型及网格划分示意图；

图2为是三角形单元节点编号及其节点力的示意图；

图3为地应力h-v坐标与x-y坐标的关系的示意图；

图4为隧道准静态开挖2D-FDEM模拟流程图；

图5为获取临界迟滞阻尼系数β_c的悬臂梁简谐振动模拟模型；

图6为迟滞阻尼系数β＝0.00106得到的P点位移-时间曲线；

图7为隧道开挖后围岩破裂碎胀大变形2D-FDEM模拟结果。

具体实施方式

以下结合附图对本发明涉及的隧道围岩破裂碎胀大变形失稳灾变过程2D-FDEM数值模拟方法进行详细说明。本发明的内容完全不限于此。

实施例

隧道开挖二维FDEM数值模拟方法如下：

步骤1.如图1所示，以圆形隧道及圆形模型边界为例，隧道直径3.0m，模型直径80m，网格细化区直径16m，网格细化区内网格尺寸h＝0.15m，模型边界处网格尺寸h＝8.0m。采用Gmsh软件建模并划分网格，导出Abaqus INP文件，而后将此文件导入至Y-GUI软件中，对三角形单元和四边形节理单元赋予力学参数，对节点和模型边界赋予对应属性，忽略模型的自重。经过单轴压缩、直接拉伸和三轴压缩试验标定后，材料力学参数如下表所示：

表1模拟参数表

步骤2.设定水平地应力σ_h＝-26.6MPa，倾角0度；垂直地应力σ_v＝-24.7MPa(负号代表压应力)，倾角90度。模型边界自由，采用下式获取所有节点的节点力，并将其反向施加至对应节点上：

式中：f_x0、f_y0、f_x1、f_y1、f_x2和f_y2分别为三角形三个节点在x方向和y方向的节点力(如图2所示)；x₀、y₀、x₁、y₁、x₂和y₂分别为三个节点的坐标；σ_xx、σ_yy、τ_xy分别为所需地应力根据坐标换算后的x、y方向正应力与剪应力：

式中：l₁、m₁、l₂、m₂分别为x轴与y轴在h轴、v轴的投影，即夹角余弦，如图3所示，在本实施例中：l₁＝1、m₁＝0、l₂＝0、m₂＝1。

步骤3.在模型中施加步骤2的节点力后，模型将产生大量动能，需待系统动能耗散完毕后方可进行开挖的模拟，即为地应力加载平衡阶段，采用下式迟滞阻尼模型将动能快速耗散完毕：

式中：m为节点质量，c为迟滞阻尼，β为迟滞阻尼系数，F_in为节点内力，F_ex为外荷载。在地应力加载阶段，迟滞阻尼系数β可取为临界阻尼β_c以快速耗散模型动能。临界阻尼系数β_c可通过悬臂梁简谐振动模拟试验得到。建立如图5所示的悬臂梁计算模型，采用任意阻尼系数β得到P点振幅-时间曲线，临界阻尼系数β_c根据下式得到：

式中：δ为在迟滞阻尼β作用下的相邻振幅比，即：

式中：u为第i个振幅，u_i+1为下一个振幅。在本实施例中，β取为0.00106，通过图5的悬臂梁简谐振动模拟得到相邻振幅比δ＝1.135，如图6所示，因此，根据上述两公式可得到临界迟滞阻尼系数β_c＝0.053。

步骤4.待步骤3中动能耗散完毕后(模型动能达到极小值，本实施设定为0.01kJ)，固定模型边界以保持住地应力。将核心材料(如图1所示)的弹性模量和粘滞阻尼μ分成至少5万次逐步软化。当上一步软化造成的围岩扰动达到平衡，即不再有新的裂纹产生后方可进行下一步的软化以确保准静态开挖卸荷过程。采用图4的流程实现核心材料准静态开挖卸荷过程。

通过上述步骤的模拟，得到了如图7所示的围岩破裂碎胀大变形模拟结果。根据本发明提出的方法，实现了具有任意模型边界和任意倾角地应力的施加，能够快速耗散地应力且实现准静态开挖卸荷软化过程，实现隧道围岩破裂碎胀大变形失稳灾变全过程的2D-FDEM数值模拟，捕获围岩裂隙场、应力场和位移场孕育演化全过程。

本发明所涉及的隧道围岩破裂碎胀大变形失稳灾变过程2D-FDEM数值模拟方法不仅仅局限在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

以上所述，仅为本发明较佳的具体实施方式，但本发明保护的范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内所做的任何修改，等同替换和改进等，均应包含在发明的保护范围之内。

Claims (8)

1.隧道围岩破裂碎胀大变形失稳灾变过程2D-FDEM数值模拟方法，其特征在于，包括如下步骤：

步骤1.建立计算模型并进行网格划分，将模型分为三个区域，包括核心材料区、网格细化区和远场边界区；模型中包括三角形单元和四边形节理单元，分别对三角形单元和四边形节理单元赋予力学参数，对节点和模型边界赋予对应属性；

步骤2.根据所需地应力，计算出所有节点的节点力，将节点力反向施加至对应节点上，此时模型边界自由；

步骤3.在模型中施加步骤2的节点力后，模型将产生大量动能，需待系统动能耗散完毕后方可进行开挖的模拟，即为地应力加载平衡阶段；采用迟滞阻尼模型和临界阻尼系数将动能快速耗散掉；

步骤4.待步骤3中动能耗散完毕后，固定模型边界以保持住地应力，迟滞阻尼系数设定为实际阻尼值；将核心材料的弹性模量和粘滞阻尼分成至少5万次逐步软化；当上一步软化造成的围岩扰动达到平衡，即不再有新的裂纹产生后方可进行下一步的软化以确保准静态开挖卸荷过程。

2.根据权利要求1所述的方法，其特征在于：所述步骤1中力学参数通过单轴压缩、直接拉伸和三轴压缩室内试验和相应的模拟试验进行标定得到。

3.根据权利要求1所述的方法，其特征在于：所述步骤2中，节点力包括三角形三个节点的力，三个节点的节点力按逆时针排序依次为：

式中：f_x0、f_y0、f_x1、f_y1、f_x2和f_y2分别为三角形三个节点在x方向和y方向的节点力；x₀、y₀、x₁、y₁、x₂和y₂分别为三个节点的坐标；σ_xx、σ_yy、τ_xy分别为所需地应力根据坐标换算后的x、y方向正应力与剪应力：

式中：σ_h、σ_v分别为水平和垂直地应力，以压应力为负；l₁、m₁、l₂、m₂分别为x轴与y轴在h轴、v轴的投影，即夹角余弦。

4.根据权利要求1所述的方法，其特征在于：所述步骤3中，所述迟滞阻尼模型如下：

式中：m为节点质量，c为迟滞阻尼，β为迟滞阻尼系数，E为弹性模量，ρ为材料密度，F_in为节点内力，F_ex为外荷载，x为节点位移。

5.根据权利要求1所述的方法，其特征在于：所述步骤3中，地应力加载带来的系统动能通过将迟滞阻尼系数β取为临界阻尼系数β_c来快速耗散，临界阻尼系数β_c可通过悬臂梁简谐振动模拟试验得到。

6.根据利要求5所述的方法，其特征在于：采用任意阻尼系数β得到振幅-时间曲线，临界阻尼系数β_c根据下式得到：

式中：δ为在迟滞阻尼β作用下的相邻振幅比，即：

式中：u为第i个振幅，u_i+1为下一个振幅。

7.根据权利要求1所述的方法，其特征在于：所述步骤4中，所述动能耗散完毕的判断方法为模型动能达到极小值，该极小值的范围为0.001-0.01kJ。

8.根据权利要求1所述的方法，其特征在于：所述步骤4中，利用下述方法以实现准静态开挖卸荷软化过程并反映掌子面径向支撑效应逐渐消失的过程：①计算模型总动能KE；②判定模型总动能KE是否小于临界动能CKE且当前软化次数n是否小于总软化次数N；③当满足KE<CKE且n<N条件时，n＝n+1；若不满足条件，n不变；④采用下式计算当前隧道材料的弹性模量和粘滞阻尼：

式中：E_n、E₀为隧道材料当前弹性模量和初始弹性模量；μ_n、μ₀为隧道材料当前粘滞阻尼和初始粘滞阻尼；power代表幂函数。

Images