# Modelling the Crack Initiation and Propagation of Notched Rock Beam under the Three-Point Bending Load

Authors: Xun Xi; Xu Wu; Qifeng Guo; Meifeng Cai
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Prediction of rock fracture is essential to understand the rock failure mechanism. The three-point bending test has been one of the most popular experiments for the determination of rock fracture parameters. However, the crack initiation and propagation of rock beam with the center notch and offset notch have not been fully understood. This paper develops a numerical method for modelling the notched beam cracking based on nonlocal extended finite element method (i.e., XFEM) and mixed mode rock fracture model. An example is worked out to demonstrate the application of the numerical method and verified with experimental results. The crack length development, crack pattern, crack opening and slipping displacements, and the load-crack mouth of displacement (P-CMOD) curve are obtained. The effects of offset notch location and mechanical properties on the crack length development, P-CMOD curve, and crack pattern are investigated and discussed. It has been found that the peak load of the notched beam nearly linearly increases with the increase of the notch offset ratio. The cracking of rock beam with offset notch is dominated by mode I fracture, but mode II fracture contributes more when crack deflection occurs. The fracture energy significantly affects the peak load, while it has little effect on the prepeak and postpeak slopes in the P-CMOD curve.


---

## Body

## 1. Introduction

The failure of rock or rock mass in deep rock engineering is closely related to the crack initiation, propagation, coalescences, and nucleation from micro- to macroscales [1–3]. The energy storage and dissipation during rock cracking significantly affect deep rock failure behaviour [4]. Prediction of rock fracture is important to understand the rock mechanical behaviours and rock failure mechanism [5]. Due to the difficulties in controlling the rock cracking for direct tensile tests, some indirect tensile tests (e.g., Brazilian test, three-point bending beam test, and semicircular bending test) have been widely used to obtain the tensile failure properties of rock [4, 6–8]. Through controlling the crack mouth of displacement (CMOD) of the notch, the whole load-CMOD curve including the postpeak softening part can be obtained by three-point bending beam tests [8]. Further, an improved understanding of rock fracture can be achieved through investigating the crack initiation and propagation of rock beam with the center notch or offset notch under the three-point bending load.As a standard or suggested method for determining the strength, fracture energy, and fracture toughness of rock and concrete from RILEM (International Union of Laboratories and Experts in Construction Materials, Systems, and Structures) [8], ASTM (American Society for Testing and Materials) [9], and ISRM (International Society for Rock Mechanics) [10], the central notched beams have been widely employed to investigate the fracture properties of rock and concrete [2, 7, 11–13]. Gong et al. [4] investigated the storage and dissipation evolution process during rock tensile failure based on the Brazilian test, point load test, and semicircular bending test. They found that there were strongly linear relationships between the elastic energy, dissipated energy, and total input energy [4]. Fakhimi et al. [2, 12] compared the characterization of fracture process zone around the crack tip by granite beams with single center notch for four different sizes and found that the length and width of the process zone increased with the increase of the specimen size. Dong et al. [6] investigated the fracture properties of rock-concrete interfaces by making a notch at the geometric center of the rock-concrete beam and found the fracture energy of the interface was much smaller than that for concrete or rock. Furthermore, the fracture of rock beam with offset notches attracted more attention for understanding the nominal strength and crack patterns for offset notched rock under the three-point bending load [14–17]. Zuo et al. [14] employed the scanning electron microscope to investigate the effect of the offset notch on the cracking of basalt beam with dimensions of 25 mm × 10 mm × 5 mm. They found that the distance between the offset notch and the beam centerline significantly affected the crack patterns [14]. Furthermore, Zuo et al. [18] investigated the effects of shale bedding on crack behaviour and fracturing mechanism by scanning electron microscope (SEM) with a loading system and found that the peak load, fracture toughness, and fracture energy all increased with the bedding angle from 0° to 90°. Luo and Gong [19] carried out a series of single cyclic loading-unloading flexural fracture tests on marble, red sandstone, and granite and proposed flexural energy storage coefficient and flexural energy dissipation coefficient to represent rock capacities for storing energy and dissipating energy, respectively. Guo et al. [15] made granite beams with single center notch and offset notch with dimensions of 145 mm × 50 mm × 18 mm and carried out three-point bending tests. They found that the larger the offset ratio was, the larger the peak load and nominal strength of the specimen were [15]. Xu and Cao [17] carried out three-point bending tests by rock-like material beams with single offset notch at different locations and found that the peak load increased linearly with the increase of the offset ratio. However, it is not easy to make an intact large-scale rock beam without preexisting cracks because lots of natural fractures or joints exist in rock mass. Moreover, experimental observations of crack initiation and propagation of rock beam, especially for the postpeak softening stages, are time-consuming, technically challenging, and expensive. Therefore, the effect of offset notch on the crack initiation and propagation of rock beam has not been fully understood.Computational simulations of rock fracture have been an active research field, which enhanced the understanding of failure mechanism of rock [16–18, 20–30]. Fakhimi et al. [2] employed a bonded particle model to investigate the fracture process of rock beam with a center notch and found that a cohesionless crack can develop at peak load. Xu and Cao [17] modelled the cracking of cemented tailings backfill beam with offset notches and found that the offset ratio affected the crack patterns of the beam by PFC software. Huang et al. [28] employed cohesive elements to model the cracking of layered rock under the semicircular bending test and found that the fracture load of sandstone increased with the increase of inclination angle between loading direction and bedding plane. Wei et al. [30] modelled the cracking of rock with multiple flaws by continuum-based discrete element method. Dong et al. [23] modelled the interface cracking of the rock-concrete composite beam with offset notch by remeshing finite element method and found that shear fracture contributed to the fracture process of beam. Jia and Zhang [26] modelled the three-point bending beam tests by the damage model in RFPA (rock failure process analysis) software and found that the location of the offset notch affected the failure mode of the rock beam. Among the existing pieces of literature on numerical modelling of three-point bending tests of rock beam, most were focused on the beam with the central notch and others discussed more about the crack patterns and failure mode for the beam with offset notches. It is of significance to predict rock cracking and understand rock fracture behaviour for notched three-point beam tests. The extended finite element method (XFEM) provided a very helpful tool for modelling rock arbitrary cracking without the requirement for remeshing [31]. Moreover, the nonlocal stress averaging for cracking can make the prediction of materials cracking more accurate [32]. Therefore, it is worth modelling the cracking of notched rock beam under the three-point bending load by the nonlocal XFEM and mixed mode fracture model.This paper attempts to develop a numerical method based on nonlocal XFEM for modelling notched rock beam cracking under the three-point bending load. Cohesive crack behaviour is used to describe the mixed mode fracture of rock. An example of rock beams with the notch at different locations is worked out to demonstrate the numerical method. The crack length development, crack pattern, crack opening and slipping displacements, and load-CMOD (P-CMOD) curve are obtained. The numerical results are compared with those from experiments. Finally, the effects of offset notch location and mechanical parameters on the crack initiation, prorogation, crack pattern, and P-CMOD curves are investigated and discussed.


## 2. Constitutive Model for Rock Fracture

Rock exhibits the tensile strain-softening behaviour due to an inelastic zone being developed ahead of the crack tip, often referred to as fracture process zone (FPZ) [7, 25]. The cohesive crack model first proposed by Jia and Zhang [26] has been employed to simulate discrete cracking in the fracture process zone of rock and concrete. Before crack initiation, the stress of rock linearly develops as a linear function:(1)σn=Knδn,σs=Ksδs,where σn and σs are the normal and shear stresses, respectively; δn and δs are corresponding normal and shear strains, respectively; Kn and Ks are the normal and tangential stiffness, respectively.When the maximum principal stress reaches the criterion value (i.e., cohesive strength), crack initiation will occur and a damage value is introduced to reduce the stiffness for stress softening, i.e.,(2)σn=1−DKnδn,σs=1−DKsδs,where D is the overall damage variable which is 0 before damage initiation and 1 after complete failure. As shown in Figure 1, the damage evolution between initiation of damage and final failure can follow a linear or exponential function [33, 34]. The areas under the normal stress-displacement curve and shear stress-displacement curve are mode I and mode II fracture energy, respectively.Figure 1
Constitutive model for rock fracture.
(a)(b)The mixed mode fracture is considered and the mixed mode fracture energy is defined as follows [34]:(3)Gc=1mI/GI+mII/GII,where GC is the mixed mode fracture energy; GI and GII are the pure mode I and pure mode II fracture energies, respectively. mI and mII are the work ratios done by the normal and shear forces, respectively.The relationship between the fracture energy and the fracture toughness can be established by Irwin’s formula [35]:(4)KC=GCE1−ν2,where KC is the fracture toughness; E is the elastic modulus; and ν is Poisson’s ratio.Furthermore, the damage value can be calculated by the following equation:(5)D=δm,fδm,max−δm,0δm,maxδm,f−δm,0,where δm,max is the maximum effective relative displacement during the loading history; δm,0 is the critical effective relative displacement when the damage starts; δm,f is the effective relative displacement when complete failure occurs; δm,f can be calculated by the mixed mode fracture energy as follows:(6)δm,f=2GcKδm,0.Once the damage value is determined, the residual stresses can be obtained. By the above equations, the stress-strain relationship of rock under mixed mode fracture is established.


## 3. XFEM Simulation

### 3.1. Basic Principles of XFEM

To overcome the problem associated with matching the geometry of the discontinuity as the crack propagation, the extended finite element method was first introduced by Belytschko and Black [36]. The presence of discontinuities is ensured by the special enriched functions in conjunction with additional degrees of freedom while the continuous displacements are derived from the traditional FEM with the retained sparsity and symmetry. As shown in Figure 2, in XFEM, the Heaviside enrichment function and the crack tip enrichment functions are introduced to represent discontinuities and crack tip fields, respectively. Furthermore, the total displacement can be expressed as follows [37, 38]:(7)u=∑I=1NNIxuIcon+HxaI+∑α=14FαxbIα,where µ is the displacement in the computational domain; N is the number of the Gauss integral points in the domain; NIx is the continuous nodal shape functions; µIcon is the continuous nodal displacement for the traditional finite element solution; Hx is the Heaviside function for achieving the displacement jump across the crack surface; aI is the nodal enriched degree of freedom; Fαx is the asymptotic crack tip function; bIα is the nodal enriched degree of freedom.Figure 2
Illustration of normal and tangential coordinates for the crack.The Heaviside function is expressed as follows:(8)H=1,ifx−x∗⋅n>0,−1,otherwise,where x is a Gauss point; x∗ is the closest point to x on the crack face; n is the unit outward normal to the crack at x∗.In a polar coordinate system (r, θ), the asymptotic crack tip function is expressed as follows:(9)Fα=rsinθ2,rcosθ2,rsinθsinθ2,rsinθcosθ2.


### 3.2. Nonlocal Averaging of Stress Field near the Crack Tip

Once the crack initiation criteria are satisfied, a new crack will be created in the enriched element. The newly introduced crack is always orthogonal to the maximum principal stress direction. However, the direction will be affected by the local element in the mesh. To reduce the mesh dependency and improve the accuracy of crack direction, a nonlocal calculation technique is used as illustrated in Figure3. The elements in the defined radius will be used to calculate the crack direction. Furthermore, a Gaussian function is used to define the weighting schemes for nonlocal averaging as follows [38]:(10)wr=12π3/2rc3e−r2/2rc2,where rc is the radius around the crack tip for averaging and r is the location of the element in the averaging region. Through this nonlocal and nonuniform calculation, a higher weighting is given to elements close to the crack tip.Figure 3
Illustration of nonlocal averaging of the stress around the crack tip.


## 3.1. Basic Principles of XFEM

To overcome the problem associated with matching the geometry of the discontinuity as the crack propagation, the extended finite element method was first introduced by Belytschko and Black [36]. The presence of discontinuities is ensured by the special enriched functions in conjunction with additional degrees of freedom while the continuous displacements are derived from the traditional FEM with the retained sparsity and symmetry. As shown in Figure 2, in XFEM, the Heaviside enrichment function and the crack tip enrichment functions are introduced to represent discontinuities and crack tip fields, respectively. Furthermore, the total displacement can be expressed as follows [37, 38]:(7)u=∑I=1NNIxuIcon+HxaI+∑α=14FαxbIα,where µ is the displacement in the computational domain; N is the number of the Gauss integral points in the domain; NIx is the continuous nodal shape functions; µIcon is the continuous nodal displacement for the traditional finite element solution; Hx is the Heaviside function for achieving the displacement jump across the crack surface; aI is the nodal enriched degree of freedom; Fαx is the asymptotic crack tip function; bIα is the nodal enriched degree of freedom.Figure 2
Illustration of normal and tangential coordinates for the crack.The Heaviside function is expressed as follows:(8)H=1,ifx−x∗⋅n>0,−1,otherwise,where x is a Gauss point; x∗ is the closest point to x on the crack face; n is the unit outward normal to the crack at x∗.In a polar coordinate system (r, θ), the asymptotic crack tip function is expressed as follows:(9)Fα=rsinθ2,rcosθ2,rsinθsinθ2,rsinθcosθ2.


## 3.2. Nonlocal Averaging of Stress Field near the Crack Tip

Once the crack initiation criteria are satisfied, a new crack will be created in the enriched element. The newly introduced crack is always orthogonal to the maximum principal stress direction. However, the direction will be affected by the local element in the mesh. To reduce the mesh dependency and improve the accuracy of crack direction, a nonlocal calculation technique is used as illustrated in Figure3. The elements in the defined radius will be used to calculate the crack direction. Furthermore, a Gaussian function is used to define the weighting schemes for nonlocal averaging as follows [38]:(10)wr=12π3/2rc3e−r2/2rc2,where rc is the radius around the crack tip for averaging and r is the location of the element in the averaging region. Through this nonlocal and nonuniform calculation, a higher weighting is given to elements close to the crack tip.Figure 3
Illustration of nonlocal averaging of the stress around the crack tip.


## 4. Worked Example and Verification

A rock beam model with single notch at different locations is used to demonstrate the application of the developed numerical method. As shown in Figure4, the dimensions of the beam are 145 mm in length, 50 mm in height, and 18 mm in width. A notch of 10 mm in height is created in the rock beam. The span length is 127 mm. The geometric parameters of the beam keep the same as the granite beam from the experiments [15]. The offset ratio for describing the location of offset notch is expressed as follows:(11)r=2cS,where r is the offset ratio; c is the distance between the notch and the centerline of the beam; S is the length of the span. Five cases for the beam with offset ratios of the notch 0, 10%, 20%, 30%, and 40% are considered in this paper.Figure 4
Geometry of the worked example.Figure5 shows the numerical model for the beam with the notch of 20% offset ratio. The bending loads are applied by the rigid bodies which frictionally contact with the rock beam. There are 19,404 elements in the model and the size of element in the interest area is 0.5 mm which is sufficiently fine. By moving the location of the preexisting notch, the numerical models for the five cases with different offset ratio can be obtained. Therefore, the mesh for the rock beam with the center and offset notch is the same, which avoids the effect of meshing on rock cracking. The mode I and mode II fracture energies are set at the same value. All the basic parameters and their sources in the model are listed in Table 1.Figure 5
Numerical model for the beam with the notch of 20% offset ratio: (a) the whole model and (b) enlarged part around the notch.Table 1
The basic parameters for the numerical simulations.
SymbolsDescriptionValueEYoung’s modulus30 GPa [15]vPoisson’s ratio0.25 [15]σmaxCohesive strength9 MPa [15]GIMode I fracture energy100 N/m [15]GIIMode II fracture energy100 N/m [15]rcRadius of nonlocal averaging0.5 mmFigure6 illustrates the cracking process and the maximum principal stress of the rock beam with 20% offset ratio notch with the loading displacement increasing. It can be seen that, with the loading displacement increasing, the stress firstly concentrates around the crack tip. After reaching the maximum allowable principal stress, a new crack is initialized and propagated towards the top of the beam. The crack is inclined to the loading point from the notch while it becomes straight at the final stage. Figure 7 shows the opening and slipping displacements of the final crack. It can be found that the opening displacement closer to preexisting notch is larger. However, the large values for the slipping displacement discretely occur in the middle part of the inclined crack. Moreover, the slipping displacement is much smaller than the opening displacement. Therefore, the cracking of rock beam with offset notch is dominated by mode I fracture but mode II fracture contributes more when crack deflection occurs. And the cracking of rock beam with offset notch transfers from mixed mode fracture to mode I fracture.Figure 6
The crack processes and maximum principal stress under the loading displacement: (a) 0.018 mm; (b) 0.036 mm; (c) 0.073 mm; (d) 0.38 mm.
(a)(b)(c)(d)Figure 7
The displacements of the final crack: (a) opening displacement; (b) slipping displacement.
(a)(b)To investigate the effect of offset notch location on the cracking, numerical simulations for the beams with offset notch ratios 0, 10%, 20%, 30%, and 40% are carried out. Figure8 illustrates the crack length development with the bending load (displacement) increasing. The crack length is obtained by calculating the overall length of the damaged elements, which is achieved by an in-house Python script. It can be seen that the loading displacement to crack initially increasing is almost the same for all the beams. However, the smaller the notch offset ratio is, the faster the crack propagation is. The larger the notch offset ratio is, the larger the final crack length is. Figure 9 shows the final crack patterns for the rock beams with different notch of offset ratios. The crack patterns are obtained from every single numerical model and then combined together. It can be seen that, for every offset notched beam, the new crack obliquely propagates from the preexisting notch to the loading point and then deflects to a straight crack. The larger the offset ratio, the farther the final crack away from the loading point, which has a good agreement with crack patterns from previous experiments [14].Figure 8
The crack length developments with the loading displacement increasing for different notch offset ratios.Figure 9
Crack patterns for the rock beams with different notch offset ratios.The load-CMOD curves from the numerical simulations are obtained as Figure10. The values of CMOD are calculated by the enriched nodes at the bottom of notch. It can be found that, with the increase of crack mouth of displacement, the load gradually increases to the peak and then decreases to zero. The postpeak curves are in bilinear shapes. The larger the offset ratio is, the greater the prepeak slope is and the larger the peak load is.Figure 10
P-CMOD curves for the rock beams with different notch offset ratios.Figure11 illustrates the comparison of peak loads from the numerical simulations and experiments [15]. It can be seen that the peak load is close to linearly increase with the increase of the notch offset ratio for both experimental and numerical results. The numerical results have a very good agreement with those from experiments [15].Figure 11
Experimental verification of the peak load.


## 5. Parametric Studies and Discussion

The effects of Young’s modulus on the crack length development and P-CMOD curve are investigated. All the other parameters keep the same as listed in Table1. The notch offset ratio of the rock beam is 20% for the parametric studies. Figure 12 shows the crack length developments with the loading displacement increasing for the different values of Young’s modulus. It can be found that the larger Young’s modulus is, the smaller the loading displacement to crack initiation is and the faster the crack propagation is. The final crack length is the same for all the models. Therefore, Young’s modulus has little effect on the final crack length. Figure 13 illustrates the P-CMOD curves for the rock beams with different Young’s modules. It can be seen that the larger Young’s modulus is, the greater the prepeak slope is and the larger the peak load is. It is interesting to find that the smaller Young’s modulus is, the slower the decreasing stages are. Therefore, Young’s modulus significantly affects the crack propagation, peak load, prepeak slopes, and postpeak slopes of the P-CMOD curves. It is necessary to accurately determine Young’s modulus of rock for prediction of the rock fracture.Figure 12
Effect of Young’s modulus on the crack length development.Figure 13
Effect of Young’s modulus on the P-CMOD curve.Figure14 illustrates the crack length development for different cohesive strengths. It can be found that the larger the strength is, the larger the loading displacement to crack initiation is. The strength has little effect on the final crack length. Figure 15 shows the P-CMOD curves for the rock beam with different cohesive strength. It can be seen that the larger the strength is, the larger the peak load is. However, the strength has little effect on the prepeak slope and postpeak slope in the P-CMOD curves.Figure 14
Effect of cohesive strength on the crack length development.Figure 15
Effect of cohesive strength on the P-CMOD curve.Figure16 illustrates the crack length development for the rock beam with different values of fracture energy. It can be seen that the loading displacements to crack initiation are almost the same. The fracture energy has very little effect on the final crack length. However, the crack propagation speed is slightly smaller for larger fracture energy. Figure 17 shows P-CMOD curves for the rock beam with different values of fracture energy. It can be seen that the fracture energy significantly affects the peak load. The larger the fracture energy is, the larger the peak load is. However, the fracture energy has little effect on the increasing and postpeak slope in the P-CMOD curves.Figure 16
Effect of fracture energy on the crack length development.Figure 17
Effect of fracture energy on the P-CMOD curve.Figure18 shows the final crack patterns for the rock beams discussed previously. It can be seen that the material mechanical parameters, i.e., Young’s modulus, strength, and fracture energy, have almost no effect on the crack pattern.Figure 18
Crack patterns for the rock beams with the notch of 20% offset ratio for different mechanical parameters: (a)E = 20 GPa; σmax = 9 MPa; Gf = 100 N/m; (b) E = 30 GPa; σmax = 9 MPa; Gf = 100 N/m; (c) E = 40 GPa; σmax = 9 MPa; Gf = 100 N/m; (d) E = 30 GPa; σmax = 6 MPa; Gf = 100 N/m; (e) E = 30 GPa; σmax = 12 MPa; Gf = 100 N/m; (f) E = 30 GPa; σmax = 9 MPa; Gf = 80 N/m; (g) E = 30 GPa; σmax = 9 MPa; Gf = 120 N/m.
(a)(b)(c)(d)(e)(f)(g)


## 6. Conclusions

In this paper, a numerical method for notched rock beam cracking under the three-point bending load has been developed based on the nonlocal XFEM and mixed mode fracture model. The cohesive crack behaviour in the fracture process zone was employed to describe the rock cracking. A worked example for rock beams with notch offset ratio of 0, 10%, 20%, 30%, and 40% has been presented to demonstrate the application of the derived method and then verified with experimental results. The effects of the offset ratio, Young’s modulus, strength, and fracture energy on the crack length development, crack pattern, and P-CMOD curves were investigated and discussed. Conclusions can be given as follows:(1)
The derived method based on nonlocal XFEM and rock mixed mode fracture is suitable for modelling the crack initiation and propagation of notched rock beam under the three-point bending load. The arbitrary crack is produced without the limitations of the mesh.(2)
The peak load of the notched rock beam is close to linearly increase with the increase of notch offset ratio. The numerical results have a good agreement with those from experiments.(3)
The cracking of rock beam with offset notch is dominated by mode I fracture but mode II fracture contributes more when crack deflection occurs. The cracking of rock beam with offset notch transfers from mixed mode fracture to mode I fracture.(4)
The material mechanical parameters, that is, Young’s modulus, strength, and fracture energy, have no effect on the crack pattern. The peak load of notched rock beam increases with the increase of Young’s modulus, strength, and fracture energy. The fracture energy has little effect on the prepeak and postpeak slopes in the P-CMOD curves.



---

*Source: 8810023-2020-09-07.xml*