[发明专利]一种FDEM-Voronoi颗粒模型的能量数值计算方法

权利要求书

1.一种FDEM-Voronoi颗粒模型的能量数值计算方法，该方法仅适用于二维平面问题，其特征在于：包括如下计算步骤：

步骤1：FDEM-Voronoi数值模型建立：构建FDEM-Voronoi数值模型的网格拓扑关系，模型采用双尺度Voronoi颗粒模型与两套粘结节理模型耦合方法构建，即岩体先采用Voronoi颗粒模型做一级尺度上的剖分，进而在Voronoi颗粒内部划分二级尺度的三角形网格结构，Voronoi颗粒间节理，采用第一套节理参数，Voronoi颗粒内部三角形单元间节理，采用第二套节理参数，上述双尺度网格结构和双节理模型共同表征完整岩体力学行为；应力诱发破裂将依据双尺度网格结构构型和节理破裂准则进行破裂演化；

步骤2：两套粘结节理模型设定与参数标定：节理力学参数包括：内摩擦角、粘聚力、抗拉强度、I型断裂能和II型断裂能，两者关系要求Voronoi颗粒间节理单元力学参数要低于Voronoi 颗粒内部节理单元 ；

步骤3：计算结果初始数据库建立：运行FDEM-Voronoi数值模型，保存每个时间步的节点、三角形单元和节理单元计算结果到数据库，节点计算结果包括节点力、节点位移以及节点坐标；三角形单元计算结果包括应力和应变；节理单元计算结果包括张开位移和滑动位移；

步骤4：能量成分计算：分析岩体破裂过程中主要能量成分，建立能量平衡方程，并分别计算各能量成分，FDEM-Voronoi颗粒模型中，岩体破裂过程中主要能量成分包括岩体储存的应变能、岩体破裂能、岩块动能以及岩块间摩擦耗能。

2.如权利要求1所述的一种FDEM-Voronoi颗粒模型的能量数值计算方法，其特征在于，步骤(4)中各能量成分计算步骤与方法如下：

步骤4.1：分析能量成分：FDEM-Voronoi颗粒模型能量平衡方程可表示为式(1)：

W_r＝W_k+W_p+W_f+W_d (1)

式中，W_r为应变能，W_k为动能，W_p为破裂能，W_f为摩擦耗能，W_d为阻尼耗能；

其中，阻尼耗能(W_d)是在FDEM理论框架下为了保证数值计算稳定性而导致的一种能量耗散，不存在物理意义，其在数值上近似等于动能的2倍，进而式(1)也可以写成式(2)，

W_r＝3W_k+W_p+W_f (2)

步骤4.2：计算应变能(W_r)：首先计算应变能增量(ΔU)，然后累加得到应变能(W_r)，应变能增量(ΔU)的计算公式可表示为式(3)，

式中，ΔU为应变能增量，i为三角形单元编号；A_i为编号为i的三角形单元面积；σ为当前时间步编号为i的三角形单元应力，σ'为前一时间步编号为i的三角形单元应力；Δε为编号为i的三角形单元当前时间步相较于前一时间步的应变增量，n为三角形单元总数；

应变能的计算公式可以表示为式(4)，

式中，t为当前时间步，T为总时间步，ΔU(t)为当前时间步的应变能增量，

步骤4.3：计算动能(W_k)：先以平均弹射速度划分区间后采用式(5)计算单个速度分区p内的所有三角形单元的动能速度划分区间规则为每单位速度为一个速度区间，

式中，p为弹射速度区间编号，m_p、v_p分别为某一弹射速度区间p内的所有三角形单元的动能、质量和平均速度，m_p、v_p采用公式(6)计算，

式中，i为三角形单元编号，A_p为某一弹射速度区间p内所有三角形单元的面积，ρ为三角形单元密度，分别为三角形单元节点x和y方向上的速度，n为某一弹射速度区间p内三角形单元总数，

动能的计算公式可以表示为式(7)，

式中，p为速度区间编号，m为速度区间划分总数；

步骤4.4：计算破裂能(W_p)：破裂能包括张拉破裂能(U_I)、剪切破裂能(U_II)和拉剪混合破裂能(U_III)，分别用公式(8)、(9)和(10)进行计算，

式中，U_I为张拉破裂能，U_II为剪切破裂能，U_III为拉剪混合破裂能，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；Gf_I、Gf_II分别为Voronoi颗粒内部节理单元I、II型断裂能；Gf_I'、Gf'_II分别为Voronoi颗粒间节理单元I、II型断裂能；h_i、h'_j分别为编号为i的Voronoi颗粒内部节理单元边长和编号为j的Voronoi颗粒间节理单元边长，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数，Gf_I、Gf_II、Gf_I'和Gf'_II是在FDEM理论框架下根据岩石力学试验确定的常数；

破裂能的计算公式可以表示为式(11)：

W_p＝U_I+U_II+U_III (11)

步骤4.4：计算摩擦耗能(W_f)，首先计算摩擦耗能增量(ΔU_fric)，然后累加得到摩擦耗能(W_f)，摩擦耗能增量(ΔU_fric)计算为式(12)；

式中，ΔU_fric为摩擦耗能增量，i、j分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元编号；为当前时间步Voronoi颗粒内部节理单元法向力，为前一时间步Voronoi颗粒内部节理单元法向力；为当前时间步Voronoi颗粒间节理单元法向力、为前一时间步Voronoi颗粒间节理单元法向力；ΔS_i为编号为i的Voronoi颗粒内部节理单元当前时间步与前一时间步的滑动位移增量，ΔS_j为编号为j的Voronoi颗粒间节理单元当前时间步与前一时间步的滑动位移增量；μ₁、μ₂分别为Voronoi颗粒内部节理单元和Voronoi颗粒间节理单元摩擦系数，n为Voronoi颗粒内部节理单元总数，m为Voronoi颗粒间节理单元总数；

摩擦耗能的计算公式可以表示为式(13)，

式中，t为时间步，T为总时间步，ΔU_fric(t)为当前时间步的摩擦耗能增量；

节理单元的法向力(F_n)计算公式为：

F_n＝F_ycosθ-F_xsinθ (14)

式中，θ为理想接触面倾角，-90°＜θ＜900；F_x、F_y分别为理想接触面的x和y方向上的节点等效力，方向规定为x轴正向为正，y轴正向为正；理想接触面为随着数值计算时间进程变化的假想接触面，F_n的计算应考虑产生加速度的影响，这与节理面两侧的F_x、F_y的相对大小相关，采用以下判别条件分别计算：

当时，

当时，

式中，0、1、2、3为节理单元四个节点编号；为编号为0的节点x方向上节点力，为编号为0的节点y方向上节点力，其他节点以此类推。