场富集有限元方法模拟多裂纹起裂、扩展和连接过程

周小平; 贾志明

doi:10.11779/CJGE202206002

场富集有限元方法模拟多裂纹起裂、扩展和连接过程 English Version

重庆大学土木工程学院, 重庆 400045

基金项目:

国家自然科学基金项目 51839009

详细信息

作者简介:
周小平(1970—)，男，博士，教授，博士生导师，主要从事岩石力学与工程方面的研究工作。E-mail：xiao_ping_zhou@126.com

中图分类号: TU451
计量
- 文章访问数: 222
- HTML全文浏览量: 34
- PDF下载量: 200
出版历程
- 收稿日期: 2021-07-11
- 网络出版日期: 2022-09-22
- 刊出日期: 2022-05-31

Field-enriched finite element method for numerical simulation of initiation, propagation and coalescence of multiple cracks

School of Civil Engineering, Chongqing University, Chongqing 400045, China

摘要

摘要: 研究含多裂纹岩石的力学响应和裂纹行为对岩体结构的设计与稳定性分析具有重要的指导意义。提出了场富集有限元方法研究脆性岩石材料中多裂纹的演化规律，包括裂纹的起裂、扩展和连接过程。提出了多裂纹在扩展过程中出现的各种交汇情况的解决方案。场富集有限元方法可以直接处理复杂的交叉裂纹，而不需要像扩展有限元一样引入额外的富集函数。通过数值算例，充分说明了场富集有限元方法在处理各种复杂裂纹系统扩展演化方面的能力。
- 脆性岩石材料 /
- 场富集有限元方法 /
- 裂纹扩展和连接 /
- 复杂多裂纹
Abstract: The study on the mechanical response and cracking behaviors of brittle rock materials with multiple cracks is of vital significance for the design and stability analysis of rock engineering structures. A field-enriched finite element method (FE-FEM) is proposed to study the evolution behaviors of multiple cracks in rock materials, including crack initiation, propagation and coalescence. The solutions to the crack coalescence problem during simulation are proposed. The field-enriched finite element method can directly deal with the complex multiple crack problem, while the extra enriched function needs to be introduced in the extended finite element method (XFEM). The analytical results of the present numerical examples demonstrate that the proposed numerical method has the capability to handle complex multiple crack propagation and coalescence.
- brittle rock material /
- field-enriched finite element method /
- crack propagation and coalescence /
- multiple cracks

HTML全文

0. 引言

岩体中存在的各种原生裂隙和节理等不连续面将岩体切割为不连续体，同时这些不连续面控制着岩体的强度和变形。岩体在开挖卸荷条件下的变形破坏往往是这些不连续面演化的结果，它们的扩展和连接直接影响着岩体工程的稳定性。因此，研究多裂纹的扩展演化机制对岩体稳定性控制具有重要的指导意义。

为理解多裂纹扩展和连接的机理，国内外学者已经做了大量的理论分析和试验研究。在理论方面，Freji-Ayoub等^[1]对含随机裂纹结构的理论解进行了推导。Rybaczuk等^[2]将分形方法应用到多裂纹扩展演化的理论研究。Chen^[3]提出了针对规则和随机分布裂纹的奇异积分方程方法。在试验方面，Chen等^[4]采用实时监测手段对含多裂纹的岩石类材料进行了单轴压缩试验，总结了5种搭接模式。Zhang等^[5]采用声-光结合的监测手段对单轴压缩下含多裂纹岩石的裂纹扩展进行了研究。

这些理论和试验的成果为数值模拟多裂纹演化规律提供了参考。由于多裂纹扩展涉及复杂的力学和几何问题，数值模拟方法成为研究多裂纹演化机制的有力工具。大量的数值模拟方法已经被用于分析多裂纹问题。例如，Carpinteri等^[6]采用边界元方法(BEM)模拟多裂纹的扩展问题。Denda等^[7]开发了一种复杂变量方法，利用边界元法对中心裂纹和边裂纹问题进行了数值模拟。Budyn等^[8]利用扩展有限元方法(XFEM)模拟复杂多裂纹的演化过程，追踪裂纹的扩展路径。Zhou等^[9]采用了一种细化的交汇策略处理多裂纹之间的交汇问题。Wang等^[10]提出了共轭键基近场动力学模型处理多裂纹扩展和交汇，并且模拟结果与试验结果匹配良好。Azadi等^[11]采用一种自适应重划分网格方法模拟任意形状的多裂纹扩展，比较了各种裂纹扩展准则和各自的裂纹路径。Zhou等^[12]提出了非局部晶格模型模拟脆性岩石中多裂纹的扩展演化。Zhou等^[13]提出了广义粒子动力学方法模拟多裂纹的扩展。Jia等^[14]在传统相场方法框架下提出了基于Hoek–Brown准则的相场演化方程，研究了压缩条件下的多裂纹演化以及搭接模式。Xu等^[15]提出了一种线性无关的高阶扩展数值流形方法用于模拟多裂纹的演化。另外，无网格方法^[16-18]由于插值不依赖于网格信息，被用于研究多裂纹扩展问题。

本文提出了场富集有限元方法研究多裂纹的扩展演化模式和相互作用机理，采用预测–修正的方法来处理裂尖交汇问题，包括裂尖与裂尖交汇、裂尖与裂纹段交汇、裂尖与自由边界交汇。通过引入加载乘子处理竞争裂纹扩展问题。通过典型数值算例说明该算法的有效性。与相场法相比，场富集有限元方法计算效率大大提高，即使在加载点或支撑端存在应力集中，也可以降低损伤演化的敏感性，有效避免非物理裂纹；同时，由于场演化方程与力平衡方程之间没有耦合关系，数值收敛性大大提高，网格的依赖性相对较弱；另外，场富集有限元方法可以灵活使用各种断裂准则来模拟问题的实际情况。与扩展有限元相比，场富集方法无需引入水平集和富集函数来描述裂纹，这大大简化了繁琐的编程，并减少了为处理各种断裂问题而引入的各种富集函数的数量。

1. 场富集有限元法理论与数值离散

场富集有限元方法将断裂力学理论与连续损伤力学理论相联系，通过引入一个损伤场区分材料的破坏与完整，将裂纹处理为损伤带，通过数值离散，该损伤场与位移场共同控制裂纹的演化。损伤场对本构模型中的刚度影响是通过衰减函数来实现的，从而影响位移场的分布。在计算的位移场和应力场的基础上，采用相互作用积分的域形式求解裂尖的应力强度因子，通过采用最常见的最大周向应力准则判断裂纹的扩展方向，给定一个小的裂纹扩展增量，引入加载乘子保证力学平衡和断裂准则条件同时满足。扩展后的裂纹提供损伤场的富集依据，损伤场采用指数距离加权方法进行富集。

1.1 基本理论

考虑一个弹性体 $\varOmega$ 由外边界 $\partial \varOmega$ 包围，位移和力边界分别表示为 $\partial {\varOmega _{\text{u}}}$ 和 $\partial {\varOmega _{\text{p}}}$ ，它们满足条件： $\partial {\varOmega _{\text{u}}} \cup \partial {\varOmega _{\text{p}}} = \partial \varOmega$ 和 $\partial {\varOmega _{\text{u}}} \cap \partial {\varOmega _{\text{p}}} =\varnothing$ 。内部裂纹表示为 $\varGamma$ 。

场富集方法的静态平衡方程和边界条件为

$\nabla \cdot \boldsymbol{\sigma}+\boldsymbol{b}=0(\text { in } \varOmega),$

(1a)

$\boldsymbol{\sigma} \cdot n=\boldsymbol{f} \left(\text { on } \partial \varOmega_{\mathrm{p}}\right),$

(1b)

$\boldsymbol{u}=\overline{\boldsymbol{u}}\left(\text { on } \partial \varOmega_{\mathrm{u}}\right),$

(1c)

式中， $\boldsymbol{b}$ ，f为体力和面力，u为物体位移，n为边界的外法向单位矢量，包含损伤的柯西应力 $\boldsymbol{\sigma }$ 表示为

$\boldsymbol{\sigma}=g(d) \boldsymbol{\sigma}_{0}=g(d) C_{0}: \boldsymbol{\varepsilon},$

(2)

式中， ${\boldsymbol{\sigma} _{\text{0}}}$ 为未损伤的柯西应力张量， ${C_0}$ 为线弹性刚度张量，小应变张量为 $\boldsymbol{\varepsilon} = 1/2(\nabla \boldsymbol{u} + \nabla {\boldsymbol{u}^{\text{T}}})$ 。衰减函数 $g(d)$ 通常定义为 $g(d) = {(1 - d)^2}$ ，它将损伤场 $d$ 与材料的力学性质耦合起来，通过改变损伤场来控制刚度衰减。 $d = 1$ 对应完全破坏的材料，而 $d = 0$ 对应完整的材料。衰减函数的形式可以任意选择，但需要满足以下条件：①非负性；②在区间[0, 1]上单调递减；③ $g(0) = 1$ ，且 $g(1) = g'(1) = 0$ 。

1.2 数值离散

式（1a）~（1c）可以通过传统的有限元进行离散求解。控制方程的弱形式表示为

$\int_{\varOmega} \boldsymbol{\sigma}: \delta \boldsymbol{\varepsilon} \mathrm{d} \varOmega-\int_{\varOmega} \boldsymbol{b} \cdot \delta \boldsymbol{u} \mathrm{d} \varOmega-\int_{\partial \varOmega_{\mathrm{p}}} \boldsymbol{f} \cdot \delta \boldsymbol{u} \mathrm{d} S=0。$

(3)

在二维背景下，将数值模型离散为4节点的四边形单元(m=4)，采用Voigt记法，位移场和损伤场的离散近似表示为

$\left.\begin{array}{l} \boldsymbol{u}=\sum\limits_{i=1}^{m} \boldsymbol{N}_{i}^{\boldsymbol{u}} \boldsymbol{u}_{i} , \\ d=\sum\limits_{i=1}^{m} N_{i} d_{i}, \end{array}\right\}$

(4)

式中，

$\left.\begin{array}{l} \boldsymbol{N}_{i}^{\mathrm{u}}=\left[\begin{array}{cc} N_{i} & 0 \\ 0 & N_{i} \end{array}\right], \\ \boldsymbol{u}_{i}=\left[\begin{array}{c} u_{x} \\ u_{y} \end{array}\right]_{i}, \end{array}\right\}$

(5)

其中，N_i为节点i处的形函数，u_i，d_i为节点处的位移和损伤。

位移场和损伤场对应的导数( $\boldsymbol{\varepsilon}$ ， $\nabla d$ )表示为

$\left.\begin{array}{l} \boldsymbol{\varepsilon}=\sum\limits_{i=1}^{m} \boldsymbol{B}_{i}^{\mathrm{u}} \boldsymbol{u}_{i}, \\ \nabla d=\sum\limits_{i=1}^{m} \boldsymbol{B}_{i}^{\mathrm{d}} d_{i}, \end{array}\right\}$

(6)

式中，

$\left. \begin{array}{l} &\boldsymbol{B}_{i}^{\mathrm{u}}=\left[\begin{array}{cc} N_{i, x} & 0 \\ 0 & N_{i, y} \\ N_{i, y} & N_{i, x} \end{array}\right], \\ &\boldsymbol{B}_{i}^{\mathrm{d}}=\left[\begin{array}{c} N_{i, x} \\ N_{i, y} \end{array}\right] \quad 。 \end{array} \right\}$

(7)

将上面的离散形式引入到控制方程的弱形式中，并且由于检验函数的任意性，离散的弱形式表示为

$\boldsymbol{f}_{i}^{\mathrm{int}}=\boldsymbol{f}_{i}^{\mathrm{ext}},$

(8)

式中，

$\boldsymbol{f}_{i}^{\mathrm{int}}=\int_{\varOmega}(1-d)^{2}\left(\boldsymbol{B}_{i}^{\mathrm{u}}\right)^{\mathrm{T}} \boldsymbol{\sigma}_{0} \mathrm{~d} \varOmega=K_{i j} \boldsymbol{u}_{j},$

(9)

$\boldsymbol{f}_{i}^{\mathrm{ext}}=\int_{\varOmega}\left(\boldsymbol{N}_{i}^{\mathrm{u}}\right)^{\mathrm{T}} \boldsymbol{b }\mathrm { d } \varOmega+\int_{\partial \varOmega_{\mathrm{p}}}\left(\boldsymbol{N}_{i}^{\mathrm{u}}\right)^{\mathrm{T}} \boldsymbol{f} \mathrm{d} S,$

(10)

其中，

$\boldsymbol{K}_{i j}=\frac{\partial \boldsymbol{f}_{i}^{\mathrm{int}}}{\partial \boldsymbol{u}_{j}}=\int_{\varOmega}(1-d)^{2}\left(\boldsymbol{B}_{i}^{\mathrm{u}}\right)^{\mathrm{T}} C_{0} \boldsymbol{B}_{j}^{\mathrm{u}} \mathrm{d} \varOmega 。$

(11)

因此，控制方程的离散形式可以写为

$\boldsymbol{K}_{i j} \boldsymbol{u}_{j}=\boldsymbol{f}_{i}^{\mathrm{ext}}。$

(12)

1.3 应力强度因子计算

本文采用相互作用积分的域形式^[19-20]求解应力强度因子。考虑一个含裂纹岩石的两个状态，状态1( $u_i^{(1)}$ ， $\varepsilon _{ij}^{(1)}$ ， $\sigma _{ij}^{(1)}$ )对应于当前状态，而状态2( $u_i^{(2)}$ ， $\varepsilon _{ij}^{(2)}$ ， $\sigma _{ij}^{(2)}$ )是一个附加状态。

相互作用积分的域形式为

$I^{(1, 2)}=\int_{A}\left(-W^{(1, 2)} \delta_{1 j}+\sigma_{i j}^{(1)} \frac{\partial u_{i}^{(2)}}{\partial x_{1}}+\sigma_{i j}^{(2)} \frac{\partial u_{i}^{(1)}}{\partial x_{1}}\right) \frac{\partial q}{\partial x_{j}} \mathrm{~d} A,$

(13)

式中，A为积分面积，如图 1所示。q为一个定义在积分域上的光滑权函数，在围线S₁上取值为1，在围线S₂上取值为0。

图 1 相互作用积分的域形式

Figure 1. Domain form of interaction integral

下载: 全尺寸图片幻灯片

W^{(1, 2)}是相互作用应变能，

$W^{(1, 2)}=\sigma_{i j}^{(1)} \varepsilon_{i j}^{(2)}=\sigma_{i j}^{(2)} \varepsilon_{i j}^{(1)}。$

(14)

对一般的混合模式问题，相互作用积分I^{(1, 2)}和应力强度因子的关系：

${I^{(1, 2)}} = \frac{2}{{{E^{\text{*}}}}}(K_{\text{I}}^{(1)}K_{\text{I}}^{(2)} + K_{{\text{II}}}^{(1)}K_{{\text{II}}}^{(2)}),$

(15)

式中，E*由弹性模量E和泊松比 $\nu$ 定义为

$E^{*}=\left\{\begin{array}{ll} E & (\text { 平面应力 }) \\ E /\left(1-v^{2}\right) & \text { (平面应变 }) \end{array}\right. 。$

(16)

通过合理选取附加状态2的近似解，可以分别得到Ⅰ型或Ⅱ型情况下的应力强度因子和相互作用积分的关系。取纯Ⅰ型的近似解， $K_{\text{I}}^{(2)} = 1$ 和 $K_{{\text{II}}}^{(2)} = 0$ ，得到：

$K_{\text{I}}^{(1)} = \frac{2}{{{E^{\text{*}}}}}{I^{(1, {\text{ModeI}})}}。$

(17)

同样地可得

$K_{{\text{II}}}^{(1)} = \frac{2}{{{E^{\text{*}}}}}{I^{(1, {\text{ModeII}})}} 。$

(18)

1.4 裂纹扩展准则

本文采用最大周向应力准则^[21]，该准则认为裂纹将沿着使周向应力 ${\sigma _{\theta \theta }}$ 最大的方向扩展。该周向应力是一个主应力，因此通过使剪应力为0可以确定裂纹扩展角 ${\theta _c}$ 。

在一般混合模式载荷下，靠近裂尖的剪应力有如下形式：

${\sigma _{r\theta }} = \frac{1}{{\sqrt {2{{\rm{ \mathsf{ π} }}}r} }}\cos \left( {\frac{\theta }{2}} \right)\left( {\frac{1}{2}{K_{\text{I}}}{\text{sin}}\theta + \frac{1}{2}{K_{{\text{II}}}}(3{\text{cos}}\theta - 1)} \right) = 0。$

(19)

求解这个方程得到：

$\theta_{\mathrm{c}}= \begin{cases}2 \arctan \frac{1}{4}\left(\frac{K_{\mathrm{I}}}{K_{\mathrm{II}}} \pm \sqrt{\frac{K_{\mathrm{I}}^{2}}{K_{\mathrm{II}}}+8}\right) & \left(K_{\mathrm{II}} \neq 0\right) \\ 0 & \left(K_{\mathrm{II}}=0\right)\end{cases} 。$

(20)

当裂纹开始扩展时，等效应力强度因子 ${K_{{\text{eq}}}}$ 和材料的断裂韧度满足下面关系：

${K_{{\text{eq}}}} \leqslant {K_{{\text{IC}}}},$

(21)

${K_{{\text{eq}}}} = {K_{\text{I}}}{\cos ^3}\frac{{{\theta _{\text{c}}}}}{2} - \frac{3}{2}{K_{{\text{II}}}}\cos \frac{{{\theta _{\text{c}}}}}{2}\sin {\theta _{\text{c}}}。$

(22)

2. 场富集方法及富集节点的搜索方法

2.1 场富集方法

本文通过“指数距离加权方法”定量富集节点处的损伤值。指数距离加权法可以有效地将损伤变量的分布限定在一个非常小的窄带内，并且可以实现损伤值随着偏离裂纹的距离的增加而快速衰减，满足损伤变量衰减分布的要求。距离指数函数的形式为 $d(s) = {{\text{e}}^{ - \left| {s/l} \right|}}$ ，s是某个节点与裂纹的距离，l是一个改变裂纹损伤程度的控制因子^[22]。

2.2 富集节点的搜索方法

如图 2所示，采用盒形和圆形搜索策略追踪完全损伤窄带。具体地说，采用盒形搜索策略找出位于完全损伤窄带内的节点，这些节点紧邻真实裂纹段。同时，通过圆形搜索策略补充两个盒形区域交界处的节点，从而形成一个连续变化的完全损伤带。

图 2 富集节点检索方式

Figure 2. Searching strategy of enriched nodes

下载: 全尺寸图片幻灯片

3. 裂纹扩展

3.1 裂纹交汇处理方法

场富集方法处理裂纹交汇问题不需要像XFEM那样额外增加交汇富集来处理交汇时的不连续。下面对各种情况的处理方式展开叙述。

（1）裂尖与裂尖交汇

如所示，两条裂纹各自独立计算扩展方向，裂纹1的裂尖从 ${O_1}$ 扩展到 ${O'_1}$ ，同时裂纹2的裂尖从 ${O_2}$ 扩展到 ${O'_2}$ ，计算 ${O'_1}{O'_2}$ 之间的距离 $\left| {\overline {{{O'}_1}{{O'}_2}} } \right|$ ，当 $\left| {\overline {{{O'}_1}{{O'}_2}} } \right|$ 小于定义的容差时，两个裂尖搭接，并将这两个裂尖冻结，之后不再继续扩展。

图 3 裂尖与裂尖交汇方式

Figure 3. Intersection strategy of two crack tips

下载: 全尺寸图片幻灯片

（2）裂尖与裂纹段交汇

在裂尖靠近另一条裂纹时，在不符合裂尖与裂尖交汇的情况时，需要考虑裂尖与另一条裂纹的裂纹段交汇。如所示，需要根据当前裂尖与裂纹段之间的相对位置分情况处理：①过 ${O'_1}$ 作裂纹段的垂线，交于点P，若 $\angle {O_1}{O'_1}P$ 为钝角，则认为裂纹按 $\overrightarrow {{{O'}_1}P}$ 扩展到裂纹段，此时冻结该裂尖。②过 ${O'_1}$ 作裂纹段的垂线，交于点P，由于 $\angle {O_1}{O'_1}P$ 较小，认为真实裂纹扩展不会按这种情况交汇，因此作出修正，过 ${O'_1}$ 作 $\overrightarrow {{O_1}{{O'}_1}}$ 的垂线，与裂纹段交于点Q，认为裂尖按 $\overrightarrow {{{O'}_1}Q}$ 扩展，并与裂纹段交汇，此时冻结该裂尖。

图 4 裂尖与裂纹段交汇方式

Figure 4. Intersection strategy of crack tip of a crack with crack segment of another crack

下载: 全尺寸图片幻灯片

（3）裂尖与自由边界交汇

如所示，当裂尖扩展到自由边界附近时，需要判断裂尖与自由边界是否相交，以及交点位置。在中，O为当前扩展步的裂尖， $O{'_1}$ 为上一扩展步的裂尖，r为相互作用积分半径，一般取值为3倍的单元特征长度^[23-24]，也为裂纹扩展长度，OA和OB为裂尖到自由边界的垂线，P和Q分别为垂足，这里假设OP < OQ，即裂尖到自由边界的最短垂直距离是OP。

图 5 裂尖到达自由边界

Figure 5. Crack tip approaching free boundary

下载: 全尺寸图片幻灯片

如果相互作用积分半径r大于裂尖到自由边界的最短距离OP，则认为裂尖与自由边界交汇，交点为P，此时冻结该裂尖。

3.2 多裂纹扩展处理方法

在多裂纹扩展模拟中，引入加载乘子^[25]控制裂纹扩展始终处于准静态条件下，即所有裂尖的等效应力强度因子均小于等于材料的断裂韧度。如果采用 ${K_{{\text{eq}}}} > {K_{{\text{IC}}}}$ 作为断裂判据，则裂纹出现动态扩展，不符合准静态裂纹扩展条件^[8]。本文采用尝试–校正的方式检索竞争性裂纹，具体步骤如下：

给定一个初始的加载乘子 $\lambda = 1$ ，假设有n条裂纹，2n个裂尖。每个增量步加载一个 $\Delta p$ ，在每个加载步施加后，计算所有裂尖处的等效应力强度因子 ${K_{{\text{eq}}}}$ ，找到最大的等效应力强度因子 $K_{{\text{eq}}}^{\max }$ ，设定一个容差值 ${K_{{\text{tol}}}}$ 使得 ${K_{{\text{cal}}}} = K_{{\text{eq}}}^{{\text{max}}} - {K_{{\text{tol}}}}$ ，比较 ${K_{{\text{cal}}}}$ 与 ${K_{{\text{IC}}}}$ 的大小。若 ${K_{{\text{cal}}}} < {K_{{\text{IC}}}}$ ，则表明没有裂纹起裂，载荷更新 ${p_t} = {p_{t - 1}} + \Delta p$ ，继续计算下一个载荷。若 ${K_{{\text{cal}}}}$ $\geqslant$ ${K_{{\text{IC}}}}$ ，则表明至少有一个裂尖达到了扩展条件，为了表述的普遍性，假设有两个裂尖达到扩展要求，分别记为 ${T_{\text{a}}}$ 和 ${T_{\text{b}}}$ ，有 $K_{{\text{eq}}}^{{T_{\text{a}}}} \approx K_{{\text{eq}}}^{{T_{\text{b}}}} > {K_{{\text{IC}}}}$ ，这两个裂尖在当前载荷步下一定会扩展，将这两个裂尖扩展一段长度，新裂尖分别为 ${T_{{\text{a'}}}}$ 和 ${T_{{\text{b'}}}}$ 。在当前时间步下，更新裂纹构型后，计算 ${T_{{\text{a'}}}}$ 和 ${T_{{\text{b'}}}}$ 的等效应力强度因子 $K_{{\text{eq}}}^{{T_{{\text{a'}}}}}$ 和 $K_{{\text{eq}}}^{{T_{{\text{b'}}}}}$ ，计算加载乘子 $\lambda = {K_{{\text{IC}}}}/{\text{min}}\left\{ {K_{{\text{eq}}}^{{T_{{\text{a'}}}}}, K_{{\text{eq}}}^{{T_{{\text{b'}}}}}} \right\}$ ，将当前时间步载荷调整为 $\lambda {\text{*}}{p_t}$ ，更新裂纹构型后重新计算所有裂尖的等效应力强度因子，这使得载荷大小调整到使 ${T_{\text{a}}}$ 和 ${T_{\text{b}}}$ 刚好满足扩展要求，此时比较其它裂尖的等效应力强度因子与断裂韧度的关系。如果其它裂尖的 ${K_{{\text{eq}}}}$ 均小于 ${K_{{\text{IC}}}}$ ，则表明在 ${T_{\text{a}}}$ 和 ${T_{\text{b}}}$ 分别扩展到 ${T_{{\text{a'}}}}$ 和 ${T_{{\text{b'}}}}$ 的过程中，没有其它裂纹扩展，在调整后的载荷基础上再增加一个 $\Delta p$ ，进入下一个载荷步；如果有其它裂尖的 ${K_{{\text{eq}}}}$ 大于断裂韧度，假设是裂尖 ${T_{\text{c}}}$ ，则说明在 ${T_{\text{a}}}$ 和 ${T_{\text{b}}}$ 分别扩展到 ${T_{{\text{a'}}}}$ 和 ${T_{{\text{b'}}}}$ 的过程中， ${T_{\text{c}}}$ 也发生了扩展，即 ${T_{\text{a}}}$ ， ${T_{\text{b}}}$ 和 ${T_{\text{c}}}$ 的扩展属于竞争裂纹扩展，同样，将新扩展后的裂尖 ${T_{\text{c}}}$ 记为 ${T_{{\text{c'}}}}$ ，恢复到未调整时的载荷，更新裂纹构型后，重新计算 ${T_{{\text{a'}}}}$ ， ${T_{{\text{b'}}}}$ 和 ${{\text{T}}_{{\text{c'}}}}$ 的等效应力强度因子，计算加载乘子 $\lambda = {K_{{\text{IC}}}}/\max \left\{ {K_{{\text{eq}}}^{{T_{{\text{c'}}}}}, {\text{min}}\left\{ {K_{{\text{eq}}}^{{T_{{\text{a'}}}}}, K_{{\text{eq}}}^{{T_{{\text{b'}}}}}} \right\}} \right\}$ ，调整载荷为 $\lambda {\text{*}}{p_t}$ ，重新计算其它裂尖的等效应力强度因子，与断裂韧度比较，重复上面的步骤，直到找出所有扩展的裂尖，并更新加载乘子使得所有裂尖的等效应力强度均小于等于断裂韧度，然后进入下一个加载步。

4. 数值算例

场富集方法提供了处理多裂纹、分叉裂纹和交叉裂纹的数值手段，下面选取了一些典型的算例验证该方法的能力，同时通过数值模拟，分析这些算例中的一些裂纹演化特征以及断裂机理。下列算例均假定为弹脆性、各向同性的均质材料。

（1）含3条平行裂纹的3点弯曲梁

在这个算例中，为了检验裂尖应力强度因子和裂纹扩展路径，选取两种裂纹水平间距分析裂纹间的相互作用。如图 6所示，3点弯曲梁的底边包含3根平行的边裂纹，梁的长度L=100 mm，宽W=10 mm，预制裂纹长度均为a₁=a₂=a=2 mm，两种裂纹之间的水平间距：B=0.2W和B=2W。材料的力学参数：弹性模量E=20 GPa，泊松比 $\nu =$ 0.2，断裂韧度 ${K_{{\text{IC}}}} =$ ${\text{47}}{\text{.43 MPa}} \cdot {\text{m}}{{\text{m}}^{1/2}}$ ^[26]。梁的上边界中点施加集中力，下边界由两个间距为8W的支撑点支撑，左端支撑点完全约束，而右边支撑点仅约束y方向位移。将整个数值模型离散为10000个4节点四边形单元，在平面应力条件下进行计算。

图 6 含3条平行裂纹的3点弯曲梁几何尺寸与加载条件

Figure 6. Geometry and loading condition of three-point bending beam containing three flaws

下载: 全尺寸图片幻灯片

绘制了无量纲应力强度因子 ${K^{\text{*}}}$ 随裂纹间距比值B/W的变化，其中无量纲应力强度因子 ${K^{\text{*}}} = K/$ $\sigma \sqrt a$ ， $\sigma = 24P/H$ ，P为施加的集中载荷。获得的结果与Civelek等^[27]的解析解结果一致，这表明提出的数值方法能够获得满意的数值解。从图中还可以看出，中间裂纹的 $K_{{\text{II}}}^*$ 不随裂纹间距的改变而发生变化，始终为零，表明中心裂纹始终沿直线扩展。从中也可以注意到，在较小的B/W值下，中间裂纹的 $K_{\text{I}}^*$ 小于两侧裂纹的 $K_{\text{I}}^*$ ；随着B/W的增大，中间裂纹的 $K_{\text{I}}^*$ 开始大于两侧裂纹的 $K_{\text{I}}^*$ 。

图 7 无量纲应力强度因子

${K^{\text{*}}}$ 的数值解和理论解^[27]对比

Figure 7. Comparison of analytical and numerical normalized stress intensity factors at different B/W

下载: 全尺寸图片幻灯片

图 8展示了B=0.2W条件下3点弯曲梁在裂纹扩展过程中的应力分布情况及损伤分布情况。在3条裂纹扩展过程中，两侧裂纹逐渐远离中间裂纹，中间裂纹沿直线扩展，说明在较小的裂纹间距情况下，裂纹间的相互作用表现为排斥。

图 8 B=0.2W条件下的Von Mises云图和损伤图

Figure 8. Von Mises contour and damage field contour at B=0.2W

下载: 全尺寸图片幻灯片

图 9展示了B=2W条件下，3点弯曲梁在裂纹扩展过程中的应力分布情况，以及损伤分布情况。中间裂纹依然沿直线扩展，裂尖处有较高的应力集中现象，而两侧裂纹逐渐偏向中间裂纹，且裂尖处的应力集中区域较小。这说明中间裂纹屏蔽了周边裂纹的应力，使得它们的应力强度因子变小。裂纹间的相互作用表现为吸引。

图 9 B=2W条件下的Von Mises云图和损伤图

Figure 9. Von Mises contour and damage field contour at B=2W

下载: 全尺寸图片幻灯片

两种裂纹间距（B=0.2W和B=2W）情况下的载荷–位移曲线如图 10所示，同时，图 8，9对应的加载步也在曲线上高亮标记出。

图 10 三点弯曲梁的载荷–位移曲线图

Figure 10. Load-displacement curves of three-point bending beam containing three flaws

下载: 全尺寸图片幻灯片

（2）周期性分布多裂纹问题

如图 11所示，含周期性分布多裂纹的矩形板受拉伸载荷作用，矩形板长L=50 mm，宽W=30 mm，裂纹长度a=4 mm，裂纹水平间距s=4 mm，垂直间距2d=4 mm，r=7 mm，c=9 mm。物理力学参数与上一算例相同。为表述方便，对裂纹和裂尖进行编号，裂纹编号为“Cn”，n为裂纹数，裂尖编号为“Tm”，m取值为1或2。

图 11 周期性裂纹几何和加载条件

Figure 11. Periodic crack geometry and loading conditions

下载: 全尺寸图片幻灯片

裂纹总共扩展了41步，裂纹的破坏形态如图 12所示，为了清晰地说明裂纹扩展过程和交汇情况，在载荷–位移曲线（图 13）中标注出裂纹的扩展和交汇顺序，并在表 1中予以说明。

图 12 数值模型最终破坏图

Figure 12. Final failure of numerical model containing periodic cracks

下载: 全尺寸图片幻灯片

图 13 载荷–位移曲线

Figure 13. Load-displacement curves of plate containing periodic cracks

下载: 全尺寸图片幻灯片

表 1 裂纹扩展和交汇顺序

Table 1. Sequence of periodic crack propagation and coalescence

标记点	裂纹交汇顺序
点a	裂纹开始起裂
点b	C2的裂尖与C4的T2和C5的T1相互吸引，C17的裂尖与C14的T2和C15的T1相互吸引
点c	C1的T1，C3的T2，C16的T1，C18的T2同时扩展到模型边界
点d	裂纹贯穿整个模型试样

下载: 导出CSV

| 显示表格

（3）随机分布的多裂纹问题

如图 14所示，含多条随机分布的裂纹的矩形板受拉伸载荷作用，板的尺寸和加载条件同上面算例。物理力学参数也与上面相同。该算例用来证明提出的数值方法处理随机分布的多裂纹的能力。

图 14 随机分布裂纹的几何和加载条件

Figure 14. Random crack geometry and loading conditions

下载: 全尺寸图片幻灯片

如图 15所示，裂纹总共扩展了89步，在裂纹起裂到完全贯穿整个试样的过程中，经历了多次交汇，通过表 2简单阐述裂纹之间的交汇顺序。同样，这些交汇发生的时间步也在载荷–位移曲线（图 16）上标明。

图 15 数值模型最终破坏图

Figure 15. Final failure of numerical model containing random cracks

下载: 全尺寸图片幻灯片

表 2 裂纹扩展和交汇顺序

Table 2. Sequence of random crack propagation and coalescence

标记点	裂纹交汇顺序
点a	裂纹开始起裂
点b	C11的T2扩展到数值模型边界
点c	C11的T1和C10的T2交汇
点d	C10的T1与C14的T1交汇
点e	C9的T1与C8交汇
点f	C8的T1与C7交汇
点g	C7的T1扩展到数值模型边界
点h	C9的T2与C10的T1交汇

下载: 导出CSV

| 显示表格

图 16 载荷–位移曲线

Figure 16. Load-displacement curves of plate containing random cracks

下载: 全尺寸图片幻灯片

（4）随机分布的复杂裂纹问题

如图 17所示，含随机分布的分叉裂纹、十字交叉裂纹和星形裂纹的方板受双轴拉伸作用。方板的边长为10 mm，裂纹的每个分支长度a=0.5 mm，分支角度是45°的整数倍。该算例通过随机分布的复杂裂纹构型验证提出的数值方法处理多裂纹扩展的能力。

图 17 裂纹分布和几何加载条件

Figure 17. Distribution of multiple cracks and geometry and loading conditions

下载: 全尺寸图片幻灯片

图 18中展示了裂纹的扩展演化以及交汇形式，同时该模型的载荷–位移曲线也绘制在图 19中。

图 18 复杂裂纹最终破坏图

Figure 18. Final failure of complex multiple cracks

下载: 全尺寸图片幻灯片

图 19 含复杂多裂纹板的载荷–位移曲线

Figure 19. Load-displacement curves of plate containing complex cracks

下载: 全尺寸图片幻灯片

5. 结论

（1）本文提出的场富集有限元方法（FE-FEM）能够有效地模拟多裂纹扩展和连接问题，包括分叉裂纹、交叉裂纹等各种复杂裂纹的组合。

（2）针对多裂纹在扩展过程中出现的各种交汇情况提出了解决方案，分别处理了裂尖与另一条裂纹的裂尖交汇，裂尖与另一条裂纹的裂纹段交汇，以及裂尖到达数值模型的自由边界处三种情形，并通过算例证明了处理策略的合理性。

（3）针对在每个加载步下多裂纹可能出现竞争性扩展，即在满足准静态扩展的条件下，在一个加载步下出现多个裂尖扩展，通过引入加载乘子，调节载荷，通过不断地调整验算，最终使得每个载荷步下裂纹尖端的等效应力均小于等于材料的断裂韧度。

（4）场富集有限元方法可以直接处理交叉裂纹，而不需要像扩展有限元一样引入额外的富集函数。通过数值算例，充分地说明了场富集有限元方法在处理各种交叉裂纹组成的复杂裂纹系统扩展演化方面的能力。

图 1 相互作用积分的域形式

Figure 1. Domain form of interaction integral

下载: 全尺寸图片幻灯片

图 2 富集节点检索方式

Figure 2. Searching strategy of enriched nodes

下载: 全尺寸图片幻灯片

图 3 裂尖与裂尖交汇方式

Figure 3. Intersection strategy of two crack tips

下载: 全尺寸图片幻灯片

图 4 裂尖与裂纹段交汇方式

Figure 4. Intersection strategy of crack tip of a crack with crack segment of another crack

下载: 全尺寸图片幻灯片

图 5 裂尖到达自由边界

Figure 5. Crack tip approaching free boundary

下载: 全尺寸图片幻灯片

图 6 含3条平行裂纹的3点弯曲梁几何尺寸与加载条件

Figure 6. Geometry and loading condition of three-point bending beam containing three flaws

下载: 全尺寸图片幻灯片

无量纲应力强度因子 ${K^{\text{*}}}$ 的数值解和理论解^[27]对比

Figure 7. Comparison of analytical and numerical normalized stress intensity factors at different B/W

下载: 全尺寸图片幻灯片

图 8 B=0.2W条件下的Von Mises云图和损伤图

Figure 8. Von Mises contour and damage field contour at B=0.2W

下载: 全尺寸图片幻灯片

图 9 B=2W条件下的Von Mises云图和损伤图

Figure 9. Von Mises contour and damage field contour at B=2W

下载: 全尺寸图片幻灯片

图 10 三点弯曲梁的载荷–位移曲线图

Figure 10. Load-displacement curves of three-point bending beam containing three flaws

下载: 全尺寸图片幻灯片

图 11 周期性裂纹几何和加载条件

Figure 11. Periodic crack geometry and loading conditions

下载: 全尺寸图片幻灯片

图 12 数值模型最终破坏图

Figure 12. Final failure of numerical model containing periodic cracks

下载: 全尺寸图片幻灯片

图 13 载荷–位移曲线

Figure 13. Load-displacement curves of plate containing periodic cracks

下载: 全尺寸图片幻灯片

图 14 随机分布裂纹的几何和加载条件

Figure 14. Random crack geometry and loading conditions

下载: 全尺寸图片幻灯片

图 15 数值模型最终破坏图

Figure 15. Final failure of numerical model containing random cracks

下载: 全尺寸图片幻灯片

图 16 载荷–位移曲线

Figure 16. Load-displacement curves of plate containing random cracks

下载: 全尺寸图片幻灯片

图 17 裂纹分布和几何加载条件

Figure 17. Distribution of multiple cracks and geometry and loading conditions

下载: 全尺寸图片幻灯片

图 18 复杂裂纹最终破坏图

Figure 18. Final failure of complex multiple cracks

下载: 全尺寸图片幻灯片

图 19 含复杂多裂纹板的载荷–位移曲线

Figure 19. Load-displacement curves of plate containing complex cracks

下载: 全尺寸图片幻灯片

表 1 裂纹扩展和交汇顺序

Table 1 Sequence of periodic crack propagation and coalescence

标记点	裂纹交汇顺序
点a	裂纹开始起裂
点b	C2的裂尖与C4的T2和C5的T1相互吸引，C17的裂尖与C14的T2和C15的T1相互吸引
点c	C1的T1，C3的T2，C16的T1，C18的T2同时扩展到模型边界
点d	裂纹贯穿整个模型试样

下载: 导出CSV

表 2 裂纹扩展和交汇顺序

Table 2 Sequence of random crack propagation and coalescence

标记点	裂纹交汇顺序
点a	裂纹开始起裂
点b	C11的T2扩展到数值模型边界
点c	C11的T1和C10的T2交汇
点d	C10的T1与C14的T1交汇
点e	C9的T1与C8交汇
点f	C8的T1与C7交汇
点g	C7的T1扩展到数值模型边界
点h	C9的T2与C10的T1交汇

下载: 导出CSV

参考文献(27)

[1]	FREIJ-AYOUB R, DYSKIN A V, GALYBIN A N. The dislocation approximation for calculating crack interaction[J]. International Journal of Fracture, 1997, 86(4): 57–62.
[2]	RYBACZUK M, STOPPEL P. The fractal growth of fatigue defects in materials[J]. International Journal of Fracture, 2000, 103(1): 71–94. doi: 10.1023/A:1007635717332
[3]	CHEN Y Z. General case of multiple crack problems in an infinite plate[J]. Engineering Fracture Mechanics, 1984, 20(4): 591–597. doi: 10.1016/0013-7944(84)90034-1
[4]	CHENG H, ZHOU X P, ZHU J, et al. The effects of crack openings on crack initiation, propagation and coalescence behavior in rock-like materials under uniaxial compression[J]. Rock Mechanics and Rock Engineering, 2016, 49(9): 3481–3494. doi: 10.1007/s00603-016-0998-9
[5]	ZHANG J Z, ZHOU X P. AE event rate characteristics of flawed granite: from damage stress to ultimate failure[J]. Geophysical Journal International, 2020, 222(2): 795–814. doi: 10.1093/gji/ggaa207
[6]	CARPINTERI A, MONETTO I. Snap-back analysis of fracture evolution in multi-cracked solids using boundary element method[J]. International Journal of Fracture, 1999, 98(3/4): 225–241. doi: 10.1023/A:1018660600546
[7]	DENDA M, DONG Y F. Complex variable approach to the BEM for multiple crack problems[J]. Computer Methods in Applied Mechanics and Engineering, 1997, 141(3/4): 247–264.
[8]	BUDYN E, ZI G, MOËS N, et al. A method for multiple crack growth in brittle materials without remeshing[J]. International Journal for Numerical Methods in Engineering, 2004, 61(10): 1741–1770. doi: 10.1002/nme.1130
[9]	ZHOU X P, CHEN J W. Extended finite element simulation of step-path brittle failure in rock slopes with non-persistent en-echelon joints[J]. Engineering Geology, 2019, 250: 65–88. doi: 10.1016/j.enggeo.2019.01.012
[10]	WANG Y T, ZHOU X P, WANG Y, et al. A 3-D conjugated bond-pair-based peridynamic formulation for initiation and propagation of cracks in brittle solids[J]. International Journal of Solids and Structures, 2018, 134: 89–115. doi: 10.1016/j.ijsolstr.2017.10.022
[11]	AZADI H, KHOEI A R. Numerical simulation of multiple crack growth in brittle materials with adaptive remeshing[J]. International Journal for Numerical Methods in Engineering, 2011, 85(8): 1017–1048. doi: 10.1002/nme.3002
[12]	ZHOU X P, FU L, QIAN Q H. A 2D novel non-local lattice bond model for initiation and propagation of cracks in rock materials[J]. Engineering Analysis with Boundary Elements, 2021, 126: 181–199. doi: 10.1016/j.enganabound.2021.03.002
[13]	ZHOU X P, BI J, QIAN Q H. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws[J]. Rock Mechanics and Rock Engineering, 2015, 48(3): 1097–1114. doi: 10.1007/s00603-014-0627-4
[14]	JIA Z M, ZHOU X P, BERTO F. Compressive-shear fracture model of the phase-field method coupled with a modified Hoek–Brown criterion[J]. International Journal of Fracture, 2021, 229(2): 161–184. doi: 10.1007/s10704-021-00546-7
[15]	XU D D, WU A Q, LI C. A linearly-independent higher-order extended numerical manifold method and its application to multiple crack growth simulation[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2019, 11(6): 1256–1263. doi: 10.1016/j.jrmge.2019.02.007
[16]	石路杨, 余天堂. 多裂纹扩展的扩展有限元法分析[J]. 岩土力学, 2014, 35(1): 263–272. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201401040.htm SHI Lu-yang, YU Tian-tang. Analysis of multiple crack growth using extended finite element method[J]. Rock and Soil Mechanics, 2014, 35(1): 263–272. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201401040.htm
[17]	BARBIERI E, PETRINIC N, MEO M, et al. A new weight-function enrichment in meshless methods for multiple cracks in linear elasticity[J]. International Journal for Numerical Methods in Engineering, 2012, 90(2): 177–195. doi: 10.1002/nme.3313
[18]	RABCZUK T, BORDAS S, ZI G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics[J]. Computational Mechanics, 2007, 40(3): 473–495. doi: 10.1007/s00466-006-0122-1
[19]	YAU J F, WANG S S, CORTEN H T. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[J]. Journal of Applied Mechanics, 1980, 47(2): 335–341. doi: 10.1115/1.3153665
[20]	SHIH C F, ASARO R J. Elastic-plastic analysis of cracks on bimaterial interfaces: part Ⅰ—small scale yielding[J]. Journal of Applied Mechanics, 1988, 55(2): 299–316. doi: 10.1115/1.3173676
[21]	ERDOGAN F, SIH G C. Closure to "discussion of 'on the crack extension in plates under plane loading and transverse shear'"[J]. Journal of Basic Engineering, 1963, 85(4): 527. doi: 10.1115/1.3656899
[22]	ZHOU X P, JIA Z M, WANG L F. A field-enriched finite element method for brittle fracture in rocks subjected to mixed mode loading[J]. Engineering Analysis with Boundary Elements, 2021, 129: 105–124. doi: 10.1016/j.enganabound.2021.04.023
[23]	SUKUMAR N, PRÉVOST J H. Modeling quasi-static crack growth with the extended finite element method Part Ⅰ: computer implementation[J]. International Journal of Solids and Structures, 2003, 40(26): 7513–7537. doi: 10.1016/j.ijsolstr.2003.08.002
[24]	MOËS N, DOLBOW J, BELYTSCHKO T. A finite element method for crack growth without remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150. doi: 10.1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J
[25]	刘丰, 郑宏, 李春光. 基于NMM的EFG方法及其裂纹扩展模拟[J]. 力学学报, 2014, 46(4): 582–590. https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201404012.htm LIU Feng, ZHENG Hong, LI Chun-guang. The nmm-based efg method and simulation of crack propagation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(4): 582–590. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201404012.htm
[26]	PALUSZNY A, MATTHÄI S K. Numerical modeling of discrete multi-crack growth applied to pattern formation in geological brittle media[J]. International Journal of Solids and Structures, 2009, 46(18/19): 3383–3397.
[27]	CIVELEK M B, ERDOGAN F. Crack problems for a rectangular plate and an infinite strip[J]. International Journal of Fracture, 1982, 19(2): 139–159. doi: 10.1007/BF00016570

施引文献(30)

期刊类型引用(10)

1.	刘毅，王海啸，王斌，周傲，位伟. 基于逐工况反演的基坑围护结构变形分析与探讨. 武汉大学学报(工学版). 2024(09): 1221-1231 . 百度学术
2.	薛秀丽，廖欢，曾超峰，刘运思，曾兴. 既有地下结构水-土阻隔效应对基坑抽水引发地层变形影响机制. 岩土工程学报. 2023(01): 103-111 . 本站查看
3.	朱宾. TRD技术在民用建筑基坑施工中的应用. 四川水泥. 2023(02): 169-171 . 百度学术
4.	陈心春. 超大深基坑TRD连续墙厚度对渗透系数的影响研究. 四川水泥. 2023(03): 42-44+59 . 百度学术
5.	陈建军，刘延昭，吴高创，李俊，徐超，王金昌. 杭州某国道明挖下穿对既有市政高架桥桥墩影响分析. 地基处理. 2023(06): 527-533 . 百度学术
6.	王铭，郭豪，胡光辉，张壮，杜金亮，刘塬熙. 基于尖点突变分析和变形预测的基坑沉降稳定性. 黑龙江科技大学学报. 2023(06): 889-895 . 百度学术
7.	仉文岗，胡鑫芸，章润红，陈春霞，黎泳钦，叶文煜，张智超，陈荣淋. 各向异性黏土中的深基坑开挖单撑失效数值模拟研究（英文）. Journal of Central South University. 2023(12): 4168-4181 . 百度学术
8.	董海君，马榕键，陈骞，张士龙. 深厚淤泥质土层深基坑变形特性研究. 铁道建筑技术. 2022(10): 87-91 . 百度学术
9.	夏曾银，潘军，盛鲁腾，程雪松，郑刚，冀叶涛，孟灵波. 注浆和隔离墙对基坑引发隧道变形的联合控制作用研究. 水利水电技术(中英文). 2022(09): 175-185 . 百度学术
10.	陈心春. 深基坑TRD连续墙稳定性研究分析. 江西建材. 2022(09): 108-110 . 百度学术