Numerical Simulation of Crack Propagation and Branching Behaviors in Heterogeneous Rock-like Materials

Abstract

The characterization and understanding of crack evolution in non-uniform geological structures are crucial for predicting the mechanical response of rock-like materials or structures under varying loading conditions. In this study, an improved Peridynamic model with a degree of heterogeneity characterized by random pre-breaking “bonds” coefficients is introduced to capture the intricacies of crack initiation, propagation, and branching behaviors in heterogeneous rock-like materials. MATLAB discrete programs for heterogeneous material models and PD simulation programs based on the FORTRAN language were developed. The effectiveness of the heterogeneous PD model in simulating crack propagation and branching patterns in heterogeneous materials has been verified through dynamic and static (quasi-static) loading cases with pre-notch. The different levels of heterogeneity not only affect the direction of crack propagation but also determine the crack deflection direction and branching patterns. The crack propagation path appears to possess obvious asymmetry in the crack propagation direction. As the load applied continues to increase, the asymmetric multi-crack branching phenomenon will occur. The higher the level of heterogeneity, the more complex the behaviors of crack propagation and branching become. This research provides valuable insights into the interplay of material heterogeneity and crack evolution, offering a foundation for improved numerical simulations and contributing to the broader field of geomechanics.

1. Introduction

Crack propagation and branching behaviors of rock-like materials, such as rock mass, cement, concrete, and other brittle materials, significantly affect the stability and safety of explorations of deep resources, building structures, and underground space. In addition, they are involved in various deep geotechnical engineering applications such as hydrogen storage, shale gas extraction, and coal mining [1,2,3,4]. A good understanding of crack propagation and the branching of brittle materials with pre-existing defects, such as rock mass or rock-like material, is significant for preventing the catastrophic failure of deep geotechnical engineering. The numerical simulation of fracture and failure in rock-like materials stands as a cornerstone in geomechanics, offering invaluable insights into the behavior of rock-like material formations under diverse loading conditions [5,6,7]. As the integrity and stability of geological structures profoundly impact engineering endeavors, understanding the mechanisms of fracture in rock-like materials becomes paramount.

Rock-like materials, such as rock mass or concrete, have complex physical and mechanical properties that differ from metal materials, such as heterogeneity, discontinuity, and anisotropy [8,9,10,11]. In engineering, on-site experiments, physical laboratory simulations, and numerical simulations are generally used to analyze the laws of material or structure movement, providing economic and safe suggestions and measures for engineering. However, using on-site and laboratory physical simulation experiments, only a limited amount of strain data can be measured, and the dynamic failure characteristics and microscopic evolution characteristics of the material or structure cannot be accurately observed. Existing mechanics experimental instruments, such as MTS, can load specimens with quasi-static or static loads, but there are certain difficulties in applying dynamic and static loads. Model experiments have varying degrees of conditional limitations in terms of cost, duration, and operating conditions. The state of destruction cannot be reproduced and preserved for a long time. On-site experiments are costly and greatly affected by terrain, construction conditions, geological conditions, and other factors [12,13,14]. Numerical experiments have significant advantages in terms of cost, dynamics, and microscopic failure, making numerical simulation an important means of studying the mechanical properties of rock-like materials.

Crack propagation and branching behavior have important significance in brittle materials. Many computational methods have been helpful in simulating the propagation and branching of cracks in brittle materials. The numerical simulation of crack propagation is a crucial area of research in geomechanics and has profound implications for the understanding of material behavior under various loading conditions. Various numerical methods have been developed to model crack propagation in rock-like materials, each with its own set of assumptions, advantages, and limitations. Finite Element Method (FEM), Discrete Element Method (DEM), Extended Finite Element Method (XFEM), and Messfree Methods (MF) are among the widely employed techniques [15,16,17,18,19]. These methods enable researchers to simulate the initiation, propagation, and interaction of cracks under different loading conditions. FEM has been extensively used in simulating crack propagation by discretizing the rock mass into finite elements and solving for the stress and strain fields. In FEM analyses, the material constitutive model assumption plays an important role for the correct estimation of the response. For continuum mechanics, the most prominent material laws for the analysis of geomaterials are the modified Cam-Clay, Hoek–Brown, and Barton models [20,21,22]. However, FEM may face challenges in accurately capturing crack paths and complexities associated with multiple interacting cracks [23]. On the other hand, DEM focuses on the discrete representation of individual particles, providing a more granular understanding of crack initiation and propagation. Nevertheless, DEM may encounter challenges in representing the continuum behavior of material at smaller scales [24]. The XFEM is designed to address some of the limitations of traditional FEM by allowing cracks to be modeled independently of the underlying mesh. While XFEM enhances crack representation, challenges remain in its application, particularly in handling complex geological structures and material interfaces [25]. As an integral non-local continuum mechanics method, Peridynamics (PD) discretizes the material or structure into a series of points, uses the integral equation to describe the motion, and obtains the position, velocity, and force of the points by solving the non-local integral equation of motion, so as to obtain the motion trajectory and material response of structures subject to external force [26,27]. The PD equation of motion contains the differential of time and the integral of space, does not contain the spatial derivative of any displacement field, and is defined in both continuous and discontinuous regions. It overcomes the difficulty of analyzing discontinuous problems by continuum mechanics and realizes the solution of multi-scale problems according to a unified framework. The PD method does not need to preset the crack initiation position and crack propagation path when dealing with fracture and failure problems. In PD, using spatial integral equations instead of PDEs to describe the force distribution on objects avoids the singularity of differential simulations in traditional continuum mechanics when solving discontinuous problems. It has shown great potential in dealing with discontinuity, fracture and damage, interactions between multi-cracks, spontaneous nucleation, and propagation and branching of cracks [28,29].

The numerical modeling of heterogeneous geomaterials is a critical domain within geomechanics, essential for comprehending the intricate behavior of geological formations characterized by non-uniform material properties [30]. The inherent complexity arising from heterogeneity in materials or structures necessitates sophisticated numerical approaches to accurately capture their mechanical response under diverse loading conditions. The key to conducting mechanical analysis of heterogeneous geomaterial is to establish a computational mechanical model that reflects the heterogeneity characteristics. Currently, numerical modeling methods for describing the heterogeneity of geomaterials include probability statistics [31,32], spatial variation [33], digital image technology [34,35], and machine learning [36,37]. To some extent, the modeling method of probability statistics can characterize the heterogeneity of geomaterials but the dual characteristics of the overall structural and local randomness of geomaterial parameters are not reflected in this method. In the spatial variation-based modeling method, the purpose of reducing estimation variance is achieved by smoothing the discreteness of measured data, although sometimes there may be some outliers that are smoothed out. Therefore, this method can only ensure the local optimal estimate rather than the overall optimal estimate. The modeling method based on digital image technology relies on the processing of digital images of geomaterials, which can accurately and objectively reflect the micro-structure of geomaterials when analyzing relatively small specimens or structures. For engineering projects, digital image technology needs to be processed and analyzed, which also limits the application of digital image technology in large-scale geotechnical engineering projects [38], along with the enormous amount of image information.

This study aims to explore a spectrum of numerical modeling methods employed for simulating the behavior of heterogeneous rock-like materials, providing an encompassing view of both their strengths and limitations. To overcome the cumbersome modeling, computational complexity, and limitations of traditional numerical simulation methods in dealing with discontinuous problems, the sub-heterogeneous PD model is employed in this paper, which is based on the random pre-breaking “bonds” [39]. In this model, only the degree of heterogeneity represented by the pre-breaking “bonds” coefficient needs to be set to establish a heterogeneous model of materials or structures, without the need for an additional model pre-process, and the model is able to reflect the characteristics of the random distribution of material heterogeneity. The size of the simulation model has no effect on the efficiency of heterogeneous modeling. The sub-homogeneous PD model is particularly suitable for applications in large-scale discontinuous and heterogeneous geotechnical engineering. The accurate simulation of fractures is pivotal for assessing the safety and performance of infrastructures such as tunnels, dams, and slopes.

2. Theory of Bone-Based Peridynamics

2.1. Basic Theory

As shown in Figure 1, at a certain time $t$, the PD equation of motion at a point $x$ is [26]:

$$\rho \ddot{u}\left(x,t\right)={\displaystyle {\int }_{H_x}f\left(u\left(x^{\prime },t\right)-u\left(x,t\right),x^{\prime }-x\right)\mathrm{d}{V}_{x^{\prime }}}+b\left(x,t\right)$$

(1)

where $u$ is the displacement vector field, $b$ is the body face vector, and $f$ is a pairwise force vector function in a PD bond that connects material point $x^{\prime }$ with point $x$.

The radius $\delta$ is called the horizon range, $H_x$ is the collection of all the material points within the range $\delta$, specifically:

$$H_x=H\left(x,\delta \right)=\left\{\left.x^{\prime }\in R:‖x^{\prime }-x‖\le \delta \right\}\right.$$

(2)

By introducing the relative distance $\xi$ in reference to the configuration and relative displacement $\eta$, it can help to discretize Equation (1). $\xi =x^{\prime }-x$ is the relative position vector and $\eta =u^{\prime }-u$ is the relative displacement vector.

The mechanical properties of the material are related to the micro-modulus $c$ and critical elongation $s_0$. The deformation energy density of the material point $x$ in its horizon $\delta$ is:

$$W=\frac{1}{2}{\displaystyle {\int }_{H_x}w\left(\eta ,\xi \right)}\mathrm{d}{V}_{\xi }=\frac{1}{2}{\displaystyle {\int }_0^{\delta }\left(\frac{c{s}^2\xi }{2}\right)}4\pi {\xi }^{2}d\xi$$

(3)

where $W$ is deformation energy density, $w\left(\eta ,\xi \right)$ is the scalar micro-potential function, and $1/2$ indicates that each point occupies half of that of the related two material points’ deformation energy density.

The pairwise force magnitude for micro-elastic materials with damage is derivable from the scalar micro-potential $w\left(\eta ,\xi \right)$:

$$f\left(\eta ,\xi \right)=\frac{𝜕w\left(\eta ,\xi \right)}{𝜕\eta }=\left\{\begin{array}{ll}\frac{\xi +\eta }{‖\xi +\eta ‖}g\left(s\right)\mu \left(t,\xi \right)& ‖\xi ‖\le \delta \\ \hfill 0\hfill & \hfill ‖\xi ‖>\delta \hfill \end{array}\right.$$

(4)

$$g\left(s\right)=cs, \forall s$$

(5)

$$\mu \left(t,\xi \right)=\left\{\begin{array}{cc}1\hfill & \hfill \hspace{1em}if s\left(\xi \right)<s_0\hfill \\ 0\hfill & \hfill \hspace{1em}otherwise\hfill \end{array}\right.$$

(6)

$$s=\frac{‖\xi +\eta ‖-‖\xi ‖}{‖\xi ‖}$$

(7)

where $\mu \left(t,\xi \right)$ is a time-related scalar function, $s$ is the value of “bond” elongation between two material points, and $s_0$ is the critical value of “bond elongation” for breakage.

In PD, the local damage value ranges from 0 to 1 and is represented by the ratio of the number of bonds broken to the total number of bonds in the horizon of point $x$.

$$\phi \left(x, t\right)=1-\frac{{\displaystyle {\int }_{H_x}\mu \left(x\prime -\mathit{x}, t\right)\mathrm{d}V\prime }}{{\displaystyle {\int }_{H_x}\mathrm{d}V\prime }}$$

(8)

Considering the variation in the spatial decline of long-range force, in the improved model, the constitutive force function $f\left(\eta ,\xi \right)$ is given as:

$$f\left(\xi ,\eta \right)=\left\{\begin{array}{cc}\frac{\xi +\eta }{‖\xi +\eta ‖}{c}_0\kappa \left(\xi ,\delta \right)s& \hspace{1em}‖\xi ‖\le \delta \\ 0& \hspace{1em}‖\xi ‖>\delta \end{array}\right.$$

(9)

where $\kappa \left(\xi ,\delta \right)$ is the kernel function that can reflect the characteristics of the long-range force spatial decline between points. In the present study, kernel function $\kappa \left(\xi ,\delta \right)$ is as follows:

$$\kappa \left(\xi ,\delta \right)={\left(1-\frac{‖\xi ‖}{\delta }\right)}^2$$

(10)

With the classical elastic energy density theory $c\left(\xi ,\delta \right)$ equal to:

$$c\left(\xi ,\delta \right)=\left\{\begin{array}{cc}\frac{10E}{\pi {\delta }^4\left(1-2\nu \right)}{\left(1-\frac{‖\xi ‖}{\delta }\right)}^2& 3\mathrm{D}\\ \frac{60E}{\pi {\delta }^3\left(1-\nu \right)}{\left(1-\frac{‖\xi ‖}{\delta }\right)}^2& plane stress\\ \frac{60E}{\pi {\delta }^3\left(1+\nu \right)\left(1-2\nu \right)}{\left(1-\frac{‖\xi ‖}{\delta }\right)}^2& plane strain\\ \frac{12E}{A{\delta }^2}{\left(1-\frac{‖\xi ‖}{\delta }\right)}^2& 1\mathrm{D}\end{array}\right.$$

(11)

where $E$ is Young’s modulus and $\nu$ is Poisson’s ratio.

For 2D problems as shown in Figure 2, the fracture energy $G_0$ in PD is given by:

$$G_0=2{\displaystyle {\int }_0^{\delta }{\displaystyle {\int }_{\mathrm{z}}^{\delta }{\displaystyle {\int }_0^{{\cos }^{-1}\left(z/\xi \right)}\left[\frac{c_0\kappa \left(\xi ,\delta \right)s_0^2\xi }{2}\right]\xi \mathrm{d}\phi \mathrm{d}\xi \mathrm{d}z}}}$$

(12)

where $z$ is the distance from material point $x$ to the fracture surface and $\theta$ is the rotation angle.

Substituting Equation (12) into Equation (13), the approximate expression of energy release rate $G_0$ in PD can be obtained as:

$$G_0=\frac{\left(\pi -1\right){\delta }^4{c}_0{s}_0^2}{120}$$

(13)

Approximate expression of energy release rate $G_0$ and critical elongation $s_0$ can be obtained. For 2D problems, the critical elongation $s_0$ can be obtained as:

$$s_0=\left\{\begin{array}{cc}\sqrt{\frac{4\pi G_0}{3\left(\pi -1\right)E\delta }}& \hspace{1em}plane stress\\ \sqrt{\frac{5\pi G_0}{4\left(\pi -1\right)E\delta }}& \hspace{1em}plane strain\end{array}\right.$$

(14)

2.2. Numerical Simulation Method

Solving the motion equation of the PD computing model involves time integration, spatial discretization, and integration. The unit volume force generated by the point $x$ in the reference configuration can be obtained by summing, and the spatial integral equation can be discretized as follows:

$${\rho }_i{\ddot{u}}_i^n+{C\dot{u}}_i^n={\displaystyle \sum _{j=1}^{N}f\left(u_j^n-u_i^n,x_j-x_i\right){ V}_j+b\left(x_i^n\right)}$$

(15)

where ${\rho }_i$ is the material density, ${\ddot{u}}_i^n$ is the acceleration, ${\displaystyle \sum _{j=1}^N}$ is the combined internal force of all other points in the horizon $\delta$ of point $x$, $C$ is the artificial damping coefficient, $C=0$ in dynamic analysis, $u_j^n=u\left(x_j,t\right)$, $n$ is the time-step number, and subscripts denote the point number.

During simulations, material point $x_j$ located in the horizon range of point $x_i$ may not be fully contained, as shown in Figure 3. The green points represent that they have complete integrated volume, but the red points are not within the horizon size of point $x_i$. The blue points mean they are out of the horizon of point $x_i$. In order to ensure the accuracy of the integration, the integral volume needs to be corrected; for example, the red points.

By introducing a volume correction factor $v_j^c$, the integral volume of each point can be expressed as [40]:

$$V_j=\left\{\begin{array}{cc}{\left|\Delta x\right|}^2& \hspace{1em}\xi <\left(\delta -r_j\right)\\ v_j^c{\left|\Delta x\right|}^2& \hspace{1em}\left(\delta -r_j\right)<\xi <\delta \end{array}\right.$$

(16)

where $r_j$ is half the length of the discretized subdomain, $r_j=\frac{\left|\Delta x\right|}{2}$, and volume correction factor $v_j^c$ can be determined as:

$$v_j^c=\left\{\begin{array}{cc}1& \hspace{1em}{\xi }_{ij}<\left(\delta -r_j\right)\\ \frac{\delta + r_j - {\xi }_{ij}}{2r_j}& \hspace{1em}\left(\delta -r_j\right)<{\xi }_{ij}<\delta \end{array}\right.$$

(17)

When $\left(\delta -r_j\right)<{\xi }_{ij}<\delta$, the volume correction factor $v_j^c$ varies linearly between $\left(\frac{1}{2},1\right)$. After introducing the volume correction factor $v_j^c$, Equation (15) can be rewritten as:

$$\rho {\ddot{u}}_i^n+C{\dot{u}}_i^n={\displaystyle \sum _{j=1}^{N}f\left(u_j^n-u_i^n,x_j-x_i\right)}\left(v_j^c{V}_j\right)+b\left(x_i^n\right)$$

(18)

The numerical integration of time can be obtained by the velocity Verlet difference scheme. The displacement $u\left(x,t\right)$ and velocity $\dot{u}\left(x,t\right)$ of the material point are:

$$\left\{\begin{array}{l}{\dot{u}}^{mid}=\frac{\rho - 0.5C\Delta t}{\rho }{\dot{u}}_i^n+\frac{\Delta t}{2\rho }{(L_u+b)}^n\\ {\dot{u}}_i^{n+1}=\frac{\rho }{\rho + 0.5C\Delta t}{\dot{u}}^{mid}+\frac{\Delta t}{2(\rho + 0.5C\Delta t)}{(L_u+b)}^{n+1}\\ u_i^{n+1}=u_i^n+\frac{\rho - 0.5C\Delta t}{\rho }{\dot{u}}_i^n\Delta t+\frac{{(\Delta t)}^2}{2\rho }{(L_u+b)}^n\end{array}\right.$$

(19)

where $\Delta t$ is the time step, which is far less than the critical time step $\Delta t_c$ and needs to meet $\Delta t\ll \Delta t_c=\frac{\Delta \chi }{C_L}$, $\Delta \chi$ is the minimum length of the “bond”, and $C_L$ is the longitudinal wave velocity in the material.

3. Heterogeneity Characterization in the PD Model

3.1. PD Modeling of Heterogeneous Materials

When defects in materials or structures have obvious geometric shapes, the following two methods can usually be used in PD simulation models to process the points within the defect range. One is to uniformly discretize points within the entire PD simulation model, delete the points within the defect range, and renumber them. Similarly, in the other method, points in the entire PD simulation model are uniformly discretized, without the need to delete points within the defect range. Only the “bonds” that pass through the defect range need to be cut down in advance. To determine whether the “bonds” are broken, the method shown in Figure 4 can be used. AB represents defects and CD represents “bonds” that pass through the defect range.

For pre-breaking criterion Figure 4a, the following equation needs to be met:

$$\left(CA\times AB\right)\cdot \left(DA\times AB\right)<0$$

(20)

For pre-breaking criterion Figure 4b, the following equation needs to be met:

$$\left(AC\times AB\right)\cdot \left(BC\times AB\right)<0$$

(21)

In existing homogeneous PD models, it is generally believed that the physical and mechanical parameters related to the mechanical response of materials, such as elastic modulus, density, and fracture energy, are uniform and constant. Due to the varying number of natural defects in materials or structures, they do not match the actual situation. When the defect size in a heterogeneous material is much smaller than the size of the computational model, it can be considered a locally homogeneous material [41,42]. However, when there are a large number of micro defects in the material that cannot be ignored and have no specific geometric shape, existing PD models cannot accurately describe the mechanical response of the material or structure. As shown in Figure 5, by introducing a random pre-breaking “bond” coefficient into a uniform discrete computational model to characterize the heterogeneity of rock-like materials, the details of our work are present in the literature [39].

Figure 6 shows the degree of material heterogeneity and the corresponding damage index represented by the random pre-breaking of “bonds”. It should be pointed out that the damage index shown in the nephogram is not the damage caused by the external load on the material, but represents the proportion of random pre-breaking “bonds” and the level of material heterogeneity [43,44] as the initial condition of the model in the simulation process.

3.2. Numerical Solution Process

Based on the above theories and simulation methods, this paper uses the FORTRAN language to compile the simulation program of the heterogeneous PD model, and the flow chart of the PD model simulation program is shown in Figure 7. In the pre-processing module, the discrete model program is compiled based on MATLAB 2018a, which can obtain the required information such as the coordinates of material points. In the post-processing module, the simulation file information is read with the help of Ensight, which also outputs the simulation results.

4. Numerical Examples of Specimens with Pre-Notch

4.1. Load Case 1: Quasi-Static Loading

In order to study the influence of material heterogeneity on crack growth under quasi-static loading, numerical simulation of crack growth is carried out for three-point bending beams with pre-notch. Considering the three-point bending beam shown in Figure 8, a pre-notch with length h = 0.04 m and width w = 0.004 m is set in the middle of the beam bottom, and the geometric dimensions of the beam are L = 0.6 m, H = 0.1 m, l = 0.5 m. The three-point bending beam model is orthogonally uniformly dispersed to a 600 × 100 × 1 lattice, and the distance between material points is |Δx| = 0.001 m. The heterogeneity represented by the pre-breaking “bond” coefficient is also considered in order to avoid premature failure at the loading boundary, and the horizon size is selected as δ = 6|Δx|. The 160 material points at the pre-notch are deleted. The total number of material points in the simulation model is 59,840.

Table 1. Physical and mechanical parameters of the three-point bending beam.

Density	Young’s Modulus	Fracture Energy Density	Poisson’s Ratio
2483 kg/m3	7.2 GPa	23 N/m	1/3

presents the physical and mechanical parameters of the tested cement three-point bending beam with pre-notch.

Figure 9 shows the crack growth process of the homogeneous three-point bending beam model with pre-notch. From the simulation results, it can be seen that the micro cracks start to aggregate at the tip of the pre-notch and form macro cracks, which propagate along a straight line in the direction of the applied concentrated load until the beam is completely fractured. The cracks generated in this model are typical I-open cracks. However, most materials in nature have typical heterogeneous characteristics. For the heterogeneous three-point bending beam model with pre-notch under quasi-static loading, the crack growth path will be significantly different from that of the homogeneous model.

The fracture and failure of naturally heterogeneous materials have strong nonlinear characteristics and are gradual accumulation processes. When the external load reaches a certain critical value, the crack will develop in a chain method, and the damage will evolve from a random distribution to the crack surface. Figure 10 shows the PD simulation results of the crack growth and fracture failure of the three-point bending beam with pre-notch when heterogeneity is M = 2.5%. Similar to the homogeneous model, micro cracks aggregate from the tip of the pre-notch and form macroscopic cracks. It can be seen from the damage nephogram and horizontal displacement nephogram during crack growth that the crack develops along the direction of the applied concentrated load but deflects about 7° to the right.

The degree of heterogeneity of the material has a significant impact on its mechanical response behaviors such as fracture failure. Figure 11 shows the PD simulation results of the crack propagation path and fracture patterns of the three-point bending beam with pre-notch when the degree of heterogeneity M = 5%. Similar to the two previous simulation models, micro cracks start to aggregate at the tip of the pre-notch and form macro cracks. It can be seen from the damage nephogram and horizontal displacement nephogram during crack growth that the crack develops along the direction of the applied concentrated load but deflects approximately 10° to the left.

4.2. Load Case 2: Dynamic Loading

In order to study the influence of material heterogeneity on crack growth under dynamic loading, numerical simulation of crack growth is carried out for a rectangular plate with pre-existing notch, as shown in Figure 12. The rectangular plate has a side length of L = 0.1 m and width of 2H = 0.05 m, as well as pre-notch with a length of L/2. The physical and mechanical parameters of the tested concrete rectangular plate with pre-existing notch are shown in

Table 2. Physical and mechanical parameters of the rectangular plate simulation model.

Density	Young’s Modulus	Fracture Energy Density	Poisson’s Ratio
2460 kg/m3	13.3 GPa	85 N/m	1/3

. The model adopts an orthogonal uniform dispersion of a 400 × 200 × 1 lattice, the distance between material points is |Δx| = 0.00025 m, and the total number of discrete material points is 80,000. Because the heterogeneity represented by the pre-breaking “bond” coefficient is considered, in order to avoid premature failure at the loading boundary, a large horizon size is selected, with δ = 6|Δx|. For the model, the boundary condition of force loading is adopted. The instantaneous force load is, respectively, applied to the three-layer material points at the upper and lower ends of the rectangular plate. The left and right sides of the rectangular plate are free boundaries without any constraint.

Figure 13, Figure 14 and Figure 15 show the crack growth process of the homogeneous rectangular plate model with pre-notch under different transient loads. From the numerical simulation results, it can be seen that under the action of an external load, micro cracks all initiate, expand, and form macro cracks from the tip of the pre-notch, and different external load sizes have an important impact on the crack growth path: when the load is small, the crack grows along a straight line in the direction of the pre-notch. With the increase in the external load, the crack propagates along a straight line as a pre-notch, and at the same time, the crack branches. When the external load continues to increase, the crack branching behaviors will occur earlier. Because the simulation model is uniformly discrete and completely homogeneous, the crack growth path is completely symmetrical in the direction of crack growth under different instantaneous loads.

The heterogeneity of the material has an important influence on their mechanical response mechanism under an external load. Figure 16 and Figure 17 show the crack growth process of rectangular plates with pre-notch under different levels of heterogeneity when the instantaneous load is 1.0 MPa and 2.0 MPa, respectively. At this time, the micro crack will still initiate from the tip of the pre-notch and propagate along the pre-notch direction. Different from the homogeneous model, the crack no longer propagates along a straight line, and different levels of heterogeneity will have a greater impact on the crack growth path. The crack propagation path appears to have obvious asymmetry in the crack propagation direction. As the load applied continues to increase, the asymmetric multi-crack branching phenomenon will occur.

5. Discussion

Rock-like materials are naturally generated porous media, composed of mineral particles of different scales and discontinuous structures such as pores and micro-cracks of various scales, exhibiting the typical characteristics of heterogeneous materials. Due to the presence of micro defects such as pores and voids, the complexity of rock-like materials is greatly increased. The material properties of porous media are one of the important sources of heterogeneity in rock-like materials and have a significant impact on physical and mechanical properties such as density and strength, as well as mechanical response mechanisms such as damage accumulation and progressive failure under external loads. By constructing an improved PD model that reflects heterogeneity, the mechanical mechanisms of damage accumulation and progressive failure in rock-like materials can be better explained.

Case studies are conducted to analyze the crack propagation and fracture processes of three-point bending beams with pre-notch under quasi-static loading. The different levels of heterogeneity not only affect the crack propagation direction but also determine the crack deflection direction and deflection angle. It can be seen from the simulation results of the above example based on the heterogeneous PD model that the degree of material heterogeneity has an important influence on fracture and failure characteristics. The different degrees of heterogeneity not only affect the crack propagation direction but also determine the crack deflection direction and deflection angle. The higher the degree of material heterogeneity, the larger the deflection angle.

The heterogeneity of materials has an important influence on their mechanical response mechanisms under external loads. The crack propagation and branching behaviors in heterogeneous rock-like materials have been investigated by improved PD simulation, through the PD simulation analysis of rectangular plates with pre-notch under dynamic loading. The greater the applied instantaneous load, the earlier the crack branching behaviors will appear. The higher the level of material heterogeneity, the more complex the crack branching phenomenon will appear in the process of crack growth. Sufficient loading magnitude is an essential factor for crack propagation, specifically for crack branching.

6. Conclusions

The study presented in this paper focuses on the influence of material heterogeneity on its mechanical response. The improved PD model with degrees of heterogeneity characterized by random pre-breaking “bonds” coefficients is introduced to capture the intricacies of crack initiation, propagation, and branching in heterogeneous rock-like materials, which are helpful to directly investigate the internal microscopic failure process of rock-like materials in various geotechnical engineering applications.

(1)

In the improved PD model, uniform discretization is also performed, and the physical and mechanical parameters of the material are uniformly valued based on the heterogeneity coefficient. During simulation, only the heterogeneity coefficient represented by pre-breaking “bonds” is needed to establish the heterogeneous PD model, without additional pre-processing of the model, and can reflect the random distribution characteristics of the heterogeneity of rock-like materials. The size of the computational model and the number of material points have no effect on the efficiency of heterogeneous modeling, making it particularly suitable for applications in large-scale discontinuous and heterogeneous geotechnical engineering.

(2)

For heterogeneous materials, the crack propagation path appears to have obvious asymmetry in the crack propagation direction. As the load applied continues to increase, the asymmetric multi-crack branching phenomenon will occur. The higher the level of heterogeneity, the more complex the behaviors of crack propagation and branching become. The simulation results verify the effectiveness of the heterogeneous PD model in simulating the fracture and failure of heterogeneous materials, in which the heterogeneity of microscopic materials can be well matched, providing a high-precision model with realistic heterogeneity for numerical investigations.

(3)

The results indicate that the non-uniform characteristics of the material have a significant impact on its mechanical response mechanism under external loading. These findings further validate the effectiveness of the sub-homogeneous PD model in addressing fracture and failure issues in heterogeneous brittle rock-like materials. This research provides valuable insights into the interplay of material heterogeneity and crack evolution, offering a foundation for improved numerical simulations and contributing to the broader field of geomechanics.