# 3D numerical manifold method for crack propagation in rock materials using a local tracking algorithm

Boyi Su, Tao Xu, Genhua Shi, Michael Heap, Xianyang Yu, Guanglei Zhou

► To cite this version:

Boyi Su, Tao Xu, Genhua Shi, Michael Heap, Xianyang Yu, et al.. 3D numerical manifold method for crack propagation in rock materials using a local tracking algorithm. Journal of Rock Mechanics and Geotechnical Engineering, 2024, 129 (8), 10.1016/j.jrmge.2024.04.038. hal-04744028

HAL Id: hal-04744028

https://hal.science/hal-04744028v1

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Journal Pre-proof

3D numerical manifold method for crack propagation in rock materials using a local
tracking algorithm

Boyi Su, Tao Xu, Genhua Shi, Michael J. Heap, Xianyang Yu, Guanglei Zhou

PII: S1674-7755(24)00460-8

DOI: https://doi.org/10.1016/j.jrmge.2024.04.038

Reference: JRMGE 1763

To appear in: *Journal of Rock Mechanics and Geotechnical Engineering*

Received Date: 3 November 2023

Revised Date: 3 February 2024

Accepted Date: 14 April 2024

<figure><img src="image_2.png" /></figure>

Please cite this article as: Su B, Xu T, Shi G, Heap MJ, Yu X, Zhou G, 3D numerical manifold method for crack propagation in rock materials using a local tracking algorithm, *Journal of Rock Mechanics and Geotechnical Engineering*, https://doi.org/10.1016/j.jrmge.2024.04.038.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition
of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of
record. This version will undergo additional copyediting, typesetting and review before it is published
in its final form, but we are providing this version to give early visibility of the article. Please note that,
during the production process, errors may be discovered which could affect the content, and all legal
disclaimers that apply to the journal pertain.

© 2024 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Published by Elsevier B.V.
3D numerical manifold method for crack propagation in rock materials using a local tracking algorithm

Boyi Su<sup>a</sup>, Tao Xu<sup>a,*</sup>, Genhua Shi<sup>b</sup>, Michael J. Heap<sup>c,d</sup>, Xianyang Yu<sup>c</sup>, Guanglei Zhou<sup>f</sup>

<sup>a</sup>Center for Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang, 110819, China

<sup>b</sup>University of Chinese Academy of Sciences, Beijing, 100049, China

<sup>c</sup>Université de Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, 5 rue René Descartes, Strasbourg, F-67084, France

<sup>d</sup>Institut Universitaire de France (IUF), Paris, France

<sup>e</sup>State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China

<sup>f</sup>Key Laboratory of Mining Disaster Prevention and Control, School of Energy and Mining Engineering, Shandong University of Science and Technology, Qingdao, 266590, China

*Corresponding author. E-mail address: xutao@mail.neu.edu.cn

**Abstract:** The modeling of crack growth in three-dimensional (3D) space poses significant challenges in rock mechanics due to the complex numerical computation involved in simulating crack propagation and interaction in rock materials. In this study, we present a novel approach that introduces a 3D numerical manifold method (3D-NMM) with a geometric kernel to enhance computational efficiency. Specifically, the maximum tensile stress criterion is adopted as a crack growth criterion to achieve strong discontinuous crack growth, and a local crack tracking algorithm and an angle correction technique are incorporated to address minor limitations of the algorithm in a 3D model. The implementation of the program is carried out in Python, using object-oriented programming in two independent modules: a calculation module and a crack module. Furthermore, we propose feasible improvements to enhance the performance of the algorithm. Finally, we demonstrate the feasibility and effectiveness of the enhanced algorithm in the 3D-NMM using four numerical examples. This study establishes the potential of the 3D-NMM, combined with the local tracking algorithm, for accurately modeling 3D crack propagation in brittle rock materials.

**Keywords:** 3D numerical manifold method (3D NMM), crack propagation, local tracking algorithm, brittle materials

# 1. Introduction

Discontinuities such as joints and cracks are abundant in natural rocks (Hudson and Priest, 1983), and the failure process of rock masses is commonly associated with their propagation, interaction, and coalescence (Hoek et al., 2002). An acceleration in the rate of crack propagation is often an indicator of impending macroscopic failure (Li et al., 2024). In rock engineering, crack propagation is a crucial factor that determines the stability of rock masses and geoengineering structures and, as a result, it is a fundamental topic in rock mechanics.

In order to simulate discontinuities in rock masses, Shi (1992a, b) developed the numerical manifold method (NMM). In comparison to the finite element method (FEM), NMM uses two cover systems, providing it with the capability to effectively address problems involving the continuous-discontinuous behavior of rock. The generation of a mathematical cover in the NMM is not limited to triangular or tetrahedral FEM meshes, and other types of mesh can be used that have specific problem-solving advantages. Terada and co-workers (Terada et al., 2003; Terada and Kurumatani, 2004) developed the NMM by using a square quadrilateral mathematical cover, compared the accuracy of the NMM with that of the FEM, and found that the NMM achieves a comparable or slightly superior accuracy to the FEM. The moving least squares numerical manifold method (MLS-NMM) combines a moving least squares interpolation with the NMM (Zheng et al., 2014; Li et al., 2022; Zhang et al., 2022; Chen et al., 2023). The MLS-NMM can handle problems with complex geometries and the large displacements caused by cracking. The MLS-NMM was used to analyze various physical phenomena, such as crack propagation, porous media seepage, and nonlinear heat transfer (Zheng et al., 2014; Li et al., 2022; Zhang et al., 2022; Chen et al., 2023). MLS-NMM is an effective numerical method for rock engineering. Compared with MLS-NMM, the developed 3D NMM has an innovation on crack propagation in three dimensions. We applied a geometric kernel in 3D NMM so that we can easily access geometric calculations. Common geometric calculations, such as the crack surface generation algorithm and the block cutting algorithm can be carried out by the geometric kernel.

The particle-based NMM (PNMM) is a numerical method that combines the discrete element method (DEM) for particle flow with the NMM (Li et al., 2018; Li and Zhao, 2019). The PNMM introduces the concept of particles into the NMM, providing advantages in terms of simplex geometry
boundaries and contact operations with particles. The PNMM enables a more accurate simulation of dynamic problems, such as rock impacts and explosions. The NMM also possesses advantages when dealing with discontinuities caused by crack propagation. For example, the NMM only requires that the physical cover is cut and therefore does not need to generate new nodes along the crack boundaries, eliminating the need for remeshing during the computation process. An et al. (2014) and Ning et al. (2023) used the NMM to simulate crack propagation problems such as pre-cracked Brazilian discs, brittle materials under dynamic loading, and the progressive failure of rock slopes. Furthermore, several NMM models were developed to model the time-dependent deformation of rock, such as the virtual crack and grain models (Yu et al., 2018, 2021, 2022, 2023; Zhou et al., 2022). Hu and Rutqvist (2020) developed a comprehensive modeling approach for coupled processes based on NMM. In this approach, the fractures were divided into three parts depending on their multiple scales. The model developed by Hu and Rutqvist (2020) is able to accurately represent fracture intersection and shear on discrete fracture scales and rigorously treat contact along discontinuous fracture surfaces.

Three-dimensional (3D) NMM involves many fields, such as crack propagation, contact theory, and slope stability. For example, He and Ma (2010) developed a 3D NMM that incorporated the process of block cutting using a description of the topological structure of the blocks. Tong et al. (2023) developed a 3D NMM to simulate steady-state and transient heat conduction problems. The results of Tong et al. (2023) showed that the 3D NMM can solve the heat conduction problem with a high accuracy. Liu et al. (2019) developed a 3D MLS-NMM, which greatly simplified the process of preprocessing by using MLS nodes to generate a mathematical cover. For crack propagation, a non-local tracking algorithm and the phase field method can be integrated with the NMM to calculate the progress of crack propagation (Yang et al., 2018; Yang and Chen, 2023; Yang et al., 2016). Kang et al. (2023) developed the contact theory for 3D NMM based on the cover contact theory and used it to better understand complex structural deformation and dynamic corrosion problems.

The process of rock failure can be explained by damage theory or fracture theory. In numerical calculations of damage theory, the mechanical parameters are degraded to simulate the failure of arbitrary elements. These elements are called damage elements. Compared with damage theory, fracture theory has the ability to determine the state of crack propagation. Annavarapu et al. (2016) classified crack tracking algorithms into two categories: the explicit and implicit approaches. The level set method is one of the implicit approaches for tracking the crack surface. Each crack is represented by two level set functions, one that tracks the crack surface and another that tracks the crack front, so that the intersection can locate the entire crack. Stolarska et al. (2001) proposed a hybrid approach that combines the level set method and the extended FEM (x-FEM) to study crack propagation. The x-FEM eliminates the need for remeshing as the crack propagates. Gravouil et al. (2002) and Moës et al. (2002) developed a 3D crack tracking algorithm using the level set method. For explicit approaches, a one-dimensional crack occurs in a two-dimensional (2D) numerical calculation, while a 2D crack occurs in a 3D numerical calculation. Since the fractures can be arbitrary with respect to the underlying mesh, it is necessary to develop tracking algorithms to track crack surfaces. Jäger et al. (2008) classified crack tracking algorithms applied to numerical models into four types: fixed, local, non-local, and global algorithms. The local tracking algorithm is a simple and easy-to-understand method derived from 2D crack tracking algorithms. Liu et al. (2014) applied a 2D local tracking algorithm to the x-FEM, generating highly discontinuous and scalable cracks within the solution domain, resulting in improved solution accuracy. Building upon this work, Naderi et al. (2016) extended the local tracking algorithm to 3D x-FEM, enabling the initiation and propagation of arbitrary cracks in 3D solids, and derived equilibrium equations to enforce a strong discontinuity between the elements on either side of the crack. The global tracking algorithm is a crack tracking algorithm that is not constrained by the mesh and typically relies on complete information of the solution domain to determine the location and direction of crack propagation. The global crack tracking algorithm was proposed based on the continuous strong discontinuity method, which uses a viscous perturbation method to ensure the uniqueness of the crack (Oliver and Huespe, 2004; Oliver et al., 2004). Jäger et al. (2009) applied a global tracking algorithm to 3D FEM, where the crack initiation and propagation are controlled by three sets of equations. Riccardi et al. (2017) integrated the global tracking algorithm into the FEM to simulate the 2D fracture processes of quasi-brittle materials, and proposed a novel crack tracking algorithm that effectively prevents the occurrence of discontinuities during crack propagation. Saloustros et al. (2019) analyzed the different types of crack analysis methods applicable to the FEM and summarized the failure criteria that take the crack propagation direction into account. They also discuss the potential problems and challenges of tracking algorithms.

Crack propagation along the edge of elements is also an effective crack tracking method. In the finite-discrete element method (FDEM), cohesive elements fill the boundary of the element, and so the crack propagates along the edge of the element. Yahaghi et al. (2023) developed grain-based model based on FDEM that considers the actual microstructure of rocks with grain-growth tessellations to model crack propagation. Zheng et al. (2023) developed a modified joint element constitutive model to simulate the nonlinear failure behavior of rocks. Compared to crack propagation along the edges of elements, the method of crack propagation passing through the elements does not restrict the propagation direction by the grid and can simulate crack propagation with more accuracy.

In this study, we implemented a crack local tracking algorithm using a 3D NMM. This algorithm has previously been used in 2D and 3D FEMs with
numerous numerical examples that demonstrate its efficiency (Liu et al., 2014; Naderi et al., 2016). This paper is structured as follows. Section 2 discusses the basic concepts of the 3D NMM and the local tracking algorithm. Section 3 describes the geometry kernel used in the 3D NMM and explains the procedure for crack propagation in the numerical calculations. In Section 4, we present and discuss four numerical examples to verify the validity of the 3D-NMM and the local tracking algorithm. Section 5 discusses the limitations of and provides some perspectives for the local tracking algorithm. Finally, in Section 6, we offer some concluding remarks and present an outlook for future work.

## 2. 3D NMM based on local tracking algorithm

### 2.1 Basic principle of 3D NMM

The NMM is a numerical method that combines continuous and discontinuous media to solve differential equations. NMM has developed two cover systems, viz. the mathematical cover (MC) and the physical cover (PC). These two cover systems enable NMM to accurately calculate problems that involve discontinuities.

The MC is an essential component in NMM, determining the precision of the solution. It must fully encapsulate the complete solution domain, and the size and shape of the MCs are not strictly restricted. Common ways to generate the MC include regularly arranged squares or cubes and an FEM mesh.

The PC serves as the boundary of the solution domain and defines the integral domain. Each MC cuts the solution domain, resulting in one or more corresponding PCs. In some cases, one MC may generate multiple PCs. For example, if a crack cuts one MC into two parts, then the MC will generate two PCs, as shown in Fig. 1. $M_1$ and $M_2$ correspond to the MCs, the trapezoidal with a crack corresponds to the solution domain, the MCs and the discontinuous crack cut the solution domain to generate physical covers (PCs). $P_{1-1}$, $P_{2-1}$ and $P_{2-2}$ correspond to the PCs. The PCs cut each other to generate manifold elements, and $E_{1-1}$, $E_{1-1,2-1}$, $E_{1-1,2-2}$ and $E_{2-2}$ correspond to the manifold elements.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Mathematical cover (MC), physical cover (PC), and manifold element.</figcaption></figure>

The displacement of any point within the solution domain is determined by the cover function and the weight function. Cover functions are displacement functions that are defined on each PC. Although PCs on either side of the same crack may be generated by the same MC, the different cover functions can lead to discontinuities between the two sides of the crack. The cover function used in 3D NMM is shown below:

$$ \left\{ \begin{array}{l} u_{e(i)x}(x, y, z) \quad (x, y, z) \in U_i \\ u_{e(i)y}(x, y, z) \quad (x, y, z) \in U_i \\ u_{e(i)z}(x, y, z) \quad (x, y, z) \in U_i \end{array} \right\} \qquad (1) $$

where $u_{e(i)x}$, $u_{e(i)y}$ and $u_{e(i)z}$ are the cover functions; and $U_i$ is the region of PC. Cover functions are usually continuous and smooth in their domain. Hence, the constant function ($u_{e(i)j}(x, y, z) = d_0, j = x, y, z$) and the first-order polynomial ($u_{e(i)j}(x, y, z) = d_0 + d_1x + d_2y + d_3z, j = x, y, z$) can be selected as a cover function. In 3D NMM, the cover function can be expressed in terms of the constant function:

$$ \left\{ \begin{array}{l} u_{e(i)x}(x, y, z) \\ u_{e(i)y}(x, y, z) \\ u_{e(i)z}(x, y, z) \end{array} \right\} = \left\{ \begin{array}{l} 1 \quad 0 \quad 0 \\ 0 \quad 1 \quad 0 \\ 0 \quad 0 \quad 1 \end{array} \right\} \left\{ \begin{array}{l} d_{e(i)x} \\ d_{e(i)y} \\ d_{e(i)z} \end{array} \right\} \qquad (2) $$

where $d_{e(i)x}$, $d_{e(i)y}$ and $d_{e(i)z}$ are the coefficients of the cover functions. A weight function is used as a weighted average to associate the cover functions with each other and to generate a global displacement function. The global displacement function can be represented as
$$ \begin{bmatrix} u_x(x, y, z) \\ u_y(x, y, z) \\ u_z(x, y, z) \end{bmatrix} = \sum_{i}^{n} \omega_{e(i)}(x, y, z) \begin{bmatrix} u_{e(i)x}(x, y, z) \\ u_{e(i)y}(x, y, z) \\ u_{e(i)z}(x, y, z) \end{bmatrix} \quad (3) $$

where $\omega_{e(i)}$ is the weight function, and $n$ is the number of PCs shared by the same manifold element. At a given point, the weight function of each PC must be greater than or equal to 0. If the PC contains the point, the weight function must be greater than 0, and if the PC does not contain the point, the weight function must be equal to 0. The sum of all weight functions of PC at a certain point must be equal to 1, as shown below:

$$ \left. \begin{array}{l} 0 \le \omega_{e(i)}(x, y, z) \le 1 \quad (x, y, z) \in U_i \\ \omega_{e(i)}(x, y, z) = 0 \quad (x, y, z) \notin U_i \\ \sum_{i}^{n} \omega_{e(i)}(x, y, z) = 1 \end{array} \right\} \quad (4) $$

In a given manifold element, the global displacement is influenced by the adjacent PCs. In recent research, tetrahedra are used as manifold elements. The MC is generated by a 3D FEM mesh, where the boundary of the solution domain is aligned with the FEM mesh. Each element has four vertices that correspond to four PCs, making $n = 4$. PC is a region that is formed by the aggregation of surrounding manifold elements, as shown in Fig. 2. In 3D NMM, we use an FEM mesh as PC. An element is the intersection zone of four different PCs, and each PC consists of all the elements that share a specific vertex.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Mesh, physical covers (PC) and manifold elements in the 3D NMM.</figcaption></figure>

The coordinate of vertices $e(i)$ is $(x_{e(i)}, y_{e(i)}, z_{e(i)})$. The expression of weight function can be written as

$$ \begin{pmatrix} \omega_{e(1)}(x, y, z) \\ \omega_{e(2)}(x, y, z) \\ \omega_{e(3)}(x, y, z) \\ \omega_{e(4)}(x, y, z) \end{pmatrix} = \begin{pmatrix} f_{11} & f_{12} & f_{13} & f_{14} \\ f_{21} & f_{22} & f_{23} & f_{24} \\ f_{31} & f_{32} & f_{33} & f_{34} \\ f_{41} & f_{42} & f_{43} & f_{44} \end{pmatrix} \begin{pmatrix} 1 \\ x \\ y \\ z \end{pmatrix} \quad (5) $$

where

$$ \begin{pmatrix} f_{11} & f_{12} & f_{13} & f_{14} \\ f_{21} & f_{22} & f_{23} & f_{24} \\ f_{31} & f_{32} & f_{33} & f_{34} \\ f_{41} & f_{42} & f_{43} & f_{44} \end{pmatrix}^T = \begin{pmatrix} 1 & x_{e(1)} & y_{e(1)} & z_{e(1)} \\ 1 & x_{e(2)} & y_{e(2)} & z_{e(2)} \\ 1 & x_{e(3)} & y_{e(3)} & z_{e(3)} \\ 1 & x_{e(4)} & y_{e(4)} & z_{e(4)} \end{pmatrix}^{-1} \quad (6) $$

## 2.2 Crack growth criterion based on maximum tensile stress

The selection of an appropriate crack propagation criterion is crucial for accurately calculating crack propagation. In this study, the maximum tensile stress criterion is chosen as the crack growth criterion to simulate the tensile failure of rock, as specified below:

$$ \sigma_{\max} = \sigma_{T} \quad (7) $$

where $\sigma_{\max}$ is the maximum principal stress of the potential crack element, and $\sigma_T$ is the tensile strength of the element. Each model starts with an initial global crack surface, and subsequent crack surfaces propagate along this initial surface. An element is classified as a potential crack element when its boundary surface contains one or more crack tips. If the maximum principal stress of a potential crack element reaches the tensile strength, the global crack surface will propagate, and a new crack surface will form passing through the element. To calculate strain based on the global displacement function,
Eq. (3) can be expressed in matrix form as follows:

$$ \begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix} = \sum_{i=1}^{4} \begin{pmatrix} \omega_{e(i)} & & \\ & \omega_{e(i)} & \\ & & \omega_{e(i)} \end{pmatrix} \begin{pmatrix} d_{e(i)x} \\ d_{e(i)y} \\ d_{e(i)z} \end{pmatrix} = Nd \qquad (8) $$

We can calculate strain from geometric equation as follows:

$$ \begin{pmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{yz} \end{pmatrix} = \begin{pmatrix} \frac{\partial u_x}{\partial x} & & & & & \\ & \frac{\partial u_y}{\partial y} & & & & \\ & & \frac{\partial u_z}{\partial z} & & & \\ & & & \frac{\partial u_x}{\partial x} & & \\ & & & & \frac{\partial u_y}{\partial y} & \\ & & & & & \frac{\partial u_z}{\partial z} \end{pmatrix} = \begin{pmatrix} \frac{\partial}{\partial x} & & & & & \\ & \frac{\partial}{\partial y} & & & & \\ & & \frac{\partial}{\partial x} & & & \\ & & & \frac{\partial}{\partial z} & & \\ & & & & \frac{\partial}{\partial x} & \\ & & & & & \frac{\partial}{\partial y} \end{pmatrix} \begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix} \qquad (9) $$

Substituting Eq. (8) into Eq. (9) results in

$$ \begin{pmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{yz} \end{pmatrix} = \sum_{i=1}^{4} \begin{pmatrix} \frac{\partial}{\partial x} & & & & & \\ & \frac{\partial}{\partial y} & & & & \\ & & \frac{\partial}{\partial x} & & & \\ & & & \frac{\partial}{\partial z} & & \\ & & & & \frac{\partial}{\partial x} & \\ & & & & & \frac{\partial}{\partial z} \end{pmatrix} \begin{pmatrix} \omega_{e(i)} & & & & \\ & \omega_{e(i)} & & & \\ & & \omega_{e(i)} & & \\ & & & \omega_{e(i)} & & \\ & & & & \omega_{e(i)} & \\ & & & & & \omega_{e(i)} \end{pmatrix} \begin{pmatrix} d_{e(i)x} \\ d_{e(i)y} \\ d_{e(i)z} \end{pmatrix} = LNd \qquad (10) $$

where L is the matrix form of differential operator of geometric equation, and Nd is the global displacement vector of specific point. Then perform the dot product of matrix L and matrix N, we obtain the relationship between strain and cover function:

$$ \begin{pmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{yz} \end{pmatrix} = \sum_{i=1}^{4} \begin{pmatrix} \frac{\partial \omega_{e(i)}}{\partial x} & & & & & \\ & \frac{\partial \omega_{e(i)}}{\partial y} & & & & \\ & & \frac{\partial \omega_{e(i)}}{\partial z} & & & \\ & & & \frac{\partial \omega_{e(i)}}{\partial x} & & \\ & & & & \frac{\partial \omega_{e(i)}}{\partial y} & \\ & & & & & \frac{\partial \omega_{e(i)}}{\partial z} \end{pmatrix} \begin{pmatrix} d_{e(i)x} \\ d_{e(i)y} \\ d_{e(i)z} \end{pmatrix} = Bd \qquad (11) $$

where B is the transform matrix between strain and cover function. Then substituting Eq. (5) into Eq. (11), B can be written into

$$ B = \begin{pmatrix} f_{12} & f_{22} & f_{32} & f_{42} & & \\ & f_{13} & f_{23} & f_{33} & f_{43} & \\ & & f_{14} & f_{24} & f_{34} & f_{44} \\ f_{13} & f_{12} & f_{23} & f_{32} & f_{43} & f_{42} \\ f_{14} & f_{12} & f_{24} & f_{23} & f_{34} & f_{44} \\ & & f_{14} & f_{24} & f_{34} & f_{43} \end{pmatrix} \qquad (12) $$

The stress is given by 3D linear elastic constitutive equation:
$$ \begin{pmatrix} \sigma_x \\ \sigma_y \\ \sigma_z \\ \tau_{xy} \\ \tau_{xz} \\ \tau_{yz} \end{pmatrix} = \frac{E}{1+\mu} \begin{pmatrix} \frac{1-\mu}{1-2\mu} & \frac{\mu}{1-2\mu} & \frac{\mu}{1-2\mu} \\ \frac{\mu}{1-2\mu} & \frac{1-\mu}{1-2\mu} & \frac{\mu}{1-2\mu} \\ \frac{\mu}{1-2\mu} & \frac{\mu}{1-2\mu} & \frac{1-\mu}{1-2\mu} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{pmatrix} \begin{pmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{xy} \\ \gamma_{xz} \\ \gamma_{yz} \end{pmatrix} \quad (13) $$

where E is the elastic modulus, and μ is the Poisson's ratio. In addition to the complex procedure of cutting elements and PCs, we employed continuum damage mechanics to temporarily describe the behavior of the failed elements. The computational mesh remains unchanged during the calculations. The cutting algorithm of the failed elements and PCs, followed by the calculation of stress applied to the crack surface, will be studied in future research. Continuum damage mechanics is used to characterize the discontinuity after crack propagation. In the present work, continuum damage mechanics has a better programming and calculation efficiency. The elastic modulus of an element can be written as

$$ E = (1 - D) E_0 \quad (14) $$

$$ D = \begin{cases} 0 & (\text{uncracked}) \\ 1 & (\text{cracked}) \end{cases} \quad (15) $$

where E represents the actual elastic modulus, D denotes the damage factor, and E₀ signifies the original elastic modulus.

When the element is not failed, D is set to 0. However, when the element fails, D increases to almost 1, which is determined by the specific model, as shown in Fig. 3. When crack propagation is complete, D will be set to a maximum value. The relationship between D and stress is binary rather than linear. The damage of an element results in discontinuous regions of displacement on both sides of the crack.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Relationship between damage factor D and element states.</figcaption></figure>

## 2.3 Crack propagation based on local tracking algorithm

Various explicit crack tracking algorithms have been developed to model strong discontinuities in numerical methods, such as the fixed crack tracking algorithm, the local crack tracking algorithm, and the global crack tracking algorithm. The local tracking algorithm is considered to be an efficient method to interpret and implement when compared to other tracking algorithms. It has been developed in the context of xFEM to model crack propagation in two or three dimensions. In this section, we aim to incorporate the local tracking algorithm into 3D NMM. In the present work, due to the complexity of geometric calculation in 3D NMM, we only involve the process of crack propagation. All examples only have one initiation global crack surface, and the local crack surface propagates from the global crack surface.

Elements that are susceptible to cracking in two dimensions can be categorized into two types, as shown in Fig. 4: elements with one crack tip, and elements with more than one crack tip. To simplify the process of crack propagation, it is common to describe the global crack surface as a single line, which means that most elements only have one or two crack tips. This approach facilitates an easier modeling and understanding of the crack and reduces the complexity in the numerical calculations. The element with one crack tip typically influences the direction of crack propagation. The local crack surface is determined by two factors: the crack tip and the propagation angle. The crack tip is determined by the adjacent element, while the propagation angle is determined by crack growth criteria such as the maximum tensile stress criterion. Elements with two crack tips are used to smoothen out the
global crack surface. When an element reaches its ultimate strength, the direction of the maximum tensile stress cannot influence the direction of the crack surface. In this case, we only need to connect the two crack tips to generate crack surface. Elements with three crack tips is more complex in a 2D local tracking algorithm. According to Yu et al. (2021), the propagation of elements with three crack tips follows the tip and boundary model (TB model, Fig. 4c) and occurs only when two global crack surfaces intersect.

<figure><img src="image_6.png" /><figcaption>Fig. 4. The local tracking algorithm in 2D: (a) The propagation process of crack in element with one crack tip, (b) The propagation process of crack in element with two crack tips, and (c) The propagation process of a crack in an element with three crack tips.</figcaption></figure>

Although the local tracking algorithm has the same principle in 2D and 3D models, they differ in the flatness of the global crack surface. In 2D models, when an element is cracked, the crack surface propagates into only one other element and generates only one new crack tip. The final formation of the global crack surface is a series of line segments from head to tail. However, in 3D models, when an element is cracked, the local crack surface intersects with the element and generates at least two new crack tips, which then propagate into at least two other elements. If not constrained, the final formation of the global crack surface may become a broken surface instead of a smooth surface. To ensure the smoothness of the global crack surface, we have adopted a conservative local crack generation strategy, which we will illustrate next.

Before describing the specifics of the 3D local tracking algorithm, it is essential to clarify certain concepts. The 3D NMM employs tetrahedral elements, and it is assumed that the model contains only one global crack surface. In order to simplify the local tracking algorithm, several hypotheses were formulated. Firstly, it was assumed that each element edge contains only one crack tip, representing the intersection of the crack surface with the edge. Secondly, it was assumed that each surface of the element features only one crack tip, which corresponds to the intersection of the crack surface with the element surface.

The potential failure element in the recent time step must satisfy two conditions. Firstly, it should be adjacent to a surface that belongs to the failure element and has a crack tip. Secondly, the element must meet the criterion of maximum tensile stress. The classification of potential failure elements is based on the number of crack tips present on their surfaces.

In the simplest scenario, an element possesses only one crack tip, as shown in Fig. 5. The propagation of crack surfaces in these potential elements is primarily influenced by the location of the crack tip and the direction of maximum tensile stress. To determine the local crack surface, it is essential for the crack surface to intersect the crack tip of the adjacent element surface while aligning its normal vector direction closest to the direction of maximum tensile stress, as described in Naderi et al. (2016), and can be seen in Fig. 6. When we only consider the normal vector of maximum tensile stress, the global crack surface will kink because of the new crack surface (imaginary crack surface in Fig. 6). As a result, we apply a rotation to the new crack surface to keep the global crack surface continuous (actual cracked surface in Fig. 6). Subsequently, a plane called the cutting plane is derived by intersecting the crack tip and the normal vector. The normal vector of the cutting plane can be calculated by the Gram-Schmidt orthogonalization method
from linear algebra (Leon et al. 2013).

<figure><img src="image_7.png" /><figcaption>Fig. 5. The local tracking algorithm in 3D: (a) Element with one crack tip, (b) Element with two crack tips, and (c) Element with three crack tips.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Flipping the local crack surface keeps the global crack surface continuous.</figcaption></figure>

The normal vector can be written as

$$ \vec{n}' = \vec{n} - \frac{(\vec{m}, \vec{n})}{(\vec{m}, \vec{m})} \vec{m} \qquad (16) $$

where $\vec{m}$ is the vector of the relate crack tip, $\vec{n}$ is the direction vector of maximum tensile stress, and $\vec{n}'$ is the normal vector of the cutting plane. The cutting plane is then used to partition the potential failure element, resulting in the formation of the local crack surface.

The second situation pertains to potential failure elements with two crack tips. In these elements, the local crack surface must pass through both crack tips. The cutting plane is determined by these two crack tips and is not influenced by the direction of maximum tensile stress. The normal vector of the cutting plane is perpendicular to the vector formed by the two crack tips. The normal vector can be written as

$$ \vec{n}' = \vec{m}_1 \times \vec{m}_2 \qquad (17) $$

where $\vec{m}_1$ and $\vec{m}_2$ are the vectors of the related crack tips. By cutting the element with this plane, the local crack surface is generated. The resulting cross-section of the element can take the form of a quadrilateral or a triangle, and the presence of a new crack tip on the cross-section will impact the adjacent uncracked element.

The final case involves potential failure elements with three crack tips. The local crack surface of these elements must encompass all three crack
tips. To accommodate these crack tips, we construct a tetrahedron by connecting all the vertices of the crack tips, resulting in a crack tetrahedron. In a 3D
crack tracking program, the crack tetrahedron serves the same purpose as the crack surface. Each edge of the crack tetrahedron has the potential to generate
a new crack tip, which can affect adjacent non-cracked elements.

In the 2D local tracking algorithm, the global crack surface is represented as a polyline composed of multiple line segments. For a given element, the local crack surface leads to only one new crack tip where it intersects with other non-cracked elements. However, in three dimensions, the local crack surfaces can take the form of triangles or quadrilaterals, resulting in multiple new crack tips. Representing the global crack surface as a smooth plane is challenging in 3D local tracking algorithms. In our research, we made the assumption that there exists only one continuous crack surface. Therefore, potential failure elements that were not mentioned earlier (such as elements without any crack tip) should be avoided to prevent the conservative local crack generation strategy from prematurely stopping crack surface propagation.

## 2.4 Crack propagation angle correction

In the local tracking algorithm, the direction of crack propagation is controlled by the element with only one crack tip. Elements with more than one
crack tip ensure that the whole crack surface is smooth. However, due to the uncertainty of the mesh, the propagation of the crack surface can result in
actual errors. The crack propagation angle correction is used in the local tracking algorithm. According to Xue et al. (2023), a complex fracture state leads
to many cracks, such as wing cracks and anti-wing cracks. As a result, only the main crack propagation is considered in the present work. According to
Naderi et al. (2016) and Yang et al. (2016), the main crack propagation angles are always larger than 90°. However, due to the influence of the grid, the
local tracking algorithm does not exclude the possibility of the propagation of acute-angle cracks. Hence, an angle correction is applied to avoid the
possibility of acute angle crack propagation in the calculations.

Fig. 7 shows an element $E_2$ with only one crack tip. When the global crack surface propagates from $E_1$ to $E_2$, the direction of the crack surface is determined by the direction of the crack tip and the maximum tensile stress. However, due to the surface shared by the two elements, the angle between the two crack surfaces is implicitly limited to a fan area of 180° ahead of the crack tip. The new crack surface is perpendicular to the direction of the maximum principal stress. However, the plane perpendicular to the maximum principal stress consists of two-half planes (solid red line S1 and dotted red line S2 in Fig. 7c and d). Due to the surface shared by two elements, the angle of crack surface propagation is limited to 180°. Different meshes may lead to different results, as seen when comparing Fig. 7c with Fig. 7d. In fracture mechanics, the main crack surface should propagate along the original plane. To achieve this, we need to determine the half-plane with the obtuse angle between the crack surfaces. However, due to the shared surface between the two elements, there is only one half-plane available for generating the new crack surface, and it is not always the half-plane required to ensure crack propagation in the original plane.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Angle limitation in the local tracking algorithm: (a) Crack propagation situation in 3D space; (b) The angle of crack propagation is limited by the mesh; (c) Result 1: Crack surface propagate to the top of the mesh; and (d) Result 2: Crack surface propagate to the bottom of the mesh.</figcaption></figure>

To avoid errors in crack propagation, we have forwarded a new and accurate method for the local tracking algorithm. To ensure accurate crack propagation using the local tracking algorithm, it is essential to clarify two basic principles. Firstly, the computational mesh must remain unchanged during crack propagation. Secondly, one element should have at most one crack surface. Based on these principles, we have proposed a modified method for crack propagation.

When crack surface tends to propagate from element $E_1$ to $E_2$, $E_2$ has one crack tip that determines the direction of the crack propagation. Next, the
angle between the two crack surfaces needs to be checked. If the angle is acute (Fig. 7d), the desired half-plane is $S_2$, and due to the angle limitation of the mesh, the crack surface propagation is stopped to avoid errors. Other methods to modify the local tracking algorithm will be studied in forthcoming contributions.

# 3 3D NMM Program implementation

## 3.1 Basic data structure

To increase the extensibility of the program and to reduce the code coupling, we have modularized the 3D NMM program. Currently, it consists of two independent modules: the calculation module and the crack module. The calculation module is responsible for assembling the element matrix and calculating the displacement of each PC. On the other hand, the crack module is designed to facilitate crack surface propagation calculations in 3D space. These two modules can work independently, including modules that may be added in the future, as shown in Fig. 8.

<figure><img src="image_10.png" /><figcaption>Fig. 8. Flowchart for 3D-NMM.</figcaption></figure>

The 3D NMM program is implemented in Python using object-oriented programming, enabling seamless manipulation of scattered data by representing them as concrete objects. The program consists of algorithms and data. In the recent implementation, the module acts as a wrapper for an assembly line factory. The algorithm functions as a production machine, referred to as the module component, while the data serve as the product on the assembly line, referred to as the module object, the module object of calculation and crack as shown in Tables 1 and 2.

<table><caption>Table 1. Module object in the calculation module for the 3D NMM program.</caption><thead><tr><th>Module object</th><th>Property</th><th>Data type</th><th>Meaning</th></tr></thead><tbody><tr><td rowspan="10">Element3D</td><td>ID</td><td>Int</td><td>Element ID</td></tr><tr><td>Cracked</td><td>Int</td><td>Whether element cracked</td></tr><tr><td>Material_dict</td><td>Dict</td><td>Detail material parameters</td></tr><tr><td>PC_ID</td><td>List[Int]</td><td>IDs of the related PCs</td></tr><tr><td>PC_displacement</td><td>List[List]</td><td>Displacement of the related PCs</td></tr><tr><td>Initial_stress</td><td>np.ndarray (size = (6, 1))</td><td>Initial stress or strain at the begin of one time step</td></tr><tr><td>Initial_strain</td><td></td><td></td></tr><tr><td>Stiff_matrix</td><td>np.ndarray (size = (12, 12))</td><td>Element matrix, used to assemble global matrix (Matrix <i>a</i> in program)</td></tr><tr><td>Mass_matrix</td><td></td><td></td></tr><tr><td>Fixed_matrix</td><td></td><td></td></tr><tr><td>Initial_force</td><td>np.ndarray (size = (12, 1))</td><td>Element force, used to assemble global force (Vector <i>f</i> in program)</td></tr></tbody></table>
<table><caption>Table 2. Module object in the crack module for the 3D NMM program.</caption><tr><td></td><td>Loading_force</td><td></td><td></td></tr><tr><td></td><td>Body_force</td><td></td><td></td></tr><tr><td></td><td>Fixed_force</td><td></td><td></td></tr><tr><td></td><td>Mass_force</td><td></td><td></td></tr><tr><td></td><td>Special_point</td><td>List[Pointer]</td><td>Pointer to the special points (Point3D)</td></tr><tr><td></td><td>Geometric_entity</td><td>vecCell</td><td>Geometric entity used to geometric calculation</td></tr><tr><td rowspan="4">Point3D</td><td>ID</td><td>Int</td><td>Point ID and the Element ID which it belongs to</td></tr><tr><td>Element_id</td><td></td><td></td></tr><tr><td>Type</td><td>Int</td><td>Loading point or fixed point</td></tr><tr><td>Coordinate</td><td>List[-float]</td><td>Coordinate of the point</td></tr><tr><td></td><td>Displacement</td><td>np.ndarray (size = (6, 1))</td><td>Displacement increase of the point in each time step</td></tr></table>

<table><thead><tr><th>Module object</th><th>Property</th><th>Data type</th><th>Meaning</th></tr></thead><tbody><tr><td rowspan="2">All module object<br>in crack module</td><td>ID</td><td>Int</td><td>ID of module object</td></tr><tr><td>Geometric_entity</td><td>vecCell</td><td>Geometric entity used to geometric calculation</td></tr><tr><td rowspan="5">CrackElement3D</td><td>Cracked</td><td>Int</td><td>Whether crack element cracked</td></tr><tr><td>Elementsurface</td><td>List[Pointer]</td><td>Pointer to element surfaces</td></tr><tr><td>Crack Cole</td><td>List[Pointer]</td><td>Pointer to crack surface (at most one)</td></tr><tr><td>Stress</td><td>np.ndarray (size = (6, 1))</td><td>Stress or strain of the crack elements</td></tr><tr><td>Strain</td><td></td><td></td></tr><tr><td rowspan="3">ElementSurface3D</td><td>Cracked</td><td>Int</td><td>Whether crack element cracked</td></tr><tr><td>Crack element</td><td>List[Pointer]</td><td>Pointer to crack elements (at most two)</td></tr><tr><td>Crack edge</td><td>List[Pointer]</td><td>Pointer to crack edge (at most one)</td></tr><tr><td rowspan="2">CrackSurface3D</td><td>Crack element</td><td>List[Pointer]</td><td>Pointer to crack elements (only one)</td></tr><tr><td>Crack edge</td><td>List[Pointer]</td><td>Pointer to crack edges</td></tr><tr><td rowspan="2">CrackEdge3D</td><td>Element surface</td><td>List[Pointer]</td><td>Pointer to element surfaces (only one)</td></tr><tr><td>Crack surface</td><td>List[Pointer]</td><td>Pointer to crack surface (at most two)</td></tr></tbody></table>

A module component is an abstract object that represents a data handling process. For example, the calculation module consists of three module components: an element creator, a matrix assembler, and an element refresher (The first letter of the abstract object in programming will be capitalized). At the beginning of the calculation, the element creator wraps scatter data into elements. The matrix assembler extracts matrices from each element and assembles them into a global matrix. When the displacement of PCs is solved, the element refresher updates the properties of the PCs and elements.

A module object is an abstract object that wraps a series of scatter data, it usually has its own physical entities. Module objects are handled by different module components in turn. The fundamental module objects we work with are elements, which correspond to the manifold elements in the physical system. In addition to elements, other objects are used to assist in modular calculations, like crack surfaces in the crack module. Properties is also an abstract conception in programming. Discrete data is organized in abstract object and accessed as properties. Different types of elements may have the same properties, such as an element ID, material parameters, a stress tensor, a strain tensor, and more. In different modules, elements may have unique properties specific to each module. For example, elements in the calculation module have the property to calculate the element matrix, while elements in the crack module have a pointer to access adjacent elements.

Different types of elements serve a common and important function: representing geometric entities. In the program, a geometric entity is a property within an element that can be accessed at any time as shown in Fig. 9. To represent a geometric entity, we introduce the concept of a geometric kernel. Geometric kernel is a concept from computer aided design (CAD), a geometric entity encompasses both geometry and topology. A geometric entity consists of vertexes, edges and surfaces. We can access the information of these substructures from the selected geometric entity.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Geometry kernel of the 3D NMM in every module.</figcaption></figure>

The geometric kernel improves computational efficiency in two ways. First, the geometric kernel can be applied to every part involved in the geometric calculation in 3D NMM. We must test every piece of code to ensure the reliability of the program. Different pieces of code usually require different test examples. Hence, we wrapped the pieces of code that involve the geometric entity and geometric calculation. We only need to test this part of the packaged code to ensure its reliability, so that it can be applied to every part in 3D NMM. Second, in coding the geometric kernel, we can apply some open-source algorithms such as VTK. Open-source algorithms also have good computational efficiency and robustness. Geometry pertains to the coordinate information of each point within the entity. For instance, in the case of a tetrahedron element, geometry is represented by the coordinates of its four vertices. Topology refers to the set of properties that remain unchanged under specific geometric transformations. The method for modeling topology is boundary representation (B-rep) (Lienhardt, 1991). In the case of a tetrahedron element, topology consists of a tetrahedron composed of four triangles. Each triangle is comprised of three edges, and each edge consists of two vertices. All of these elements represent the topological information of the tetrahedron element and possess their own distinct geometry and topology.

## 3.2 Implementation of calculation and crack module

The calculation module encompasses several functions, such as computing the element matrix, assembling the global matrix, solving for PC displacement, and calculating the interpolation of element displacement, as shown in Fig. 10a. These functions are implemented through the module component. In 3D NMM, there are six types of element matrices: stiff matrix, initial stress matrix, point load matrix, body load matrix, fixed point matrix, and inertia matrix. After determining the PC displacement, the subsequent step involves calculating the displacement of the manifold element vertices, which is achieved through shape function interpolation.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Flowcharts in each module in the proposed 3D NMM: (a) Calculation module, and (b) Crack module.</figcaption></figure>
The crack module is specifically designed to calculate 3D crack surface propagation, as shown in Fig. 10b. Due to the involvement of computational geometry in this process, numerous processes and geometric objects, such as elements, need to be concretely defined to generate module components and module objects. Within the crack module, there are three types of module components: crack element creator, element cracker, and crack element refresher. Additionally, four types of module objects are used to facilitate the calculation of crack propagation: crack element, element surface, crack surface, and crack edge.

Within the crack module, each module object possesses its own geometric entity, which becomes a property of that object. The crack element serves as the representation of a manifold element within the physical system and acts as the foundational unit handle in the crack module. Additionally, there are four link states that exist within the crack module: Crack element with element surface, crack element with crack surface, crack surface with crack edge, and element surface with crack edge, as shown in Fig. 11. Fig. 11a shows the crack element in the program. We can obtain an arbitrary crack element by ID, and we can obtain another module object by pointers in the crack element.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Detailed description of the geometric kernel in the crack module in the 3D NMM.</figcaption></figure>

A crack element is associated with four pointers that indicate its connection to element surfaces. These element surfaces represent the geometric surfaces of a manifold element within the physical system. Each element surface also maintains two pointers that allow access to its adjacent crack elements. Consequently, accessing adjacent elements through pointers is the most common manipulation during crack propagation.

Furthermore, each crack element contains a pointer to a crack surface, implying that a manifold element can have at most one local crack surface, which is dependent on its crack state. The crack surface represents the local crack surface within the physical system and includes a pointer to its superior crack element. The crack state of an element is indicated by one of four possibilities (Fig. 12): (1) Element with initial cracking, (2) Undamaged element, (3) Potential damaged element, and (4) Damaged element. The four states are only used to indicate the element states in the program. The damage factor *D* does not have any relationship with the four element states.

<figure><img src="image_14.png" /><figcaption>Fig. 12. The four states of elements in programming.</figcaption></figure>

Each crack surface is associated with multiple pointers to crack edges, with the number of pointers determined by the shape of the cross-section, which can be either a quadrilateral or a triangle. A crack edge is formed at the intersection of a crack surface and an element surface, representing the boundary of the local crack surface. Additionally, a crack edge maintains two pointers to its adjacent crack surfaces.

An element surface contains only one pointer, which points to the corresponding crack edge, depending on the state of the associated crack element. Similarly, it is possible to obtain the element surface from the crack edge using a similar approach. The calculation of the number of crack edges intended
to crack within a crack element is crucial in the local tracking algorithm, and it is concurrently updated with the crack state of the elements. It is worth noting that a tetrahedron crack surface is special, as it possesses six crack edges. However, we only consider the edges that intersect with the element surface.

The process of an element transitioning from a state of uncracked to cracked involves multiple steps. First, the neighboring crack element is cracked and generates crack edge. The state of the potentially cracked element is marked to indicate that it can be cracked. When the stress on the crack element exceeds a predefined threshold, the corresponding crack plane is generated using the employed local tracking algorithm, and the state of the crack element is marked as a cracked element. The program determines the count of crack edges on the element surface and calculates the normal vector of the local crack surface. Then, it identifies all points where the crack surface intersects with the element surface and connects them to generate crack edges and crack surface. Finally, the information within the model is updated by incorporating local crack surfaces into the global crack surface, modifying the crack flag of the crack element, and establishing pointers from the crack element to the crack surface and from the element surface to the crack edge. Additionally, the propagated crack in the model is the dominant one, and other cracks have only a little effect on the strength of the model. If a crack element reaches its strength limit but does not have a crack edge on its element surface, it will not be cracked.

# 4. Numerical benchmark examples

## 4.1 Stress concentration around a circular hole

Stress concentration around a circular hole is a classic example in elastic mechanics (Jaeger et al., 2009) and has been used, for example, to formulate the widely-used pore-emanating crack model of Sammis and Ashby (1986). We use this example here to validate 3D NMM. The study of stress concentration marks the beginning of fracture mechanics, but it is beyond the scope of our discussion.

In 3D NMM, all boundary conditions are applied according to special points. We use densely regularly arranged special points to simulate the applied surface boundary conditions, where the loading points and fixed points are also special points. In the present work, loading points are used to simulate the stress boundary condition. The fixed points are also used to simulate the constant velocity boundary condition, where the coordinates of these fixed points change with time steps.

Detailed information about the model can be found in Fig. 13a. The numerical model is computed only using the calculation module and consists of 30,419 manifold elements and 6654 PCs. The dimensions of the model are 10 m × 2 m × 10 m. An equivalent force with an opposite direction is applied to both the upper and lower surfaces. In 3D NMM, all boundary conditions are applied at specific points. Each surface is subjected to 961 loading points (31 × 31), with a force value of 10 N applied to each loading point. The model is subjected to an equivalent tensile stress ($q_0$) of 480.5 Pa.

The results of the example are presented in Fig. 13b and c, showing a maximum stress component of approximately 1450 Pa around the hole, which is roughly three times the value of $q_0$. Fig. 13c displays the stress component in the z direction extracted from the measuring line. The two dotted lines represent the maximum and minimum values along the thickness direction. Additionally, the curve of the analytic solution exhibits three other significant points, confirming the validity of the 3D NMM.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Results of stress concentration around the circular hole using 3D NMM: (a) Boundary condition and mesh, (b) Stress component in the z direction, and (c) Curve of stress component in the z direction.</figcaption></figure>

## 4.2 Uniaxial stress on pre-crack disc

All the models to validate 3D NMM are made of linear elastic material, the mechanical properties of which are provided in Table 3. Because of the behavior of the material is elastic prior to failure, scaling the size up or down does not significantly affect the results, and the default unit of size is meters. All the model uses tetrahedrons as manifold elements.

<table><caption>Table 3. Material parameters of numerical examples for the Brazilian disc simulations using 3D NMM.</caption><thead><tr><td>Young's modulus, E (GPa)</td><td>Poisson's ratio, μ</td><td>Density (kg/m³)</td><td>Young's modulus of fixed point spring (GPa)</td><td>Tensile strength (MPa)</td></tr></thead><tbody><tr><td>10</td><td>0.2</td><td>1</td><td>1000</td><td>0.1</td></tr></tbody></table>

The second numerical simulation involves a pre-cracked Brazilian disc test. The pre-cracked Brazilian disc has a diameter of 0.2 m and a thickness of 0.04 m. The boundary condition is such that a fixed velocity of 0.0001 m/step is applied to the upper surface, while the lower surface remains fixed without displacement. The model consists of 3599 elements and 1077 PCs. The length of the initial global crack surface is 0.02 m. There are two crack angles: 90° and 60°. The initial crack surface propagates all the generated crack surfaces.
The failure areas of the two models are located along the central axis of the disc and perpendicular to the direction of maximum tensile stress. This numerical example demonstrates that the local tracking algorithm can propagate the crack in its original direction or cause it to change direction. Furthermore, it validates the reasonableness of using the maximum tensile stress criterion for crack propagation. The process of crack propagation is illustrated in Figs. 14 and 15. Our results compare well with those of Yang et al. (2016), who also modeled crack propagation using a 3D-NMM code.

<figure><img src="image_16.png" /><figcaption>Fig. 14. Model results for a Brazilian disc simulation containing a 90° angle pre-crack using 3D NMM: (a) Model setup and crack angle, (b) Step = 0, (c) Step = 12, and (d) Step = 24.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Model results for a Brazilian disc simulation containing a 60° angle pre-crack using 3D NMM: (a) Model setup and crack angle, (b) Step = 0, (c) Step = 12, and (d) Step = 24.</figcaption></figure>

## 4.3 Fracture of L-shape concrete specimens

The third numerical example examines the fracture behavior of an L-shaped concrete specimen (Fig. 16). Fig. 16a illustrates the geometry and boundary conditions of the specimen. A displacement-based loading scheme is employed, with fixed points and a defined time step. Specifically, this model consists of 5642 elements and 1711 PCs. We apply a constant velocity in a narrow rectangular area (width of 1 cm) on the surface of model. As shown in Fig. 16, the primary mode of failure observed in the L-shaped specimen is tensile failure, and the resulting pattern of the global crack surface compares well with the results of Naderi et al. (2016), who used a 3D augmented FEM code (3D A-FEM). It is important to note that the deformation of the global crack surface is expected to deviate from a planar shape, allowing for the evaluation of the ability of the crack surface to navigate using the local tracking algorithm.

## 4.4 Rectangular plate with single hole

The final example investigates crack propagation in a rectangular plate with a single hole, and it is compared with the results from Yang and Chen (2023). Displacement is applied to the upper surface of the plate at each time step, while the bottom surface remains fixed. The model comprises a total of 3712 elements and 1173 PCs. The initial global crack surface is located at heights of 10 mm (Fig. 17) or 14 mm (Fig. 18) within the model. This initial global crack surface leads to the generation of crack surfaces in 8 manifold elements.
<figure><img src="image_18.png" /><figcaption>Fig. 16. Model results for an L-shaped specimen simulation using 3D NMM: (a) Model setup illustrating the geometry and boundary conditions of the specimen, (b) Step = 10, (c) Step = 20, (d) Step = 30, (e) Final result of numerical example showing the global crack surface, and (f) Result from Naderi et al. (2016) for comparison (all unit in mm).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Model results for a rectangular plate simulation with a single hole ($h = 10$) using 3D NMM: (a) Model setup illustrating the geometry and boundary conditions of the specimen, (b) Mesh of the model, (c) Step = 15, (d) Step = 30, (e) Step = 60, and (f) Result from Yang and Chen (2023) for comparison (all unit in mm).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 18. Model results for a rectangular plate simulation with a single hole ($h = 14$) using the 3D NMM: (a) Step = 15, (b) Step = 30, (c) Step = 60, and (d) Result from Yang and Chen (2023) for comparison.</figcaption></figure>

The objective of these last two examples is to assess the capability of the local tracking algorithm to simulate the pattern of the global crack surface. When the initial crack surface is at a height of 10 mm, the global crack surface propagates directly across the entire specimen. However, when the height is set to 14 mm, the global crack surface propagates and then turns toward the center hole of the specimen. Our model results are consistent with the results presented in Yang and Chen (2023), who used a 3D-NMM code with tetrahedral meshes (Figs. 17 and 18).

## 5 Discussion and limitation

The MC system and PC system offer a unique approach to integrating continuous and discontinuous problems within a single governing equation, distinguishing it from the FEM. In the case of dealing with problems caused by crack propagation, it is only necessary to cut the PC and the corresponding manifold elements. This approach enables the computation of strong discontinuous problems without the need for remeshing during the calculation process. Within the 3D NMM program, the accurate implementation of geometric calculations, such as polyhedron cutting, holds critical importance (Yang et al., 2016; Yang and Chen, 2023). The program in our research relies on a robust geometry kernel, which handles all geometric representations and operations. Recent studies have highlighted two primary functions of the geometric kernel: storing geometric information and executing geometric operations. We have adopted the boundary representation (B-rep) method for storing 3D entities (Requicha and Voelcker, 1983), which originates from computer graphics. Additionally, we have utilized or intend to utilize certain open-source algorithm libraries (Schroeder et al., 1998). Looking ahead, future work will involve a variety of geometric calculations. These include the generation of manifold elements using MCs and PCs during pretreatment, polyhedron cutting with arbitrary planes, and contact determination between two distinct blocks.

There are various numerical methods available for modeling crack propagation, including xFEM (Zi and Belytschko, 2003) and the phase field method (Borden et al., 2012; Xu et al., 2022). The xFEM utilizes the sign function to approximate displacement across the crack, while the phase field method considers the evolution of potential energy in the control equation. However, neither of these methods physically cuts the mesh. In contrast, the NMM physically cuts the manifold elements and the PCs to model strong discontinuous cracks (Yu et al., 2022). This approach leads to more accurate crack tracing and calculations following fracture. Although the crack propagation in this study exhibits similarities to real-world scenarios, it does not perfectly align in terms of both time and space. Moving forward, the future plan for crack calculation involves separating the PCs and manifold elements during the calculation process. This separation aims to provide a quantitative description of cracks in both time and space.

Tracking cracks in three dimensions presents a substantial computational challenge. Several tracking algorithms have been developed specifically for crack tracking in 3D space, such as global tracking algorithms (Riccardi et al., 2017) and the level set method (Gravouil et al., 2002; Moës et al., 2002). Each algorithm has its own unique advantages. Among these algorithms, the local tracking algorithm stands out as the simplest approach to track cracks. While local tracking algorithms have been developed successfully in FEM (Liu et al., 2014; Naderi et al., 2016). The local track algorithm in three dimension comes with certain limitations. The presence of crack surfaces within individual elements leads to a global crack surface that is not always smooth. In cases where there is a significant crack angle, the global crack surface can exhibit kinking and bifurcation phenomena. To mitigate these issues, it is important to ensure that the crack tracking path is not excessively long and that the crack angles are not too large. However, further research is necessary to explore alternative methods and enhance the local tracking algorithm's performance in three dimension.
## 6 Conclusions

In the present paper, we introduce the development of a 3D NMM by using a geometric kernel to facilitate the multiple geometric calculations involved and to enhance computational efficiency. This geometric kernel is designed to store the geometric and topological information of various entities, such as MCs and manifold elements. It also facilitates geometric calculations, including block cutting operations. The maximum tensile stress criterion is adopted as a crack growth criterion to achieve strong discontinuous crack growth, and a local crack tracking algorithm and an angle correction technique are incorporated into the 3D-NMM to address minor limitations of the algorithm in a 3D model. This local tracking algorithm represents the effect of cracks on the displacement field through the damage (i.e. reduction) of mechanical parameters. Subsequently, the program is implemented in Python, using object-oriented programming in two independent modules: a calculation module and a crack module to enhance the performance of the algorithm. Finally, we demonstrate the feasibility, accuracy, and effectiveness of the enhanced algorithm in the 3D-NMM using four numerical examples. Our results demonstrate that the 3D-NMM, incorporating crack local tracking algorithms, can effectively handle 3D cracking problems. This study establishes the potential of the 3D-NMM, combined with the local tracking algorithm, for accurately modeling 3D crack propagation in brittle rock materials. As a result, the 3D-NMM emerges as an exciting prospect to model deformation and failure scenarios in natural and geoengineered rock masses, improving hazard monitoring and mitigation strategies.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 42172312 and 52211540395). M. Heap acknowledges support from the Institut Universitaire de France (IUF).

## References

An, X., Ning, Y., Ma, G., He, L., 2014. Modeling progressive failures in rock slopes with non-persistent joints using the numerical manifold method. Int. J. Numer. Anal. Methods Geomech. 38, 679-701.

Annavarapu, C., Settgast, R.R., Vitali, E., Morris, J.P., 2016. A local crack-tracking strategy to model three-dimensional crack propagation with embedded methods. Comput. Meth. Appl. Mech. Eng. 311, 815-837.

Borden, M.J., Verhoosel, C.V., Scott, M.A., Hughes, T.J.R., Landis, C.M., 2012. A phase-field description of dynamic brittle fracture. Comput. Meth. Appl. Mech. Eng. 217-220, 77-95.

Chen, Y., Zheng, H., Yin, B., Li, W., 2023. The MLS-based numerical manifold method for Darcy flow in heterogeneous porous media. Eng. Anal. Bound. Elem. 148, 220-242.

Gravouil, A., Moës, N., Belytschko, T., 2002. Non-planar 3D crack growth by the extended finite element and level sets-Part II: Level set update: NON-PLANAR 3D CRACK GROWTH-PART II. Int. J. Numer. Methods Eng. 53, 2569-2586.

He, L., Ma, G., 2010. DEVELOPMENT OF 3D NUMERICAL MANIFOLD METHOD. Int. J. Comput. Methods. 07, 107-129.

Hoek, E., Carranza-Torres, C., Corkum, B., 2002. Hoek-Brown failure criterion-2002 edition. Proceedings of NARMS-Tac, 1, 267-273.

Hu, M., Rutqvist, J., 2020. Numerical manifold method modeling of coupled processes in fractured geological media at multiple scales. J. Rock Mech. Geotech. Eng. 12, 667-681.

Hudson, J.A., Priest, S.D., 1983. Discontinuity frequency in rock masses. Int. J. Rock Mech. Min. Sci. Abstracts. Elsevier, 73-89.

Jaeger, J.C., Cook, N.G.W., Zimmerman, R. 2009. Fundamentals of rock mechanics. John Wiley & Sons.

Jäger, P., Steinmann, P., Kuhl, E., 2008. Modeling three-dimensional crack propagation-A comparison of crack path tracking strategies. Int. J. Numer. Methods Eng. 76, 1328-1352.

Jäger, P., Steinmann, P., Kuhl, E., 2009. Towards the treatment of boundary conditions for global crack path tracking in three-dimensional brittle fracture.
Comput. Mech. 45, 91-107.

Kang, G., Ouyang, Q.M., Ning, Y.J., Chen, P.W., 2023. Development of three-dimensional numerical manifold method with cover-based contact theory. Eng. Anal. Bound. Elem. 155, 182-196.

Leon, S.J., Björck, Å., Gander, W., 2013. Gram-Schmidt orthogonalization: 100 years and more. Numer. Linear Algebr. Appl. 20, 492-532.

Li, W., Zheng, H., Yu, X., Jia, C., Sun, X., 2022. MLS-Based Numerical Manifold Method for Modeling the Cracked Rock Considering the Contact of the Crack Surface. Front. Earth Sci. 9:825508.

Li, X., Zhang, Q., Li, J., Zhao, J., 2018. A numerical study of rock scratch tests using the particle-based numerical manifold method. Tunn. Undergr. Space Technol. 78, 106-114.

Li, X., Zhao, J., 2019. An overview of particle-based numerical manifold method and its application to dynamic rock fracturing. J. Rock Mech. Geotech. Eng. 11, 684-700.

Li, Z., Xu, T., Zhao, L., Liu, Y., Xu, Y., Heap, M.J., Utili, S., Liu, B., Su, B., 2024. Enhancing stability analysis of open-pit slopes via integrated 3D numerical modeling and data monitoring. Eng. Fail. Anal. 163, 108495.

Lienhardt, P., 1991. Topological models for boundary representation: a comparison with n- generalized maps. Computer-Aided Design, 23, 59-82.

Liu, F., Zhang, K., Liu, Z., 2019. Three-dimensional MLS-based numerical manifold method for static and dynamic analysis. Eng. Anal. Bound. Elem. 109, 43-56.

Liu, W., Yang, Q.D., Mohammadizadeh, S., Su, X.Y., 2014. An efficient augmented finite element method for arbitrary cracking and crack interaction in solids: EFFICIENT A-FEM FOR ARBITRARY CRACKING AND CRACK INTERACTION IN SOLIDS. Int. J. Numer. Methods Eng., 99, 438-468.

Moës, N., Gravouil, A., Belytschko, T., 2002. Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model: NON-PLANAR 3D CRACK GROWTH-PART I. Int. J. Numer. Methods Eng. 53, 2549-2568.

Naderi, M., Jung, J., Yang, Q.D., 2016. A three dimensional augmented finite element for modeling arbitrary cracking in solids. Int. J. Fract., 197, 147-168.

Ning, Y., Lu, Q., Liu, X., 2023. Fracturing failure simulations of rock discs with pre-existing cracks by numerical manifold method. Eng. Anal. Bound. Elem. 148, 389-400.

Oliver, J., Huespe, A.E., 2004. Continuum approach to material failure in strong discontinuity settings. Comput. Meth. Appl. Mech. Eng., 193, 3195-3220.

Oliver, J., Huespe, A.E., Samaniego, E., Chaves, E.W.V., 2004. Continuum approach to the numerical simulation of material failure in concrete. Int. J. Numer. Anal. Methods Geomech., 28, 609-632.

Requicha, A., Voelcker, H., 1983. Solid Modeling: Current Status and Research Directions. IEEE Comput. Graph. Appl., 3, 25-37.

Riccardi, F., Kishta, E., Richard, B., 2017. A step-by-step global crack-tracking approach in E-FEM simulations of quasi-brittle materials. Eng. Fract. Mech. 170, 44-58.

Saloustros, S., Cervera, M., Pelà, L., 2019. Challenges, Tools and Applications of Tracking Algorithms in the Numerical Modelling of Cracks in Concrete and Masonry Structures. Arch. Comput. Method Eng. 26, 961-1005.

Sammis, C.G., Ashby, M.F., 1986. The failure of brittle porous solids under compressive stress states. Acta Metallurgica, 34, 511-526.

Schroeder, W., Martin, K.M., Lorensen, W.E., 1998. The visualization toolkit an object-oriented approach to 3D graphics. Prentice-Hall, Inc.

Shi, G.H. 1992a. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of deformable block structures. Eng. Comput. 9, 157-168

Shi, G.H. 1992b. Manifold method of material analysis. Army Research Office Research Triangle Park NC.

Stolarska, M., Chopp, D.L., Moës, N., Belytschko, T., 2001. Modelling crack growth by level sets in the extended finite element method. Int. J. Numer. Methods Eng. 51, 943-960.

Szwedzicki, T., 2003. Rock mass behaviour prior to failure. Int. J. Rock Mech. Min. Sci. 40, 573-584.

Terada, K., Asai, M., Yamagishi, M., 2003. Finite cover method for linear and non-linear analyses of heterogeneous solids. Int. J. Numer. Methods Eng., 58, 1321-1346.

Terada, K., Kurumatani, M., 2004. Performance assessment of generalized elements in the finite cover method. Finite Elem. Anal. Des. 41, 111-132.

Tong, D., Yi, X., Tan, F., Jiao, Y., Liang, J., 2023. Three-dimensional numerical manifold method for heat conduction problems with a simplex integral on the boundary. Sci. China Technol. Sci. 67, 1007-1022.
Xu, B., Xu, T., Xue, Y., Heap, M.J., Ranjith, P.G., Wasantha, P.L.P., Li, Z., 2022. Phase-field modeling of crack growth and interaction in rock. Geomech. Geophys. Geo-Energy Geo-Resour. 8, 180.

Xue, Y., Xu, T., Heap, M.J., Meredith, P.G., Mitchell, T.M., Wasantha, P.L.P. 2023. Time-dependent cracking and brittle creep in macrofractured sandstone. Int. J. Rock Mech. Min. Sci. 162, 105305.

Yahaghi, J., Liu, H., Chan, A., Fukuda, D., 2023. Development of a three-dimensional grain-based combined finite-discrete element method to model the failure process of fine-grained sandstones. Comput. Geotech., 153, 105065.

Yang, S., Cao, M., Ren, X., Ma, G., Zhang, J., Wang, H., 2018. 3D crack propagation by the numerical manifold method. Comput. Struct. 194, 116-129.

Yang, S., Chen, R., 2023. A new strategy for 3D non-persistent crack propagation by the numerical manifold method with tetrahedral meshes. Eng. Anal. Bound. Elem. 148, 190-204.

Yang, Y., Tang, X., Zheng, H., Liu, Q., He, L., 2016. Three-dimensional fracture propagation with numerical manifold method. Eng. Anal. Bound. Elem. 72, 65-77.

Yu, X.Y., Xu, T., Heap, M., Zhou, G.L., Baud, P., 2018. Numerical Approach to Creep of Rock Based on the Numerical Manifold Method. Int. J. Geomech., 18, 04018153.

Yu, X.Y., Xu, T., Heap, M.J., Baud, P., Reuschlé, T., Heng, Z., Zhu, W.C., Wang, X.W., 2021. Time-dependent deformation and failure of granite based on the virtual crack incorporated numerical manifold method. Comput. Geotech., 133, 104070.

Yu, X.Y., Xu, T., Heap, M.J., Heng, Z., Ranjith, P.G., Su, B., Wasantha, P.L.P., Sun, G., 2023. Time-dependent virtual crack model of rock with application to slope stability. Eng. Anal. Bound. Elem. 154, 172-185.

Yu, X.Y., Xu, T., Heap, M.J., Heng, Z., Zhu, W., Zhou, G., Su, B., 2022. A Virtual Crack-Based Numerical Manifold Approach to Crack Initiation, Propagation and Coalescence in Granite. Rock Mech. Rock Eng. 55, 7791-7816.

Zhang, L., Guo, F., Zheng, H., 2022. The MLS-based numerical manifold method for nonlinear transient heat conduction problems in functionally graded materials. Int. Commun. Heat Mass Transf. 139, 106428.

Zheng, H., Liu, F., Li, C., 2014. The MLS-based numerical manifold method with applications to crack analysis. Int. J. Fract., 190, 147-166.

Zheng, Y., Yan, C., Zheng, H., 2023. Modified joint element constitutive model for FDEM to simulate the nonlinear mechanical behavior of rocks. Comput. Geotech., 164, 105831.

Zhou, G.L., Xu, T., Konietzky, H., Zhu, W., Heng, Z., Yu, X.Y., Zhao, Y., 2022. An improved grain-based numerical manifold method to simulate deformation, damage and fracturing of rocks at the grain size level. Eng. Anal. Bound. Elem. 134, 107-116.

Zi, G., Belytschko, T. 2003. New crack-tip elements for XFEM and applications to cohesive cracks. Int. J. Numer. Methods Eng. 57, 2221-2240.

<figure><img src="image_21.png" /></figure>

Boyi Su is a PhD candidate in Center for Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang, China. His research interests mainly focus on numerical simulations of crack propagation and interaction in brittle materials, especially with three-dimensional numerical manifold method (3D NMM).
<figure><img src="image_22.png" /></figure>

**Tao Xu** received his PhD from Northeastern University, Shenyang, China in 2005. He is a professor and a doctoral supervisor in the School of Resources and Civil Engineering at Northeastern University, China. His research focuses on the time-dependent deformation and fracturing of brittle rocks, associated long-term stability of rock slope, and thermo-hydro-mechanical coupling process of rock mass. Tao is actively participating in professional activities, acting as an assessor of National Natural Science Foundation of China, China Scholarship Council (CSC), Australian Research Council, Czech Science Foundation, some provincial and ministerial science and technological awards, and he has ever academically visited City University of Hong Kong, EPFL (École Polytechnique Fédérale de Lausanne) in Switzerland, Monash University in Australia, University of Melbourne in Australia, University of Strasbourg in France, University College London in UK, and TU Graz in Austria. He has authored over 100 international journal papers and co-authored three monographs.
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>