A Numerical Research on Crack Process of Gypsum Containing Single Flaw with Different Angle and Length in Uniaxial Loading

To investigate the crack behaviour of rock or rock-likematerial in uniaxial loading, a series of numerical simulationswere conducted on gypsum specimens containing a single flaw with different inclination angle (0–90) and length (10mm–30mm). Based on the numerical simulations results, the effect of flaw length on peak strength, the CI stress, and the CD stress were investigated with different inclination angles.The results show that the peak strength decreased initially and then increasedwith increasing of the flaw angle. Meanwhile, the peak strength decreased gradually when the length of the preexisting flaw increased. When the inclination angle was 30, 45, and 60, the reduction degree of peak strength increased with increasing of the flaw length. The CI stress and CD stress not only depend on the inclination angle but also depend on flaw length. Four types of crack were observed in numerical simulations. The present research facilitates increased understanding of crack behaviour of rock under different conditions.


Introduction
The cracking behaviour of rock and rock-like specimens containing artificial flaws under compressive loading has been experimentally investigated for decades. For example, gypsum [1][2][3][4][5], glass [6], polymethyl methacrylate (PMMA) [7][8][9], diastone [1-3, 7, 10-13], sandstone [8,14], granite [7], and limestone [9] were tested. Similar observations have been made in their experiments using different terms-two types of cracks are observed in compression: tensile or wing cracks and shear or secondary cracks ( Figure 1). Wing or tensile cracks were usually found to be the first cracks to appear from the preexisting flaws. Secondary cracks are shear cracks and always initiate from the tips of the flaws after the first cracks. They are initially stable, but they may become unstable during the loading process near crack coalescence or near specimen failure. Some studies revealed that crack inclination angles have significant effect on the crack propagation and coalescence process. Bobet [15] analysed the initiation of secondary cracks in gypsum with two flaws (angle from 30 ∘ to 60 ∘ ) under compression. It observed all flaw inclination angles, wing cracks, and secondary cracks initiate at the same stress level. Park and Bobet [16] investigated crack coalescence in gypsum specimens with open and closed flaws under compression and one of their results showed that the mean wing and secondary crack initiation stresses increase with flaw inclination angle and coalescence stress slightly increases as the flaw angle increases. Yang and Jing [14] analysed strength failure and crack coalescence behaviour of brittle sandstone samples containing a single fissure (angle from 0 ∘ to 75 ∘ ) under uniaxial compression. They found that the fissure angle has a key effect on the strength and deformation behaviour of sandstone samples under uniaxial compression, and the mechanical parameters (uniaxial compressive strength Young's modulus and peak axial strain) of flawed sandstone firstly decrease and then increase with increasing fissure angle. Zhou et al. [17,18] conducted a series of uniaxial compression experiments and numerical simulation on rock-like materials (formed from a mixture of sand, plaster, limestone, and water at mass ratio of 126 : 9 : 9 : 16) containing multiple flaws (angle from 30 ∘ to 60 ∘ ) and concluded that the crack initiation modes depend on the flaw angle and the nonoverlapping length.
A number of experimental and numerical investigations have been conducted for fractured rock under compression since Nemat-Nasser and Horii [20] investigated the fracturing mechanisms of flaw plates (model material) and obtained that flaw length is one of the parameters controlling the failure mode of sample. Wong et al. [21], Tang and Kou [22], and Tang et al. [23] performed uniaxial and biaxial compression on numerical model samples containing a number of large, preexisting flaws and a row of suitably oriented smaller flaws by RFPA2D software, which shows that the crack growth is stable and stops at some finite crack length with the increase of lateral pressure, while a lateral tensile stress with even a small value will result in an unstable crack growth after a critical crack length is attained. Yang and Jing [14] performed a series of uniaxial compressive tests on brittle sandstone which contained a single fissure (angle from 30 ∘ to 60 ∘ and length from 0 to 25 mm). Zhang  It was found that the uniaxial compressive strength, Young's modulus, and peak axial strain of sandstone samples with preexisting single fissure are all lower than those of intact sandstone sample, which is closely related to the fissure length and fissure angle.
The above experimental and numerical results for precracked rock or rock-like material have distinctly shown that inclination angle and length of cracks in rock or rocklike material have an important influence on the crack behaviours. However, it should be noted that some previous investigations focused mainly on a fixed length of precrack with varying inclination angle for rock and the length of crack was not taken in account. Some previous literatures analysed the peak strength and deformation modulus with different crack lengths but lower inclination angle (30 ∘ -60 ∘ ), due to hard fabrication of crack. Less experimental studies were performed on rock with different crack lengths and different inclination angles (0-90 ∘ ) to investigate the initiation and propagation process of precrack rocks. Therefore, in this paper, to better understand the effects of inclination angle and length of crack on the peak strength, crack initial stress, crack secondary stress, and failure behaviour of precrack rock, a series of numerical simulations with distinct element code (UDEC) on the uniaxial compression tests were conducted for precrack gypsum specimen with varying inclination angle (0-90 ∘ ) and length (10-30 mm) of cracks. Moreover, the Voronoi block in UDEC is adopted to better describe crack initiation, propagation, and coalescence process.

Introduction about Distinct Element
Method. UDEC is a two-dimensional numerical program based on the distinct element method for discontinuity modelling. It has been widely used in the literature to explore damage of rock [24][25][26][27][28][29], as well as Brazilian test analysis [30,31]. UDEC especially joint model used in UDEC captures several features that are representative of the physical response of joints. For joint model in UDEC, all blocks have connection with contacts (corner-to-corner contacts, edge-to-corner contacts, and edge-to-edge contact). However, edge-to-edge contact is important, because it corresponds to the case of a rock joint closed along its entire length ( Figure 2). If the maximum tensile stress exceeds the tensile strength of the contact or the maximum shear stress exceeds the shear strength of the contact, contact will break. And a microtensile crack or a microshear crack, respectively, will form. The contact behaviour is described by where and are the normal and shear stresses, and are the normal and shear displacements, and are the tensile and residual tensile strengths, max is the shear strength, and are the cohesive and the residual cohesive strengths, and are the friction and the residual friction angles, and Δ is the incremental contact shear displacement.

Procedure of the Simulation.
In this study, the UDEC model is the gypsum specimen which was selected from previous literature [1]. The reason is that the experiment tests about gypsum specimen containing single crack under different conditions were conducted by Wong [1], and the physical experiments were also repeated numerically by PFC in previous literatures [32,33]. The UDEC model was set up with elastic but indivisible Voronoi blocks. Model size was chosen to replicate those of the diametrical cross section of specimens used for laboratory testing (i.e., 50 mm width  and 100 mm height) (see Figure 2). Average edge length of Voronoi blocks was selected as 0.002 m which yielded about 1320 Voronoi blocks in total within the model. Every Voronoi block is divided into a plurality of triangular-shaped block. The maximum edge length of the triangular zones is 0.00133 m. Figure 2 shows a typical UDEC model and how individual Voronoi blocks are numerically interconnected.
For a distinct element model such as the UDEC that synthesizes macroscale material behaviour from the interactions of microscale components, the input properties of the microscopic constituents are usually not known. Such intrinsic characteristics imply that it is nearly impossible to calibrate a set of macro output properties from the assembly such that they are exactly equal to that of the physical tests. Only a set of comparable properties can be obtained from the calibrated model. The parameters of gypsum in experimental tests and numerical simulations are shown in Table 1. Figure 3 shows the stress-strain curves of the intact gypsum specimen of UDEC and experiment test under uniaxial compression.
The second stage of simulation involved a series of numerical simulation tests with single crack under different conditions: the crack length of gypsum specimen is fixed at 10 mm, 15 mm, 20 mm, 25 mm, and 30 mm, respectively, while the crack inclination angle measured horizontally is varied from 0 ∘ to 90 ∘ (see Figure 2).  In the present numerical study, specimens are loaded under a uniaxial vertical compression in a velocity controlled manner. A sufficiently low loading rate of 0.005 m/s is applied to ensure that the specimen remains in a quasistatic equilibrium throughout the test. It takes approximately 1000000 steps to load one specimen to complete failure.

Numerical Simulation Results and Analysis
The results of the simulation work performed as described in the previous section were analysed in different ways to meaningfully describe crack behaviour of gypsum in different conditions under uniaxial compression. Figure 4 shows the axial stress-axial strain curves and failed contact number-axial strain curves of UDEC model containing a single flaw 4 Shock and Vibration with different inclination angles and lengths under uniaxial loading, respectively. As shown in Figure 4, all axial stressaxial strain behaviour can be divided into four stages. To show this part clearly, an axial stress-axial strain curve for model containing a single flaw with inclination angle 30 ∘ and length 10 mm under uniaxial compression is depicted, and four pictures of UDEC model at different stages were included in Figure 5.

Stress-Strain Characteristics.
Before the initiation of macrocrack from the preexisting crack, only a few microcracks (white colour) initiated at the axial stress of 9 MPa. In this stage, the curve of failed contact number is almost unchanged. As the axial stress reached 17.6 MPa, macrotensile crack (white colour) initiated from both ends of preexisting crack. The number of failed contacts increased steadily with increasing axial strain in this stage. When the axial stress increased to about 25.5 MPa, the secondary crack initiated from left top of preexisting crack. The number of failed contacts increased rapidly when the axial strain is changed a little. The peak value of stress-strain curve is 26 MPa. After the peak point, a large macrocracks' coalescence has happened and the model failed.
The failed contact numbers under different conditions are shown in Figure 6. It can be observed from Figure 6 that the failed contact number at first increased and then decreased with increasing of inclination angle under the same length of preexisting crack. When the inclination angle was 0 ∘ and 90 ∘ , respectively, the failed contact number is almost unchanged with length of preexisting crack from 10 mm to 30 mm. However, when the inclination angle was 15 ∘ and 60 ∘ , respectively, the failed contact number increased gradually when length of preexisting crack increased. But when the inclination angle was 30 ∘ and 45 ∘ , respectively, the failed contact number increased gradually when length of preexisting crack decreased. And this decrease trend is most obvious when the inclination angle was 45 ∘ , as shown in Figure 6.
The peak strengths of model containing single flaw are shown in Figure 7. It is shown from Figure 7 that the peak strength at first decreased and then increased with increasing of the flaw angle. The behaviours are akin to those observed by Zhou et al. [17], Yang and Jing [14], Xu et al. [34], and Wasantha et al. [35]. Meanwhile, the peak strength decreased gradually when length of preexisting crack increased. To show the extent of the reduction in peak strength of the model 6 Shock and Vibration where cm is the peak strength of intact model (the value is 33.5 MPa). cn is the peak strength of model containing single flaw. The reduction degrees for peak strengths of different models are listed in Table 2 and Figure 8.
For the model with flaw length of 10 mm and 15 mm, increased from 2.6% to 24.2% when the inclination angle was increased from 0 ∘ to 90 ∘ . For the model with flaw length of 20 mm, 25 mm, and 30 mm, increased from 2.9% to 50.5% when the inclination angle was increased from 0 ∘ to 90 ∘ . For the model with flaw length > 15 mm, increased when the inclination angle was 30 ∘ , 45 ∘ , and 60 ∘ , respectively. For the model with flaw length < 15 mm, almost unchanged when the inclination angle was 0 ∘ to 90 ∘ , respectively. When the inclination angle was 90 ∘ , 75 ∘ , and 0 ∘ , the flaw length has no influence on . When the inclination angle was 30 ∘ , 45 ∘ , and 60 ∘ , increased with increasing of the flaw length .  Figure 9: The methods used to establish crack initial stress and crack damage stress [19]. discussed. Nicksiar and Martin (2012) [19] reviewed the various methods about crack initiation stress in previous literatures. In this study, Poisson's ratio method was used to obtain the crack initiation stress (see Figure 9(b)). The ratio remains constant until cracking occurs. The stress associated with the change in Poisson's ratio is referred to as the CI stress [19,36]. And the volumetric strain method was used to obtain the damage stress (see Figure 9(a)). When the unstable crack begins to grow, the volumetric strain reaches maximum. The stress associated with maximum volumetric strain is referred to as the CD stress [19,37,38]. In general, Poisson's ratio and volumetric strain are calculated based on the lateral and axial strains. In UDEC model, the lateral and axial strains are calculated using the average displacement in the and directions at dozens of locations in some sampling area to monitor the specimen response ( Figure 10).

Crack Initiation Stress and
The CI stresses ci and ci / cn are shown in Table 3 and Figure 11 by a series of calculations. As shown from Table 3, a general increase of the CI stress ci with the increase of flaw inclination angle from 30 ∘ to 90 ∘ is observed. The CI stress ci decreased mainly when flaw inclination angle < 30 ∘ . The CI stress ci decreased gradually with increasing of flaw length . The effect is relatively more prominent for 30 ∘ , 45 ∘ , and 60 ∘ . ci / cn displayed the similar law in Figure 11. ci / cn decreased mainly when flaw inclination angle < 45 ∘ . ci / cn increased with the increase of flaw inclination angle from 45 ∘ to 90 ∘ . ci / cn increased with the decrease of flaw length . In addition, for the model with flaw inclination angle of 0 ∘ , 15 ∘ , 75 ∘ , and 90 ∘ , ci / cn changed from about 0.6 to 0.8. For the model with flaw inclination angle of 30 ∘ , 45 ∘ , and 60 ∘ , ci / cn changed from about 0.1 to 0.3.
The CD stresses cd and cd / cn are also shown in Table 4 and Figure 12 by a series of calculations. The CD stresses cd and cd / cn decreased with the increase of flaw length . And the CD stress cd at first decreased and then increased with the increase of flaw inclination angle from 0 ∘ to 90 ∘ . cd / cn decreased mainly when flaw inclination angle < 30 ∘ . cd / cn   increased with the increase of flaw inclination angle from 30 ∘ to 90 ∘ . In addition, when flaw inclination angle < 15 ∘ or > 75 ∘ , the CD stress cd is substantially equal to the peak strength. With the decrease of flaw length , the effect is relatively more obvious.

Description of Crack Behaviour.
In a flaw-containing specimen, the block sizes have effect on the crack processes [33]. Wong results showed seven crack types with different trajectories and initiation mechanism ( Figure 13). Based on Wong results, this section will analyse the crack types of model containing single flaw with different inclination angle and length. The crack behaviours of model containing single flaw with different inclination angle and length were analysed in Figures 14, 15, and 16. Four pictures of each model show the crack initiation, propagation, crack type at peak stress, and sketch of cracks.
For length of 10 mm and flaw angle of 30 ∘ , 45 ∘ , and 60 ∘ , the crack initiation and crack propagation processes are shown in Figure 14. As can be seen from Figure 14, the first macroscopic cracks which are wing cracks (T1) are initiated from the tips of flaws. Cracks initiating later than the first cracks were usually referred to as secondary cracks. For the flaw angle of 30 ∘ , a portion of the secondary cracks initiated upwards near the left tip, which are antitensile cracks (T3). A portion of the secondary cracks initiated downwards from the left tip, which are mixed tensile-shear cracks (S-T) ( Figure 14(30-4)). For the flaw angle of 45 ∘ , the secondary cracks consisted of some microcracks clusters initiated at a distance away from the preexisting flaw, which are tensile cracks (T) (Figure 14(45-4)). For the flaw angle of 60 ∘ , three secondary cracks initiated from the left tip of flaw, two of them initiated downwards are mixed tensile-shear cracks (S-T), the third crack initiated upwards is antitensile crack (T3). One secondary crack initiated downwards from the left tip of flaw is antitensile crack (T3) (Figure 14(60-4)). Figure 15 shows the crack initiation and crack propagation processes with length of 20 mm and flaw angle of which are substantially parallel to the first crack. The one initiated upwards is antitensile crack (T3). And other cracks are mixed tensile-shear cracks (S-T). Three secondary cracks initiated from right flaw tip, which are type 2 tensile cracks (T2) (Figure 15(30-4)). For the flaw angle of 45 ∘ , two first cracks initiated from the flaw tip, which are wing cracks (T1). Two secondary cracks initiated upwards from right flaw tip are mixed tensile-shear cracks (S-T). Two secondary cracks  initiated downwards from left flaw tip are type 2 tensile cracks (T2) (Figure 15(45-4)). For the flaw angle of 60 ∘ , two macroscopic wing cracks (T1) initiated from the flaw tip as the first cracks and then four secondary cracks initiated from the right flaw tip, one of which is type 2 tensile crack (T2). Other secondary cracks are tensile cracks (T) (Figure 15(60-4)). Figure 16 shows the crack initiation and crack propagation processes with length of 30 mm and flaw angle of    = 10 mm  T1, T3, S-T  T1, T  T, T1, T3, S-T  = 20 mm  T1, T2, T3, S-T  T1, T2, S-T  T, T1, T2  = 30 mm  T,T1, T2, T3  T1, T2, S-T  T1, S, S-T of macroscopic cracks initiated from the flaw tip. The one initiated from the left flaw tip is mixed tensile-shear crack (S-T). While the one initiated from the right flaw tip is type 2 tensile crack (T2) (Figure 16(45-4)). For the flaw angle of 60 ∘ , two macroscopic tensile wing cracks (T1) initiated from the flaw tip as the first cracks. One secondary crack as shear crack (S) initiated at distance from the right flaw tip. Two secondary cracks initiated near the left flaw tip are mixed tensile-shear cracks (S-T) ( Figure 16(60-4)).
Five crack patterns are observed in single flaw with different geometries shown in Table 5: tensile crack (T), type 1 tensile crack (T1), type 2 tensile crack (T2), type 3 tensile crack (T3), mixed tensile-shear crack (S-T), and shear crack (S). Meanwhile, the crack types changed with different geometries. This indicated that the geometries have great effect on the crack type.

Discussion
From the loading of UDEC model, the peak strength, the CI stress, and the CD stress were obtained. The CI stress ratios ( ci / cn ) at first decreased and then increased as the flaw inclination angle increased, which was also observed in some references [7,32,39]. With the increasing of flaw length, the variation shape of the CI stress ratios changed from "V" into "U"; note that some points are not within the scope of the law. One plausible explanation was due to definition of the CI stress. In this paper, Poisson's ratio method was used to obtain the crack initiation stress [19]. But in actual operation, the CI stress is defined with 10% change of Poisson's ratio in some models. Another plausible explanation was due to random distribution of block sizes.
In the present numerical study, tensile crack (T), type 1 tensile crack (T1), type 2 tensile crack (T2), type 3 tensile crack (T3), mixed tensile-shear crack (S-T), and shear crack (S) were obtained. Compared with Wong literature [2,3], some of the paths of cracking were different in the present study. One plausible explanation was due to influence of grain size. The grain size in UDEC influences the onset of cracking and possible initiating mechanisms [19]. Table 3 shows the crack process with length of 10 mm and flaw angle of 30 ∘ , 45 ∘ , and 60 ∘ under different grain sizes. The number of microcracks increased with the increase of grain size. The first crack which is wing crack (T1) has the same trend. With the increase of grain size, the trend becomes better. Table 6 also shows difference about the type of secondary crack. For example, for the flaw angle of 45 ∘ , the secondary cracks for = 853 appear to be tensile crack (T). But the secondary cracks for = 2295 appear to be antitensile crack (T3), mixed tensile-shear crack (S-T), and shear crack (S) together.
So, the grain size has effect on type of crack in numerical simulation. In particular, the smaller the grain sizes, the more the time required for calculations. In the future research, to obtain good results, a more reasonable block size should be adopted.

Conclusion
A series of numerical simulations on rock-like specimens containing single flaw with different inclination angle and flaw length were subjected to uniaxial compressive loads. The main conclusions can be summarized as follows: (1) The peak strength at first decreased and then increased with increasing of the flaw angle. Moreover, the peak strength decreased gradually when length of preexisting crack increased. When the inclination angle was 90 ∘ and 0 ∘ , the flaw length has no influence on . When the inclination angle was 30 ∘ , 45 ∘ , and 60 ∘ , increased with increasing of the flaw length .
(2) The CI stress and the CI stress ratios ( ci / cn ) depend on the inclination angle and flaw length. The CI stress, the CI stress ratios ( ci / cn ), and the CD stress at first decreased and then increased as the flaw inclination angle increased. With the increasing of flaw length, the variation shape of the CI stress and the CI stress ratios changed from "V" into "U"; the inclination angle and flaw length almost have no effect on cd / cn .