QUASI-STATE-BASED PERIDYNAMICS METHOD FOR THE WHOLE PROCESS OF ROCK BRITTLE FAILURE
-
摘要: 近场动力学(peridynamics, PD)作为一种新兴方法, 正展现出成为一种分析复杂岩体力学响应方法的巨大潜力. 然而, 键基近场动力学与连续介质力学联系薄弱, 缺乏对岩石多类型断裂模式和峰后行为的处理能力, 而且接触-摩擦模型严重简化. 为此, 提出了一种准态基近场动力学方法. 首先, 通过基于键型区分的力密度计算方法, 实现了键的应力张量计算, 加深了其与连续介质力学的联系; 然后, 提出了一种通用式的断裂准则引入思路, 推导了一种适于当前理论的单直线型虚拟裂缝模型, 实现了对岩石多类型脆性破坏的全过程模拟; 最后, 通过粒子网格孪生算法, 实现了不同介质状态的物质点和孪生网格的同步处理, 通过势函数和库伦摩擦定律, 实现了对各类接触情形和接触程度的准确判断和表征, 并增加了考虑摩擦系数的摩擦力密度计算. 已有算例表明, 模拟结果与实验结果的吻合程度很高. 这使得完善后的近场动力学的近场动力学处理非连续介质形成-演化过程的能力得到了有效提升, 扩充了原有的理论框架和适用范围.Abstract: Peridynamics (PD), as an emerging methodology, is demonstrating substantial potential as a means to analyze the mechanical responses of complex rock formations. However, the bond-based peridynamics method exhibits weak connections with continuum mechanics (CM). It cannot handle multiple fracture modes and post-peak behavior in rocks and significantly oversimplifies the contact-friction model. The calculation system of state-based peridynamics is more complex and less stable. Therefore, this paper proposed a quasi-state-based peridynamics method. This paper first achieves a simultaneous calculation of bond force density and stress tensor by employing a force density calculation method that distinguishes based on bond types, thereby deepening its association with continuum mechanics. Subsequently, a general criterion introduction approach (e.g., maximum tensile stress criterion and Mohr-Coulomb criterion) is provided, leading to the derivation of a single linear fictitious crack model suitable for the current theoretical framework. This enables a comprehensive simulation of multiple types of brittle failure in rock materials. Finally, using a particle mesh twin algorithm, binding fracture criterion and the discretization process of twin meshes, simultaneous calculations of material points and twin meshes in different medium states are accomplished. Through potential function and Coulomb's friction law, accurate assessments and characterizations of contact situations and degrees are achieved, incorporating friction force density calculations considering the friction coefficient. Case studies demonstrate a good agreement between simulation results and experimental outcomes. This signifies a significantly improved capability of the refined peridynamics in handling the formation-evolution processes of discontinuous media, thereby substantially broadening its original theoretical framework and applicability.
-
引 言
在土木、建筑和采矿等岩石工程中, 固体材料和结构的变形-损伤-破坏问题始终是固体力学和断裂力学等学科研究的经典问题[1-2], 尤其是分布范围较为广泛的岩石材料, 更是各个领域的关注重点. 对于岩石材料而言, 裂缝萌生是其初始破坏的表现形式[3]. 随着计算机软/硬件技术的发展, 数值模拟方法的仿真能力也取得了极大进展, 已经成为研究岩石材料裂缝萌生-扩展-分岔-交汇过程的有效手段[4].
对岩石材料裂缝发育过程的研究是一个复杂且多方面的过程, 包括了从连续介质到非连续介质、从弹塑性力学到损伤断裂力学和从宏观理论模型到微观理论模型等[5]. 目前, 依据是否将岩石视为连续介质, 数值模拟方法可以被大致分为3类: 连续方法、非连续方法和连续-非连续方法[6-7]. 然而, 这些方法受限于理论基础(例如局部模型和偏微分方程等)与拟模拟物理现象(例如裂缝尖端扩展和波场的传播等)之间的不匹配性, 从而在计算精度、效率和正确性等方面表现出不适应性[8]. 此时, 一种基于非局部理论构建数值框架和基于空间积分方程描述物质力学行为的新兴方法——近场动力学迅速发展[9-11].
近场动力学不仅可以自然地模拟连续介质向非连续介质的转化过程, 还可以模拟非连续介质的演化过程. 近些年来, 该理论愈发受到关注. 最初的近场动力学被称为键基近场动力学, 它重新描述了物质点内力与变形的关系, 但是不涉及应力和应变等概念, 这使其与连续介质力学的联系较为薄弱[12]. 然而, Silling等[13]众多学者认为键基近场动力学模拟脆性材料断裂过程的潜力尚未得到充分开发.
面对连续介质的转化过程, 近场动力学作为一种无网格方法, 在推导之初就引入标量形式的损伤统计函数, 不需要任何外部辅助性的裂缝扩展方向判别方法[14], 这使得近场动力学具有较低的网格依赖性, 即可以避免裂缝扩展方向的近似判别(例如, XFEM等)和较高的计算成本(例如PF等). 键基近场动力学大多采用临界键伸长率或改进的相关断裂准则. 然而, 宏观尺度的断裂类型往往不局限于拉伸破坏, 剪切或混合破坏往往更具有普遍性和突发性, 这意味着, 应该避免单一断裂准则在宏/微观尺度上引起的断裂类型判别矛盾的争议[15-16]. 不仅如此, 鉴于岩石材料复杂受力条件和峰后行为, 需要考虑近场动力学本构模型的适配性问题, 即多种断裂准则和相应弹塑性本构模型的综合性研究亟待开展.
面对非连续介质的演化过程, 近场动力学接触-摩擦模型的研究也在取得进展, 但比较缓慢[17]. 最初, Silling[18]建立了刚性碰撞接触模型和可变形碰撞接触模型, 并加以实现. 随后, 美国Sandia国家实验室开发了Peridigm软件, 同时定义了物质点的接触力计算表达式[19]. 然而, 这些开源的近场动力学分析软件, 只关注如何避免介质间的相互渗透. 此后, Silling等[20]提出了一种改进键基近场动力学接触-摩擦模型; Kamensky等[21]提出了原生接触-摩擦模型和非常规态基近场动力学接触-摩擦模型, 等等. 然而, 已有近场动力学接触-摩擦模型都是通过圆形区域表征物质点接触范围, 通过距离表征物质点接触程度. 通过圆域表征物质点的接触范围会导致接触情形“漏判”和接触程度表征失真, 通过距离表征多维问题的接触程度时具有局限性, 这些接触行为处理方式缺乏坚实的物理意义[18]. 最终, 这两种处理方式加剧了接触行为和摩擦行为的模拟失真, 极大地削弱了近场动力学处理非连续介质演化过程的能力.
虽然近场动力学是一种新兴方法, 但是随着相关研究的不断推进, 它正展现出成为一种评价复杂岩体力学响应的研究手段的巨大潜力. 在键基近场动力学的理论框架下, 本文提出了一种系统化的准态基近场动力学方法, 弥补了该理论处理岩石介质状态转化全过程的现有不足, 以期完善和丰富理论框架, 扩展其在岩石工程领域的应用范围.
1. 基本理论
1.1 键基近场动力学的基本理论简述
如图1所示, 在键基近场动力学的理论框架中[9], 假设某材料构型占据一空间域$ \varPhi $, 在某一时刻t, 该空间域内任一物质点x与其周围一定范围Hx(近场范围)内的其他物质点存在相互作用, 近场范围内其他物质点$ {\boldsymbol{x}}' $为该物质点x (简称中心物质点)的族内物质点, 并且满足以下关系
![]()
图 1 键基近场动力学的材料构型和力密度作用方式Figure 1. The material configuration and force density in the bond-based peridynamics$$ {{\boldsymbol{x}}}'\in {H}_{x}:{d}_{{\boldsymbol{x}}{{\boldsymbol{x}}}'}\leqslant\delta $$ (1) 式中, $ d_{{\boldsymbol{x x}}'} $为欧几里得距离; $ \delta $为近场范围半径.
族内物质点与中心物质点之间的作用力被称为力密度, 它的大小主要取决于物质点的初始相对位置矢量、相对位移矢量和材料参数. 通过力密度函数f综合反映物质点相互作用, 该函数描述内力与变形的关系. 基于此, 构建满足牛顿第二定律形式的运动方程
$$ \rho ({\boldsymbol{x}}) \ddot {\boldsymbol{u}} ({\boldsymbol{x}},t) = \int_{{H_x}} {\boldsymbol{f}} \left( {{\boldsymbol{u}}({\boldsymbol{x}},t),{{\boldsymbol{u}}' }\left( {{{\boldsymbol{x}}' },t} \right),{\boldsymbol{x}},{{\boldsymbol{x}}' },t} \right){\text{d}}{V_{{{\boldsymbol{x}}' }}} + {\boldsymbol{b}}({\boldsymbol{x}},t) $$ (2) 式中, $ \rho $为材料密度; $ \boldsymbol{u} $为物质点x的位移, 加速度$ \ddot {\boldsymbol{u}} = {\partial ^2}{\boldsymbol{u}}({\boldsymbol{x}},t)/{\partial ^2}t $; $ V_{{\boldsymbol{x}}'} $为中心物质点x的族内物质点$ {\boldsymbol{x}}' $的体积; b为体力密度.
近场动力学通过“键”的概念描述物质点之间的相互作用关系. 物质点的相对位置矢量$ {\boldsymbol{\xi}} = {\boldsymbol{x}}' - {\boldsymbol{x}} $和相对位移矢量$ {\boldsymbol{\eta}} = {\boldsymbol{u}}'({\boldsymbol{x}}',t) - {\boldsymbol{u}}({\boldsymbol{x}},t) $的定义可以被用以表征键的长度变化程度, 使力密度函数得以简化
$$ {\boldsymbol{f}}\left( {{\boldsymbol{u}}({\boldsymbol{x}},t),{{\boldsymbol{u}}' }\left( {{{\boldsymbol{x}}' },t} \right),{\boldsymbol{x}},{{\boldsymbol{x}}' },t} \right) = {\boldsymbol{f}}({\boldsymbol{\eta}} ,{\boldsymbol{\xi}} )\frac{{{\boldsymbol{\eta}} + {\boldsymbol{\xi}} }}{{||{\boldsymbol{\eta}} +{\boldsymbol{ \xi}} ||}} $$ (3) 1.2 准态基近场动力学的基本理论
1.2.1 材料变形分析方法
(1)基于键型区分的键变形计算方式
如图2所示, 以二维情形为例, 初始时刻, 在任一近场范围内, 键被划分成平行键(P-bond)和倾斜键(I-bond). 其中, 平行于任一坐标轴的键被称为平行键(平行于横轴的为水平键, 即H-bond; 平行于纵轴的为垂直键, 即V-bond); 不平行于任何坐标轴的键被称为倾斜键.
在变形时刻, 假设物质点在两个方向上产生的位移分量分别为u和v, 通过建立近场动力学微分算子与上述键型之间的联系, 可以获取键的应变分量[22-24]. 如图2所示, 此时, 键的变形计算不仅会考虑端部物质点, 还会考虑周围物质点(贡献物质点), 因此本文所提方法是一种准态基近场动力学方法, 具体如下.
(i) 平行键和倾斜键的应变分量计算
水平键的应变分量
$$ \left. \begin{split} & {\varepsilon _{\text{H}}^x = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}X}} \approx \frac{{{u_i} - {u_j}}}{{{X_i} - {X_j}}}} \\ & {\varepsilon _{\text{H}}^y = \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}Y}} \approx \frac{1}{2}\left( {\frac{{{v_k} - {v_l}}}{{{Y_k} - {Y_l}}} + \frac{{{v_m} - {v_n}}}{{{Y_m} - {Y_n}}}} \right)} \\ & \varepsilon _{\text{H}}^{xy} = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}Y}} + \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}X}} \approx \frac{1}{2}\left[ \left( {\frac{{{u_k} - {u_l}}}{{{Y_k} - {Y_l}}} + \frac{{{v_i} - {v_j}}}{{{X_i} - {X_j}}}} \right) +\right.\\ &\qquad \left.\left( {\frac{{{u_m} - {u_n}}}{{{Y_m} - {Y_n}}} + \frac{{{v_i} - {v_j}}}{{{X_i} - {X_j}}}} \right) \right] \end{split} \right\} $$ (4) 垂直键的应变分量
$$ \left. \begin{split} & {\varepsilon _{\text{V}}^x = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}X}} \approx \frac{1}{2}\left( {\frac{{{u_k} - {u_l}}}{{{X_k} - {X_l}}} + \frac{{{u_m} - {u_n}}}{{{X_m} - {X_n}}}} \right)} \\ & {\varepsilon _{\text{V}}^y = \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}Y}} \approx \frac{{{v_i} - {v_j}}}{{{Y_i} - {Y_j}}}} \\ & \varepsilon _{\text{V}}^{xy} = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}Y}} + \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}X}} \approx \frac{1}{2}\left[ \left( \frac{{{u_i} - {u_j}}}{{{Y_i} - {Y_j}}} + \frac{{{v_k} - {v_l}}}{{{X_k} - {X_l}}} \right) +\right.\\ &\qquad \left.\left( {\frac{{{u_i} - {u_j}}}{{{Y_i} - {Y_j}}} + \frac{{{v_m} - {v_n}}}{{{X_m} - {X_n}}}} \right) \right] \end{split} \right\} $$ (5) 倾斜键的应变分量
$$ \left. \begin{split} & {\varepsilon _{\text{I}}^x = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}X}} \approx \frac{1}{2}\left( {\frac{{{u_i} - {u_k}}}{{{X_i} - {X_k}}} + \frac{{{u_j} - {u_l}}}{{{X_j} - {X_l}}}} \right)} \\ & {\varepsilon _{\text{I}}^y = \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}Y}} \approx \frac{1}{2}\left( {\frac{{{v_i} - {v_l}}}{{{Y_i} - {Y_l}}} + \frac{{{v_j} - {v_k}}}{{{Y_j} - {Y_k}}}} \right)} \\ & \varepsilon _{\text{I}}^{xy} = \frac{{{\mathrm{d}}u}}{{{\mathrm{d}}Y}} + \frac{{{\mathrm{d}}v}}{{{\mathrm{d}}X}} \approx \frac{1}{2}\left[ \left( {\frac{{{u_i} - {u_l}}}{{{Y_i} - {Y_l}}} + \frac{{{v_i} - {v_k}}}{{{X_i} - {X_k}}}} \right) +\right.\\ &\qquad \left.\left( {\frac{{{u_j} - {u_k}}}{{{Y_j} - {Y_k}}} + \frac{{{v_j} - {v_l}}}{{{X_j} - {X_l}}}} \right) \right] \end{split} \right\} $$ (6) 式中, ε为键的应变分量, 其上下角标分别为应变分量标记和键型标记; i, j, k, l, m和n均为物质点序号; X和Y为物质点的初始坐标.
(ii) 平行键和倾斜键的力密度计算
水平键、垂直键和倾斜键的力密度分量
$$\left.\begin{split} & \begin{aligned} & {f_{\text{H}}^x = {{\boldsymbol{T}}^{{\text{PD}}}}\left[ {\begin{array}{*{20}{l}} {\varepsilon _{\text{H}}^x} \\ {\varepsilon _{\text{H}}^y} \end{array}} \right]} \\ & {f_{\text{H}}^y = 0} \end{aligned} \\ &\begin{aligned} & {f_{\text{V}}^x = 0} \\ & {f_{\text{V}}^y = {{\boldsymbol{T}}^{{\text{PD}}}}\left[ {\begin{array}{*{20}{l}} {\varepsilon _{\text{V}}^x} \\ {\varepsilon _{\text{V}}^y} \end{array}} \right]} \end{aligned} \\ & \begin{aligned} & {f_{\text{I}}^x = \frac{{\left| {{X_j} - {X_i}} \right|}}{{\Delta x}}{{\boldsymbol{T}}^{{\text{PD}}}}\varepsilon _{\text{I}}^{xy}} \\ & {f_{\text{I}}^y = \frac{{\left| {{Y_j} - {Y_i}} \right|}}{{\Delta y}}{{\boldsymbol{T}}^{{\text{PD}}}}\varepsilon _{\text{I}}^{xy}} \end{aligned} \end{split} \right\}$$ (7) 式中, $ f $为力密度分量, 其上下角标分别为力密度分量标记和键型标记; $ \Delta x $和$ \Delta y $为离散间隔; TPD为PD刚度矩阵, 其与弹性系数矩阵D存在转化关系
$$ {{\boldsymbol{T}}^{{\text{PD}}}} = {\boldsymbol{D}}{{\boldsymbol{\lambda}} _{{\text{matr }}}}{\delta _{ij}} $$ (8) 具体展开为
$${{\boldsymbol{T}}^{{\text{PD}}}} = \left( \left[ {\begin{array}{*{20}{c}} {\dfrac{E}{{1 - v_{\text{p}}^2}}}&{\dfrac{{Ev}}{{1 - v_{\text{p}}^2}}}&0 \\ {\dfrac{{Ev}}{{1 - v_{\text{p}}^2}}}&{\dfrac{E}{{1 - v_{\text{p}}^2}}}&0 \\ 0&0&{\dfrac{E}{{2\left( {1 + {v_{\text{p}}}} \right)}}} \end{array}} \right]{{\mathrm{or}}}\right.$$ $$\begin{split} &\left.\left[ {\begin{array}{*{20}{c}} {\dfrac{{E\left( {1 - {v_{\text{p}}}} \right)}}{{\left( {1 + {v_{\text{p}}}} \right)\left( {1 - 2{v_{\text{p}}}} \right)}}}&{\dfrac{{E{v_{\text{p}}}}}{{\left( {1 + {v_{\text{p}}}} \right)\left( {1 - 2{v_{\text{p}}}} \right)}}}&0 \\ {\dfrac{{E{v_{\text{p}}}}}{{\left( {1 + {v_{\text{p}}}} \right)\left( {1 - 2{v_{\text{p}}}} \right)}}}&{\dfrac{{E\left( {1 - {v_{\text{p}}}} \right)}}{{\left( {1 + {v_{\text{p}}}} \right)\left( {1 - 2{v_{\text{p}}}} \right)}}}&0 \\ 0&0&{\dfrac{E}{{2\left( {1 + {v_{\text{p}}}} \right)}}} \end{array}} \right] \right)\cdot\\ &\qquad \left[ {\begin{array}{*{20}{c}} {\dfrac{1}{{[m(m + 1)/2]{V_x}\Delta x}}} \\ {\dfrac{1}{{[m(m + 1)/2]{V_x}\Delta x}}} \\ {\displaystyle\sum {\sum {\left| {{Y' } - Y} \right|} } {V_x}} \end{array}} \right]{\delta _{ij}}\\[-1pt] \end{split}$$ (9) 式中, m为$ \delta $与$ \Delta x $的比值; E为弹性模量; vp为泊松比; $ {\boldsymbol{\lambda}}_{\text {matr }} $为转换系数矩阵[25]; $ \delta_{ij} $为克罗内克积.
(2)自适应力密度分量方向分配算法
对于当前方法而言, 力密度分量方向需要重新分配. 以二维情形为例, 首先, 如图3(a)所示, 将物质点i的近场范围划分为4个相同的扇形区域I, II, III和IV; 然后, 如图3(b)所示, 以区域Ⅰ为例, 物质点i与物质点j1, j2和j3分别为垂直键、倾斜键和水平键的端部物质点. 那么, 各类键型的力密度分量分别为
![]()
图 3 自适应力密度分量方向分配算法Figure 3. An adaptive algorithm for force density component direction$$\left. \begin{split} &\left[ {\begin{array}{*{20}{c}} {{{f}}_{\text{H}}^{x(i)}} \\ {{{f}}_{\text{H}}^{y(i)}} \\ {{{f}}_{\text{H}}^{x\left( {{j_{_\ell }}} \right)}} \\ {{{f}}_{\text{H}}^{y\left( {{j_\ell }} \right)}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathrm{sign}}} \left( {{i^x}} \right)} \\ {{{\mathrm{sign}}} \left( {{i^y}} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^x} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^y} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {f_{\text{H}}^{x(i)\left( {{j_\ell }} \right)}} \\ {f_{\text{H}}^{y(i)\left( {{j_\ell }} \right)}} \\ {f_{\text{H}}^{x(i)\left( {{j_\ell }} \right)}} \\ {f_{\text{H}}^{y(i)\left( {{j_{_\ell }}} \right)}} \end{array}} \right]\\ & \left[ {\begin{array}{*{20}{c}} {{{f}}_{\text{V}}^{x(i)}} \\ {{{f}}_{\text{V}}^{y(i)}} \\ {{{f}}_{\text{V}}^{x\left( {{j_{_\ell }}} \right)}} \\ {{{f}}_{\text{V}}^{y\left( {{j_{_\ell }}} \right)}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathrm{sign}}} \left( {{i^x}} \right)} \\ {{{\mathrm{sign}}} \left( {{i^y}} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^x} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^y} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f_{\text{V}}^{x(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{V}}^{y(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{V}}^{x(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{V}}^{y(i)\left( {{j_{_\ell }}} \right)}} \end{array}} \right]\\ &\left[ {\begin{array}{*{20}{c}} {{{f}}_{\text{I}}^{x(i)}} \\ {{{f}}_{\text{I}}^{y(i)}} \\ {{{f}}_{\text{I}}^{x\left( {{j_{_\ell }}} \right)}} \\ {{{f}}_{\text{I}}^{y\left( {{j_{_\ell }}} \right)}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathrm{sign}}} \left( {{i^x}} \right)} \\ {{{\mathrm{sign}}} \left( {{i^y}} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^x} \right)} \\ {{{\mathrm{sign}}} \left( {j_{_\ell }^y} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {f_{\text{I}}^{x(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{I}}^{y(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{I}}^{x(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{\text{I}}^{y(i)\left( {{j_{_\ell }}} \right)}} \end{array}} \right]\end{split} \right\}$$ (10) 式中, f为区域Ⅰ完成方向分配的力密度分量, 其上下角标分别为力密度分量标记和键型标记; $\ell $为键的端部物质点的序号, $ \ell $ = 1, 2, 3; $ f_{\mathrm{H}}^{x(i)\left(j_{\ell}\right)} = f_{\mathrm{H}}^{x} $, 其余与此类似, 不再赘述; i和$j_\ell $为键的端部物质点; sign(·)矩阵为力密度分量方向分配矩阵
$$ \left[ {\begin{array}{*{20}{c}} {{{\mathrm{sign}}} \left( {{i^x}} \right)} \\ {{{\mathrm{sign}}} \left( {{i^y}} \right)} \\ {{{\mathrm{sign}}} \left( {j_\ell ^x} \right)} \\ {{{\mathrm{sign}}} \left( {j_\ell ^y} \right)} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} \alpha \\ \beta \\ { - \alpha } \\ { - \beta } \end{array}} \right]^{\mathrm{T}}} $$ (11) 式中, α 和 β分别为X方向和Y方向的具体力密度分量方向, 初始值可以指定为 +1或−1.
最后, 如图3(c)和图3(d)所示, 以区域Ⅰ为参考, 通过轴对称转换原则和中心对称转换原则, 分别完成其余3个区域的力密度分量方向的分配, 具体为
$$ \left[ {\begin{array}{*{20}{c}} {{{f}}_{{\text{Areal }}}^{x(i)}} \\ {{{f}}_{{\text{Areal }}}^{y(i)}} \\ {{{f}}_{{\text{Areal }}}^{x\left( {{j_{_\ell }}} \right)}} \\ {{{f}}_{{\text{Areal }}}^{y\left( {{j_\ell }} \right)}} \end{array}} \right] = {{{{\boldsymbol{Dire}}}}} \left[ {\begin{array}{*{20}{c}} {f_{{\text{Type}}}^{x(i)\left( {{j_\ell }} \right)}} \\ {f_{{\text{Type }}}^{y(i)\left( {{j_{_\ell }}} \right)}} \\ {f_{{\text{Type}}}^{x(i)\left( {{j_\ell }} \right)}} \\ {f_{{\text{Type }}}^{y(i)\left( {{j_{_\ell }}} \right)}} \end{array}} \right] $$ (12) 式中, f为完成方向分配的力密度分量, 其下角标为区域标记, Areal 为 Area I, Area II, Area III 或Area IV; Type为键型标记, 可为H, V和I; Dire矩阵为4个区域的力密度分量方向转化矩阵, 由α 和 β组成.
(3)键的应力张量
水平键、垂直键和倾斜键的应力分量分别为
$$\left.\begin{split} &\left[ {\begin{array}{*{20}{c}} {\sigma _{\text{H}}^x} \\ {\sigma _{\text{H}}^y} \\ {\tau _{\text{H}}^{xy}} \end{array}} \right] = {\boldsymbol{D}}\left[ {\begin{array}{*{20}{l}} {\varepsilon _{\text{H}}^x} \\ {\varepsilon _{\text{H}}^y} \\ {\varepsilon _{\text{H}}^{xy}} \end{array}} \right]\\ &\left[ {\begin{array}{*{20}{l}} {\sigma _{\text{V}}^x} \\ {\sigma _{\text{V}}^y} \\ {\tau _{\text{V}}^{xy}} \end{array}} \right] = {\boldsymbol{D}}\left[ {\begin{array}{*{20}{c}} {\varepsilon _{\text{V}}^x} \\ {\varepsilon _{\text{V}}^y} \\ {\varepsilon _{\text{V}}^{xy}} \end{array}} \right]\\ &\left[ {\begin{array}{*{20}{l}} {\sigma _{\text{I}}^x} \\ {\sigma _{\text{I}}^y} \\ {\tau _{\text{I}}^{xy}} \end{array}} \right] = {\boldsymbol{D}}\left[ {\begin{array}{*{20}{l}} {\varepsilon _{\text{I}}^x} \\ {\varepsilon _{\text{I}}^y} \\ {\varepsilon _{\text{I}}^{xy}} \end{array}} \right] \end{split}\right\}$$ (13) 式中, $\sigma $和$\tau $均为键的应力分量.
1.2.2 材料状态转化方法
(1)键的断裂准则
在第1.2.1节中, 键的应变分量计算方式依据于几何方程. 键的应变张量不仅可以配合TPD计算键的力密度分量, 而且可以配合D计算键的应力张量. 如此一来, 可以为经典应力断裂准则的引入提供理论基础. 以二维情形为例, 本文详细介绍拉/剪两类准则的引入思路, 具体如下.
(i) 键的主应力计算
$$ {\begin{array}{*{20}{l}} {\sigma _{{\text{Type }}}^1} \\ {\sigma _{{\text{Type }}}^3} \end{array} = \frac{{\sigma _{{\text{Type }}}^x + \sigma _{{\text{Type }}}^y}}{2} \mp \sqrt {{{\left( {\frac{{\sigma _{{\text{Type }}}^x - \sigma _{{\text{Type }}}^y}}{2}} \right)}^2} + {{\left( {\tau _{{\text{Type }}}^{xy}} \right)}^2}} } $$ (14) 式中, $ \sigma_{\text {Type }}^{3} $为最大主应力; $ \sigma_{\text {Type }}^{1} $为最小主应力.
(ii) 最大拉应力准则和莫尔-库仑准则
当键可能发生拉伸断裂时, 利用最大拉应力准则进行判断
$$ f_{{\text{Type }}}^{{\text{Tensile }}} = {\sigma _{\text{t}}} - \sigma _{{\text{Type }}}^3 < 0 $$ (15) 式中, $ f_{\text {Type }}^{\text {Tensile }} $为拉伸屈服函数; $ \sigma_{\mathrm{t}} $为抗拉强度.
当键可能发生剪切断裂时, 利用莫尔-库仑准则进行判断
$$ f_{{\text{Type}}}^{{\text{Shear }}} = \sigma _{{\text{type}}}^1 - \sigma _{{\text{type}}}^3{N_\varphi } + 2c\sqrt {{N_\varphi }} < 0 $$ (16) $$ {N_\varphi } = \frac{{1 + \sin \varphi }}{{1 - \sin \varphi }} $$ (17) 式中, $ f_{\text {Type }}^{\text {Shear }} $为剪切屈服函数; $ c $为黏聚力; φ为内摩擦角.
这两种准则被分别用以判断微观尺度上键的拉伸和剪切断裂, 从而避免近场动力学在宏/微观尺度上具有断裂类型判别矛盾的争议. 应当指出, 因为岩石材料具有抗压而不抗拉的特点, 所以, 先判断键是否发生拉伸断裂, 再判断键是否发生剪切断裂.
(2)虚拟裂缝模型
在断裂力学中, 虚拟裂缝模型的本质是通过模拟黏聚力衰减过程, 复现介质状态和承载能力逐渐变化的过程. 目前, 因为近场动力学的现有弹塑性本构模型不再适用于当前理论框架, 本文提出一种适于近场动力学的单直线型虚拟裂缝模型, 下面介绍具体内容.
(i) 初始阶段
这个阶段是虚拟裂缝的萌生阶段. 当键发生断裂时, 虚拟裂缝随之产生, 这意味着虚拟裂缝的萌生过程是不可逆的. 根据键的断裂类型, 虚拟裂缝可被分为Ⅰ型(拉伸虚拟裂缝)和Ⅱ型虚拟裂缝(剪切虚拟裂缝). 应当指出, 虚拟裂缝具有传递应力的能力, 它与真实裂缝都是真实存在的.
(ii) 过渡阶段
这个阶段是虚拟裂缝的发育阶段. 首先, 如图4所示, 利用法向张开度w和法向黏聚力密度$ f_{{\mathrm{n}}}^{\text {coh }} $(具有黏聚力效应, 量纲为FL−6, 切向黏聚力密度$ f_{{\mathrm{s}}}^{\text {coh }} $同样如此)表征拉伸虚拟裂缝的发育程度, 利用切向滑移量s和$ f_{{\mathrm{s}}}^{\text {coh }} $表征剪切虚拟裂缝的发育程度
![]()
图 4 适于准态基近场动力学的单直线虚拟裂缝模型和弹塑性本构模型Figure 4. Single-linear fictitious crack models and elastic-plastic constitutive models suitable for quasi-state-based peridynamics$$ \left.\begin{split} &{{f}}_{{{\mathrm{n}}}}^{\text{coh}} = \left\{\begin{aligned} &{{f}}_{{{\mathrm{n}}}}^{\mathrm{lim}}\left(1-\frac{w}{{w}_{\text{n}}}\right),\quad 0\leqslant w \leqslant {w}_{\text{n}}\\ &0,\quad w > {w}_{\text{n}}\end{aligned}\right. \\ &{{f}}_{\text{s}}^{\text{coh}} = \left\{\begin{aligned} &{{f}}_{\text{s}}^{\mathrm{lim}}\left(1-\frac{s}{{s}_{\text{p}}}\right),\quad 0\leqslant s \leqslant {s}_{\text{p}}\\ &0,\quad s > {s}_{\text{p}}\end{aligned}\right.\end{split}\right\} $$ (18) 式中, $ {f}_{\rm{n}}^{\lim } $和$ {f}_{\rm{s}}^{\lim } $分别为键发生拉伸/剪切断裂时的临界力密度, 这两种临界力密度的计算方式与第1.2.1节的内容相同; $ w_{{\mathrm{n}}} $和$ s_{\mathrm{p}} $为临界的法向张开度和切向滑移量.
然后, 因为虚拟裂缝的发育阶段必然消耗能量, 在断裂能未被完全消耗之前, 虚拟裂缝具有传递应力的能力. 此时, 利用Ⅰ型断裂能$ G_{{\mathrm{f}}}^{{\mathrm{I}}} $、Ⅱ型断裂能$ G_{\mathrm{f}}^{\mathrm{II}} $和材料强度参数可获取拉伸和剪切虚拟裂缝的$ w_{\mathrm{n}} $和$ s_{{\mathrm{p}}} $.
最后, 考虑到虚拟裂缝的发育阶段往往不是简单的单调扩展过程, 通过引入卸载/重新加载阶段模拟虚拟裂缝的间歇性闭合现象.
(iii) 结束阶段
这个阶段是虚拟裂缝的消失阶段. 当拉伸和剪切虚拟裂缝的w和s达到$ w_{\mathrm{n}} $和$ s_{{\mathrm{p}}} $时, 虚拟裂缝将转化为真实裂缝, 即不再具有传递应力的能力.
1.2.3 材料接触行为分析方法
(1)空间网格的反衍和离散过程
近场动力学作为一种无网格法, 求解体系不会考虑空间网格的影响, 这也是现有接触-摩擦模型忽视空间网格的根本原因. 因此, 本文提出粒子网格孪生算法反衍物质点的空间网格, 实现两者的同步计算. 然后, 通过绑定断裂准则和空间网格的离散过程, 实现空间网格的状态转化过程, 从而为接触情形的准确判断和接触程度的准确表征提供坚实的理论基础.
粒子网格孪生算法是通过等效网格节点, 实现物质点的空间网格构建. 在这里, 网格节点状态完全依赖于物质点状态, 即孪生网格. 孪生网格可以等效地反映材料构型的变形行为. 如图5(a)所示, 物质点位于孪生网格的中心, 孪生网格的形状并不固定, 可以是任意多边形, 只需要满足多数量的连续性即可. 如图5(b)和图5(c)所示, 孪生网格可以被分为节点共享型和节点独立型. 节点共享型孪生网格指的是物质点的孪生网格具有相同的节点. 例如, 图5(b)中的物质点1和2的孪生网格节点都包括节点5和8. 节点独立型孪生网格指的是物质点的孪生网格具有不相同的节点. 这些网格节点会存在空间重合的情形. 例如, 图5(c)中的物质点1的孪生网格节点包括节点9 ~ 12, 物质点2的孪生网格节点包括节点13 ~ 16. 其中, 节点12和15存在空间重合, 它们会同时考虑物质点1和2的影响. 此时, 物质点1被称为节点15的相关物质点, 物质点2被称为节点15的原生物质点. 应当指出, 因为节点独立型孪生网格在离散化过程中无需新增节点, 易于数值实现, 所以本文将其作为反衍目标.
孪生网格的节点构造过程分为初始阶段和实时阶段. 在初始阶段, 它的坐标信息取决于原生物质点
$$ \left( {X_{{\text{node }}}^i,Y_{{\text{node }}}^i} \right) = f\left( {X_{{\text{N}} - {\text{point }}}^{[{\text{m}}]} \pm \lambda \Delta X,Y_{{\text{N}} - {\text{point }}}^{[{\text{m}}]} \pm \lambda \Delta Y} \right) $$ (19) 式中, $ X^{i}_{{\mathrm{node}}} $和$ Y^{i}_{{\mathrm{node}}} $分别为网格节点在X和Y方向的初始坐标, i为网格节点序号; $ X_{\mathrm{N}-\mathrm{point}}^{[m]} $和$ Y_{\mathrm{N}-\mathrm{point} }^{[\mathrm{m}]} $分别为原生物质点在X和Y方向的初始坐标, [m]为原生物质点序号; $ \lambda $为坐标折减系数.
在实时阶段, 它的坐标信息会考虑两类物质点
$$ \left( {\dot X_{{\text{node }}}^i,\dot Y_{{\text{node }}}^i} \right) = \frac{1}{{[n]}}\sum\limits_{{\text{[m]}} = 1}^{[n]} {\left( {\dot X_{{\text{M}} - {\text{point}}}^{[{\text{m}}]},\dot Y_{{\text{M}} - {\text{point}}}^{[{\text{m}}]}} \right)} $$ (20) 式中, [n]为两类物质点的数量之和; $ \dot{X}_{{\mathrm{M-point}}}^{[{\mathrm{m}}]} $和$ \dot{Y}_{\mathrm{M-point}}^{[\mathrm{m}]} $分别分别为上述两类物质点在X和Y方向的实时坐标. 针对可能发生的键断裂情况, 将断裂准则与孪生网格离散过程进行同步绑定. 孪生网格离散过程也是网格节点的相关物质点数量减少的过程. 随着相关物质点数量的逐渐减少, 原有空间重合程度就会减弱, 最终表现为“连续状态”的孪生网格发生相互脱离的效果, 如图6所示.
![]()
图 6 物质点的孪生网格的状态转化过程Figure 6. Twin mesh state transformation process of material point(2)物质点对的接触检索
一般地, 接触检索是接触-摩擦模型建立的首要任务. 因此, 在上述基础上, 本文通过设计“粗判”和“细判”检索方法确定发生接触的物质点对, 从而触发对接触/摩擦力密度的计算.
如图7(a)所示, “粗判”检索是判断物质点对的直线距离是否低于标准距离
$$ {d}_{x{x}'}\leqslant{d}_{\text{STD}} $$ (21) 式中, $ d_{x x'} $为欧几里得距离; $ d_{\mathrm{STD}} $为预设标准距离.
因为近场动力学基于非局部理论, 当物质点失去相互作用时, 两者的实时位置可能并不相邻. “粗判”检索可以迅速筛选这类物质点对.
针对上述结果, 如图7(b)所示, “细判”检索是判断物质点的某一网格节点是否侵入至对方的孪生网格中, 具体为
$$ \left.\begin{aligned} &\boldsymbol{AB}\cdot {{{\boldsymbol{Q}}}_{1}{\boldsymbol{B}}} \geqslant 0 \\ &\boldsymbol{CB}\cdot {{{\boldsymbol{Q}}}_{1}{\boldsymbol{B}}}\geqslant 0 \\ &{{{\boldsymbol{CQ}}}_{2}}\cdot {{{\boldsymbol{Q}}}_{1}{{\boldsymbol{Q}}}_{2}} \geqslant 0 \\ &{{{\boldsymbol{AQ}}}_{2}}\cdot {{{\boldsymbol{Q}}}_{1}{{\boldsymbol{Q}}}_{2}} \geqslant 0\end{aligned}\right\} $$ (22) (3)基于势函数的物质点接触力密度算法
经典的基于势函数的接触算法往往通过距离相关量表征接触程度[26]. 显然, 距离作为一维物理量, 容易受到限制. 因此, 本文提出了基于面积比例的势函数的接触算法, 具体内容如下.
针对“细判”检索的结果, 将发生接触的物质点分别称为接触物质点和目标物质点, 如图8所示. 允许此时两个物质点的孪生网格存在互嵌区Ω, 该区域内的某一微元对接触物质点的微小接触力密度
$$ {\text{d}}{{\boldsymbol{f}}_{{\text{cf}}}} = {\text{d}}{{\boldsymbol{f}}_{\text{C}}} - {\text{d}}{{\boldsymbol{f}}_{\text{T}}} $$ (23) 式中, dfC为接触物质点的微元嵌入目标物质点受到的微小接触力密度, dfT同理.
假设互嵌区内的微元利用$ \omega (x,y) $进行表示, 通过势函数$ \varTheta (\cdot) $表征物质点对侵入物质点的排斥作用. 那么, 上述微小接触力密度都可以被表征为
$$ {\text{d}}{{\boldsymbol{f}}_{{\text{C/T}}}} = - \nabla \varTheta (\omega ){\text{d}}A $$ (24) 式中, $ \nabla \varTheta (\omega ) $为势函数梯度, dA为微元面积.
如此一来, 微小接触力密度可以进一步表征为
$$ \left.\begin{aligned} &\text{d}{{\boldsymbol{f}}}_{\text{C}} = -\nabla {\varTheta }_{\text{T}}\left({\omega }_{1}\right)\text{d}A \\ &\text{d}{{\boldsymbol{f}}}_{\text{T}} = -\nabla {\varTheta }_{\text{C}}\left({\omega }_{2}\right)\text{d}A\end{aligned}\right\} $$ (25) 由于空间位置重叠, 对互嵌区进行积分, 即将式(25)代入到式(23)中, 整合后得到总接触力密度
$$ {{\boldsymbol{f}}}_{\text{cf}} = {\displaystyle {\int }_{S = {S}_{\text{C}}\cap {S}_{\text{T}}} \text{d}}{{\boldsymbol{f}}}_{\text{cf}}\text{d}A = {\displaystyle {\int }_{S = {S}_{\text{C}}\cap {S}_{\text{T}}}\left[ \nabla {\varTheta }_{\text{C}}\left({\omega }_{2}\right)-\nabla {\varTheta }_{\text{T}}\left({\omega }_{1}\right) \right]}\text{d}A $$ (26) 利用格林公式进行积分转换
$$ {{\boldsymbol{f}}_{{\text{cf}}}} = \int_{\varGamma = {\varGamma _{\text{C}}} \cap {\varGamma _{\text{T}}}} {{{\boldsymbol{n}}_\varGamma }} \left( {{\varTheta _{\text{C}}} - {\varTheta _{\text{T}}}} \right){\text{d}}\varGamma $$ (27) 式中, $ {\boldsymbol{n}}_{\varGamma} $是孪生网格边界$ \varGamma $的单位外法向量.
此时, 互嵌区的接触力密度是分布在互嵌区边界上的力密度. 通过集度与力的关系, 确定集度$ {{\boldsymbol{P}}_\varGamma } $为
$$ {{\boldsymbol{P}}_\varGamma } = {{\boldsymbol{n}}_\varGamma }\left( {{\varTheta _{\text{C}}} - {\varTheta _{\text{T}}}} \right) $$ (28) 然后, 通过建立适宜的势函数, 使其和接触程度建立联系, 本文提出了一种基于面积比例的势函数
$$ \varTheta (\omega ) = \left\{ \begin{aligned} & {{K_{\text{n}}}{S_{{\text{r }}}},}\quad {\omega \in \varOmega } \\ & {0,}\quad {\omega \notin \varOmega } \end{aligned} \right. $$ (29) 式中, Kn为接触刚度; Sr为内嵌区面积和孪生网格的面积比例, 无量纲.
然后, 对该势函数进行统一标定. 这是因为, 非均匀离散或网格变形都会引起孪生网格的差异过大, 势函数会受到孪生网格变化的影响, 易造成嵌入程度相同而接触力密度不同的情形, 如图9所示.
为了使势函数基于统一的标准对接触程度进行表征, 同时保证势函数具有坚实的物理意义. 本文引入一个标尺面积Sstd, 该标尺面积可视具体情况而定. 例如, 最小的网格尺寸. 此时的势函数表达式可进一步写为
$$ \varTheta (\omega ) = \left\{\begin{aligned} &{K}_{\text{n}}\frac{{S}_{\text{ini}}}{{S}_{}},\quad \text{small-strain}\text{, }\omega \in \varOmega \\ &{K}_{\text{n}}\frac{{S}_{\text{ini}}}{{S}_{\text{std}}},\quad \text{large-strain}\text{, }\omega \in \varOmega \\ &0,\quad \omega\notin \varOmega \end{aligned}\right. $$ (30) 式中, Sini为统一标定后的互嵌区面积; S为均匀离散的孪生网格面积.
此时, 通过式(27)得到的接触力密度分布在互嵌区边界上. 由于物质点对的接触行为需要通过集中式力密度进行表征, 因此需要转化. 如图10所示, 接触物质点内嵌区边界上的力密度fC21和fC22, 配合单位向量e12和$ {\boldsymbol{e}}_{11}' $, 可以获得等效集中力密度$ \boldsymbol{f}_{\mathrm{C} 21}' $和$ \boldsymbol{f}_{\mathrm{C} 22}' $. 此时, 目标物质点的集中力密度分量为
$$ \left. \begin{aligned} & {{\boldsymbol{f}}_{\text{T}}' = {{\boldsymbol{f}}_{{\text{C}}21}}{\boldsymbol{e}}_{12}' + {{\boldsymbol{f}}_{{\text{C}}22}}{\boldsymbol{e}}_{12}' } \\ & {{{\boldsymbol{f}}_{\text{T}}} = {{\boldsymbol{f}}_{{\text{C}}21}}{{\boldsymbol{e}}_{12}} + {{\boldsymbol{f}}_{{\text{C}}22}}{{\boldsymbol{e}}_{12}}} \end{aligned} \right\} $$ (31) 1.2.4 材料摩擦行为分析方法
本文在上述接触算法的基础上, 根据库伦摩擦定律, 将摩擦力密度计算过程分为静摩擦阶段和动摩擦阶段
$$ \left.\begin{aligned} &{{\boldsymbol{f}}}_{\text{s}} = {{\boldsymbol{f}}}_{\text{fri\_s}},\quad \text{static friction}\text{ } \\ &{{\boldsymbol{f}}}_{\text{s}} = {{\boldsymbol{f}}}_{\text{fri\_n}},\quad \text{sliding friction}\text{ }\end{aligned} \right\}$$ (32) 式中, fs为物质点的摩擦力密度; ffri_s为静摩擦阶段的摩擦力密度; ffri_n为滑动摩擦阶段的摩擦力密度.
因为摩擦力和法向接触力相关, 如图11所示, 将物质点的接触力密度沿滑动界面的平行方向和垂直方向进行分解, 获得切向和法向接触力密度, $ \boldsymbol{f}_{\text {cfx }}^{1} $和$ \boldsymbol{f}_{\text {cfx }}^{2} $为法向接触力密度, $ \boldsymbol{f}_{\text {cfy }}^{1} $和$\boldsymbol{f}_{\text {cfy }}^{2} $为切向接触力密度.
在静摩擦阶段, 即$ {\boldsymbol{f}}_{s}\leqslant \left|{\boldsymbol{f}}_{{\mathrm{fti}}\_{\mathrm{s}}}^{\max }\right| $(最大静摩擦力密度), 物质点的静摩擦力密度为
$$ {{\boldsymbol{f}}_{{\text{s }}}} = {{\boldsymbol{f}}_{{\text{fri\_s }}}} = {K_{\text{s}}}\frac{{\Delta \eta }}{\delta }{{\mathrm{sign}}} \left( {{v^{[{\text{m}}]}}} \right) $$ (33) 式中, Ks为摩擦刚度; $ \Delta \eta $为相对位移; sign(v[m])为速度相关符号函数.
在滑动摩擦阶段, 即$ {\boldsymbol{f}}_{s}>\left|{\boldsymbol{f}}_{ {{\mathrm{fri}}\_{\mathrm{s}} }}^{\max }\right| $, 物质点的动摩擦力密度为
$$ {{\boldsymbol{f}}_{\text{s}}} = {{\boldsymbol{f}}_{{{{\mathrm{fri}}\_{\mathrm{n}}}}}} = {\mu _{}}\left| {{{\boldsymbol{f}}_{{{{\mathrm{cf}}\_{\mathrm{n}}}}}}} \right|{ {\boldsymbol{S}}_{{v^{[{\text{m}}]}}}} $$ (34) 式中, $ \mu $为摩擦系数; $ {\boldsymbol{f}}_{ {{\mathrm{cf}}\_{\mathrm{n}} }} $为法向接触力密度; $ {\boldsymbol{S}}_{{v^{[{\text{m}}]}}} $为物质点速度相关的单位向量.
此时, 在近场动力学的理论框架下, 将物质点对的接触力密度和摩擦力密度纳入统一表征
$$ \begin{split} & \rho \ddot {\boldsymbol{u}} ({\boldsymbol{x}},t) = \int {\boldsymbol{f}} \left( {{\boldsymbol{u}}({\boldsymbol{x}},t),{{\boldsymbol{u}}' }\left( {{{\boldsymbol{x}}' },t} \right),{\boldsymbol{x}},{{\boldsymbol{x}}' },t} \right){v_{c\left( {{{\boldsymbol{x}}' }} \right)}}{\text{d}}{V_{{{\boldsymbol{x}}' }}} + \\ &\quad \int {{{\boldsymbol{f}}_{{\text{cf}}}}} \left( {{\boldsymbol{u}}({\boldsymbol{x}},t),{{\boldsymbol{u}}''}\left( {{{\boldsymbol{x}}''},t} \right),{\boldsymbol{x}},{{\boldsymbol{x}}''},t} \right){v_{c\left( {{{\boldsymbol{x}}''}} \right)}}{\text{d}}{V_{{{\boldsymbol{x}}''}}} + \\ &\quad \int {{{\boldsymbol{f}}_{\text{s}}}} \left( {{\boldsymbol{u}}({\boldsymbol{x}},t),{{\boldsymbol{u}}''}\left( {{{\boldsymbol{x}}''},t} \right),{\boldsymbol{x}},{{\boldsymbol{x}}''},t} \right){v_{c\left( {{{\boldsymbol{x}}''}} \right)}}{\text{d}}{V_{{{\boldsymbol{x}}''}}} + {\boldsymbol{b}}({\boldsymbol{x}},t) \end{split} $$ (35) 式中, $ v_{{c}\left({\boldsymbol{x}}'\right)} $为族内物质点体积修正因子; $ v_{c\left(\boldsymbol{x}''\right)} $为与中心物质点发生接触/摩擦的物质点体积修正因子, $ V_{{\boldsymbol{x}}''} $为其体积; $ \boldsymbol{f}_{{\mathrm{cf}}}\left({{\boldsymbol{u}}}({{\boldsymbol{x}}}, t), {{\boldsymbol{u}}}''\left(\boldsymbol{x}'', t\right), \boldsymbol{x}_{,}, \boldsymbol{x}'', t\right) $为接触力密度函数; $ \boldsymbol{f}_{{\mathrm{s}}}\left(\boldsymbol{u}(\boldsymbol{x}, t), \boldsymbol{u}^{\prime \prime}\left(\boldsymbol{x}^{\prime \prime}, t\right), \boldsymbol{x}, \boldsymbol{x}^{\prime \prime}, t\right) $为摩擦力密度函数.
2. 数值实现
2.1 体积求和方式
作为一种无网格法, 近场动力学将材料构型视为具有体积信息和质量信息的物质点集, 将近场动力学运动方程转化为适于时间和空间离散化的控制方程后, 可以将每个物质点视为研究对象进行分析, 从而描述材料的力学响应. 一般地, 将计算模型转化为离散化系统包括两个过程: 首先, 将标准有限元网格等效为离散网格; 然后, 将位于网格中心的一阶高斯积分点等效为物质点, 并对其各类力学行为进行统一表征处理(见式(35)). 在面对上述离散化系统时, 本文自主开发的数值计算程序将物质点体积积分求解过程转化为族内物质体积求和求解过程
$$ \begin{split} &\rho \ddot {\boldsymbol{u}}\left({{\boldsymbol{x}}}_{i},t\right) = {\displaystyle \sum _{i = 1}^{{N}_{i}}{\displaystyle \sum _{j = 1}^{{N}_{j}}{\boldsymbol{f}}}}(\cdot){v}_{c\left({{\boldsymbol{x}}}_{j}'\right)}{V}_{\left({{\boldsymbol{x}}}_{j}'\right)} + \\ &\qquad {\displaystyle \sum _{i = 1}^{{N}_{i}}{\displaystyle \sum _{l = 1}^{{N}_{l}}{{\boldsymbol{f}}}_{\text{cf}}}}(\cdot){v}_{c\left({{\boldsymbol{x}}}''_{l}\right)}{V}_{\left({{\boldsymbol{x}}}''_{l}\right)} + \\ &\qquad {\displaystyle \sum _{i = 1}^{{N}_{i}}{\displaystyle \sum _{m = 1}^{{N}_{m}}{{\boldsymbol{f}}}_{\text{s}}}}(\cdot){v}_{c\left({{\boldsymbol{x}}}''_{m}\right)}{V}_{\left({x}''_{m}\right)} + {\boldsymbol{b}}\left({{\boldsymbol{x}}}_{i},t\right)\end{split} $$ (36) 式中, Ni为中心物质点数量; Nj是第i个中心物质点的族内物质点数量; $ V_{\left({\boldsymbol{x}}_{j}'\right)} $为第i个中心物质点的第j个族内物质点的体积, $ {v}_{c\left({{\boldsymbol{x}}}_{j}'\right)} $为其体积修正因子; Nl是与第i个中心物质点发生接触的物质点数量; $ {V}_{\left({{\boldsymbol{x}}}_{j}'\right)} $为第l个发生接触的物质点的体积, $ v_{{{c}}\left({\boldsymbol{x}}_{l}''\right)} $为其体积修正因子; Nm是与第i个中心物质点发生摩擦的物质点数量; $ V_{\left({\boldsymbol{x}}_{m}''\right)} $为第m个发生摩擦的物质点的体积, $ v_{c\left({\boldsymbol{x}}''_{m}\right)} $为其体积修正因子.
2.2 迭代求解方式
本文通过引入基于人工阻尼的动态松弛方法, 使离散化系统的动态物质点可以在人工阻尼的作用下达到平衡状态
$$ \begin{split} & \rho \ddot {\boldsymbol{u}} \left( {{{\boldsymbol{x}}_i},t} \right) + C\mathop {\boldsymbol{u}}\limits^. ({\boldsymbol{x}},t) = \sum\limits_{i = 1}^{{N_i}} {\sum\limits_{j = 1}^{{N_j}} {\boldsymbol{f}} } ( \cdot ){v_{c\left( {{\boldsymbol{x}}_j' } \right)}}{V_{\left( {{\boldsymbol{x}}_j' } \right)}} + \\ &\qquad \sum\limits_{i = 1}^{{N_i}} {\sum\limits_{l = 1}^{{N_l}} {{{\boldsymbol{f}}_{{\text{cf}}}}} } ( \cdot ){v_{c\left( {{\boldsymbol{x}}_l^{\prime \prime }} \right)}}{V_{\left( {{\boldsymbol{x}}_l^{\prime \prime }} \right)}} + \\ &\qquad \sum\limits_{i = 1}^{{N_i}} {\sum\limits_{m = 1}^{{N_m}} {{{\boldsymbol{f}}_{\text{s}}}} } ( \cdot ){v_{c\left( {{\boldsymbol{x}}_m^{\prime \prime }} \right)}}{V_{\left( {{\boldsymbol{x}}_m^{\prime \prime }} \right)}} + {\boldsymbol{b}}\left( {{{\boldsymbol{x}}_i},t} \right) \end{split} $$ (37) 式中, $ C \dot{\boldsymbol{u}}(\boldsymbol{x}, t) $为局部阻尼项; C为局部阻尼系数.
然后, 通过显示中心差分方法实现上述迭代过程.
3. 数值验证
3.1 基准算例
基准算例1: 岩样尺寸为0.1 m × 0.4 m × 0.001 m, 在岩样上端面施加大小为0.01 m/s的垂直向上速度, 在岩样下端面施加法向约束, 并岩样上下对称面设置允许开裂位置. 该试验的各类材料参数取值和加载方式均依据于文献[27]. 该算例的$ G_{\mathrm{f}}^{\mathrm{I}} $的取值为500 N/m.
基准算例2: 岩样尺寸为0.2 m × 0.2 m × 0.001 m, 在岩样上端面施加大小为0.01 m/s的水平向右速度, 在岩样上/下端面施加法向约束, 仍在上下对称面设置允许开裂位置. 该试验的各类材料参数取值和加载方式均依据文献[28]. 该算例共设计了3组试验, 第1组 ~ 第3组试验的区别仅在于$ G_{{\mathrm{f}}}^{{\mathrm{II}}} $的取值, 分别为1000, 1250和1500 N/m.
图12和图13分别为基准算例1和2的载荷-位移曲线. 由图12可以发现, 峰值载荷为195.3 N, 对应的峰值应力为1.954 MPa, 这与岩样的理论抗拉强度(2 MPa)极为接近. 此外, $ G_{\mathrm{f}}^{{\mathrm{I}}} $对峰后载荷刚刚接近于残余阶段(由局部阻尼项引起的阻尼力造成)的临界位移有明显影响, 临界位移为5.127 × 10−4 m. 应当指出, 该临界位移与拉伸虚拟裂缝的理论$ w_{{\mathrm{n}}} $十分接近. 这些内容均可以证明当前本构模型计算结果的正确性.
![]()
图 13 基准算例2的第1组 ~ 第3组试验的载荷-位移曲线Figure 13. Loading-displacement curves of tests 1-3 in fundamental test 2由图13可以发现, 这3组试验的载荷-位移曲线可以被划分为4个阶段, 分别为弹性阶段、硬化阶段、软化阶段和残余阶段. 峰后斜率与$ G_{{{\mathrm{f}}}}^{{\mathrm{II}}} $具有正向相关关系, 峰后斜率随着$ G_{{{\mathrm{f}}}}^{{\mathrm{II}}} $的增加而升高, 曲线均呈现近似线性变化的趋势, 岩样主要表现为塑性, 最终剪应力均维持在0.1743 MPa附近. 从理论上讲, 将3组试验的载荷-位移曲线等效转化为应力-应变曲线后, 应力-应变曲线与坐标轴所围面积即为$ G_{{{\mathrm{f}}}}^{{\mathrm{II}}} $, 各组试验的$ G_{{{\mathrm{f}}}}^{{\mathrm{II}}} $为1076.20, 1341.86和1608.15 N/m, 相对误差分别约为7.62%, 7.35%和7.21%.
3.2 预制单倾斜裂缝砂岩压缩实验
如图14所示, 砂岩试样的尺寸为0.06 m × 0.12 m × 0.0001 m, Kn = 5.777 × 1022 N·m−6, E = 28.12 GPa, vp = 0.2, $ \rho $ = 2620 kg/m3, $ \mu $ = 0.8, $ \sigma_{{\mathrm{t}}} $ = 12 MPa, φ = 35°, c = 81.9 MPa, 物质点总数为22149, 在平面应变条件下进行计算. 在岩样中部预制一条长度为0.01 m并与X轴夹角为45°的倾斜裂缝(方案1). 岩样上/下端面的加载速度v的大小为0.005 m/s. 然后, 仅将裂缝长度分别增加至0.015 m (方案2)和0.02 m (方案3). 在3个方案中均布置了如图14所示的4条测线, 分别为T1, T2, T3和T4. 应当指出, 该试验的计算参数和方案设计均与文献[29]相同.
![]()
图 14 中心预制45°倾斜裂缝的矩形砂岩试样计算示意图Figure 14. Calculation schematic diagram of a rectangular sandstone specimen with centrally prefabricated 45° inclined fracture图15分别给出了室内试验结果[30]和本文的数值模拟结果. 如图15(a)所示, 预制裂缝尖端沿着岩样的横向和竖向进行了延伸, 最终形成了两条翼型裂缝. 在整个裂缝扩展路径上, 分布大量全部键断裂的“橙色”物质点(全破坏物质点), 这表明岩样在破坏过程中产生了细小碎屑. 这是FEM, XFEM, PF等数值方法难以模拟的效果.
如图15(b)和图15(c)所示, 随着预制裂缝长度的增加, 预制裂缝尖端的延伸程度和裂缝类型都发生了变化. 岩样产生了4种类型裂缝, 分别为翼型裂缝、抗张拉裂缝、倾斜次生裂缝和共面次生剪切裂缝, 这些裂缝与预制裂缝之间形成了贯通. 此外, 岩样损伤程度与预制裂缝长度呈正相关关系, 即随着预制裂缝长度的增加, 岩样的损伤程度也增加.
在经典近场动力学的接触-摩擦模型中, 通过圆域式接触范围表征方式和距离式接触程度表征方式, 模拟非连续介质的相互排斥情形. 然而, 这是一种简化式的处理方式, 模拟得到的非连续介质力学行为的误差较为明显.
为了准确分析物质点的接触行为, 本文对监测物质点60, 525, 820的接触情况进行了观察和统计, 筛选预分析的时步范围(物质点对的直线距离不变, 而接触程度发生变化). 最终, 在时步数目$ N_{\text {step }} = 37\;600 \sim 37\;800 $, $ N_{\text {step}} = 18\;200 \sim 18\;400 $, $ N_{\text {step}} = 57\;600 \sim 57\;800 $内, 观察到接触力密度分别呈现“下降”、“下降”和“上升”的趋势, 如图16所示. 然而, 经典模型的接触力密度在该时段内不会产生变化. 显然, 距离式接触程度表征方式的分析和计算存在失真的情况, 这种情形在整个计算过程中会普遍存在, 这是影响上述裂缝扩展过程的重要因素之一.
![]()
图 16 监测物质点的接触力密度的变化曲线Figure 16. The curves of the contact force density of monitoring points与此同时, 本文将物质点的接触形式分成平行式接触(P-contact)和非平行式接触(NP-contact). 平行式接触指的是物质点的内嵌区边界有两条完整且平行的孪生网格边界的接触情形, 其余则为非平行式接触. 其中, 非平行式接触包括一种对角式接触(D-contact), 它指的是内嵌区边长均小于0.1倍的孪生网格边长的情形. 如图17所示, 上述3个监测物质点的接触总次数分别为109561, 221561和22411. 对于3个监测物质点而言, 非平行式接触所占的比例都是最大的, 对角式接触所占比例分别为39.14%, 44.51%和33.56%, 即远超总接触次数的1/3. 然而, 经典模型会“漏判”上述占比很大的对角式接触. 显然, 圆域式接触范围表征方式存在较大的误差, 这同样是影响上述裂缝扩展过程的重要因素.
![]()
图 17 监测物质点的各类接触形式的数量和比例Figure 17. Number and proportion of various contact forms at monitoring points3.3 预制多倾斜裂缝砂岩压缩实验
如图18所示, 砂岩试样的尺寸为0.08 m × 0.16 m × 0.0001 m, 计算参数和加载条件与文献[29]的相同. 在岩样中部预制了3条长度为0.015 m的倾斜裂缝, 分别为C1, C2和C3, 它们与X轴夹角均为45°, C1和C3的内部裂缝尖端夹角$ \theta $为75°(方案1). 然后, 保持各项情况不变, 将夹角$ \theta $依次增加至90°(方案2)和105°(方案3).
![]()
图 18 多倾斜裂缝的矩形砂岩试样计算示意图Figure 18. Calculation schematic diagram of a rectangular sandstone specimen with multi-inclined fractures图19给出了方案1的预制倾斜裂缝扩展过程云图、最大主应力σ3云图、剪应力云图、位移场云图和应力-应变曲线. 首先, 在位移载荷的作用下, 岩样内部会进行应力调整. 然后, 当Nstep = 28541时, 岩样应力达到点I, 岩样的最大主应力和剪应力(绝对值)主要分布在3条预制倾斜裂缝的尖端附近, 并且应力场呈现对称性. 由于裂缝尖端存在应力集中现象, 裂缝尖端将会率先产生损伤, 并逐渐向外扩散. 此外, 虽然倾斜裂缝使得两个方向的位移场的连续性出现不同程度的波动, 但是它们的分布状态仍比较规则. 接下来, 当Nstep = 76356时, 岩样应力达到点II 翼型裂缝初步形成, 预制倾斜裂缝尖端附近仍然存在应力集中现象, 两个位移场的连续性波动程度进一步加深. 随后, 当Nstep = 91590时, 岩样应力达到点III全破坏物质点数量增加, 岩样的应力场和位移场的分布变得不再规则, C1和C2已经完成联通, C1和C3产生相互联通的趋势, 裂缝类型变得较为复杂, 包括翼型裂缝、抗翼型裂缝、准共面次生裂缝和倾斜次生裂缝. 最后, 当Nstep = 102985时, 岩样应力达到点IV 岩样发生宏观失效.
![]()
图 19 方案1的砂岩试样的压缩破坏形态Figure 19. Compressive failure process of sandstone specimen in scheme 1图20给出了室内试验结果[30]、已有数值模拟结果[29]和本文的数值模拟结果. 如图20所示, 预制裂缝的倾斜角度对裂缝类型几乎没有影响, 主要包括翼型裂缝、抗翼型裂缝、准共面次生裂缝和倾斜次生裂缝. 预制裂缝的贯通趋势是裂缝扩展的主要方向, 翼型裂缝和抗翼型裂缝是连接预制裂缝的主要裂缝.
![]()
图 20 不同倾角的预制裂缝扩展过程Figure 20. Propagation processes of prefabricated cracks with different inclination angles图21分别给出了多种物理量的变化趋势. 通过图21(a)可以发现, 砂岩试样的峰值强度先由85.9 MPa降为72.5 MPa, 然后逐渐增加至79.2 MPa, 这与室内实验结果的峰值强度规律变化是一致的(87.5, 74和79.9 MPa). 裂缝起裂应力的数值远低于峰值强度, 但呈现出与峰值强度相同的变化趋势, 先由24.4 MPa降为20.9 MPa, 然后逐渐增加至26.9 MPa, 这也与室内实验结果的峰值强度变化规律是一致的(19.1, 17.1和24 MPa). 通过图21(b)可以发现, 随着倾斜角度的增加, 预制裂缝首次贯通耗时的变化幅度并不明显. 方案3的耗时最长, 这与预制裂缝间距最远有关. 然而, 贯通式裂缝的长度(不包括预制裂缝长度)先由17.23 mm增加至21.51 mm, 然后增加至25.67 mm. 在方案3中, 两条抗翼型裂缝间还存在准共面次生裂缝的影响, 因此方案3的贯通式裂缝的长度是最大的.
4. 结论
本文提出了一种系统化的准态基近场动力学方法. 首先, 通过区分键型和搜索贡献物质点, 得到了基于键型区分的力密度计算方法, 实现了键的力密度和应力张量的同步计算; 然后, 引入最大拉应力准则和莫尔-库仑准则, 推导了两类适于当前近场动力学理论框架的单直线型虚拟裂缝模型, 实现了对材料多类型断裂模式的判断和对峰后阶段承载能力减弱的模拟; 接下来, 通过粒子网格孪生算法, 实现了不同介质状态下的物质点和孪生网格的同步计算, 通过孪生网格和势函数的组合, 对不同的接触情形和接触程度进行准确判断和表征; 最后, 根据库伦摩擦定律, 实现了考虑摩擦系数的摩擦力密度计算.
通过基准算例和复杂预制裂缝扩展算例, 充分证明了本文所提方法对非连续介质形成-演化过程的分析和处理能力. 与经典近场动力学相比, 本文所提方法与连续介质力学的联系更加紧密, 对断裂模式和峰后行为的处理能力更强, 对接触情形和摩擦情形的求解效果更好, 可以更加突出近场动力学的天然优势, 使其处理复杂受载条件下非连续介质形成-演化过程的能力得到了提升, 扩充了近场动力学的理论框架和适用范围.
-
![]()
图 1 键基近场动力学的材料构型和力密度作用方式
Figure 1. The material configuration and force density in the bond-based peridynamics
![]()
图 3 自适应力密度分量方向分配算法
Figure 3. An adaptive algorithm for force density component direction
![]()
图 4 适于准态基近场动力学的单直线虚拟裂缝模型和弹塑性本构模型
Figure 4. Single-linear fictitious crack models and elastic-plastic constitutive models suitable for quasi-state-based peridynamics
![]()
图 6 物质点的孪生网格的状态转化过程
Figure 6. Twin mesh state transformation process of material point
![]()
图 13 基准算例2的第1组 ~ 第3组试验的载荷-位移曲线
Figure 13. Loading-displacement curves of tests 1-3 in fundamental test 2
![]()
图 14 中心预制45°倾斜裂缝的矩形砂岩试样计算示意图
Figure 14. Calculation schematic diagram of a rectangular sandstone specimen with centrally prefabricated 45° inclined fracture
![]()
图 16 监测物质点的接触力密度的变化曲线
Figure 16. The curves of the contact force density of monitoring points
![]()
图 17 监测物质点的各类接触形式的数量和比例
Figure 17. Number and proportion of various contact forms at monitoring points
![]()
图 18 多倾斜裂缝的矩形砂岩试样计算示意图
Figure 18. Calculation schematic diagram of a rectangular sandstone specimen with multi-inclined fractures
![]()
图 19 方案1的砂岩试样的压缩破坏形态
Figure 19. Compressive failure process of sandstone specimen in scheme 1
![]()
图 20 不同倾角的预制裂缝扩展过程
Figure 20. Propagation processes of prefabricated cracks with different inclination angles
-
[1] 赵伦洋, 赖远明, 牛富俊等. 基于增量变分理论的多孔岩石多尺度弹塑性本构模型. 岩石力学与工程学报, 2022, 41(S2): 3163-3173 (Zhao Lunyang, Lai Yuanming, Niu Fujun, et al. A multi-scale elastoplastic constitutive model of porous rock based on incremental variational theory. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(S2): 3163-3173 (in Chinese) Zhao Lunyang, Lai Yuanming, Niu Fujun, et al. A multi-scale elastoplastic constitutive model of porous rock based on incremental variational theory. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(S2): 3163-3173 (in Chinese)
[2] 朱其志, 陈毓栋, 余健等. 不同状态灰砂岩应力松弛特性试验及本构模型研究. 河北工程大学学报(自然科学版), 2022, 39(1): 38-45 (Zhu Qizhi, Chen Yudong, Yu jian, et al. Experimental and constitutive modeling study on stress relaxation properties of sandstone in different states. Journal of Hebei University of Engineering (Natural Science Edition), 2022, 39(1): 38-45 (in Chinese) Zhu Qizhi, Chen Yudong, Yu jian, et al. Experimental and constitutive modeling study on stress relaxation properties of sandstone in different states. Journal of Hebei University of Engineering (Natural Science Edition), 2022, 39(1): 38-45 (in Chinese)
[3] 刘小宇, 杨政, 张慧梅. 准脆性材料抗压强度能量平衡尺寸效应模型. 力学学报, 2022, 54(6): 1613-1629 (Liu Xiaoyu, Yang Zheng, Zhang Huimei. Energy balance size effect model of compressive strength for quasi-brittle materials. Chinses Journal of Theoretical and Applied Mechanics, 2022, 54(6): 1613-1629 (in Chinese) Liu Xiaoyu, Yang Zheng, Zhang Huimei. Energy balance size effect model of compressive strength for quasi-brittle materials. Chinses Journal of Theoretical and Applied Mechanics, 2022, 54(6): 1613-1629 (in Chinese)
[4] 周小平, 贾志明. 场富集有限元方法模拟多裂纹起裂、扩展和连接过程. 岩土工程学报, 2022, 44(6): 988-996 (Zhou Xiaoping, Jia Zhiming. Field-enriched finite element method for numerical simulation of initiation, propagation and coalescence of multiple cracks. Chinese Journal of Geotechnical Engineering, 2022, 44(6): 988-996 (in Chinese) Zhou Xiaoping, Jia Zhiming. Field-enriched finite element method for numerical simulation of initiation, propagation and coalescence of multiple cracks. Chinese Journal of Geotechnical Engineering, 2022, 44(6): 988-996 (in Chinese)
[5] 赵伦洋, 赖远明, 牛富俊等. 硬脆性岩石多尺度损伤蠕变模型及长期强度研究. 中南大学学报(自然科学版), 2022, 53(8): 3071-3080 (Zhao Lunyang, Lai Yuanming, Niu Fujun, et al. Multi-scale damage creep model and long-term strength for hard brittle rocks. Journal of Central South University (Science and Technology), 2022, 53(8): 3071-3080 (in Chinese) Zhao Lunyang, Lai Yuanming, Niu Fujun, et al. Multi-scale damage creep model and long-term strength for hard brittle rocks. Journal of Central South University (Science and Technology), 2022, 53(8): 3071-3080 (in Chinese)
[6] 楚锡华. 颗粒材料的离散颗粒模型与离散——连续耦合模型及数值方法. [博士论文]. 大连: 大连理工大学, 2007 (Chu Xihua. The discrete particle and coupled discrete-continuum models and numerical methods for granular materials. [PhD Thesis]. Dalian: Dalian University of Technology, 2007 (in Chinese) Chu Xihua. The discrete particle and coupled discrete-continuum models and numerical methods for granular materials. [PhD Thesis]. Dalian: Dalian University of Technology, 2007 (in Chinese)
[7] 田锋. 基于体的三维块体接触算法研究及其在洞室围岩破坏过程中应用. [硕士论文]. 阜新: 辽宁工程技术大学, 2022 (Tian Feng. Study on the 3d block contact algorithm based on the volume and its application during the failure processes of cavern surrounding rock. [Master Thesis]. Fuxin: Liaoning Technical University, 2022 (in Chinese) Tian Feng. Study on the 3d block contact algorithm based on the volume and its application during the failure processes of cavern surrounding rock. [Master Thesis]. Fuxin: Liaoning Technical University, 2022 (in Chinese)
[8] 周小平, 王允腾, 钱七虎. 爆破荷载作用下岩石破坏特性的“共轭键”基近场动力学数值模拟研究. 中国科学: 物理学 力学 天文学, 2020, 50(2): 52-64 (Zhou Xiaoping, Wang Yunteng, Qian Qihu. Numerical simulations of failure characteristics of rock materials under blasting loads using the conjugated bond-pair-based peridynamics. Sci Sin-Phys Mech Astron, 2020, 50(2): 52-64 (in Chinese) Zhou Xiaoping, Wang Yunteng, Qian Qihu. Numerical simulations of failure characteristics of rock materials under blasting loads using the conjugated bond-pair-based peridynamics. Sci Sin-Phys Mech Astron, 2020, 50(2): 52-64 (in Chinese)
[9] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48: 175-209 doi: 10.1016/S0022-5096(99)00029-0
[10] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Computers and Structures, 2005, 83(17): 1526-1535
[11] Silling SA, Epton M, Weckner O, et al. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88: 151-184
[12] 朱其志, 李惟简, 尤涛等. 扩展键基近场动力学的几何非线性扩展. 同济大学学报(自然科学版), 2022, 50(4): 455-462 (Zhu Qizhi, Li Weijian, You Tao, et al. Geometrical nonlinear extension of extended bond-based peridynamic model. Journal of Tongji University (Natural Science), 2022, 50(4): 455-462 (in Chinese) Zhu Qizhi, Li Weijian, You Tao, et al. Geometrical nonlinear extension of extended bond-based peridynamic model. Journal of Tongji University (Natural Science), 2022, 50(4): 455-462 (in Chinese)
[13] Silling SA, Müge FC. Peridynamic model for microballistic perforation of multilayer graphene. Theoretical and Applied Fracture Mechanics, 2021, 113: 102947 doi: 10.1016/j.tafmec.2021.102947
[14] Hobbs M, Dodwell HJ. An examination of the size effect in quasi-brittle materials using a bond-based peridynamic model. Engineering Structures, 2022, 262: 1-18
[15] You T, Zhu QZ, Li PF, et al. Incorporation of tension-compression asymmetry into plastic damage phase-field modeling of quasi-brittle geomaterials. International Journal of Plasticity, 2020, 124: 71-95 doi: 10.1016/j.ijplas.2019.08.003
[16] 章青, 顾鑫, 郁杨天. 冲击载荷作用下颗粒材料动态力学响应的近场动力学模拟. 力学学报, 2016, 48(1): 56-63 (Zhang Qin, Gu Xin, Yu Yangtian. Peridynamics simulation for dynamic response of granular materials under impact loading. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56-63 (in Chinese) Zhang Qin, Gu Xin, Yu Yangtian. Peridynamics simulation for dynamic response of granular materials under impact loading. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56-63 (in Chinese)
[17] Wang LW, Sheng XY, Luo JB. A peridynamic frictional contact model for contact fatigue crack initiation and propagation. Engineering Fracture Mechanics, 2022, 264: 1-17
[18] Silling SA. EMU User’s Manual. Albuquerque: Sandia National Laboratories, 2004
[19] Parks ML, DJ, Littlewood JA, et al. Peridigm Users’ Guide. Albuquerque: Sandia National Laboratories, 2012
[20] Silling SA, Stewart A. Meshfree peridynamics for soft materials. 2016
[21] Kamensky D, Behzadinasab M, Foster JT, et al. Peridynamic modeling of frictional contact. Journal of Peridynamics and Nonlocal Modeling, 2019, 1(2): 107-121 doi: 10.1007/s42102-019-00012-y
[22] Zhang YN, Madenci E, Gu X, et al. A coupled hydro-mechanical peridynamic model for unsaturated seepage and crack propagation in unsaturated expansive soils due to moisture change. Acta Geotech, 2023, 18: 6297-6313
[23] Anicode SVK, Madenci E. Direct coupling of dual-horizon peridynamics with finite elements for irregular discretization without an overlap zone. Engineering with Computers, 2023, DOI: 10.1007/s00366-023-01800-3
[24] Hu Y, Madenci E. Bond-based peridynamics with an arbitrary poisson’s ratio//57th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2016
[25] Hu YL, Yu Y, Madenci E. Peridynamic modeling of composite laminates with material coupling and transverse shear deformation. Composite Structures, 2020, 253(1): 112760
[26] Zhou L, Esteban R, Bryan E, et al. A smooth contact algorithm for the combined finite discrete element method. Computational Particle Mechanics, 2020, 7(5): 1-15
[27] 侯艳丽, 张楚汉. 用三维离散元实现混凝土Ⅰ型断裂模拟. 工程力学, 2007, 24(1): 37-43 (Hou Yanli, Zhang Chuhan. Mode Ⅰ-fracture simulation of concrete based on 3D distinct element method. Engineering Mechanics, 2007, 24(1): 37-43 (in Chinese) Hou Yanli, Zhang Chuhan. Mode Ⅰ-fracture simulation of concrete based on 3D distinct element method. Engineering Mechanics, 2007, 24(1): 37-43 (in Chinese)
[28] 王学滨, 芦伟男, 钱帅帅等. 静水压力条件下开挖直径及卸荷时间对巷道围岩变形-开裂影响的连续-非连续方法模拟. 应用力学学报, 2020, 37(4): 1841-1848 (Wang Xuebin, Lu Weinan, Qian Shuaishuai, et al. Numerical simulation of effects of the excavation unloading time and chamber diameter on deformation-cracking processes of the chamber surrounding rock under hydrostatic pressure. Chinese Journal of Applied Mechanics, 2020, 37(4): 1841-1848 (in Chinese) Wang Xuebin, Lu Weinan, Qian Shuaishuai, et al. Numerical simulation of effects of the excavation unloading time and chamber diameter on deformation-cracking processes of the chamber surrounding rock under hydrostatic pressure. Chinese Journal of Applied Mechanics, 2020, 37(4): 1841-1848 (in Chinese)
[29] Wang YT, Zhou XP, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics. Engineering Fracture Mechanics, 2016, 163: 248-273 doi: 10.1016/j.engfracmech.2016.06.013
[30] Yang SQ, Yang DS, Jing HW, et al. An Experimental study of the fracture coalescence behaviour of brittle sandstone specimens containing three fissures. Rock Mechanics and Rock Engineering, 2011, 45: 563-582
-
期刊类型引用(0)
其他类型引用(1)
下载:
计量
- 文章访问数: 309
- HTML全文浏览量: 103
- PDF下载量: 93
- 被引次数: 1
