# Research on the Behavior and Mechanism of Three-Dimensional Crack Growth under Uniaxial Loading

Authors: Zhibo Zhang; Shujie Li; Xuanye Qin
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

Most of the cracks in the rock masses are in a three-dimensional (3D) state, and it is always a hot topic to reveal the mechanical mechanism of 3D crack growth. In this paper, the research on the growth behavior of 3D crack is performed through laboratory experiments and numerical simulations. Cement samples with different angles of 3D crack are prepared, and the uniaxial compression experiment is carried out. The results indicate that initiation of preexisting crack with an angle of 45° is easier and shear failure characteristics of corresponding samples are obvious. Through theoretical analysis, the preexisting crack starts to grow at the end of the short axis, along the short axis end to the long axis end of the preexisting crack, the shear effect decreases gradually, and the tearing effect increases gradually. Combined with numerical simulation, the experimental and analysis results are verified, and the preexisting crack growth process is presented. The growth direction of the preexisting crack changes from perpendicular to the crack surface to parallel principal stress direction, and the maximum growth length can reach 1.2 times the minor axis radius of the preexisting crack. The research results can provide an important theoretical basis for revealing the evolution process of the cracks in rock masses.


---

## Body

## 1. Introduction

The rock is a natural and heterogenous material, possessing various defects inside and on the surface. Although these defects account for only a small part of the volume of rock, they have an important impact on the mechanical properties of the rock. When the rock is subjected to the external load, various cracks randomly distributed in the rocks will start to grow, resulting in the instability of the rock structure [1–4]. Therefore, research on the growth behavior and mechanical mechanism of the rock has important meaning for ensuring the stability of rock mass.Generally, the cracks in the rock masses are constrained by forces in three directions (3D state). However, most research simplifies the crack to 2D state due to the complexity of experimental observation and theoretical analysis of a 3D crack. On this basis, a lot of experiments and research work have been carried out, and many meaningful results have been obtained [5–12].Although the 2D crack model reduces the complexity of research, it is still different from the real 3D crack. Research on the mechanical behavior of 3D crack is an unavoidable topic. Therefore, scholars have carried out further research on the 3D crack. Zhou et al. and Yang et al. used the PMMA sample and mortar sample with the preexisting crack to conduct an experiment. The research results showed that geometric characteristics of the preexisting crack have an important effect on peak stress and failure behavior of the corresponding sample [13–15]. During the uniaxial compression process, Dyskin et al. found the growth started from the 3D crack tip to form the wrapping wing cracks, and the maximum length of the wing crack could reach 1–1.5 times the initial crack radius [16]. In order to observe the growth of 3D crack intuitively, Li et al. conducted research on 3D crack through CT scanning technology [17]. Further, some scholars observed the growth phenomenon of 3D crack during the biaxial loading process, and they found that lateral load had a significant influence on the growth of 3D crack [18–22]. Some numerical simulation results were also reported, which is consistent with the previous experimental results [23, 24]. Meanwhile, the effect of physical parameters of the 3D crack on its growth has been studied, including the number of the 3D cracks and geometry shape, and corresponding research also obtained many meaningful results [25–35].According to the above analysis, a lot of attempts were made to study the growth of the 3D crack. However, due to the difficulty in making the samples with 3D crack and observation, the research on the growth mechanism of crack with different inclination angles was relatively less. Therefore, in this paper, cement samples with different angles of 3D crack were made to conduct a uniaxial compression experiment. Combined with theoretical analysis and numerical simulation, the growth mechanism of crack with different inclination angles and the effect on the corresponding sample were revealed. The research results have important theoretical value and practical significance for ensuring the stability of rock mass.


## 2. Experiment Introduction

### 2.1. The Sample Preparation

In this experiment, the cement samples were cast using a unified plastic mold with a size of Φ50 mm × 100 mm. An aluminum foil with a thickness of 0.2 mm was processed into the elliptical shape for constructing the preexisting crack, whose semimajor axisa was 7.5 mm and semiminor axis b was 5 mm. The surface of the aluminum foil was coated with lubricating oil to reduce friction. Then, the cement with the water–cement ratio of 1 : 2 was poured into the mold slowly. Finally, the top of the molds was sealed with the plastic wrap and put into YH-40B standard constant temperature and humidity curing box for 28 days. In total, five types of samples with angles of 0°, 30°, 45°, 60°, and 90° were prepared, as shown in Figure 1.Figure 1
The samples preparation. (a) Making process of the samples. (b) Schematic diagram of the samples.
(a)(b)


### 2.2. Experimental System and Process

The experimental system consisted of a loading system, AE signal acquisition system, and shielding system, as shown in Figure2. The loading system was mainly YAW-600 electrohydraulic servo pressure testing machine. The DS5 full-waveform AE acquisition instrument was used to obtain the AE signal produced during the entire fracture process of the samples. The experiment was carried out in GP1A electromagnetic shielding room with a shielding effect being 75 dB, which can effectively reduce interference from the external environment.Figure 2
Schematic diagram of the experimental system.In this uniaxial loading experiment, a displacement control mode was adopted to apply load to the samples, and the loading rate was 0.5µm/s. Two AE transducers were arranged in the middle of the sample to collect the AE signal.


## 2.1. The Sample Preparation

In this experiment, the cement samples were cast using a unified plastic mold with a size of Φ50 mm × 100 mm. An aluminum foil with a thickness of 0.2 mm was processed into the elliptical shape for constructing the preexisting crack, whose semimajor axisa was 7.5 mm and semiminor axis b was 5 mm. The surface of the aluminum foil was coated with lubricating oil to reduce friction. Then, the cement with the water–cement ratio of 1 : 2 was poured into the mold slowly. Finally, the top of the molds was sealed with the plastic wrap and put into YH-40B standard constant temperature and humidity curing box for 28 days. In total, five types of samples with angles of 0°, 30°, 45°, 60°, and 90° were prepared, as shown in Figure 1.Figure 1
The samples preparation. (a) Making process of the samples. (b) Schematic diagram of the samples.
(a)(b)


## 2.2. Experimental System and Process

The experimental system consisted of a loading system, AE signal acquisition system, and shielding system, as shown in Figure2. The loading system was mainly YAW-600 electrohydraulic servo pressure testing machine. The DS5 full-waveform AE acquisition instrument was used to obtain the AE signal produced during the entire fracture process of the samples. The experiment was carried out in GP1A electromagnetic shielding room with a shielding effect being 75 dB, which can effectively reduce interference from the external environment.Figure 2
Schematic diagram of the experimental system.In this uniaxial loading experiment, a displacement control mode was adopted to apply load to the samples, and the loading rate was 0.5µm/s. Two AE transducers were arranged in the middle of the sample to collect the AE signal.


## 3. Experimental Results

### 3.1. Influence of Crack Angle on Crack Initiation Stress

Because it was impossible to directly measure the initiation stress of the internal 3D crack for cement samples, AE energy was used to determine indirectly initiation stress. According to Guo’s and Zhang’s research [36, 37], the first rapid increase of AE energy corresponds to the initiation process of preexisting crack. Based on this method, the initiation stress can be determined. Further, the initiation stress was normalized according to the following equation, and the calculation results are shown in Figure 3 and Table 1:(1)K=σciσmax×100%,where K represented initiation stress level and σmax and σci represented the peak stress and initiation stress of the preexisting crack, respectively.Figure 3
Evolution characteristics of AE energy with stress: (a)α = 0°; (b) α = 30°; (c) α = 45°; (d) α = 60°; (e) α = 90°.
(a)(b)(c)(d)(e)Table 1
Initiation stress and peak stress for different types of samples.
Crack angle (°)030456090Initiation stress (MPa)22.0412.539.8011.5512.14Peak stress (MPa)26.0117.2716.7814.7715.84Initiation stress level (%)8573587877According to Table1, the range of initiation stress was from 9.80 MPa to 22.04 MPa, and the initiation stress level ranged from 58% of the peak stress to 85% of the peak stress. The change trend of the initiation stress and initiation stress level was similar, and both presented a decrease-increase trend with the increase of the preexisting crack angle. When the angle of the preexisting crack was 45°, crack initiation stress and crack initiation stress level reached the minimum value, 9.80 MPa and 58%, respectively. This phenomenon indicated that initiation of the preexisting crack was easier when the angle of the preexisting crack was 45°, and the crack growth rate was relatively slow.


### 3.2. Influence of Crack Angle on Failure Mode of the Samples

According to Szwedzicki’s research [38], five failure modes of hard brittle cylindrical rock samples under uniaxial loading are defined, including simple extension, multiple extension, multiple fracturing, multiple shear, and simple shear, as shown in Figure 4.Figure 4
Classification of rock failure modes under uniaxial compression [38].Combined with observed results, the failure process of the sample with different angles of crack is described in Table2.Table 2
Description of the failure process of the sample with different angles of crack.
Angle (°)Failure processSchematic diagram of sample failureCrack-freeMost of the induced cracks grew along the direction parallel to loading, and failure mode was the simple extension mode.0Because the preexisting crack was parallel to the loading direction, induced cracks also grew along the loading direction. Compared with the failure mode of the crack-free sample, the failure mode for this type of samples was more complex, presenting a multiple extension mode.30With the angle of preexisting crack increasing, shear failure is gradually strengthened and the failure mode was the multiple shear mode.45When the angle of preexisting crack reached 45°, sample failure mode presented obvious shear failure characteristics (simple shear mode).60As the angle of the preexisting crack continued to increase, shear failure characteristics gradually weakened and there was a significant increase in the number of the induced cracks. The failure mode was similar to the multiple fracturing mode.90When the angle of the preexisting crack reached 90°, the preexisting crack was perpendicular to the loading direction. It is challenging to form stress concentration at the crack tip, and shear failure characteristics were further weakened. However, its failure mode was different from that of the crack-free sample (simple extension mode) due to the preexisting crack, presenting a multiple fracturing mode.


## 3.1. Influence of Crack Angle on Crack Initiation Stress

Because it was impossible to directly measure the initiation stress of the internal 3D crack for cement samples, AE energy was used to determine indirectly initiation stress. According to Guo’s and Zhang’s research [36, 37], the first rapid increase of AE energy corresponds to the initiation process of preexisting crack. Based on this method, the initiation stress can be determined. Further, the initiation stress was normalized according to the following equation, and the calculation results are shown in Figure 3 and Table 1:(1)K=σciσmax×100%,where K represented initiation stress level and σmax and σci represented the peak stress and initiation stress of the preexisting crack, respectively.Figure 3
Evolution characteristics of AE energy with stress: (a)α = 0°; (b) α = 30°; (c) α = 45°; (d) α = 60°; (e) α = 90°.
(a)(b)(c)(d)(e)Table 1
Initiation stress and peak stress for different types of samples.
Crack angle (°)030456090Initiation stress (MPa)22.0412.539.8011.5512.14Peak stress (MPa)26.0117.2716.7814.7715.84Initiation stress level (%)8573587877According to Table1, the range of initiation stress was from 9.80 MPa to 22.04 MPa, and the initiation stress level ranged from 58% of the peak stress to 85% of the peak stress. The change trend of the initiation stress and initiation stress level was similar, and both presented a decrease-increase trend with the increase of the preexisting crack angle. When the angle of the preexisting crack was 45°, crack initiation stress and crack initiation stress level reached the minimum value, 9.80 MPa and 58%, respectively. This phenomenon indicated that initiation of the preexisting crack was easier when the angle of the preexisting crack was 45°, and the crack growth rate was relatively slow.


## 3.2. Influence of Crack Angle on Failure Mode of the Samples

According to Szwedzicki’s research [38], five failure modes of hard brittle cylindrical rock samples under uniaxial loading are defined, including simple extension, multiple extension, multiple fracturing, multiple shear, and simple shear, as shown in Figure 4.Figure 4
Classification of rock failure modes under uniaxial compression [38].Combined with observed results, the failure process of the sample with different angles of crack is described in Table2.Table 2
Description of the failure process of the sample with different angles of crack.
Angle (°)Failure processSchematic diagram of sample failureCrack-freeMost of the induced cracks grew along the direction parallel to loading, and failure mode was the simple extension mode.0Because the preexisting crack was parallel to the loading direction, induced cracks also grew along the loading direction. Compared with the failure mode of the crack-free sample, the failure mode for this type of samples was more complex, presenting a multiple extension mode.30With the angle of preexisting crack increasing, shear failure is gradually strengthened and the failure mode was the multiple shear mode.45When the angle of preexisting crack reached 45°, sample failure mode presented obvious shear failure characteristics (simple shear mode).60As the angle of the preexisting crack continued to increase, shear failure characteristics gradually weakened and there was a significant increase in the number of the induced cracks. The failure mode was similar to the multiple fracturing mode.90When the angle of the preexisting crack reached 90°, the preexisting crack was perpendicular to the loading direction. It is challenging to form stress concentration at the crack tip, and shear failure characteristics were further weakened. However, its failure mode was different from that of the crack-free sample (simple extension mode) due to the preexisting crack, presenting a multiple fracturing mode.


## 4. Initiation Mechanism of Crack Tip

### 4.1. Mechanical Mechanism of 3D Crack Initiation

According to the theory of the maximum tensile stress, the crack will start to grow from the front edge point. For the elliptical crack used in this experiment, the growth of the elliptical crack can be considered as a mixed mode of fracture mode II and fracture mode III, and the corresponding 3D coordinate can be shown in Figure5 [39].Figure 5
The 3D coordinate diagram of elliptic crack.According to the elastic mechanics oblique section formula [40], the equivalent normal stress σe and equivalent shear stress τe acting on the preexisting crack surface can be calculated by using the following equations:(2)σe=σsin2α,(3)τe=σsinαcosα−μσsin2α,where σ is the far-field compressive stress, α is the angle of elliptic crack to the loading direction, and μ is friction coefficient between crack surfaces.Combined with Figure5, define point O as the center on the front edge of the preexisting crack in the local x-y-z coordinate system, and the following equation can be used to calculate the 3D stress component at the end of the crack with radius r and fracture bending angle θ:(4)σxx=−KII2πrsinθ22+cosθ2cos3θ2,σyy=KII2πrsinθ2cosθ2cos3θ2,τxy=KII2πrcosθ21−sinθ2sin3θ2,τxz=−KIII2πrsinθ2,τyz=KIII2πrcosθ2,σzz=−2νKIIπrsinθ2,where KII and KIII are the stress intensity factors of fracture mode II and fracture mode III, respectively. Parameter ν is Poisson’s ratio.The analytical solutions of stress intensity factors for fracture mode II and fracture mode III can be calculated through the following equation:(5)KII=τeπbk2sinψCsin2ψ+b2/a2cos2ψ1/4,KIII=τe1−νπbk2k′cosψCsin2ψ+b2/a2cos2ψ1/4,C=k2+νk′2Ek−νk′2Kk,k2=1−b2a2,k′2=b2a2,where a and b represent the semimajor axis and semiminor axis of the elliptical crack, respectively. ψ represents the ellipse polar angle, and K (k) and E (k) are the first and second complete elliptic integrals, respectively.According to the 3D coordinate transformation equation, the normal stressσN of crack initiation surface with the bending angle of θ and the torsion angle of φ can be expressed as the following equation:(6)σN=σxxsin2θcos2ϕ+σyycos2θcos2ϕ+σzzsin2ϕ−2τxysinθcosθcos2ϕ+2τyzcosθsinϕcosϕ−2τxzsinθsinϕcosϕ.Define the ratio of normal stressσN and far-field compressive stress σ as the stress concentration factor δ. Through solving the following partial differential equation, the point-by-point stress state of crack front and the change of bending angle θ and torsion angle φ with ellipse polar angle ψ can be obtained:(7)∂δ∂θ=0,∂δ∂ϕ=0.


### 4.2. 3D Crack Initiation Calculation

Due to the symmetry of the elliptical crack, only the calculation results on the 1/4 ellipse boundary were given. The specific calculation parameters are shown in Table3. Because the preexisting crack was subjected to pure shear stress when α was 0 and the preexisting crack was subjected to pure pressure stress when α was 90°, both cannot be considered as a mixed mode of fracture mode II and fracture mode III. The cracks with the angles of 15° and 75° were used to conduct theoretical analysis, replacing the cracks with the angles of 0 and 90°.Table 3
Specific parameters of 3D crack initiation calculation.
Semimajor axisa (mm)Semiminor axisb (mm)Ellipse polar angleψ (°)Friction coefficientμPoisson’s ratioνRadiusr (mm)7.505.000 ∼ 900.100.180.05The far-field compressive stress was set as 20 MPa to calculate the stress intensity factor (SIF) of fracture mode II and fracture mode III. According to equation (5), the variation of KII and KIII with elliptic polar angle ψ is shown in Figures 6 and 7.Figure 6
Variation of KII with elliptic polar angle ψ.Figure 7
Variation of KIII with elliptic polar angle ψ.As shown in Figures6 and 7, the following can be seen:(1)
ParameterKII increases monotonically with the increase of ellipse polar angle ψ, but KIII decreases monotonically with the increase of ellipse polar angle ψ. This phenomenon indicates that the shear effect of fracture mode II decreases gradually along the short axis end to the long axis end of the preexisting crack, but the tearing effect of fracture mode III increases gradually.(2)
When the crack angleα is 45°, the values of KII and KIII are always the largest. It indicates that the crack with an angle of 45° is subjected to the strongest shearing and tearing effects under the same loading conditions.Further, the relationship between the stress concentration factorδ and the ellipse polar angle ψ can be calculated, as shown in Figure 8. The following can be seen:(1)
For any angleα, the stress concentration factor δ increases monotonically with the increase of the ellipse polar angle ψ. When ψ ranges from 0° to 45°, the change of the parameter δ is obvious. When ψ exceeds 60°, the change trend of δ becomes flat. When ψ is equal to 90°, δ reaches the maximum value. It indicates that the preexisting crack starts to grow at the end of the short axis. Then, the crack continues to grow from the end of the short axis to the end of the long axis, and the growth rate decreases gradually.(2)
The stress concentration factorsδ at the end of the short axis for different crack angles α are shown in Figure 9. It can be found that the stress concentration factor δ at the end of the short axis presents an increase-decrease trend with the increase of crack angle. When the angle is 45°, the stress concentration factor δ reaches the maximum value (3.40). It means that the initiation process of the preexisting crack is easier to occur when angle α is 45°, which is consistent with the experimental result (Figure 3 and Table 1).Figure 8
Variation of stress concentration factorsδ with elliptic polar angle ψ.Figure 9
Variation of stress concentration factorsδ with crack angles α.


## 4.1. Mechanical Mechanism of 3D Crack Initiation

According to the theory of the maximum tensile stress, the crack will start to grow from the front edge point. For the elliptical crack used in this experiment, the growth of the elliptical crack can be considered as a mixed mode of fracture mode II and fracture mode III, and the corresponding 3D coordinate can be shown in Figure5 [39].Figure 5
The 3D coordinate diagram of elliptic crack.According to the elastic mechanics oblique section formula [40], the equivalent normal stress σe and equivalent shear stress τe acting on the preexisting crack surface can be calculated by using the following equations:(2)σe=σsin2α,(3)τe=σsinαcosα−μσsin2α,where σ is the far-field compressive stress, α is the angle of elliptic crack to the loading direction, and μ is friction coefficient between crack surfaces.Combined with Figure5, define point O as the center on the front edge of the preexisting crack in the local x-y-z coordinate system, and the following equation can be used to calculate the 3D stress component at the end of the crack with radius r and fracture bending angle θ:(4)σxx=−KII2πrsinθ22+cosθ2cos3θ2,σyy=KII2πrsinθ2cosθ2cos3θ2,τxy=KII2πrcosθ21−sinθ2sin3θ2,τxz=−KIII2πrsinθ2,τyz=KIII2πrcosθ2,σzz=−2νKIIπrsinθ2,where KII and KIII are the stress intensity factors of fracture mode II and fracture mode III, respectively. Parameter ν is Poisson’s ratio.The analytical solutions of stress intensity factors for fracture mode II and fracture mode III can be calculated through the following equation:(5)KII=τeπbk2sinψCsin2ψ+b2/a2cos2ψ1/4,KIII=τe1−νπbk2k′cosψCsin2ψ+b2/a2cos2ψ1/4,C=k2+νk′2Ek−νk′2Kk,k2=1−b2a2,k′2=b2a2,where a and b represent the semimajor axis and semiminor axis of the elliptical crack, respectively. ψ represents the ellipse polar angle, and K (k) and E (k) are the first and second complete elliptic integrals, respectively.According to the 3D coordinate transformation equation, the normal stressσN of crack initiation surface with the bending angle of θ and the torsion angle of φ can be expressed as the following equation:(6)σN=σxxsin2θcos2ϕ+σyycos2θcos2ϕ+σzzsin2ϕ−2τxysinθcosθcos2ϕ+2τyzcosθsinϕcosϕ−2τxzsinθsinϕcosϕ.Define the ratio of normal stressσN and far-field compressive stress σ as the stress concentration factor δ. Through solving the following partial differential equation, the point-by-point stress state of crack front and the change of bending angle θ and torsion angle φ with ellipse polar angle ψ can be obtained:(7)∂δ∂θ=0,∂δ∂ϕ=0.


## 4.2. 3D Crack Initiation Calculation

Due to the symmetry of the elliptical crack, only the calculation results on the 1/4 ellipse boundary were given. The specific calculation parameters are shown in Table3. Because the preexisting crack was subjected to pure shear stress when α was 0 and the preexisting crack was subjected to pure pressure stress when α was 90°, both cannot be considered as a mixed mode of fracture mode II and fracture mode III. The cracks with the angles of 15° and 75° were used to conduct theoretical analysis, replacing the cracks with the angles of 0 and 90°.Table 3
Specific parameters of 3D crack initiation calculation.
Semimajor axisa (mm)Semiminor axisb (mm)Ellipse polar angleψ (°)Friction coefficientμPoisson’s ratioνRadiusr (mm)7.505.000 ∼ 900.100.180.05The far-field compressive stress was set as 20 MPa to calculate the stress intensity factor (SIF) of fracture mode II and fracture mode III. According to equation (5), the variation of KII and KIII with elliptic polar angle ψ is shown in Figures 6 and 7.Figure 6
Variation of KII with elliptic polar angle ψ.Figure 7
Variation of KIII with elliptic polar angle ψ.As shown in Figures6 and 7, the following can be seen:(1)
ParameterKII increases monotonically with the increase of ellipse polar angle ψ, but KIII decreases monotonically with the increase of ellipse polar angle ψ. This phenomenon indicates that the shear effect of fracture mode II decreases gradually along the short axis end to the long axis end of the preexisting crack, but the tearing effect of fracture mode III increases gradually.(2)
When the crack angleα is 45°, the values of KII and KIII are always the largest. It indicates that the crack with an angle of 45° is subjected to the strongest shearing and tearing effects under the same loading conditions.Further, the relationship between the stress concentration factorδ and the ellipse polar angle ψ can be calculated, as shown in Figure 8. The following can be seen:(1)
For any angleα, the stress concentration factor δ increases monotonically with the increase of the ellipse polar angle ψ. When ψ ranges from 0° to 45°, the change of the parameter δ is obvious. When ψ exceeds 60°, the change trend of δ becomes flat. When ψ is equal to 90°, δ reaches the maximum value. It indicates that the preexisting crack starts to grow at the end of the short axis. Then, the crack continues to grow from the end of the short axis to the end of the long axis, and the growth rate decreases gradually.(2)
The stress concentration factorsδ at the end of the short axis for different crack angles α are shown in Figure 9. It can be found that the stress concentration factor δ at the end of the short axis presents an increase-decrease trend with the increase of crack angle. When the angle is 45°, the stress concentration factor δ reaches the maximum value (3.40). It means that the initiation process of the preexisting crack is easier to occur when angle α is 45°, which is consistent with the experimental result (Figure 3 and Table 1).Figure 8
Variation of stress concentration factorsδ with elliptic polar angle ψ.Figure 9
Variation of stress concentration factorsδ with crack angles α.


## 5. Numerical Simulation

Although the experiment is an important method to study the 3D crack initiation mechanism, the results of physical experiments can only explain the phenomenon of 3D crack initiation to a certain extent due to the limitations of sample material, observation accuracy, and data processing system. Therefore, in this paper, based on physical experiment and theoretical analysis, the numerical model was established to describe the whole process of 3D crack growth, which will further reveal the 3D crack initiation mechanism in the rock mass.


### 5.1. Calculation of the Stress Intensity Factors

#### 5.1.1. Theoretical Calculation Method

FRANC3−D simulation software was used to calculate stress intensity factors by M-integral. Because M-integral is derived from J-integral, it has a similar mathematical expression with J-integral [41], as expressed in the following equation:(8)J=∫Vσij∂ui∂x−W∗δ1j∂q∂xdV,where δ1j is the Kronecker notation, W∗ is the strain energy function, q is the smooth function, and V is the integral domain.J-integral and stress intensity factors also satisfy the following relation:(9)J=1−ν2EKI2+1−ν2EKII2+1+νEKIII2.Assuming that the material was linear elastic material, the stress field and displacement field can be defined as(10)σij=σij1+σij2,εij=εij1+εij2,ui=ui1+ui2.The stress intensity factors can be expressed as(11)KI=KI1+KI2,KII=KII1+KII2,KIII=KIII1+KIII2.In equations (10) and (11), superscripts “(1)” and “(2)” denote two independent linear elastic equilibrium states. Substituting equation (11) into equation (9), the following results can be obtained:(12)J=J1+J2+M1,2,J1=1−ν2EKI12+1−ν2EKII12+1+νEKIII12,J2=1−ν2EKI22+1−ν2EKII22+1+νEKIII22,M1,2=21−ν2EKI1KI2+1−ν2EKII1KII2+1+νEKIII1KIII2.According to equation (12), three equations are needed to obtain the three stress intensity factors. Assume that state (1) is the real equilibrium state to be solved, and state (2) is the virtual elastic equilibrium state. Furthermore, [KI(2a) = 1, KII(2a) = 0, KIII(2a) = 0], [KI(2b) = 0, KII(2b) = 1, KIII(2b) = 0], [KI(2c) = 0, KII(2c) = 0, KIII(2c) = 1] are used to replace state (2), respectively. Then, three equations can be built to obtain the stress intensity factors, as expressed in the following equation:(13)KI1=E21−ν2M1,2a,KII1=E21−ν2M1,2b,KIII1=E21+νM1,2c.In order to calculate the parametersM(1,2a), M(1,2b), and M(1,2c), substitute equation (10) into equation (8), as expressed in the following equation:(14)M1,2=∫Vσij1∂ui2∂x+σij2∂ui1∂x−W∗1,2δ1j∂q∂xdV,W∗1,2=σij1εij2+σij2εij1.In this paper, FRANC3−D is used to solve the stress intensity factors at the crack tip by using the above-mentioned M-integral method.


#### 5.1.2. Analysis of the Stress Intensity Factors

In this paper, FRANC3−D is responsible for the mesh generation and the calculation of the stress intensity factors, and model calculation after meshing can use ABAQUS. The finite element mesh generation of the numerical model is shown in Figure 10.Figure 10
Mesh generation of the numerical model.In this numerical simulation, the sample size is Φ50 × 100 mm, and the reference compressive stress is 20 MPa. The values of elastic modulus and Poisson’s ratio are consistent with the data in the previous study, that is,E = 6000 MPa, ν = 0.18, the elliptical polar angle ψ is 0° ∼ 90°, and the crack angle α can be 15°, 30°, 45°, 60°, and 75°, respectively.ByM-integral method, the variation of stress intensity factors KII and KIII with the elliptical polar angle ψ can be calculated, as shown in Figures 11 and 12, and numerical simulation is consistent with the previous analytical solution.Figure 11
Numerical simulation results of KII.Figure 12
Numerical simulation results of KIII.


### 5.2. Simulation of Preexisting Crack Growth Process

In this section, XFEM in ABAQUS is used to analyze the entire process of the crack growth with angle of 45°. The numerical model of the cement sample is established in the Part module. The size and material properties are consistent with the previous research. The maximum principal stress damage of the traction separation laws is used as the criterion for the crack initiation, and the maximum principal stress is 9.8 MPa. Considering the convergence of the calculation, the external load is set to be geometrically nonlinear. Because the displacement-controlling loading model was used in the experiment, the loading rate is set as 0.04 mm/step in the numerical simulation. Considering the complexity of the preexisting crack growth, the model is divided into eight equal parts by using the segmentation technology in the meshing process, obtaining 42000 structural units. The specific modeling and mesh generation are shown in Figure13.Figure 13
Calculation model and mesh generation.Figure14 shows the growth process of the 3D preexisting crack under uniaxial compression. The growth of the preexisting crack mainly generates the tension wing crack. The tension wing cracks are the clearest near the end of the short axis of the preexisting crack, presenting wrapping growth.Figure 14
Diagram of 3D preexisting crack growth process.Whenε is equal to 0.48εc, the wing cracks are generated near the end of the short axis of the preexisting crack and grow approximately perpendicular to the crack surface. In contrast, the kink zone near the end of the long axis grows parallel to the crack surface. When ε reaches 0.58εc, the wing cracks at the end of the short axis gradually grow to both ends to form wrapping wing cracks, and the upper and lower ends are symmetrical. With the loading process continuing, the wing cracks will wrap the whole preexisting crack. When ε reaches 0.72εc, the wing cracks growth direction is changed and starts to grow toward the maximum principal stress direction. Two petal-shaped cracks appear along the front edge of the wrapping wing cracks in the previous stage. The petal-shaped cracks grow gradually, and the overall growth rate is faster than the wrapping wing cracks. Then, the petal-shaped cracks grow along the axial loading direction and form vertical tension cracks. It is noticeable that the growth length of the wing crack is limited, and the maximum growth length is only 1.2 times the minor axis radius of the preexisting crack.


## 5.1. Calculation of the Stress Intensity Factors

### 5.1.1. Theoretical Calculation Method

FRANC3−D simulation software was used to calculate stress intensity factors by M-integral. Because M-integral is derived from J-integral, it has a similar mathematical expression with J-integral [41], as expressed in the following equation:(8)J=∫Vσij∂ui∂x−W∗δ1j∂q∂xdV,where δ1j is the Kronecker notation, W∗ is the strain energy function, q is the smooth function, and V is the integral domain.J-integral and stress intensity factors also satisfy the following relation:(9)J=1−ν2EKI2+1−ν2EKII2+1+νEKIII2.Assuming that the material was linear elastic material, the stress field and displacement field can be defined as(10)σij=σij1+σij2,εij=εij1+εij2,ui=ui1+ui2.The stress intensity factors can be expressed as(11)KI=KI1+KI2,KII=KII1+KII2,KIII=KIII1+KIII2.In equations (10) and (11), superscripts “(1)” and “(2)” denote two independent linear elastic equilibrium states. Substituting equation (11) into equation (9), the following results can be obtained:(12)J=J1+J2+M1,2,J1=1−ν2EKI12+1−ν2EKII12+1+νEKIII12,J2=1−ν2EKI22+1−ν2EKII22+1+νEKIII22,M1,2=21−ν2EKI1KI2+1−ν2EKII1KII2+1+νEKIII1KIII2.According to equation (12), three equations are needed to obtain the three stress intensity factors. Assume that state (1) is the real equilibrium state to be solved, and state (2) is the virtual elastic equilibrium state. Furthermore, [KI(2a) = 1, KII(2a) = 0, KIII(2a) = 0], [KI(2b) = 0, KII(2b) = 1, KIII(2b) = 0], [KI(2c) = 0, KII(2c) = 0, KIII(2c) = 1] are used to replace state (2), respectively. Then, three equations can be built to obtain the stress intensity factors, as expressed in the following equation:(13)KI1=E21−ν2M1,2a,KII1=E21−ν2M1,2b,KIII1=E21+νM1,2c.In order to calculate the parametersM(1,2a), M(1,2b), and M(1,2c), substitute equation (10) into equation (8), as expressed in the following equation:(14)M1,2=∫Vσij1∂ui2∂x+σij2∂ui1∂x−W∗1,2δ1j∂q∂xdV,W∗1,2=σij1εij2+σij2εij1.In this paper, FRANC3−D is used to solve the stress intensity factors at the crack tip by using the above-mentioned M-integral method.


### 5.1.2. Analysis of the Stress Intensity Factors

In this paper, FRANC3−D is responsible for the mesh generation and the calculation of the stress intensity factors, and model calculation after meshing can use ABAQUS. The finite element mesh generation of the numerical model is shown in Figure 10.Figure 10
Mesh generation of the numerical model.In this numerical simulation, the sample size is Φ50 × 100 mm, and the reference compressive stress is 20 MPa. The values of elastic modulus and Poisson’s ratio are consistent with the data in the previous study, that is,E = 6000 MPa, ν = 0.18, the elliptical polar angle ψ is 0° ∼ 90°, and the crack angle α can be 15°, 30°, 45°, 60°, and 75°, respectively.ByM-integral method, the variation of stress intensity factors KII and KIII with the elliptical polar angle ψ can be calculated, as shown in Figures 11 and 12, and numerical simulation is consistent with the previous analytical solution.Figure 11
Numerical simulation results of KII.Figure 12
Numerical simulation results of KIII.


## 5.1.1. Theoretical Calculation Method

FRANC3−D simulation software was used to calculate stress intensity factors by M-integral. Because M-integral is derived from J-integral, it has a similar mathematical expression with J-integral [41], as expressed in the following equation:(8)J=∫Vσij∂ui∂x−W∗δ1j∂q∂xdV,where δ1j is the Kronecker notation, W∗ is the strain energy function, q is the smooth function, and V is the integral domain.J-integral and stress intensity factors also satisfy the following relation:(9)J=1−ν2EKI2+1−ν2EKII2+1+νEKIII2.Assuming that the material was linear elastic material, the stress field and displacement field can be defined as(10)σij=σij1+σij2,εij=εij1+εij2,ui=ui1+ui2.The stress intensity factors can be expressed as(11)KI=KI1+KI2,KII=KII1+KII2,KIII=KIII1+KIII2.In equations (10) and (11), superscripts “(1)” and “(2)” denote two independent linear elastic equilibrium states. Substituting equation (11) into equation (9), the following results can be obtained:(12)J=J1+J2+M1,2,J1=1−ν2EKI12+1−ν2EKII12+1+νEKIII12,J2=1−ν2EKI22+1−ν2EKII22+1+νEKIII22,M1,2=21−ν2EKI1KI2+1−ν2EKII1KII2+1+νEKIII1KIII2.According to equation (12), three equations are needed to obtain the three stress intensity factors. Assume that state (1) is the real equilibrium state to be solved, and state (2) is the virtual elastic equilibrium state. Furthermore, [KI(2a) = 1, KII(2a) = 0, KIII(2a) = 0], [KI(2b) = 0, KII(2b) = 1, KIII(2b) = 0], [KI(2c) = 0, KII(2c) = 0, KIII(2c) = 1] are used to replace state (2), respectively. Then, three equations can be built to obtain the stress intensity factors, as expressed in the following equation:(13)KI1=E21−ν2M1,2a,KII1=E21−ν2M1,2b,KIII1=E21+νM1,2c.In order to calculate the parametersM(1,2a), M(1,2b), and M(1,2c), substitute equation (10) into equation (8), as expressed in the following equation:(14)M1,2=∫Vσij1∂ui2∂x+σij2∂ui1∂x−W∗1,2δ1j∂q∂xdV,W∗1,2=σij1εij2+σij2εij1.In this paper, FRANC3−D is used to solve the stress intensity factors at the crack tip by using the above-mentioned M-integral method.


## 5.1.2. Analysis of the Stress Intensity Factors

In this paper, FRANC3−D is responsible for the mesh generation and the calculation of the stress intensity factors, and model calculation after meshing can use ABAQUS. The finite element mesh generation of the numerical model is shown in Figure 10.Figure 10
Mesh generation of the numerical model.In this numerical simulation, the sample size is Φ50 × 100 mm, and the reference compressive stress is 20 MPa. The values of elastic modulus and Poisson’s ratio are consistent with the data in the previous study, that is,E = 6000 MPa, ν = 0.18, the elliptical polar angle ψ is 0° ∼ 90°, and the crack angle α can be 15°, 30°, 45°, 60°, and 75°, respectively.ByM-integral method, the variation of stress intensity factors KII and KIII with the elliptical polar angle ψ can be calculated, as shown in Figures 11 and 12, and numerical simulation is consistent with the previous analytical solution.Figure 11
Numerical simulation results of KII.Figure 12
Numerical simulation results of KIII.


## 5.2. Simulation of Preexisting Crack Growth Process

In this section, XFEM in ABAQUS is used to analyze the entire process of the crack growth with angle of 45°. The numerical model of the cement sample is established in the Part module. The size and material properties are consistent with the previous research. The maximum principal stress damage of the traction separation laws is used as the criterion for the crack initiation, and the maximum principal stress is 9.8 MPa. Considering the convergence of the calculation, the external load is set to be geometrically nonlinear. Because the displacement-controlling loading model was used in the experiment, the loading rate is set as 0.04 mm/step in the numerical simulation. Considering the complexity of the preexisting crack growth, the model is divided into eight equal parts by using the segmentation technology in the meshing process, obtaining 42000 structural units. The specific modeling and mesh generation are shown in Figure13.Figure 13
Calculation model and mesh generation.Figure14 shows the growth process of the 3D preexisting crack under uniaxial compression. The growth of the preexisting crack mainly generates the tension wing crack. The tension wing cracks are the clearest near the end of the short axis of the preexisting crack, presenting wrapping growth.Figure 14
Diagram of 3D preexisting crack growth process.Whenε is equal to 0.48εc, the wing cracks are generated near the end of the short axis of the preexisting crack and grow approximately perpendicular to the crack surface. In contrast, the kink zone near the end of the long axis grows parallel to the crack surface. When ε reaches 0.58εc, the wing cracks at the end of the short axis gradually grow to both ends to form wrapping wing cracks, and the upper and lower ends are symmetrical. With the loading process continuing, the wing cracks will wrap the whole preexisting crack. When ε reaches 0.72εc, the wing cracks growth direction is changed and starts to grow toward the maximum principal stress direction. Two petal-shaped cracks appear along the front edge of the wrapping wing cracks in the previous stage. The petal-shaped cracks grow gradually, and the overall growth rate is faster than the wrapping wing cracks. Then, the petal-shaped cracks grow along the axial loading direction and form vertical tension cracks. It is noticeable that the growth length of the wing crack is limited, and the maximum growth length is only 1.2 times the minor axis radius of the preexisting crack.


## 6. Conclusions

(1)
With the increase of preexisting crack angleα, both the initiation stress and initiation stress level present a decrease-increase trend. When α is equal to 45°, the initiation stress and initiation stress level reach the minimum value. Meanwhile, the samples with different preexisting crack angles correspond to a different failure mode. When the crack angle reaches 45°, shear failure characteristics are obvious.(2)
According to the analysis of the relationship between stress intensity factors (KII and KIII) and elliptic polar angle ψ, it can be found that the shear effect of fracture mode II decreases gradually along the short axis end to the long axis end of the preexisting crack, but the tearing effect of fracture mode III tearing effect increases gradually. What is more, the crack with an angle of 45° is subjected to the strongest shearing and tearing effects under the same loading conditions.(3)
Based on the variation of stress concentration factorsδ with elliptic polar angle ψ, the preexisting crack starts to grow from the end of the short axis. The path of crack growth is from the end of the short axis to the end of the long axis, and the growth rate decreases gradually.(4)
Through FRANC3−D and ABAQUS software, the growth process preexisting crack is simulated. The simulation results verify that preexisting crack grows from the end of the short axis to the end of the long axis. The growth direction changes from perpendicular to the crack surface to parallel principal stress direction and maximum growth length can reach 1.2 times the minor axis radius of the preexisting crack.



---

*Source: 5560253-2021-03-24.xml*