岩石三维裂纹扩展问题的近场动力学数值模拟

摘要

近场动力学（PD）方法基于键的形式处理断裂问题，对裂纹尖端可自动追踪，因而在求解三维裂纹扩展问题时具备极大的优势，但该方法中也存在数值震荡与边界误差等问题。为解决上述问题，首先应用无网格伽辽金弱形式框架对非常规态基近场动力学方法进行了探讨；随后引入近场动力学微分算子（PDDO）近似，并对比分析了该近似与重构核粒子（RKPM）近似间的差异，提出了具备更高数值精度的RKPM-PD耦合算法，并给出了该算法的隐式迭代流程。最后，通过若干数值算例验证了该算法在求解三维裂纹扩展问题上的有效性。

关键词

近场动力学; 无网格方法; 伽辽金法; 三维裂纹扩展

外部荷载作用下岩体内裂纹的扩展与贯通往往导致岩石工程发生失稳破坏，因此预测岩体中裂纹的扩展过程对于揭示岩体的变形和破坏规律以及评价岩土工程的安全性具有重要意义。目前，受困于研究手段的不足，对岩体裂纹扩展规律的研究仍集中于二维空间内。但实际工程中的岩体裂隙，无论深埋于岩体内部或只露于浅层表面，其本质都是三维裂纹扩展问题，因此，对岩体三维裂隙扩展的研究更具有实际意义^［1］。

预测三维裂纹的扩展是一项极具难度的问题。在理论研究方面，由于三维裂纹扩展问题的复杂性、多样性与求解上的巨大难度，至今仍未有理论上的相关突破。在室内试验方面， 一般采用类岩石材料或CT扫描方式对三维裂隙岩体的破裂机理进行研究^［2］，但由于岩体内部裂纹扩展观测手段的匮乏，室内试验尚无法满足工程问题的需要。

由于理论分析与室内试验的缺乏，学术界对三维裂纹扩展的研究集中于数值方法。有限元法（FEM）是最早用于模拟裂纹扩展问题的数值方法，但FEM在模拟裂纹扩展时需要不断地重新划分网格，极大地增加了计算工作量。对于复杂的三维裂纹问题，网格的切割划分更具挑战性。鉴于FEM的这些缺陷，出现了基于点近似的无网格方法^［3］。该方法不需要随裂纹扩展进行网格重构工作，降低了裂纹扩展问题的难度。但无网格方法中也引入了一些新的缺陷^［4］，如：①本质边界条件的施加相对复杂；②计算量更大；③配点型无网格计算稳定性差等。

鉴于传统数值求解方法在处理裂纹扩展问题时的不足，Silling^［5］提出了基于非局部作用思想的近场动力学（Peridynamic， PD）方法。该方法以非局部作用的积分模型代替传统理论的微分模型，避免了因在位移不连续处求导引起的奇异性问题。2007年Silling等^［6］提出了更为广泛的非常规态基PD（Non-ordinary state-based peridynamic， NOSB-PD）理论，该理论借助变形梯度张量$\mathit{F}$在PD框架下引入经典力学中的应力与应变，搭建了微观键与宏观强度准则间的关系，因而在岩土工程领域得到了广泛的应用^［7-8］。目前，NOSB-PD方法中仍存在因零能模式与边界效应等问题导致的计算精度不足的缺陷，且至今仍未有较为一致的解决方案^［9］。

为解决PD方法中的零能模式与边界误差问题，首先简介NOSB-PD模型的基本理论及应用无网格伽辽金格式对该模型的重构；随后，引入PD微分算子（PDDO）近似，并对比其与RKPM近似间的异同；基于对比结论，建立RKPM-PD耦合算法，并提出适用于该耦合算法的基于增量格式Newmark算法的隐式迭代流程；最后，采用RKPM-PD耦合算法对若干三维裂纹扩展问题进行模拟分析。

1 近场动力学理论

1.1　非常规态基近场动力学

在键基PD以及常规态基PD模型中并没有传统意义上应力、应变的概念，这为近场动力学的应用带来诸多不便。为解决这一问题，一种新的非常规态基近场动力学（Non-ordinary state-based peridynamic， NOSB-PD）方法被提出。在引入应力、应变概念后，不同材料的本构关系可以很方便地应用于近场动力学方法。在物质点$x_i$处，NOSB-PD方法的运动方程为

（1）

$\mathrm{式}\mathrm{中}:\rho$为密度，$\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-\mathrm{3}}$；下标$i\mathrm{、}j$分别代表物理量对应的物质点编号；u为位移；H为物质点$x_i$支撑域的范围；$\underline{\mathit{T}}$为力矢量状态，$\mathrm{N}·{\mathrm{m}}^{-\mathrm{6}}$；${\mathit{\xi }}_{ij}$为键，${\mathit{\xi }}_{ij}={\mathit{x}}_j-{\mathit{x}}_i$，其中x代表物质点的坐标向量，m；$\mathit{b}$为外力，$\mathrm{N}·{\mathrm{m}}^{-\mathrm{3}}$；$V$为体积，${\mathrm{m}}^{\mathrm{3}}$。在小变形情况下，$\underline{\mathit{T}}$可通过Cauchy应力张量$\mathit{\sigma }$表示，如式（2）：

（2）

式中：${\omega }_{ij}$为与键长$\left|{\mathit{\xi }}_{ij}\right|$有关的非负的影响函数；$\mathit{K}$为形状张量的逆矩阵，${\mathrm{m}}^{-\mathrm{5}}$。将式（2）代入式（1）中，小变形情况下NOSB-PD方程可写为

$$\rho {\ddot{\mathit{u}}}_i=\underset{}{\int }{}_H\left(\omega {\mathit{\sigma }}_i{\mathit{K}}_i{\mathit{\xi }}_{ij}+\omega {\mathit{\sigma }}_j{\mathit{K}}_j{\mathit{\xi }}_{ij}\right)\mathrm{d}{V}_j+{\mathit{b}}_i$$

（3）

式中：$\omega ={\omega }_{ij}={\omega }_{ji}$，${\mathit{\xi }}_{ij}=-{\mathit{\xi }}_{ji}$。

为将传统本构模型应用于近场动力学框架，NOSB-PD方法通过变形梯度张量$\mathit{F}$建立起PD理论与传统力学间的联系。在NOSB-PD模型中，变形梯度张量$\mathit{F}$被定义为

$$\mathit{F}=\underset{}{\int }{}_H{\omega }_{ij}\left({\mathit{\xi }}_{ij}+{\mathit{\eta }}_{ij}\right)\otimes {\mathit{\xi }}_{ij}\mathrm{d}{V}_j·\mathit{K}$$

（4）

式中：运算符$\otimes$代表并积运算；${\mathit{\eta }}_{ij}$代表键${\mathit{\xi }}_{ij}$的变形量，$\mathrm{m}$。可通过式（5）得到：

$$\mathit{K}={\mathit{A}}^{-\mathrm{1}}={\left(\underset{}{\int }{}_H{\omega }_{ij}{\mathit{\xi }}_{ij}\otimes {\mathit{\xi }}_{ij}\mathrm{d}\xi \right)}^{-\mathrm{1}}$$

（5）

为充分应用NOSB-PD方法与传统本构模型相耦合的优势，可在NOSB-PD方法中通过物质点的应力来定义键的损伤。定义键ξ_ij的最大主应力为键的2个端点最大主应力的平均值，当键上的最大主应力超过材料容许的最大强度时，该键发生断裂。对于岩石材料，可通过最大拉应力准则及摩尔-库仑准则来确定岩石材料的最大容许强度 ^［10］。

1.2　基于伽辽金格式的NOSB-PD方法重构

对于传统无网格伽辽金格式，在采用节点积分后，其最终方程可表示为

$$\mathit{M}\ddot{\mathit{U}}+\overline{\mathit{K}}\mathit{U}=\mathit{P}$$

（6）

$\mathbf{式}\mathbf{中}:\mathit{M}$为质量矩阵；$\overline{\mathit{K}}$为刚度矩阵；$\mathit{P}$为外荷载向量；$\mathit{U}$为所有节点的自由度向量。基于变形梯度张量$\mathit{F}$的求解公式（4），节点应变求解矩阵${\mathit{B}}_i$可表示为

$${\mathit{B}}_i=\left[\begin{array}{ccccccc}-{\displaystyle \sum _{j=\mathrm{1}}^m}{\omega }_{ij}{g}_x{V}_j& & {\omega }_{i\mathrm{1}}{g}_x{V}_{\mathrm{1}}& & \cdots & {\omega }_{in}{g}_x{V}_n& \\ & -{\displaystyle \sum _{j=\mathrm{1}}^m}{\omega }_{ij}{g}_y{V}_j& & {\omega }_{ij}{g}_y{V}_{\mathrm{1}}& \cdots & & {\omega }_{in}{g}_y{V}_n\\ -{\displaystyle \sum _{j=\mathrm{1}}^m}{\omega }_{ij}{g}_y{V}_j& -{\displaystyle \sum _{j=\mathrm{1}}^m}{\omega }_{ij}{g}_x{V}_j& {\omega }_{i\mathrm{1}}{g}_y{V}_j& {\omega }_{i\mathrm{1}}{g}_x{V}_{\mathrm{1}}& \cdots & {\omega }_{in}{g}_y{V}_n& {\omega }_{in}{g}_x{V}_n\end{array}\right]$$

（7）

式中：$\mathit{g}={\left[\begin{array}{cc}{g}_x& g_y\end{array}\right]}^{\mathrm{T}}=\mathit{K}{\mathit{\xi }}_{ij}$；$m$为$x_i$支撑域内节点$x_j$的数目。

将式（7）代入式（6），可得到式（6）与NOSB-PD方法的全局方程式（3）完全等价的结论^［11］。该结果揭示了NOSB-PD方法实质上是一种伽辽金弱形式方法。基于该结论对NOSB-PD方法中存在的零能模式与边界误差进行分析可得：①与其他无网格方法类似，NOSB-PD中的零能模式问题也是由于其在伽辽金弱形式下采用了节点积分；②NOSB-PD方法中的节点应变$\mathit{\epsilon }$是支撑域内最小二乘意义下的计算结果，而最小二乘下的近似一般不会满足δ属性，因此，NOSB-PD方法会在其边界处出现较大的计算误差。

2 近场动力学微分算子近似与重构核函数近似

在非常规态基近场动力学中变形梯度张量$\mathit{F}$求解方式的基础上，Madenci等^［12］提出了一种被称为近场动力学微分算子（PDDO）的理论，该理论同样继承了PD中以积分代替微分的思想，可通过积分形式求得高阶微分算子，可视为式（4）的扩展形式。

2.1　近场动力学微分算子理论

依据多元函数泰勒公式，在d维空间中，标量函数$u_j$在点$x_j$处展开为

$$u_j=\sum _{\left(n_{\mathrm{1}},\dots ,n_d\right)\in {\alpha }_d^n}\frac{r_{\mathrm{1}}^{n_{\mathrm{1}}}\cdots r_d^{n_d}}{n_{\mathrm{1}}!\cdots n_d!}{u}_{i,n_{\mathrm{1}}\cdots n_d}+{\rm O}\left(r^{\alpha +\mathrm{1}}\right)$$

（8）

$\mathrm{式}\mathrm{中}:$n为偏导数最高阶数；$n_d$为$d$维度的偏导数阶数；$r_d$为$d$维度的${\mathit{x}}_j$与${\mathit{x}}_i$间的距离。在式（8）中引入向量${\mathit{p}}_j、{\partial }_{\alpha }{\mathit{u}}_i$及对角矩阵$\mathit{C},\mathrm{得}\mathrm{到}$

$${\mathit{p}}_j={\left(\begin{array}{cccccc}\mathrm{1},& r_{\mathrm{1}},& \dots ,& r_{\mathrm{1}}^{n_{\mathrm{1}}}\cdots r_d^{n_d},& \dots ,& r_d^n\end{array}\right)}^{\mathrm{T}}$$

（9）

$${\partial }_{\alpha }{\mathit{u}}_i={\left(\begin{array}{cccccc}{u}_{i,\mathrm{0}\dots \mathrm{0}},& u_{i,\mathrm{0}\dots \mathrm{1}},& \dots ,& u_{i,n_{\mathrm{1}}\cdots n_d},& \dots ,& u_{i,n\dots \mathrm{0}}\end{array}\right)}^{\mathrm{T}}$$

$$\mathit{C}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\begin{array}{ccccc}\mathrm{1},& \dots ,& \frac{\mathrm{1}}{n_{\mathrm{1}}!\cdots n_d!},& \dots ,& \frac{\mathrm{1}}{n!}\end{array}\right]$$

显然，矩阵$\mathit{C}$满足$\mathit{C}={\mathit{C}}^{\mathrm{T}}$。代入式（9）后，泰勒展开式（8）可重新写为

$$u_j={\mathit{p}}_j^{\mathrm{T}}\mathit{C}{\partial }_{\alpha }{\mathit{u}}_i$$

（10）

在物质点$x_i$的整个支撑域（horizon）内，基于最小二乘原理可求得

$${\partial }_{\alpha }{u}_i={\int }_H{\omega }_{ij}\mathit{K}{\mathit{p}}_j{u}_j\mathrm{d}{V}_j$$

（11）

其中

$$\begin{array}{cc}\mathit{K}={\left(\mathit{A}\mathbf{C}\right)}^{-\mathrm{1}},& \mathit{A}={\int }_H{\omega }_{ij}{\mathit{p}}_j\otimes {\mathit{p}}_j\mathrm{d}{V}_j\end{array}$$

（12）

在保证矩阵$\mathit{A}$可逆的情况下，通过式（11）即可计算出物质点$x_i$上任意维度任意阶导数的PDDO近似。

2.2　二维空间下的PDDO近似

为便于理解，给出二维空间下一阶导数的PDDO近似形式。在二维线性基条件下，式（9）可简化为

$${\mathit{p}}_j={\left[\begin{array}{ccc}\mathrm{1}& x_j-x_i& y_j-y_i\end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}\mathrm{1}& {\xi }_x& {\xi }_y\end{array}\right]}^{\mathrm{T}}$$

（13）

$${\partial }_{\alpha }{\mathit{u}}_i={\left[\begin{array}{ccc}{u}_i& u_{i,x}& u_{i,y}\end{array}\right]}^{\mathrm{T}}$$

$$\mathit{C}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\left[\begin{array}{ccc}\mathrm{1}& \mathrm{1}& \mathrm{1}\end{array}\right]\right)$$

引入${\mathit{a}}_{\mathrm{00}}={\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}$，通过式（11）可求得位移$u_i$为

$$u_i={\int }_H{\omega }_{ij}{\mathit{a}}_{\mathrm{00}}^{\mathrm{T}}\mathit{K}{\mathit{p}}_j{u}_j\mathrm{d}{V}_j$$

（14）

对求解域进行离散，可得

$$\begin{array}{l}{u}_i={\displaystyle \sum _{j=\mathrm{1}}^m}{\omega }_{ij}{\mathit{a}}_{\mathrm{00}}^{\mathrm{T}}\mathit{K}{\mathit{p}}_j{V}_j{u}_j={\displaystyle \sum _{j=\mathrm{1}}^m}{\mathit{N}}_j{u}_j=\\ \mathit{N}\mathit{d}={\mathit{a}}_{\mathrm{00}}^{\mathrm{T}}\mathit{K}\mathit{H}\mathit{d}\end{array}$$

（15）

式中：$\mathit{N}$向量即为PDDO近似下的形函数，式（15）各分量为

$$\mathit{N}=\left[\begin{array}{cccc}{N}_{\mathrm{1}}& N_{\mathrm{2}}& \cdots & N_m\end{array}\right]$$

（16）

$$\mathit{d}={\left[\begin{array}{cccc}{u}_{\mathrm{1}}& u_{\mathrm{2}}& \cdots & u_m\end{array}\right]}^{\mathrm{T}}$$

$$\mathit{H}=\left[\begin{array}{cccc}{\omega }_{\mathrm{1}}{\mathit{p}}_{\mathrm{1}}{V}_{\mathrm{1}}& {\omega }_{\mathrm{2}}{\mathit{p}}_{\mathrm{2}}{V}_{\mathrm{2}}& \cdots & {\omega }_m{\mathit{p}}_m{V}_m\end{array}\right]$$

基于同样方式引入${\mathit{a}}_{\mathrm{10}}={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]}^{\mathrm{T}},{\mathit{a}}_{\mathrm{01}}={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}$，位移导数的离散形式为

$$\begin{array}{l}{u}_{i,x}={\mathit{N}}_{,x}\mathit{d}={\mathit{a}}_{\mathrm{10}}^{\mathrm{T}}\mathit{K}\mathit{H}\mathit{d}\\ u_{i,y}={\mathit{N}}_{,y}\mathit{d}={\mathit{a}}_{\mathrm{01}}^{\mathrm{T}}\mathit{K}\mathit{H}\mathit{d}\end{array}$$

（17）

式中： ${\mathit{N}}_{,x}$，${\mathit{N}}_{,y}$分别为位移一阶导数在PDDO方法下的形函数。

2.3　RKPM近似函数

基于Liu等^［13］研究成果，RKPM近似位移$u_i$可写为

$$u_i={\int }_H{\omega }_{ij}{\mathit{a}}_{\mathrm{00}}^{\mathrm{T}}\mathit{K}{\mathit{p}}_j{u}_j\mathrm{d}{V}_j$$

（18）

对比PDDO方案得到的位移近似函数式（14），可以看到PDDO与RKPM两方案得到的位移近似函数完全一致。

尽管PDDO与RKPM两方案得到的位移近似函数完全一致，但两者对位移导数近似的求解并不一致。RKPM方法中位移导数的近似函数通过对位移近似函数求导得到，即

$$u_{,x}={\mathit{N}}_{,x}\mathit{d}$$

$${\mathit{N}}_{,x}={\mathbf{a}}_{\mathrm{00}}^{\mathrm{T}}\left({\mathit{K}}_{,x}\mathit{H}+\mathit{K}{\mathit{H}}_{,x}\right)$$

（19）

2.4　PDDO近似与RKPM近似的性质对比

Krongauz与Belytschko ^［14］指出，为获得最优的数值精度，伽辽金格式中的近似函数应满足一致性条件与积分约束条件。

2.4.1　一致性条件

一般而言，一致性条件总是比较容易满足的。容易证得，PDDO与RKPM均满足一致性条件。在一维算例中，11个固定节点均匀布置在区间［0，10］上，插值点间距0.01，基向量为$\mathit{p}={\left[\begin{array}{cc}\mathrm{1}& x\end{array}\right]}^{\mathrm{T}}$，影响函数取为三次样条函数。假设节点位移为$y=x$，通过插值点计算得到的位移及位移导数的图像如图1所示。从图中可以看出，当基函数一致时，PDDO与RKPM计算得到的位移近似与位移导数近似具备相同阶数的一致性。

图1  一阶导数为常数时RKPM与PDDO近似的精度比较

Fig. 1  Comparison of approximate accuracy between RKPM and PDDO at $\mathit{y}=\mathit{x}$

2.4.2　积分约束条件

积分约束条件也被称作相容性（compatibility）条件，可写为

$${\int }_{\Omega }{\mathit{\Psi }}_{I,i}\left(\mathit{x}\right)\mathrm{d}\Omega ={\int }_{\Gamma }{\mathit{\Psi }}_I\left(\mathit{x}\right)n_i\mathrm{d}\Gamma$$

（20）

式中：$\mathit{\Psi }$为伽辽金方法中的位移形函数；$\Omega \mathrm{为}\mathrm{求}\mathrm{解}\mathrm{域};\Gamma \mathrm{为}\mathrm{求}\mathrm{解}\mathrm{域}\mathrm{边}\mathrm{界}\mathrm{。}$对于全局伽辽金格式，积分约束条件是强制性满足的。在RKPM近似中，位移导数的形函数是通过对位移形函数$\mathit{N}$进行求导运算而得，因此RKPM近似天然满足积分约束条件。而PDDO近似中位移导数形函数${\mathit{N}}_{,x}$与位移形函数$\mathit{N}$并不一定具有严格的求导关系，前者可能只是形函数导数的一种近似逼近，两者间并不具备严格的导数与原函数的关系。因此，一般情况下，积分约束条件式（20）未必满足。

为充分说明上述条件，取一维杆件进行分析。杆长$L=\mathrm{1}$，节点间距$\mathrm{d}x=\mathrm{0.1},$弹性模量$E=\mathrm{1},$一维杆件内有均匀分布的应力$\sigma$，为保持平衡，杆件左右两端需施加$t=\sigma$的外力荷载。显然，模型的平衡方程等价为式（20）。模型内力记为$f_{\mathrm{i}\mathrm{n}\mathrm{t}}$，并加以上标$\mathrm{P}\mathrm{D}\mathrm{D}\mathrm{O}$及$\mathrm{R}\mathrm{K}\mathrm{P}\mathrm{M}$以示区别；模型所受外力记为$f_{\mathrm{o}\mathrm{u}\mathrm{t}}$。计算结果如图2所示，可以看到，$f_{\mathrm{o}\mathrm{u}\mathrm{t}}=f_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\mathrm{R}\mathrm{K}\mathrm{P}\mathrm{M}}\ne f_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\mathrm{P}\mathrm{D}\mathrm{D}\mathrm{O}}$，即PDDO近似不满足积分约束条件，而RKPM近似满足该条件。

图2  一维杆件内外力分布

Fig. 2  Distribution of internal and external forces of one-dimensional bar

2.5　不同积分格式下RKPM与PDDO方法计算精度对比

为探究RKPM与PDDO方法在不同积分方案下的表现，采用与2.4.2节相同的体力作用下的一维杆模型，并分别采用节点积分方案及两点高斯积分方案对两近似方法的收敛性进行分析。在不同积分方案下，两方法的位移相对误差收敛性如图3所示。对比不同积分方案，可以看到：

图3  不同积分方案下PDDO与RKPM方法误差对比

Fig. 3  Error comparison of PDDO and RKPM methods in different integration schemes

（1）对于RKPM方法而言，采用高斯积分的精度远远高于节点积分。而对于PDDO方法而言，积分阶次的提高对结果的影响很小。

（2）在不同积分方案下，RKPM方法的精度均比PDDO方法的精度高。但具体而言，在节点积分时，两者的精度比较接近；在采用高斯积分的情况下，两者精度间的差距最为明显。

3 RKPM-PD耦合算法

3.1　概述

传统的无网格方法大多基于断裂力学准则处理裂纹扩展问题。在数值计算中，裂纹的追踪是非常棘手的问题，特别是在三维裂纹扩展问题中，对于裂纹面尖端的动态追踪几乎不可能实现。在PD方法中，裂纹不再需要显式的追踪，而是通过键上的损伤来表示，当键的临界伸长率超过一定极限时，两点之间的作用键发生断裂。在NOSB-PD方法中，也可认为当键上的应力超过允许应力时，两点之间的作用键发生断裂。

基于PDDO近似与RKPM近似间的异同以及PD键形式的裂纹扩展方案，提出一种RKPM-PD耦合算法，即非裂纹计算部分采用基于RKPM近似的无网格伽辽金格式以获取更精确的计算精度，裂纹处理部分则采用PD中基于键形式的损伤模型避免复杂裂纹面的追踪问题。

在RKPM-PD中，键的断裂准则与裂纹表示方法与NOSB-PD理论一致，其不同之处在于，RKPM-PD中除了更新节点与节点间键的状态外，还需要更新积分点与节点间键的状态。具体而言，在节点积分中只需考虑节点与节点间键的状态，但对于背景网格积分（高斯积分），模型中存在积分点与节点两类不同的物质点。按照PD键的概念，算法中高斯点与节点间的键、节点与节点间的键，这2类不同的键均需考虑其断裂与否。

3.2　增量格式Newmark算法

在时间迭代中，采用增量格式的Newmark算法进行计算。t时刻的Newmark方程可写为

$$\widehat{\mathit{K}}∆\mathit{U}=∆\widehat{\mathit{Q}}$$

（21）

$$\widehat{\mathit{K}}=\left(\mathit{K}+b_{\mathrm{0}}\mathit{M}\right)$$

$$∆\widehat{\mathit{Q}}=∆\mathit{Q}-\left(b_{\mathrm{1}}{\dot{\mathit{U}}}_t+b_{\mathrm{2}}{\ddot{\mathit{U}}}_t\right)$$

（22）

式中：$b_{\mathrm{0}}\mathrm{、}{b}_{\mathrm{1}}\mathrm{、}{b}_{\mathrm{2}}$均为Newmark迭代参数。在计算出$∆\mathit{U}$后，即可得到${\mathit{U}}_{t+∆t}$的值，并通过${\mathit{U}}_{t+∆t}$计算出新的${\mathit{\sigma }}_{t+∆t}$。之后，根据不同的断裂准则，利用$t+∆t$时刻状态信息判断是否有新的键断掉。对于是否有新键断裂，可分为：

（1）如果当前步中没有新的断键出现，说明$t+∆t$时刻方程达到平衡，可进入下一步计算。在通过$∆\mathit{U}$得到$t+∆t$时刻的速度及加速度后，进行变量回代，开启下一时步的计算。

（2）如果当前步中有新的断裂出现，说明$t+∆t$状态并未达到平衡，方程（21）仍需进行迭代计算。此时，对于支撑域内出现新断键的物质点，更新其键的状态，并重新计算该点上的形函数$\left[\mathit{N}\right]$及对应的包含一阶导数的矩阵$\left[\mathit{B}\right]$。在新的断裂出现后，体系内的不平衡力$\mathit{R}$为

$$\mathit{R}=∆\widehat{\mathit{Q}}-\int {\left[\mathit{B}\right]}^{\mathrm{T}}\left[\mathit{\sigma }\right]\mathrm{d}V$$

（23）

式中：$\left[\mathit{\sigma }\right]$为断键后的应力，通过$\left[\mathit{B}\right]$与${\mathit{U}}_{t+∆t}$得到。同时，新的有效刚度矩阵［$\widehat{\mathit{K}}$］为

$$\left[\widehat{\mathit{K}}\right]=\left[\mathit{K}\right]+b_{\mathrm{0}}\mathit{M}$$

（24）

式中：$\left[\mathit{K}\right]$为新的全局刚度矩阵，可通过矩阵［B］求得。新的位移增量$\left[∆\mathit{U}\right]$为

$$\left[∆\mathit{U}\right]={\left[\mathit{K}\right]}^{-\mathrm{1}}\times \mathit{R}$$

（25）

显然，$t+∆t$时刻的位移为

$${\mathit{U}}_{t+∆t}={\mathit{U}}_t+∆\mathit{U}+\left[∆\mathit{U}\right]$$

（26）

根据式（26）的${\mathit{U}}_{t+∆t}$重新计算出${\mathit{\sigma }}_{t+∆t}$，并继续判断是否有新的断键出现，如果依然出现新的断键，则返回步骤（2）中继续迭代，直至系统内不再出现新的断键。此时，结束当前时步内的循环，转至步骤（1）中。

4 三维裂纹扩展模拟计算

4.1　含倾斜预制裂纹的三维受拉平板

为体现RKPM-PD方法在处理三维裂纹扩展问题中的优势，考虑如图4所示模型，其几何尺寸为$\mathrm{2}\mathrm{m}\times \mathrm{1}\mathrm{m}\times \mathrm{0.2}\mathrm{m}$，模型下表面（y=0m）固定，上表面（y=2m）受垂直于该平面的拉伸载荷。模型左侧被一裂纹面完全贯穿，裂纹深度为0.2m， 并与Oxy平面呈45°夹角。模型节点间距为0.03m，支撑域半径为2倍节点间距。模型材料力学参数为：弹性模量$E=\mathrm{30}\mathrm{G}\mathrm{P}\mathrm{a},$泊松比$\nu =\mathrm{0.2},$密度$\rho =\mathrm{2}\mathrm{500}$ kg·m^-3。裂纹扩展采用临界伸长率准则，键的临界伸长率$s_{\mathrm{c}}=\mathrm{0.001}$。计算采用Newmark积分，时间步长$∆t=\mathrm{0.1}\mathrm{m}\mathrm{s}$， 载荷$∆\sigma =\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{P}\mathrm{a}$。

图4  三维含倾斜裂纹受拉平板示意

Fig. 4  Sketch map of 3D tensile plate with inclined crack

图5中给出预制倾斜裂纹面在受拉平板内的扩展过程。可见，随着拉应力的增大，裂纹首先在倾斜裂纹面的边缘向前扩展，随后，裂纹面从倾斜状态逐渐扭转为竖直状态；之后，裂纹面依然保持垂直于Oxy面的竖直状态沿x方向向前延伸，最终贯穿整个平板。从图中可见，裂纹穿透平板位置约在y=1处，该值近似于倾斜裂纹面y坐标的中心值。另外，在图6中给出了对应状态下平板内y方向位移的分布，可见，位移分布与裂纹面的扩展状态相吻合。

图5  含倾斜裂纹平板在拉伸载荷作用下的裂纹扩展路径

Fig. 5  Propagations of 3D plate with inclined crack at tensile load

图6  含倾斜裂纹平板在拉伸载荷作用下的纵向位移分布

Fig. 6  Snapshots of vertical displacement of 3D plate with inclined crack

4.2　含预制裂纹三维巴西圆盘

考虑不同厚度的带预制裂纹的三维巴西圆盘试验。如图 7所示，模型在Oxy平面内的直径为100mm，z方向的厚度为10mm（厚度方向为垂直于纸面的方向）。圆盘上下两端受压，采用Newmark隐式积分方案，时间步长取为$∆t=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}\mathrm{s}$，平均加载速度为每步${\mathrm{10}}^{-\mathrm{8}}\mathrm{m}$。岩石试样力学参数为：弹性模量$E=\mathrm{10}\mathrm{G}\mathrm{P}\mathrm{a},$泊松比$v=\mathrm{0.25},$密度$\rho =\mathrm{2}\mathrm{500}$ kg·m^-3，抗拉强度$T=\mathrm{0.5}\mathrm{M}\mathrm{P}\mathrm{a},\mathrm{黏}$聚力$C=\mathrm{5}\mathrm{M}\mathrm{P}\mathrm{a}$以及内摩擦角$\varphi ={\mathrm{40}}^{\mathrm{o}}$。离散模型中平均节点间距$h=\mathrm{2}\mathrm{m}\mathrm{m}$。在本算例中，预制裂纹长度为30mm， 与水平方向呈${\mathrm{45}}^{\mathrm{o}}$夹角。为探究预制裂纹深度对结果的影响，预制裂纹深度分别取10mm、6mm及2mm。

图 7  带预制裂纹的三维巴西圆盘示意

Fig. 7  Sketch map of 3B Brazilian disk with prefabricated cracks

当预制裂纹深度为10mm时，预制裂纹在厚度方向完全贯穿平板，模型在z方向完全一致。该算例的裂纹扩展路径如图 8所示，可以看到，裂纹扩展路径与二维情况下相一致，裂纹尖端分别向圆盘上下两端的压缩载荷施加处扩展，直至最后裂纹完全贯穿整个圆盘。如图 9所示，该模型采用RKPM-PD的计算结果与在三维NMM方法下的模拟结果^［10］及试验结果^［15］均保持一致。

图 8  三维巴西圆盘含完全贯穿裂纹的裂纹扩展路径

Fig. 8  Snapshots of crack path within 3D Brazilian disk with a fully penetrating crack

图 9  含完全贯穿裂纹的三维巴西圆盘在NMM方法及室内试验下的裂纹扩展路径

Fig. 9  Crack path of 3D Brazilian disk with a fully penetrating crack obtained by 3D NMM and experiment results

通过对裂纹完全贯穿试件的分析，证明该方法在三维裂纹问题中的有效性。然后考虑裂纹未完全贯穿的情况，设置裂纹嵌入深度分别为6mm、2mm，其结果分别显示于图10与图11。为显示裂纹未完全贯穿的情况，对于同一时刻的裂纹扩展结果分别给出其正面图（z=0）与背面图（z=10mm）。其中，裂纹从正面嵌入，终止于圆盘内部，未穿透背面。

图 10  三维巴西圆盘含嵌入深度6mm未贯穿裂纹的裂纹扩展路径

Fig. 10  Snapshots of crack path within 3D Brazilian disk with a penetrating crack (insert depth=6mm)

图 11  三维巴西圆盘含嵌入深度2mm未贯穿裂纹的裂纹扩展路径

Fig. 11  Snapshots of crack path within 3D Brazilian disk with a penetrating crack (insert depth=2mm)

通过对未穿透裂纹的结果分析可知，裂纹首先沿厚度方向扩展，从背面可逐步看到损伤区域的形成；之后，正面与背面的裂纹均沿竖直方向朝两端扩展。对比2个面的裂纹最终扩展轨迹，可以看到两者大体一致，细微处有所差别。同时可以看到，裂纹切割深度越浅，启裂所需计算步越多。图 12给出了3种不同切割深度下圆盘最顶端3列节点的轴向平均应力-应变曲线，从图中可以清晰看出，巴西圆盘在表面裂纹（2mm）时所能承受的载荷远远大于深度裂纹（6mm）及完全贯穿裂纹（10mm）时的载荷。

图 12  三维巴西圆盘算例中裂纹扩展过程的应力-应变曲线

Fig. 12  Axial stress-strain curves of crack propagation in 3D Brazilian disk at different cut depths

5 结语

针对NOSB-PD方法面临的零能模式与边界误差问题，提出了一种具备更高精度的RKPM-PD耦合算法，建立了一套适用于该耦合算法的隐式求解方案，并成功将其应用于三维裂纹扩展问题中。主要结论如下：

（1）NOSB-PD方法实质上等价于采用节点积分的伽辽金弱形式方法。NOSB-PD中的零能模式产生的原因是其在弱形式下采用了节点积分形式。

（2）PDDO近似与RKPM近似均满足一致性条件，但仅有RKPM近似满足相容性条件。在提高积分计算阶次后，基于RKPM近似函数的伽辽金方法计算精度有较为明显的提升，但基于PDDO近似函数的伽辽金方法计算精度并无明显提高。

（3）RKPM-PD耦合算法不受零能模式与边界效应的困扰，具备更高的计算精度，并可有效求解三维动态裂纹扩展问题。

作者贡献声明

崔 昊：具体研究工作的开展和论文撰写。

郑 宏：论文的选题、指导和修改。

李春光：公式推导部分的指导。

韩 月：论文修改。