CN112989668A - 一种FDEM-Voronoi颗粒模型的能量数值计算方法 - Google Patents

一种FDEM-Voronoi颗粒模型的能量数值计算方法 Download PDF

Info

Publication number
CN112989668A
CN112989668A CN202110349110.9A CN202110349110A CN112989668A CN 112989668 A CN112989668 A CN 112989668A CN 202110349110 A CN202110349110 A CN 202110349110A CN 112989668 A CN112989668 A CN 112989668A
Authority
CN
China
Prior art keywords
energy
voronoi
joint
formula
calculation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110349110.9A
Other languages
English (en)
Other versions
CN112989668B (zh
Inventor
邱士利
李平
张世瑞
李邵军
江权
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan Institute of Rock and Soil Mechanics of CAS
South Central Minzu University
Original Assignee
Wuhan Institute of Rock and Soil Mechanics of CAS
South Central University for Nationalities
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan Institute of Rock and Soil Mechanics of CAS, South Central University for Nationalities filed Critical Wuhan Institute of Rock and Soil Mechanics of CAS
Priority to CN202110349110.9A priority Critical patent/CN112989668B/zh
Publication of CN112989668A publication Critical patent/CN112989668A/zh
Application granted granted Critical
Publication of CN112989668B publication Critical patent/CN112989668B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种FDEM‑Voronoi颗粒模型的能量数值计算方法，本发明目的是针对高应力环境下硬岩破裂模式表征和能量计算缺乏有效数值计算方法问题，在有限元‑离散元耦合数值方法(FDEM)框架下，提供一种FDEM‑Voronoi颗粒模型，表征深部硬岩隧洞或洞室开挖区应力诱发张拉和张剪断裂机制，并基于该模型发明了破裂过程中应变能、动能、断裂能、摩擦耗能等等能量成分的数值计算方法。可用于水电、采矿、交通和铁路等行业的深部地下隧洞(道)和硐室(群)高应力诱发破裂与岩爆能量释放过程中能量评估领域。

Description

一种FDEM-Voronoi颗粒模型的能量数值计算方法
技术领域
本发明涉及岩土工程中硬脆性岩体高应力动态破裂过程中应变能、动能、断裂能、摩擦耗能等能量成分的数值计算。
技术背景
高应力条件下硬脆性岩体发生破裂后原有应变能将被释放、消耗和转化，大部分释放能量消耗在岩体破裂过程，部分能量转化为应力波向岩体外传播，部分因摩擦作用转化为热能等。高储能释放后各能量份额的比例关系以及它们随着破裂演化过程的变化规律是决定这深部岩体力学行为的关键控制因素之一。如深部岩体高应力诱发岩爆灾害，其爆坑岩块的动力冲击性和弹射速度是与岩爆发生时应变能向动能转化的比例密切相关，同时高储能岩体断裂时辐射的应力波也是造成岩爆发生的重要因素之一。可见，高应力硬脆性岩体断裂过程中能量演化过程对研究和揭示深部高应力诱发破裂和岩爆灾害发生机理至关重要。
然而，目前工程界和学术界对高应力硬脆性岩体断裂过程中能量比例关系量化缺乏有效的计算方法和数值理论。虽然国内外学者提出了众多能量数值指标，如：能量释放率(ERR)、能量储存率(ESR)、局部能量释放密度(LERD)、模拟地层功(MGW)和局部能量释放率(LERR)、相对能量释放指数(RERI)等，但这些指标大多基于弹性、弹塑性数值分析方法，对断裂过程和能量集聚和释放过程均采用连续介质力学宏观表征的方法。然而，高应力诱发硬岩破裂本质是连续到非连续的转变过程，能量宏观表征已不足以满足这种两种状态转化过程的量化需求。这正是上述指标仅用于高应力破裂或岩爆灾害倾向性评估的重要原因，难于满足岩爆等灾害深入研究的需求，如无法对岩爆爆出岩块动能进行准确计算，无法估计动态断裂震源有效量级以及震动传播过程PPV量值及其分布等。
存在的问题主要包含三个方面：
(1)缺乏有效表征硬岩连续-非连续破裂过程表征的数值方法。目前，几个高级耦合数值计算方法，如近场动力学方法、有限元-离散元耦合方法、数值流形方法等正努力填补数值方法的不足。
(2)缺乏准确表征硬岩高应力诱发断裂机制的岩体结构表征方法。虽然已经有大量基于离散元颗粒流、三维粘结块体以及扩展有限元模型的岩体破裂结构与机制计算方法，但在表征开挖面附近低围压高法向荷载形成的张拉或张剪破裂机制时都表现出明显的不适应性。本发明中提出的FDEM-Voronoi颗粒模型正是独创的解决该问题的最新方法，是在有限元-离散元理论(FDEM)框架下双尺度Voronoi颗粒模型与节理粘结过程区模型耦合，形成了类似于实验室尺度穿晶、沿晶的断裂机制并推广到现场工程尺度，能够有效解决洞室开挖边界处低围压高法向荷载形成的张拉或张剪破裂机制问题。
(3)缺乏对硬岩破裂过程能量成分定量计算的有效算法。因此，本发明在FDEM-Voronoi颗粒模型基础上考虑双尺度Voronoi颗粒模型与节理粘结过程区模型耦合控制作用，建立了应变能、动能、断裂能、摩擦耗能等能量成分的数值计算方法。本发明的创新在于发明了一种准确表征高应力岩体破裂模式的岩体结构模型，并基于该模型在FDEM理论框架下计算能量成分的数值方法。
该发明对准确表征高应力硬岩破裂演化模式、准确计算破裂过程中各能量成分和深入揭示深部硬岩岩爆灾害发生机理等具有理论意义。对评估岩爆震动参数、块体弹射速度以及破裂损伤区范围、设计合理支护系统方案及防控措施和相关参数具有重要的工程应用价值。
发明内容
本发明目的是针对高应力环境下硬岩破裂模式表征和能量计算缺乏有效数值计算方法问题，在有限元-离散元耦合数值方法(FDEM)框架下，提供一种FDEM-Voronoi颗粒模型，表征深部硬岩隧洞或洞室开挖区应力诱发张拉和张剪断裂机制，并基于该模型发明了破裂过程中应变能、动能、断裂能、摩擦耗能等等能量成分的数值计算方法。
一种FDEM-Voronoi颗粒模型的能量数值计算方法，该方法仅适用于二维平面问题，包括如下计算步骤(以中国西部某地下硐室工程模拟为示例阐述)：
步骤1：FDEM-Voronoi数值模型建立：构建FDEM-Voronoi数值模型的网格拓扑关系，模型采用双尺度Voronoi颗粒模型与两套粘结节理模型耦合方法构建，即岩体先采用Voronoi颗粒模型做一级尺度上的剖分，进而在Voronoi颗粒内部划分二级尺度的三角形网格结构，Voronoi颗粒间节理，采用第一套节理参数，Voronoi颗粒内部三角形单元间节理，采用第二套节理参数，上述双尺度网格结构和双节理模型共同表征完整岩体力学行为；应力诱发破裂将依据双尺度网格结构构型和节理破裂准则进行破裂演化；
步骤2：两套粘结节理模型设定与参数标定：节理力学参数包括：内摩擦角、粘聚力、抗拉强度、I型断裂能和II型断裂能。两者关系要求Voronoi颗粒间节理单元力学参数，尤其是抗拉强度和粘聚力要低于Voronoi颗粒内部节理单元，从而调动张拉或张剪机制；
步骤3：计算结果初始数据库建立：运行FDEM-Voronoi数值模型，保存每个时间步的节点、三角形单元和节理单元计算结果到数据库，节点计算结果包括节点力、节点位移以及节点坐标；三角形单元计算结果包括应力和应变；节理单元计算结果包括张开位移和滑动位移；
步骤4：能量成分计算：分析岩体破裂过程中主要能量成分，建立能量平衡方程，并分别根据式(1)～(16)数值计算公式计算各能量成分，FDEM-Voronoi颗粒模型中，岩体破裂过程中主要能量成分包括岩体储存的应变能、岩体破裂能、岩块动能以及岩块间摩擦耗能，各能量成分计算步骤与方法如下：
步骤4.1：分析能量成分。FDEM-Voronoi颗粒模型能量平衡方程可表示为式(1)。
Wr＝Wk+Wp+Wf+Wd (1)
式中，Wr为应变能，Wk为动能，Wp为破裂能，Wf为摩擦耗能，Wd为阻尼耗能；
其中，阻尼耗能(Wd)是在FDEM理论框架下为了保证数值计算稳定性而导致的一种能量耗散，不存在物理意义，其在数值上近似等于动能的2倍，进而式(1)也可以写成式(2)，
Wr＝3Wk+Wp+Wf (2)
步骤4.2：计算应变能(Wr)。首先计算应变能增量(ΔU)，然后累加得到应变能(Wr)，应变能增量(ΔU)的计算公式可表示为式(3)，
Figure BDA0003001880990000031
式中，ΔU为应变能增量，i为三角形单元编号；Ai为编号为i的三角形单元面积；σ为当前时间步编号为i的三角形单元应力，σ'为前一时间步编号为i的三角形单元应力；Δε为编号为i的三角形单元当前时间步相较于前一时间步的应变增量，n为三角形单元总数；
应变能的计算公式可以表示为式(4)，
Figure BDA0003001880990000041
式中，t为当前时间步，T为总时间步，ΔU(t)为当前时间步的应变能增量，
步骤4.3：计算动能(Wk)：先以平均弹射速度划分区间后采用式(5)计算单个速度分区p内的所有三角形单元的动能
Figure BDA0003001880990000048
速度划分区间规则为每单位速度为一个速度区间，例如0～2m/s、2～4m/s为两个速度分区，
Figure BDA0003001880990000042
式中，p为弹射速度区间编号，
Figure BDA0003001880990000043
mp、vp分别为某一弹射速度区间p内的所有三角形单元的动能、质量和平均速度。mp、vp采用公式(6)计算，
Figure BDA0003001880990000044
式中，i为三角形单元编号，Ap为某一弹射速度区间p内所有三角形单元的面积，ρ为三角形单元密度，
Figure BDA0003001880990000047
分别为三角形单元节点x和y方向上的速度，n为某一弹射速度区间p内三角形单元总数，
动能的计算公式可以表示为式(7)，
Figure BDA0003001880990000045
式中，p为速度区间编号，m为速度区间划分总数；
步骤4.4：计算破裂能(Wp)：破裂能包括张拉破裂能(UI)、剪切破裂能(UII)和拉剪混合破裂能(UI&II)，分别用公式(8)、(9)和(10)进行计算，
Figure BDA0003001880990000046
Figure BDA0003001880990000051
Figure BDA0003001880990000052
式中，UI为张拉破裂能，UII为剪切破裂能，UI&II为拉剪混合破裂能，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；GfI、GfII分别为Voronoi颗粒内部节理单元I、II型断裂能；Gf'I、Gf'II分别为Voronoi颗粒间节理单元I、II型断裂能；hi、h'j分别为编号为i的Voronoi颗粒内部节理单元边长和编号为j的Voronoi颗粒间节理单元边长，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数，GfI、GfII、Gf'I和Gf'II是在FDEM理论框架下根据岩石力学试验确定的常数。
破裂能的计算公式可以表示为式(11)。
Wp＝UI+UII+UI&II (11)
步骤4.4：计算摩擦耗能(Wf)。首先计算摩擦耗能增量(ΔUfric)，然后累加得到摩擦耗能(Wf)。摩擦耗能增量(ΔUfric)计算为式(12)。
Figure BDA0003001880990000053
式中，ΔUfric为摩擦耗能增量，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；
Figure BDA0003001880990000054
为当前时间步Voronoi颗粒内部节理单元法向力，
Figure BDA0003001880990000055
为前一时间步Voronoi颗粒内部节理单元法向力；
Figure BDA0003001880990000056
为当前时间步Voronoi颗粒间节理单元法向力、
Figure BDA0003001880990000057
为前一时间步Voronoi颗粒间节理单元法向力；ΔSi为编号为i的Voronoi颗粒内部节理单元当前时间步与前一时间步的滑动位移增量，ΔSj为编号为j的Voronoi颗粒间节理单元当前时间步与前一时间步的滑动位移增量；μ1、μ2分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元摩擦系数，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数。
摩擦耗能的计算公式可以表示为式(13)。
Figure BDA0003001880990000061
式中，t为时间步，T为总时间步，ΔUfric(t)为当前时间步的摩擦耗能增量。
节理单元的法向力(Fn)计算力学模型如图8所示，计算公式为：
Fn＝Fycosθ-Fxsinθ (14)
式中，θ为理想接触面倾角(-90°<θ＜90°)；Fx、Fy分别为理想接触面的x和y方向上的节点等效力(方向规定为x轴正向为正，y轴正向为正)；理想接触面为随着数值计算时间进程变化的假想接触面。Fn的计算应考虑产生加速度的影响，这与节理面两侧的Fx、Fy的相对大小相关，采用以下判别条件分别计算。
当
Figure BDA0003001880990000062
时，
Figure BDA0003001880990000063
当
Figure BDA0003001880990000064
时，
Figure BDA0003001880990000065
式中，0、1、2、3为节理单元四个节点编号；
Figure BDA0003001880990000066
为编号为0的节点x方向上节点力，
Figure BDA0003001880990000067
为编号为0的节点y方向上节点力，其他节点以此类推。
本发明具有以下优势：
本发明的创新在于发明了一种准确表征高应力岩体破裂模式的岩体结构模型，并基于该模型在FDEM理论框架下发明了能量成分的数值计算方法。采用FDEM-Voronoi数值模型能够准确描述岩体的破裂过程和岩爆发生机制，从破裂机制上表征岩爆发生瞬时能量分配与转化过程，并能通过节点位移、块体单元应力和节理单元力学参数角度定量的计算能量成分，该方法具有严格的力学基础，计算精度远优于基于图像运动特征提取速度场而确定能量成分的方法。
附图说明
图1为本发明基于FDEM-Voronoi模型的能量数值计算流程图；
图2为本发明应变能(Wr)数值计算流程图；
图3为本发明动能(Wk)数值计算流程图；
图4为本发明破裂能(Wp)数值计算流程图；
图5为本发明摩擦耗能(Wf)数值计算流程图；
图6为本发明FDEM-Voronoi数值模型图；
图7为本发明数值模型示例应力云图；
图8为本发明FDEM粘结断裂节理单元的典型行为；
图9为本发明摩擦耗能数值计算力学模型示意图；
图10为本发明应变能-时间步曲线图；
图11为本发明每个速度分区内动能分布折线图；
图12为本发明破裂能-时间步曲线图；
图13为本发明摩擦耗能-时间步曲线；
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合实例对本发明作进一步的详细描述。此处实施示例仅用于说明和解释本发明，并不用于限定本发明。
种FDEM-Voronoi颗粒模型的能量数值计算方法，该方法仅适用于二维平面问题，包括如下计算步骤(以中国西部某地下硐室工程模拟为示例阐述)：
步骤1：FDEM-Voronoi数值模型建立。构建网格拓扑关系如图6所示的FDEM-Voronoi数值模型，模型采用双尺度Voronoi颗粒模型与两套粘结节理模型耦合方法构建，即岩体先采用Voronoi颗粒模型做一级尺度上的剖分，进而在Voronoi颗粒内部划分二级尺度的三角形网格结构。Voronoi颗粒间节理，如图6中Inter-CCE所示，采用第一套节理参数，Voronoi颗粒内部三角形单元间节理，如图6中Intra-CCE所示，采用第二套节理参数，上述双尺度网格结构和双节理模型共同表征完整岩体力学行为。应力诱发破裂将依据双尺度网格结构构型和节理破裂准则进行破裂演化。FDEM-Voronoi数值模型特有优势在于一级尺度Voronoi颗粒调动张拉或张剪机制，而另一级尺度三角形网格及节理单元刻画完整部分岩石的断裂过程。FDEM粘结节理单元的典型力学行为如图7所示。
步骤2：两套粘结节理模型(如图6中Intra-CCE和Inter-CCE)设定与参数标定。节理力学参数包括：内摩擦角、粘聚力、抗拉强度、I型断裂能和II型断裂能。两者关系要求Voronoi颗粒间节理单元力学参数(尤其是抗拉强度和粘聚力)要低于Voronoi颗粒内部节理单元，从而调动张拉或张剪机制。
步骤3：计算结果初始数据库建立。运行FDEM-Voronoi数值模型，保存每个时间步的节点、三角形单元和节理单元计算结果到数据库。节点计算结果包括节点力、节点位移以及节点坐标；三角形单元计算结果包括应力和应变；节理单元计算结果包括张开位移和滑动位移。
步骤4：能量成分计算。分析岩体破裂过程中主要能量成分，建立能量平衡方程，并分别根据式(1)～(16)数值计算公式计算各能量成分。FDEM-Voronoi颗粒模型中，岩体破裂过程中主要能量成分包括岩体储存的应变能、岩体破裂能、岩块动能以及岩块间摩擦耗能。各能量成分计算步骤与方法如下：
步骤4.1：分析能量成分。FDEM-Voronoi颗粒模型能量平衡方程可表示为式(1)。
Wr＝Wk+Wp+Wf+Wd (1)
式中，Wr为应变能，Wk为动能，Wp为破裂能，Wf为摩擦耗能，Wd为阻尼耗能。
其中，阻尼耗能(Wd)是在FDEM理论框架下为了保证数值计算稳定性而导致的一种能量耗散，不存在物理意义，其在数值上近似等于动能的2倍。进而式(1)也可以写成式(2)。
Wr＝3Wk+Wp+Wf (2)
步骤4.2：计算应变能(Wr)。首先计算应变能增量(ΔU)，然后累加得到应变能(Wr)。应变能增量(ΔU)的计算公式可表示为式(3)。
Figure BDA0003001880990000081
式中，ΔU为应变能增量，i为三角形单元编号；Ai为编号为i的三角形单元面积；σ为当前时间步编号为i的三角形单元应力，σ'为前一时间步编号为i的三角形单元应力；Δε为编号为i的三角形单元当前时间步相较于前一时间步的应变增量，n为三角形单元总数。
应变能的计算公式可以表示为式(4)。
Figure BDA0003001880990000091
式中，t为当前时间步，T为总时间步，ΔU(t)为当前时间步的应变能增量。
步骤4.3：计算动能(Wk)。先以平均弹射速度划分区间后采用式(5)计算单个速度分区p内的所有三角形单元的动能
Figure BDA0003001880990000098
速度划分区间规则为每单位速度为一个速度区间，例如0～2m/s、2～4m/s为两个速度分区。
Figure BDA0003001880990000092
式中，p为弹射速度区间编号，
Figure BDA0003001880990000093
mp、vp分别为某一弹射速度区间p内的所有三角形单元的动能、质量和平均速度。mp、vp采用公式(6)计算。
Figure BDA0003001880990000094
式中，i为三角形单元编号，Ap为某一弹射速度区间p内所有三角形单元的面积，ρ为三角形单元密度，
Figure BDA0003001880990000095
分别为三角形单元节点x和y方向上的速度，n为某一弹射速度区间p内三角形单元总数。
动能的计算公式可以表示为式(7)。
Figure BDA0003001880990000096
式中，p为速度区间编号，m为速度区间划分总数。
步骤4.4：计算破裂能(Wp)。破裂能包括张拉破裂能(UI)、剪切破裂能(UII)和拉剪混合破裂能(UI&II)，分别用公式(8)、(9)和(10)进行计算。
Figure BDA0003001880990000097
Figure BDA0003001880990000101
Figure BDA0003001880990000102
式中，UI为张拉破裂能，UII为剪切破裂能，UI&II为拉剪混合破裂能，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；GfI、GfII分别为Voronoi颗粒内部节理单元I、II型断裂能；Gf'I、Gf'II分别为Voronoi颗粒间节理单元I、II型断裂能；hi、h'j分别为编号为i的Voronoi颗粒内部节理单元边长和编号为j的Voronoi颗粒间节理单元边长，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数，GfI、GfII、Gf'I和Gf'II是在FDEM理论框架下根据岩石力学试验确定的常数。
破裂能的计算公式可以表示为式(11)。
Wp＝UI+UII+UI&II (11)
步骤4.4：计算摩擦耗能(Wf)。首先计算摩擦耗能增量(ΔUfric)，然后累加得到摩擦耗能(Wf)。摩擦耗能增量(ΔUfric)计算为式(12)。
Figure BDA0003001880990000103
式中，ΔUfric为摩擦耗能增量，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；
Figure BDA0003001880990000104
为当前时间步Voronoi颗粒内部节理单元法向力，
Figure BDA0003001880990000105
为前一时间步Voronoi颗粒内部节理单元法向力；
Figure BDA0003001880990000106
为当前时间步Voronoi颗粒间节理单元法向力、
Figure BDA0003001880990000107
为前一时间步Voronoi颗粒间节理单元法向力；ΔSi为编号为i的Voronoi颗粒内部节理单元当前时间步与前一时间步的滑动位移增量，ΔSj为编号为j的Voronoi颗粒间节理单元当前时间步与前一时间步的滑动位移增量；μ1、μ2分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元摩擦系数，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数。
摩擦耗能的计算公式可以表示为式(13)。
Figure BDA0003001880990000111
式中，t为时间步，T为总时间步，ΔUfric(t)为当前时间步的摩擦耗能增量。
节理单元的法向力(Fn)计算力学模型如图8所示，计算公式为：
Fn＝Fycosθ-Fxsinθ (14)
式中，θ为理想接触面倾角(-90°<θ＜90°)；Fx、Fy分别为理想接触面的x和y方向上的节点等效力(方向规定为x轴正向为正，y轴正向为正)；理想接触面为随着数值计算时间进程变化的假想接触面。Fn的计算应考虑产生加速度的影响，这与节理面两侧的Fx、Fy的相对大小相关，采用以下判别条件分别计算。
当
Figure BDA0003001880990000112
时，
Figure BDA0003001880990000113
当
Figure BDA0003001880990000114
时，
Figure BDA0003001880990000115
式中，0、1、2、3为节理单元四个节点编号；
Figure BDA0003001880990000116
为编号为0的节点x方向上节点力，
Figure BDA0003001880990000117
为编号为0的节点y方向上节点力，其他节点以此类推。
一种FDEM-Voronoi颗粒模型的能量数值计算方法，计算流程如图1所示，包括如下步骤：
步骤1：建立如图6所示的FDEM-Voronoi颗粒模型，以中国西部某地下硐室工程模拟为例，破坏情况如图7所示。计算并导出结果构建数据库，数据库文件为所有时间步的数据，结果数据输出频率为10000，每个时间步为4.5×10-8s。数据库文件应包括三角形单元应力应变、三角形单元节点力、三角形单元节点坐标位移、节理单元张开位移和滑动位移和节理单元节点坐标等结果数据。
步骤2：应变能(Wr)计算流程如图2所示。先根据三角形单元面积、前后时间步平均应力状态和应变增量采用式(3)计算应变能增量，随后进行应变能增量的累加，从而得到总应变能，并绘制总应变能-时间步曲线，典型计算结果如图10所示，可见峰值应变能约为500kJ。
步骤3：动能(Wk)计算流程如图3所示。首先进行破裂区单元速度分区，总共分为11个速度分区，分别为0～2、2～4、4～6、6～8、8～10、10～12、12～14、14～16、16～18、18～20、20～22m/s。针对每个速度分区内每个三角形单元，根据三角形单元质量、节点速度采用方程(5)～(6)计算每个分区内动能，最后将所有分区内的三角形单元动能进行累加求和，根据方程(7)确定总动能，典型计算结果如图11所示，总动能约为83kJ。
步骤4：破裂能(Wp)计算流程如图4所示，在识别完颗粒间和颗粒内节理的基础上，首先依据节理断裂模型区分张拉裂纹、剪切裂纹和混合裂纹，针对每种裂纹分别计算。采用式(8)计算每一条张拉裂纹的张拉破裂能并求和得到所有张拉裂纹的总张拉破裂能(UI)，采用式(9)计算每一条剪切裂纹的剪切破裂能并求和得到所有剪切裂纹的总剪切破裂能(UII)，采用式(10)计算每一条混合张剪裂纹的拉剪破裂能并求和得到所有混合张剪裂纹的总拉剪破裂能(UI&II)，最后将三类裂纹的各自总能量进行累加获取总破裂能(Wp)，采用式(11)计算。典型计算结果如图12所示，总破裂能约为285kJ。
步骤5：摩擦耗能(Wfric)计算流程如图5所示，在识别完颗粒间和颗粒内节理的基础上，基于节理法向力计算力学模型，首先根据节理单元4个节点力，采用式(14)、(15)、(16)计算节理单元法向力，然后结合相邻前后两个时间步滑动位移增量采用式(12)计算摩擦耗能增量，随后进行摩擦耗能增量的累加，根据式(13)从而得到总摩擦耗能，并绘制总摩擦耗能-时间步曲线，典型计算结果如图13所示，可见摩擦耗能约为140kJ。
至此，图1所示计算流程基本完成，岩体破裂过程的应变能、新生破裂的破裂能、岩块动能以及岩块间摩擦耗能四个能量成分计算完毕。
需要说明的是，本文中所描述的具体实施例仅是对本发明专利精神作举例说明。本发明专利所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种FDEM-Voronoi颗粒模型的能量数值计算方法，该方法仅适用于二维平面问题，其特征在于：包括如下计算步骤：
步骤1：FDEM-Voronoi数值模型建立：构建FDEM-Voronoi数值模型的网格拓扑关系，模型采用双尺度Voronoi颗粒模型与两套粘结节理模型耦合方法构建，即岩体先采用Voronoi颗粒模型做一级尺度上的剖分，进而在Voronoi颗粒内部划分二级尺度的三角形网格结构，Voronoi颗粒间节理，采用第一套节理参数，Voronoi颗粒内部三角形单元间节理，采用第二套节理参数，上述双尺度网格结构和双节理模型共同表征完整岩体力学行为；应力诱发破裂将依据双尺度网格结构构型和节理破裂准则进行破裂演化；
步骤2：两套粘结节理模型设定与参数标定：节理力学参数包括：内摩擦角、粘聚力、抗拉强度、I型断裂能和II型断裂能。两者关系要求Voronoi颗粒间节理单元力学参数，尤其是抗拉强度和粘聚力要低于Voronoi颗粒内部节理单元，从而调动张拉或张剪机制；
步骤3：计算结果初始数据库建立：运行FDEM-Voronoi数值模型，保存每个时间步的节点、三角形单元和节理单元计算结果到数据库，节点计算结果包括节点力、节点位移以及节点坐标；三角形单元计算结果包括应力和应变；节理单元计算结果包括张开位移和滑动位移；
步骤4：能量成分计算：分析岩体破裂过程中主要能量成分，建立能量平衡方程，并分别根据式(1)～(16)数值计算公式计算各能量成分，FDEM-Voronoi颗粒模型中，岩体破裂过程中主要能量成分包括岩体储存的应变能、岩体破裂能、岩块动能以及岩块间摩擦耗能。
2.如权利要求1所述的一种FDEM-Voronoi颗粒模型的能量数值计算方法，其特征在于，步骤(4)中各能量成分计算步骤与方法如下：
步骤4.1：分析能量成分。FDEM-Voronoi颗粒模型能量平衡方程可表示为式(1)。
Wr＝Wk+Wp+Wf+Wd (1)
式中，Wr为应变能，Wk为动能，Wp为破裂能，Wf为摩擦耗能，Wd为阻尼耗能；
其中，阻尼耗能(Wd)是在FDEM理论框架下为了保证数值计算稳定性而导致的一种能量耗散，不存在物理意义，其在数值上近似等于动能的2倍，进而式(1)也可以写成式(2)，
Wr＝3Wk+Wp+Wf (2)
步骤4.2：计算应变能(Wr)。首先计算应变能增量(ΔU)，然后累加得到应变能(Wr)，应变能增量(ΔU)的计算公式可表示为式(3)，
Figure FDA0003001880980000011
式中，ΔU为应变能增量，i为三角形单元编号；Ai为编号为i的三角形单元面积；σ为当前时间步编号为i的三角形单元应力，σ'为前一时间步编号为i的三角形单元应力；Δε为编号为i的三角形单元当前时间步相较于前一时间步的应变增量，n为三角形单元总数；
应变能的计算公式可以表示为式(4)，
Figure FDA0003001880980000021
式中，t为当前时间步，T为总时间步，ΔU(t)为当前时间步的应变能增量，
步骤4.3：计算动能(Wk)：先以平均弹射速度划分区间后采用式(5)计算单个速度分区p内的所有三角形单元的动能
Figure FDA0003001880980000028
速度划分区间规则为每单位速度为一个速度区间，例如0～2m/s、2～4m/s为两个速度分区，
Figure FDA0003001880980000022
式中，p为弹射速度区间编号，
Figure FDA0003001880980000023
mp、vp分别为某一弹射速度区间p内的所有三角形单元的动能、质量和平均速度。mp、vp采用公式(6)计算，
Figure FDA0003001880980000024
式中，i为三角形单元编号，Ap为某一弹射速度区间p内所有三角形单元的面积，ρ为三角形单元密度，
Figure FDA0003001880980000025
分别为三角形单元节点x和y方向上的速度，n为某一弹射速度区间p内三角形单元总数，
动能的计算公式可以表示为式(7)，
Figure FDA0003001880980000026
式中，p为速度区间编号，m为速度区间划分总数；
步骤4.4：计算破裂能(Wp)：破裂能包括张拉破裂能(UI)、剪切破裂能(UII)和拉剪混合破裂能(UI&II)，分别用公式(8)、(9)和(10)进行计算，
Figure FDA0003001880980000027
Figure FDA0003001880980000031
Figure FDA0003001880980000032
式中，UI为张拉破裂能，UII为剪切破裂能，UI&II为拉剪混合破裂能，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；GfI、GfII分别为Voronoi颗粒内部节理单元I、II型断裂能；Gf'I、Gf'II分别为Voronoi颗粒间节理单元I、II型断裂能；hi、h'j分别为编号为i的Voronoi颗粒内部节理单元边长和编号为j的Voronoi颗粒间节理单元边长，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数，GfI、GfII、Gf'I和Gf'II是在FDEM理论框架下根据岩石力学试验确定的常数。
破裂能的计算公式可以表示为式(11)。
Wp＝UI+UII+UI&II (11)
步骤4.4：计算摩擦耗能(Wf)。首先计算摩擦耗能增量(ΔUfric)，然后累加得到摩擦耗能(Wf)。摩擦耗能增量(ΔUfric)计算为式(12)。
Figure FDA0003001880980000033
式中，ΔUfric为摩擦耗能增量，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；
Figure FDA0003001880980000034
为当前时间步Voronoi颗粒内部节理单元法向力，
Figure FDA0003001880980000035
为前一时间步Voronoi颗粒内部节理单元法向力；
Figure FDA0003001880980000036
为当前时间步Voronoi颗粒间节理单元法向力、
Figure FDA0003001880980000037
为前一时间步Voronoi颗粒间节理单元法向力；ΔSi为编号为i的Voronoi颗粒内部节理单元当前时间步与前一时间步的滑动位移增量，ΔSj为编号为j的Voronoi颗粒间节理单元当前时间步与前一时间步的滑动位移增量；μ1、μ2分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元摩擦系数，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数。
摩擦耗能的计算公式可以表示为式(13)，
Figure FDA0003001880980000038
式中，t为时间步，T为总时间步，ΔUfric(t)为当前时间步的摩擦耗能增量。
节理单元的法向力(Fn)计算公式为：
Fn＝Fycosθ-Fxsinθ (14)
式中，θ为理想接触面倾角(-90°<θ＜90°)；Fx、Fy分别为理想接触面的x和y方向上的节点等效力(方向规定为x轴正向为正，y轴正向为正)；理想接触面为随着数值计算时间进程变化的假想接触面。Fn的计算应考虑产生加速度的影响，这与节理面两侧的Fx、Fy的相对大小相关，采用以下判别条件分别计算。
当
Figure FDA0003001880980000041
时，
Figure FDA0003001880980000042
当
Figure FDA0003001880980000043
时，
Figure FDA0003001880980000044
式中，0、1、2、3为节理单元四个节点编号；
Figure FDA0003001880980000045
为编号为0的节点x方向上节点力，
Figure FDA0003001880980000046
为编号为0的节点y方向上节点力，其他节点以此类推。
CN202110349110.9A 2021-03-31 2021-03-31 一种FDEM-Voronoi颗粒模型的能量数值计算方法 Active CN112989668B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110349110.9A CN112989668B (zh) 2021-03-31 2021-03-31 一种FDEM-Voronoi颗粒模型的能量数值计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110349110.9A CN112989668B (zh) 2021-03-31 2021-03-31 一种FDEM-Voronoi颗粒模型的能量数值计算方法

Publications (2)

Publication Number Publication Date
CN112989668A true CN112989668A (zh) 2021-06-18
CN112989668B CN112989668B (zh) 2022-05-17

Family

ID=76338671

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110349110.9A Active CN112989668B (zh) 2021-03-31 2021-03-31 一种FDEM-Voronoi颗粒模型的能量数值计算方法

Country Status (1)

Country Link
CN (1) CN112989668B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113927738A (zh) * 2021-10-15 2022-01-14 中国重汽集团青岛重工有限公司 一种搅拌罐螺旋叶片的优化设计方法及成型方法
CN114169180A (zh) * 2021-12-15 2022-03-11 东北大学 一种基于沿晶和穿晶接触面的脆性岩石的数值试验方法
CN114861401A (zh) * 2022-04-08 2022-08-05 武汉大学 一种层状岩体fdem数值模拟输入参数标定方法
CN115081302A (zh) * 2022-07-15 2022-09-20 中国矿业大学 支护构件与硐室围岩接触及相互作用的模拟方法和系统
CN115292990A (zh) * 2022-07-18 2022-11-04 南方科技大学 一种连续-非连续耦合的二维固体破裂模拟方法
CN115964902A (zh) * 2023-03-16 2023-04-14 太原理工大学 一种高应力大断面软岩底板破裂区确定及底鼓治理方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109063239A (zh) * 2018-06-19 2018-12-21 中国地质大学（武汉） 一种水热耦合的三维数值模拟方法
CN112362520A (zh) * 2020-10-30 2021-02-12 武汉大学 一种有限元-离散元耦合数值模拟程序（fdem）输入参数快速标定方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109063239A (zh) * 2018-06-19 2018-12-21 中国地质大学（武汉） 一种水热耦合的三维数值模拟方法
CN112362520A (zh) * 2020-10-30 2021-02-12 武汉大学 一种有限元-离散元耦合数值模拟程序（fdem）输入参数快速标定方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
ALESSANDRA INSANA等: "3D voronoi tessellation for the study of mechanical behavior of rocks at different scales", 《IACMAG 2021》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113927738A (zh) * 2021-10-15 2022-01-14 中国重汽集团青岛重工有限公司 一种搅拌罐螺旋叶片的优化设计方法及成型方法
CN113927738B (zh) * 2021-10-15 2023-02-17 中国重汽集团青岛重工有限公司 一种搅拌罐螺旋叶片的优化设计方法及成型方法
CN114169180A (zh) * 2021-12-15 2022-03-11 东北大学 一种基于沿晶和穿晶接触面的脆性岩石的数值试验方法
CN114861401A (zh) * 2022-04-08 2022-08-05 武汉大学 一种层状岩体fdem数值模拟输入参数标定方法
CN114861401B (zh) * 2022-04-08 2024-04-05 武汉大学 一种层状岩体fdem数值模拟输入参数标定方法
CN115081302A (zh) * 2022-07-15 2022-09-20 中国矿业大学 支护构件与硐室围岩接触及相互作用的模拟方法和系统
US11966672B2 (en) 2022-07-15 2024-04-23 China University Of Mining And Technology Method and system for simulating contact and interaction between support member and chamber surrounding rock mass
CN115292990A (zh) * 2022-07-18 2022-11-04 南方科技大学 一种连续-非连续耦合的二维固体破裂模拟方法
CN115964902A (zh) * 2023-03-16 2023-04-14 太原理工大学 一种高应力大断面软岩底板破裂区确定及底鼓治理方法

Also Published As

Publication number Publication date
CN112989668B (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN112989668B (zh) 一种FDEM-Voronoi颗粒模型的能量数值计算方法
Vazaios et al. Assessing fracturing mechanisms and evolution of excavation damaged zone of tunnels in interlocked rock masses at high stresses using a finite-discrete element approach
Feng et al. Excavation-induced deep hard rock fracturing: methodology and applications
Gudmundsson Emplacement and arrest of sheets and dykes in central volcanoes
Qu et al. Seismic cracking evolution for anti-seepage face slabs in concrete faced rockfill dams based on cohesive zone model in explicit SBFEM-FEM frame
Das et al. A mesh-free approach for fracture modelling of gravity dams under earthquake
Dai et al. Study on the mechanism of displacement mutation for jointed rock slopes during blasting excavation
CN114297864B (zh) 一种受陡缓倾角控制的碎裂松动岩体边坡稳定性分析方法
CN103399342B (zh) 一种基于岩体应变能的瞬态卸荷诱发振动预报方法
CN111553101B (zh) 一种隧道开挖上覆岩层开裂预报方法及围岩支护方法
Mirzabozorg et al. Dam-reservoir-massed foundation system and travelling wave along reservoir bottom
CN115859714A (zh) 一种基于fem-dem联合仿真的岩石爆破全过程模拟方法
CN111475978A (zh) 一种高位远程滑坡后破坏工程防护效果的预测方法
WO2023124664A1 (zh) 岩体稳定性极限分析方法
Ou et al. Numerical analysis of seepage flow characteristic of collapse column under the influence of mining
Koehn et al. Modelling of extension and dyking-induced collapse faults and fissures in rifts
Azarfar et al. A discussion on numerical modeling of fault for large open pit mines
CN114135288B (zh) 一种冲击地压煤层巷道高压水射流割缝卸压参数优化方法
Tang et al. Discrete element simulation for investigating fragmentation mechanism of hard rock under ultrasonic vibration loading
Sun et al. Research on rockburst proneness evaluation method of deep underground engineering based on multi-parameter criterion
Liu et al. A study of the dynamic damage characteristics of Rock-Like materials with different connectivity of concealed structural surfaces
Choi et al. Periodic monitoring of physical property changes in a concrete box-girder bridge
Zhou et al. Rock dilation and its effect on fracture transmissivity
CN103425837A (zh) 一种馆藏文物防轨道交通振动的设计方法
Wang et al. Dual bilinear cohesive zone model-based fluid-driven propagation of multiscale tensile and shear fractures in tight reservoir

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant