# A 3D Peridynamic Model to Simulate Indirect Tensile Failure in Marble

Authors: Wenwu Tan; Keping Zhou; Shitong Zhou; Weicheng Ren
Journal: Geofluids
 (2022)

---

## Abstract

Driven by the commitment of achieving peak carbon dioxide emissions before 2030 and carbon neutrality before 2060, China is promoting the optimization of its energy and industry structures. In the electricity supply industry, a large number of different types of sophisticated large and giant hydropower stations are planned or being constructed in western China. In future decades, hydropower may take the dominant position in China’s electricity supply. Rock fracturing behavior plays a vital role in geotechnical hazard prevention and geotechnical engineering design. With the construction of hydropower stations, rock mechanics investigation will be of great concern. This study presents a hybrid model that coupled the finite element method and peridynamic theory to simulate the indirect tensile strength of marble. Sensitivity analyses of the loading rate, mesh size, and fracture energy release rate are performed. The numerical results indicate that both the loading rate and mesh size significantly affect the numerical representation of rock properties. After the calibration of the fracture energy release rate, Brazilian tensile strength modeling is successfully conducted, and the numerical results are consistent with the experimental results.


---

## Body

## 1. Introduction

With the rapid development of the Chinese economy, industry faces a substantial demand for energy, particularly electricity [1]. It is known that coal, hydro, nuclear, wind, and solar are the main types of energy that are used to generate electricity in power plants in China [2]. Fossil fuel-based power, which causes serious environmental pollution problems, remains China’s largest source of electricity. In addition, China commits that it will strive to peak carbon dioxide emissions before 2030 and achieve carbon neutrality before 2060 [3]. As a consequence, the government has clarified the policy regarding new thermal power plant projects. Hydropower stations will be developed more often in future decades in China. The construction of hydropower stations is accompanied by complex geotechnical problems, such as rock deformation, damage, and failure [4, 5]. Understanding the rock fracturing process is meaningful for the prevention of geotechnical hazards and the design of hydropower stations.Laboratory tests are inevitable for the investigation of rock fracturing mechanisms [6–9]. In the early years of such testing, the rock failure process was observed by the naked eye. The monitoring approach for rock deformation and fracturing processes has been largely developed with scientific and technological improvements [10]. Yang et al. evaluated the strength and infrared radiation characteristics of rock under the influence of freeze–thaw cycle treatment [11]. Niu et al. studied the stress and strain development of rock under tension using the digital image correlation method [12]. Yao and Xia performed a dynamic fracturing behavior investigation of notched semicircle granite specimens by means of a high-speed camera system [13]. Ma et al. evaluated rock block stability using a remotely positioned laser Doppler vibrometer [14]. This advanced noncontact measurement technology provides new insight into understanding the rock failure mode and mechanism. Laboratory tests involve the procedures of coring, sampling, polishing, and the configuration of monitoring equipment, which are time- and cost-consuming [15]. Scholars are pursuing efficient and low-cost approaches for rock mechanic testing [16].In the past few decades, numerical simulation has been treated as a third promising method (in terms of theoretical analysis and experimental testing) for solving mathematical or physical equations that involve geotechnical engineering problems [17–19]. With the introduction of high-performance computing, numerical models could be a very high-resolution technique of approximation in practical engineering [20]. Scholars have developed several numerical approaches, such as the finite element method (FEM), discrete element method, boundary element method, finite difference method, smoothed-particle hydrodynamics, and extended FEM. Each numerical method has advantages for solving specific problems with proper boundary and initiation conditions. However, a single numerical approach may face constraints in solving complicated geotechnical problems. For example, to overcome the disadvantage of FEM in solving the large-deformation problem of blasting, a hybrid algorithm that couples FEM and smoothed-particle hydrodynamics was developed, in which the rock near the blasthole is built by smoothed-particle hydrodynamics and the far-field rock with small deformation is built by the FEM [21]. Another example, to efficiently simulate rock crack initiation and propagation subjected to an external load, is the finite-discrete element method [22]. In rock material simulation, researchers focus on the path of crack propagation and the resulting fracturing mode; thus, the development of a numerical method that is able to simulate the fracturing behavior of rock is considered [23]. As a new nonlocal continuum mechanics formulation, peridynamic theory has great advantages in the application of material fracturing simulation and has been widely used in the field of rock mechanics investigation [24].Aiming at accurately capturing the mechanical properties of marble numerically to improve hydropower station construction, a hybrid method combining FEM and peridynamic theory is employed to simulate the Brazilian tensile strength (BTS) of rock. Sensitivity analysis of numerical parameters such as loading rate, mesh size, and fracture energy release rate are performed. The numerical results demonstrate that the fracture mode and force versus displacement match very well with the experimental results.


## 2. Bond-Based Peridynamic Theory

FEM has been widely used in the solution of continuum mechanics. However, the governing equations of FEM include the differentiation of space. As a result, a singularity arises near the discontinuity region and influences crack propagation [25]. Unlike FEM, peridynamic theory is a meshless method based on nonlocal scope theory, which was first introduced by Dr. Stewart Silling in 2000 [26]. The peridynamic theory of mechanics brings together the mathematical descriptions of continuous fluids, fractures, and particles into one framework [27]. The theory of peridynamics was first known as the “bond-based peridynamic theory”. Subsequently, state-based peridynamic theory was proposed to overcome the disadvantage of the bond-based peridynamic method in that its bonds are not restricted to central forces, or to a Poisson’s ratio of 1/4, as with the bond-based method [28–30]. In peridynamic theory, a solid material is discretized into a series of material points that represent a specific volume and mass. Furthermore, each material point interacts with other material points within a certain range adjacent to it, and the mechanical behavior of the material is characterized by the spatial integral equation. The relationships between the mechanical displacement, relative reference position, and relative displacement are described in Figure 1. Compared with classical partial differential equation models, the peridynamic nonlocal continuum model for solid mechanics is an integrodifferential equation that does not involve spatial derivatives of the displacement field. Therefore, peridynamic theory has been successfully applied to simulate crack propagation and predict damage in a variety of situations [31–34].Figure 1
The motion of a peridynamic continuum body.The equation of motion of arbitrary pointx at time t can be described as follows:
(1)ρxu¨x,t=∫Hxfη,ξdVx′+bx,t,where ρ is the mass density, Hx is the neighborhood in a specific radius δ, fη,ξ is the force function between points x and x′, the relative position is ξ=x′−x, the relative displacement is η=ux′,t−ux,t, and bx,t is the density function of the prescribed body force.In bond-based peridynamic, the material can be treated as microelastic.fη,ξ is the bond force, which can be described as follows:
(2)fη,ξ=∂2Wη,ξ∂η∂ξ=∂2Wη,ξ∂η∂ξ=∂wη,ξ∂η∀ξ,η,where Wη,ξ is a potential function.In the microelastic material, the total potential energy of the bond at pointx is
(3)Wη,ξ=12cs2ξ,where c represents the micromodulus and is related to the elastic constant of the material, which can be obtained based on the bulk modulus,
(4)c=18Kπδ4.s in Equation (3) denotes the elongation of the bond:
(5)s=ξ+η−ξξ.In the bond-based peridynamic model, the material failure criterion is that the elongation of the bond between two material points exceeds the thresholdsc. sc is related to the fracture energy release rate G0 in classical fracture mechanics,
(6)G0=πcsc2δ510.


## 3. Numerical Model and Sensitivity Analysis

### 3.1. Laboratory Test of Marble BTS

Marble cored from western China was selected for the experimental test, and this marble is composed of calcite 94.8%, mica 3.2%, and dolomite 2.0%. The density of the marble is 2698.11 kg/m3. The line loading setup is shown in Figure 2. The rock specimens were disc samples with a diameter and thickness of 50 mm and 25 mm, respectively [35]. Due to the dimension errors that arose during coring, the real thickness and diameter of the disc sample were 48.73 mm and 25.42 mm, respectively. The tested dimensions were used in the subsequent strength calculation. The surfaces of a disc specimen were polished to avoid any irregularities in the thickness of the sample. In this study, the water content was not considered. The rock specimen was dried in a 105°C oven for 24 hours. The axial load was applied by the upper loading platen under the displacement control mode at a speed of 0.01 mm/min by using a hydraulic material testing machine. The load was applied continuously until the rock specimen failed, and the maximum load at failure was recorded. The indirect tensile strength was derived in accordance with
(7)BTS=2PπDt,where P, D, and t denote the applied axial load at failure, test specimen diameter, and test specimen thickness, respectively.Figure 2
Schematic diagram of the BTS test in the laboratory.Scanning electron microscopy (SEM) is widely used to analyze surfaces and can provide a highly magnified high-resolution image of the surface of the sample studied. A TESCAN Mira3 LMU scanning electron microscope (TESCAN, Czech Republic) was used to observe the microstructure of the marble. As Figure3 shows, the rock sample is relatively homogeneous, and the structure is compact. Marble was selected for this testing because it has a low intrinsic anisotropy. Thus, it is acceptable to simplify marble as an isometric material in numerical simulations.Figure 3
Microimage of marble.


### 3.2. Numerical Model Configuration

As Figure4 shows, the numerical model (case no. M1) configuration was consistent with the model depicted in Figure 2, in which the diameter and thickness of a disc specimen were 50 mm and 25 mm, respectively. The disc specimen was sandwiched by two flat platens, and the deformation of the two platens was not considered. To save computing time, BTS modeling was performed using a hybrid model that couples the FEM and a peridynamic model. The interaction between FEM and the peridynamic model was achieved by a contact algorithm. Two flat platens were established in the FEM mesh, and a rock model is built with peridynamic theory. In this hybrid model, the FEM part was built with the conventional method. Bond-based peridynamics was implemented into the code in the FEM framework with discontinuous Galerkin theory. Technically, the peridynamic-based part of a model is built by converting the FEM model using a detaching function. The top loading platen was applied with a constant downward velocity boundary condition until the maximum axial load was recorded in the results.Figure 4
Numerical model configuration.


### 3.3. Effect of Loading Rate

In the laboratory, a conventional static BTS test of rock may take several minutes to complete. If the loading rate remains consistent with the laboratory test in numerical simulation with an explicit solver, the long computing time is unacceptable. Thus, according to the rock response to the loading rate, a quasistatic loading method can be introduced for static problem simulation, such as unconfined compressive test modeling, BTS modeling, and fracture toughness modeling [36, 37]. The principle of the quasistatic loading method is to enhance the loading rate in the numerical simulation to reach a balance between efficiency of the computing process and accuracy of the numerical results. To examine the effect of the loading rate on the BTS modeling results in a peridynamic model, 6 sets of BTS modeling schemes with different loading rates were carried out using a coarse mesh model, as shown in Figure 5. The two rigid platens contained 16200 hexahedral meshes (25116 nodes), and a peridynamic domain (86016 nodes) was built based on 10752 FEM meshes distributed in the form of a butterfly. The corresponding loading rates were 0.5 m/s, 0.1 m/s, 0.05 m/s, 0.01 m/s, 0.005 m/s, and 0.001 m/s, and the lower rigid platen was fixed.Figure 5
The convergence of the force under different loading rates.Figure5 shows that the loading rate plays a significant role in the force response. In the laboratory test, marble shows apparent brittle properties, which are characterized by the force level dropping dramatically after the peak. At a loading rate of 0.5 m/s, the rock response to the high-speed loading tends to be ductile, which deviates from the essential character of the marble. With the loading rate decreasing, the curves of force versus displacement become closer, which indicates that the influences of these loading rates on the force are converging. According to the magnified local view in Figure 5, the three curves are quite similar when the loading rate varies from 0.001 to 0.01 m/s. In particular, the difference in the peak force at loading rates of 0.001 m/s and 0.005 m/s is only 85.39 N. Obviously, the influence of the loading rate on the force is negligible when the loading rate is less than 0.005 m/s. Therefore, a loading rate of 0.005 m/s is determined for subsequent numerical simulation under the quasistatic loading condition.


### 3.4. Mesh Independence

For a typical BTS test or modeling, the expected result is that the disc specimen is split into two approximately identical halves with gentle shear damage at the points of contact between the loading platens and rock specimen. Therefore, it is necessary that a sufficiently fine mesh should be distributed along the loading diameter. It is known that tensile cracks initiate from the center of the disc specimen and propagate along the loading diameter in a successful BTS test. In a specific peridynamic model, the butterfly mesh strategy provides the finest mesh in the center of the disc specimen, which ensures the accuracy of initial tensile crack generation.In the mesh-based numerical approach, its accuracy and efficiency largely depend on a discrete representation of the model geometry that tracks the deformation of continuum bodies. On the one hand, a suitable mesh size could represent the real geometry of the research object very well. Furthermore, the fine mesh distribution is conducive to presenting a high resolution of local cracks. Hence, the mesh refinement procedure was performed at a loading rate of 0.005 m/s, and five more numerical cases with finer meshes are shown in Figure6. All the peridynamic models are built based on the butterfly FEM mesh. For the purpose of comparison, case M1 is added to Figure 6. The finest mesh model contains 1,334,000 nodes. The convergence of the peak force with increasing mesh or node numbers is shown in Figure 7. Clearly, the numerical results of the peridynamic model are mesh-dependent. When the number of nodes increases from 86,016 to 1,334,000, the peak force increases from 7019.62 N to 12045.65 N. Meanwhile, the rate of increase in the peak force decreases when the number of nodes exceeds 451,584, which demonstrates that the peak force is reaching convergence. Compared with that in case M5, the number of nodes in case M6 is over two times that in M5, but the peak force in case M6 increases by only 334.88 N (showing less than 3% variance). Therefore, both solutions in the cases of M5 and M6 are approximate to the convergence solution, and the effect of mesh size is negligible. In the subsequent numerical verification section, numerous trial-and-error methods are conducted to determine the fracture energy release rate. Taking the computing efficiency into account, case M5 was selected for subsequent analysis.Figure 6
Mesh independence analysis.Figure 7
Convergence of peak force with a finer mesh.


## 3.1. Laboratory Test of Marble BTS

Marble cored from western China was selected for the experimental test, and this marble is composed of calcite 94.8%, mica 3.2%, and dolomite 2.0%. The density of the marble is 2698.11 kg/m3. The line loading setup is shown in Figure 2. The rock specimens were disc samples with a diameter and thickness of 50 mm and 25 mm, respectively [35]. Due to the dimension errors that arose during coring, the real thickness and diameter of the disc sample were 48.73 mm and 25.42 mm, respectively. The tested dimensions were used in the subsequent strength calculation. The surfaces of a disc specimen were polished to avoid any irregularities in the thickness of the sample. In this study, the water content was not considered. The rock specimen was dried in a 105°C oven for 24 hours. The axial load was applied by the upper loading platen under the displacement control mode at a speed of 0.01 mm/min by using a hydraulic material testing machine. The load was applied continuously until the rock specimen failed, and the maximum load at failure was recorded. The indirect tensile strength was derived in accordance with
(7)BTS=2PπDt,where P, D, and t denote the applied axial load at failure, test specimen diameter, and test specimen thickness, respectively.Figure 2
Schematic diagram of the BTS test in the laboratory.Scanning electron microscopy (SEM) is widely used to analyze surfaces and can provide a highly magnified high-resolution image of the surface of the sample studied. A TESCAN Mira3 LMU scanning electron microscope (TESCAN, Czech Republic) was used to observe the microstructure of the marble. As Figure3 shows, the rock sample is relatively homogeneous, and the structure is compact. Marble was selected for this testing because it has a low intrinsic anisotropy. Thus, it is acceptable to simplify marble as an isometric material in numerical simulations.Figure 3
Microimage of marble.


## 3.2. Numerical Model Configuration

As Figure4 shows, the numerical model (case no. M1) configuration was consistent with the model depicted in Figure 2, in which the diameter and thickness of a disc specimen were 50 mm and 25 mm, respectively. The disc specimen was sandwiched by two flat platens, and the deformation of the two platens was not considered. To save computing time, BTS modeling was performed using a hybrid model that couples the FEM and a peridynamic model. The interaction between FEM and the peridynamic model was achieved by a contact algorithm. Two flat platens were established in the FEM mesh, and a rock model is built with peridynamic theory. In this hybrid model, the FEM part was built with the conventional method. Bond-based peridynamics was implemented into the code in the FEM framework with discontinuous Galerkin theory. Technically, the peridynamic-based part of a model is built by converting the FEM model using a detaching function. The top loading platen was applied with a constant downward velocity boundary condition until the maximum axial load was recorded in the results.Figure 4
Numerical model configuration.


## 3.3. Effect of Loading Rate

In the laboratory, a conventional static BTS test of rock may take several minutes to complete. If the loading rate remains consistent with the laboratory test in numerical simulation with an explicit solver, the long computing time is unacceptable. Thus, according to the rock response to the loading rate, a quasistatic loading method can be introduced for static problem simulation, such as unconfined compressive test modeling, BTS modeling, and fracture toughness modeling [36, 37]. The principle of the quasistatic loading method is to enhance the loading rate in the numerical simulation to reach a balance between efficiency of the computing process and accuracy of the numerical results. To examine the effect of the loading rate on the BTS modeling results in a peridynamic model, 6 sets of BTS modeling schemes with different loading rates were carried out using a coarse mesh model, as shown in Figure 5. The two rigid platens contained 16200 hexahedral meshes (25116 nodes), and a peridynamic domain (86016 nodes) was built based on 10752 FEM meshes distributed in the form of a butterfly. The corresponding loading rates were 0.5 m/s, 0.1 m/s, 0.05 m/s, 0.01 m/s, 0.005 m/s, and 0.001 m/s, and the lower rigid platen was fixed.Figure 5
The convergence of the force under different loading rates.Figure5 shows that the loading rate plays a significant role in the force response. In the laboratory test, marble shows apparent brittle properties, which are characterized by the force level dropping dramatically after the peak. At a loading rate of 0.5 m/s, the rock response to the high-speed loading tends to be ductile, which deviates from the essential character of the marble. With the loading rate decreasing, the curves of force versus displacement become closer, which indicates that the influences of these loading rates on the force are converging. According to the magnified local view in Figure 5, the three curves are quite similar when the loading rate varies from 0.001 to 0.01 m/s. In particular, the difference in the peak force at loading rates of 0.001 m/s and 0.005 m/s is only 85.39 N. Obviously, the influence of the loading rate on the force is negligible when the loading rate is less than 0.005 m/s. Therefore, a loading rate of 0.005 m/s is determined for subsequent numerical simulation under the quasistatic loading condition.


## 3.4. Mesh Independence

For a typical BTS test or modeling, the expected result is that the disc specimen is split into two approximately identical halves with gentle shear damage at the points of contact between the loading platens and rock specimen. Therefore, it is necessary that a sufficiently fine mesh should be distributed along the loading diameter. It is known that tensile cracks initiate from the center of the disc specimen and propagate along the loading diameter in a successful BTS test. In a specific peridynamic model, the butterfly mesh strategy provides the finest mesh in the center of the disc specimen, which ensures the accuracy of initial tensile crack generation.In the mesh-based numerical approach, its accuracy and efficiency largely depend on a discrete representation of the model geometry that tracks the deformation of continuum bodies. On the one hand, a suitable mesh size could represent the real geometry of the research object very well. Furthermore, the fine mesh distribution is conducive to presenting a high resolution of local cracks. Hence, the mesh refinement procedure was performed at a loading rate of 0.005 m/s, and five more numerical cases with finer meshes are shown in Figure6. All the peridynamic models are built based on the butterfly FEM mesh. For the purpose of comparison, case M1 is added to Figure 6. The finest mesh model contains 1,334,000 nodes. The convergence of the peak force with increasing mesh or node numbers is shown in Figure 7. Clearly, the numerical results of the peridynamic model are mesh-dependent. When the number of nodes increases from 86,016 to 1,334,000, the peak force increases from 7019.62 N to 12045.65 N. Meanwhile, the rate of increase in the peak force decreases when the number of nodes exceeds 451,584, which demonstrates that the peak force is reaching convergence. Compared with that in case M5, the number of nodes in case M6 is over two times that in M5, but the peak force in case M6 increases by only 334.88 N (showing less than 3% variance). Therefore, both solutions in the cases of M5 and M6 are approximate to the convergence solution, and the effect of mesh size is negligible. In the subsequent numerical verification section, numerous trial-and-error methods are conducted to determine the fracture energy release rate. Taking the computing efficiency into account, case M5 was selected for subsequent analysis.Figure 6
Mesh independence analysis.Figure 7
Convergence of peak force with a finer mesh.


## 4. Numerical Results and Analysis

An appropriate determination of the model parameters is essential for accurate simulation of the indirect tensile process of marble using peridynamic theory. After determining the key numerical parameters (loading rate and mesh size) that affect the accuracy of the numerical results, calibration of the fracture energy release rate is necessary to capture the target BTS of marble. In a typical BTS test, the fracture mode of a rock specimen is a mixed mode of slight shear and tensile failure. Thus, the fracture energy release rate in BTS modeling cannot be directly determined by experimental results. Therefore, the method of trial and error was used to obtain the appropriate fracture energy release rate of the marble. The numerical result of case M5 are quite close to the target. Thus, the real fracture energy release rate of marble was obtained by slightly reducing the fracture energy release rate. It is known that Poisson’s ratio being fixed at 0.25 constitutes the limitation of bond-based peridynamic theory [38, 39]. Therefore, to accurately implement the numerical simulation of the rock material, Poisson’s ratio of the simulated object should be approximately 0.25. The marble used in this study meets the needs of peridynamic simulation. In the previous sections, the material elastic modulus and fracture energy release rate were selected to be 30 GPa and 50 J/m2, respectively. With the method of variable control and trial and error, the final elastic modulus and fracture energy release rate were determined to be 15 GPa and 46 J/m2, respectively.Figure8 shows the comparison of the force versus displacement curves between the experiment and simulation, with the experimental and numerical peak forces being 10.00 kN and 9.99 kN, respectively, which indicates that the numerical result quantitatively matches the experimental results very well. In the experimental results, the rock specimen was split into two approximate halves along the loading diameter, and slight shear damage was observed in the interactive zone between the loading platen and specimen. Figure 9 shows the comparison of the failure mode between the experiment and simulation. Accordingly, in the numerical results, the tensile cracks initiated from the center of the disc specimen are consistent with the requirement of a successful BTS test. Furthermore, a triangular shear zone in the interaction zone between the loading platen and specimen is observed. All the phenomena demonstrate that the numerical experiment that captured the mechanical properties of the marble was successful. Therefore, the peridynamic simulation results of the BTS test shown in Figures 8 and 9 are consistent with the experimental results.Figure 8
Comparison of the force versus displacement curve between the experiment and simulation.Figure 9
Comparison of fracturing mode between experiment and simulation (plastic strain).


## 5. Conclusion

From the view of computing cost savings and the demand for rock fracturing behavior monitoring, this study presented a hybrid model that coupled a FEM model and a peridynamic model for indirect tensile simulation of marble. The BTS modeling result indicates that the hybrid model is mesh- and loading rate-dependent. On the basis of the convergences of force with loading rate and peak force with mesh fineness, the calibration of the fracture energy release rate of marble was conducted. The numerical results indicate that the crack initiated from the center of the disc model and propagated along the loading diameter, which is consistent with the experimental results. The hybrid numerical model provides an efficient approach for a rock fracturing behavior study. However, the bond-based peridynamic model requires a fixed Poisson’s ratio of 0.25 in 3D modeling, which may affect the deformation of rock subjected to an external load.



---

*Source: 8232332-2022-04-21.xml*