Numerical modeling of fracture process using a new fracture constitutive model with applications to 2D and 3D engineering cases

The overwhelming majority of experimental and numerical tests show the dependence of mechanical response on discontinuities (such as joints, faults, and bedding plane). In this study, fracture process is numerically investigated using finite‐difference method. A quasi‐continuum model (assuming that the real rock mass represents the sum of intact rock and fracture) has been developed to describe the fracture propagation in a fractured rock. First, this model has been validated against the experimental results and was shown to simulate satisfactorily the fracture propagation in the tensile, shear as well as the mixed (tensile and shear) modes. This model then has been used to investigate the fracture processes during the 2D long‐wall mining and 3D pillar failure. According to the 2D and 3D simulations, it was found that: (a) this model can provide a new reasonable approach to simulating fracture processes in either 2D or 3D cases; (b) roof failure is mainly caused by shear fractures rather than tensile fractures during the 2D long‐wall mining.


| INTRODUCTION
Natural discontinuities (eg, joints, faults, and bedding planes) are generated through long-term geological processes. The mechanical response of rock mass is known to be greatly impacted by these discontinuities. The existence of these discontinuities has been the principal obstacle to estimate precisely the deformability and stability of rock masses in practical engineering design. Thus, better interpretation of such influence is essential to investigate aiming at achieving realistically mechanical response of geo-materials.
fracture process (eg, initiation, propagation, coalescence, etc) are also strongly dependent on such discontinuities. This indicates that on one hand, the experimental/numerical tests using widely used intact specimen (without discontinuities) are not sufficient to constrain the real rock properties, on the other hand, the comprehensive model involving detailed information of fracture behaviors, in most cases, is not commonly available.
Considering the information achieved, the mechanism of the fracture process within the fractured rock is yet to be fully understood. As a result, different constitutive models (according to continuum or dis-continuum theory) involving the description of fracture behavior also have been proposed, but they were practically not used due to their complexity and limitations. In continuum methodology, [17][18][19][20] it was assumed that fractures should be uniformly distributed in each numerical zone. The impact of fractures on mechanical response was regarded as stress-dependent anisotropy of corresponding constitutive parameters (such as Young's modulus, Poisson's ratio) while plastic strain is employed to quantify the displacement discontinuity. However, since plastic strain is irreversible, this is contradictory with the fact that the contact states of fracture (opening, closure, and slippage) are reversible. To solve this problem, in dis-continuum methodology, fracture constitutive models [21][22][23] have been used to describe the mechanical response (eg, aperture variation, shear dilation, and slippage) for both contact and noncontact state of the fractures. A common view in these models is that rock failure is mainly attributed to the accumulation of open fractures, which are characterized by a constant normal stiffness. These obtained results show that fracture models have the ability to reproduce the real mechanical behaviors within the discontinuous rock mass. However, as mentioned above, on account that these fracture models usually have many parameters and are complex and complicated, they have been not widely used. A fracture constitutive model with simple function is therefore required to describe the fracture behaviors in a fractured rock.
Numerical methods are particularly useful for incorporating the influence of discontinuities on fracture process. Better understanding has been achieved on fracture initiation, propagation, and coalescence in discontinuous rock mass through numerous modeling simulations. [24][25][26][27][28][29][30][31] However, most of these existing methods have their limitations in representing the true physics of fracture behavior. Finite element method 32,33 is the first numerical method adopted to investigate the fracture behaviors, but re-meshing is required while calculating each step of fracture growth. This is time-consuming and may cause possible errors when interpolating variables from the old to the new mesh. Although some authors 34,35 have proposed the extended finite element method (XFEM) to solve the problem of re-meshing, XFEM is still not able to simulate the complex cases of fracture intersections. Discrete element methods provide good results in modeling fracture propagation, but cannot reflect the real characteristics in continuum field. On another aspect, most of the numerical investigations were focused on the 2D applications and only a few deal with 3D applications of engineering cases, which are typically more complicated, time-consuming, and costly or are impossible via experimental tests. Obviously, 3D analysis is more realistic and suitable for engineering applications.
In the present study, we investigate the fracture process in a fractured rock. This work has been motivated by many recent studies focusing on the fractured rock (using specimens with discontinuities) properties. An original quasi-continuum model is introduced that meets the requirement of being simple and at the same time capturing the principal features of the fracture behaviors. This model has been implemented into the finite-difference 3-D dynamic time-marching explicit code Flac3D. The proposed model was tested against the results from two types of experiments on rocks. The verifications show that the code could simulate fracture propagation in a tensile, shear, and the mixed (tensile and shear) modes. This model then has been used to investigate the fracture processes in the 2D long-wall mining and 3D pillar failure.
According to quasi-continuum theory, fractured rock could be regarded as a mixture of two parts, namely rock matrix element and fracture element ( Figure 1). In view of this point, a composed constitutive model is proposed to describe the fracture initiation and propagation in rock matrix elements as well as the fracture behaviors in fractured elements. The theoretical parts of the models are described in detail below. where d is a non-negative scalar function and g is the plastic potential function. The yield function used is the classical Mohr-Coulomb criterion (it can be any yield criterion) that is most widely used in rock engineering because of its simplicity. The yield functions (for shear mode in Equation (7) and for tensile mode Equation (8)) read where N = (1 + sin ) ∕ (1 − sin ); is the internal friction angle (°), c is the cohesion (MPa). Both and c are defined as usual from the conventional (axisymmetric) compression data that can be presented as 1 3 envelope. t is the tensile strength that can be evaluated from direct or indirect tensile tests (MPa). We assume that g is equal to F s and F t (associated rule), respectively, for shear and tensile modes. In this study, we use the second strain tensor deviator to evaluate the shear plastic strain. 36 It reads Accordingly, for tensile plastic strain, the maximum tensile strain is employed. It reads When the maximum shear (or tensile) strain is reached, calculation element fails. In general, the orientation of fracture is related to rock failure mode. For example, fracture F I G U R E 1 Sketch of quasi-continuum method ( n is normal stress; s is shear stress; f is normal stress for the fracture; f is shear stress for the fracture; m, l axis in the local coordinate system is decomposed from s )
caused by tensile mode is perpendicular to the maximum tensile principal strain; as for shear fracture, it is generated at the position where it has the least shear strength.

| Constitutive model for the fractured element
According to quasi-continuum theory, the total incremental strain of a fractured element is decomposed on the increment of the elastic strain d el ij and the fracture elastic strain d f ij . 37,38 Its expression is as follows: For the intact part of the fractured element, d el ij is related to stresses by Hooke's equations (Equation 5). As for fracture, it includes two cases: contact and noncontact states ( Figure 2). For example, for the contact state, the fracture can yield compression and shear stresses; otherwise, the fracture cannot bear any stress for the noncontact condition. Thus, the additional strain caused by the fractures should be respectively calculated with different constitutive models for the two conditions. The brief descriptions are as follows. where n , s are normal and shear stresses (MPa); k n , k s are normal and shear stiffness (k n = k s = 10 GPa∕m in this study); l c = V∕S c , which is the characteristic length of the fracture (m), V is the volume of the fractured element (m 3 ); S c is the fracture area in the fractured element (m 2 ). In the contact state, the normal stress imposed to the fracture surface plane is compressive, which means a nonzero shear stress along the fracture defined as where tan f is the friction coefficient of the fracture surface and reflects the degree of fracture roughness. Once shear stress exceeds the maximum shear stress, shear slip (or dilation) occurs.

| Noncontact state
In the noncontact state, normal and shear stresses are both zero. In fact, fracture strain is equal to the strain of intact rock matrix. It is calculated as where e 1 = K + 4G∕3.

| Transition of the contact state
In reality, the contact state of the fracture surfaces does not remain constant but always changes. The transition of fracture state can be divided into two modes ( Figure 2): contact → noncontact; noncontact → contact. Since the normal stress and strain of contact fracture satisfy the elastic relationship in the contact state, in this study, we take strain as criterion. That is, when the tensile strain is generated on the fracture surface, the fracture varies from contact state to noncontact state. On the contrary, when the tensile strain becomes zero, the fracture changes from noncontact state to contact state.

| Solution strategy
The formulated constitutive models have then been implemented into the finite-difference 3-D dynamic time-marching explicit code Flac3D 36 for the numerical modelings in this work. The global methodology and procedure for the finitedifference method (FDM) computation is briefly introduced. Figure 3 shows the flow chart of the solving procedure. The solution starts with model initiation. Then, according to the type of calculation elements, different constitutive models are employed for intact (Section 2.2) and fractured (Section 2.3) elements, respectively. For intact elements, the basic variables of elements, such as stress and displacement are computed to obtain the shear and tensile plastic strains. When the obtained either shear or tensile plastic strain in some element exceeds the corresponding maximum value, this element fails and it is immediately changed to fractured element. As for fractured elements, the contact state of the fractures must be first determined based on the description in Section 2.3.3 and then the corresponding constitutive model in Section 2.3 is adopted for continuing computation. The mechanical calculation terminates once the maximum unbalance force meets the preset error tolerance. To verify the effectivity and applicability of the models, some simulations of traditional rock mechanic tests first are carried out, as described in the next section.  Table 1. Note that since the maximum shear and tensile plastic strain decide how easy the calculation elements are fractured, the fracture process is greatly affected by these values. For all the simulation examples in this study, they are manually adjusted to matching the corresponding experimental tests.

| Results
We first begin with testing the role of the grid size (resolution) on the fracture propagation. We use homogeneous brick grids and specimen model sizes 61 × 122, 121 × 242, 181 × 362 and 241 × 482 numerical zones in the x and z directions, respectively. Nx stands for the number of zones in x direction. Figure 5 presents fracture propagations for different grid sizes at = 20 • . It is seen that the isolated failure plane thickness decreases as grid resolution increases. To clearly observe the errors, fracture paths for different grid sizes are illustrated in Figure 6. It shows that fracture path depends on the grid size but very little. Since the objective of this work was to simulate the fracture process, the following numerical models were performed with Nx = 61 in order to save calculating time.
The curves of load force vs displacement obtained from simulations are illustrated in Figure 7A. It is clear that load force decreases steeply (to approximately zero) until the peak value is reached, which also holds rather common for experimental tests. In addition, good fittings between experimental and numerical results ( Figure 7B) are found in the failure strength (peak point at each curve in Figure 7A) as well as its variations with the increase in dip angle. Figure 8 shows the numerical results of horizontal (x) and vertical (z) displacement distributions at different dip angles ( = 0 • , 10 • , 20 • , 30 • , and 45.9 • ). It is clearly observed that the evolution of displacement accompanied with fracture propagation, which allows for the fracture propagation path to be visualized. An example of the horizontal displacement distributions for different calculation steps ( Figure 9) and its final grid deformation (added by 10 times horizontal and vertical displacement) at = 20 • (Figure 10)  Figure 11. Overall, there are no significant differences in fracture paths between the experimental and numerical results except of the case at = 45.9 • . The relatively larger discrepancy at = 45.9 • might be due to the heterogeneity of the experimental specimen. These results suggest the proposed model could simulate the fracture process in the tensile mode.

| Modeling setup
The experimental tests reported by Zhao et al 14 are taken as an example to confirm fracture initiation and propagation in the shear mode. In this series of tests, the rectangular prismatic rock-like samples (40 mm × 10 mm × 100 mm) with one central prefabricated flaw (14 mm in length) were used, as shown in Figure 4B. The model includes cuboid specimen and two stiff elastic platens corresponding to steel. The specimen and platens are separated by the frictional interfaces with zero cohesion and the interface friction angle ( = 5 • , close to the real case). The elastic-plastic properties corresponding to the constitutive model are listed in Table 2.
Note that all the presented models in this section below are run with the same parameters listed in Table 2. The only parameter that varies is dip angle , and its values are specified in the Figure 12. During the testing process, the downward velocity (V z = −1 × 10 −6 m/s) is applied to the upper platen until the model fails. Figure 12A shows the numerical stress-strain curves of the specimens with one prefabricated flaw at different dip angles. It is apparent that the stress-strain curve first approximates a straight line corresponding to the linear elasticity and then exhibits a steep drop until the strain reaches a certain value, which captures well the postpeak characteristic of the experimental result. Figure 12B shows that the peak stresses for experimental and numerical results are matched rather satisfactorily as well as its variation as the dip angle increases. For instance, the peak stress is minimum when the dip angle is 60 • , which is also true for the numerical result.  Figure 13 compares the horizontal displacement distribution for experimental and numerical results. It is shown that some discrepancies in displacement fields between experimental and numerical results exist clearly, but many experimental features are also numerically captured. For instance, the positions where the maximum displacements occur are similar. Overall, the numerical results capture well the real characteristics of experimental results. As for fracture process, it can be observed from Figure 14 that on the whole, the fracture model could simulate the fracture process in the shear mode and even mixed mode (shear and tensile) well in comparison with experimental tests. However, there are still larger errors in the fracture propagation path (than those in three-point bending tests) between experimental and numerical results. These errors might be due to heterogeneities of the specimen and pre-existing micro-fractures, which result in deviations at random for some degree along the propagation path. But in most cases, like experimental tests, the fractures numerically initiate from the vicinity of the prefabricated fracture tips and propagate along the loading direction to the upper and lower boundaries of the specimen. The evolution of the horizontal displacement distributions for different calculation steps at = 30 • is presented in Figure 15. Overall, the numerical results in Figure 14A are reasonable. Besides, the fracture propagation path was clearly marked into two modes to study the failure mode (tensile and shear modes, Figure 14A). The numerical results show that when the dip angles are small (15 • and 30 • ), the failure modes are mainly tensile; when the dip angles are beyond 30 • , the macroscopic fractures exhibit a mixture of tensile and shear characteristics and tensile fractures dominate over shear fractures. It indicates that failure mode does not remain constant but varies as dip angle increases.

| Results
In general, numerical results provide good fittings to not only three-point bending tests (specimen with one single flaw) but also uniaxial compression tests (specimen with one single flaw). This implies that the new fracture models proposed in this study are able to describe well the real characteristics related to the fracture initiation, propagation and the failure pattern in tensile mode, shear mode, and the mixed mode (tensile and shear).

| Modeling setup
A numerical model measuring 5 m× 0.4 m × 1.85 m in x,y,z directions, respectively, is created in Flac3D to simulate the physical model performed by Kang et al. 41 The model contains mainly five rock seams: coal seam, fine sandstone, sandy mudstone, mudstone, siltstone, as shown in Figure 16. In the numerical model, the coal seam is at a depth of approximately 0.275 m and has a thickness of 0.15 m. Bedding planes are introduced as the ultrathin thin interlayer with elastic properties between different rock seams. The corresponding constitutive properties are listed in Table 3.  process, shear and tensile fractures are tracked during the calculation. Note that only the zones above coal seam are monitored. Blue lines stand for shear fractures while red lines represent tensile fractures. It can be found that, new shear fractures are generated in the overburden (in particular around the corner of roof) due to the stress concentration caused by the extraction process and then extend deeper into the roof as the face continues to advance. However, tensile fractures mainly appear in the form of bedding separation (the horizontal straight lines shown in Figure 17) because of the enhanced vertical subsidence (vertical displacement). This is so in the literature on long-wall mining researches. 42 Moreover, a comparison of the final fracture characteristics between the physical and numerical model is made in Figure 18. We can see that the numerically captured failure patterns and the experimentally captured failure patterns are rather similar on many aspects. Open fractures (characterized as the fractures in rock matrixes) zones located horizontally some distance above the immediate roof and vertically on the left side of the crushed roof, which is the similar case for the physical model. The left and right dip angles of the crushed zones in numerical model are 61° and 57°, respectively, in comparison with 55° and 60° in the physical model. The larger discrepancy in left dip angle is attributed to the deviation in the center of subsidence area ( Figure 19) and in the larger abutment pressure along the floor of coal seam near the left boundary of the model (Figure 20). Note that the abutment pressure distribution along the floor of coal seam ( Figure 20) produced by simulation (dash line) approximates physical model (solid line) rather satisfactorily.

| Results
In addition, the initial fracture interval is numerically identified (dash line shown in Figure 17F). The main difference in failure pattern is that on the right side of this line the immediate roof shows crushed failure characteristics while bedding separation is the main failure type on the left side of this line. However, in physical modeling, the crushed failure is also observed on the left side of this line. The reason might be that in reality the roof strata on the left side of this line falls down and then broke into pieces due to the gravity effect, which is not captured in the simulation. After statistics for rock matrixes, shear fractures dominate over tensile fractures by a ratio of approximately 10:1, which is in good agreement with numerical results (9:1, Gao et al 42 ) and field observation. 43 Furthermore, herein, we define damage as the ratio of the number of fractured zones to that of total zones above the coal seam. Figure 21 displays the progressive damage during the mining process. All the three series of points (shear, tensile, and total) show the same trend but different slopes. The function D p (d) (D p stands for progressive damage, and d is the mining distance) must meet the following requirements: (a) have qualitatively the same trend as the points in Figure 21; (b) the slope remains more or less constant at larger mining distance. The numerous trials have led to the following simple functions for shear Equation (17) and total damage Equation (18) Note that d in the equations is limited due to the special experimental condition. D s p and D tot p represent shear and total damage, respectively. Shear damage shows a hyperbolic relationship with mining distance. And this is the case for total damage, which is reasonable on account that shear failure mode is dominant for rock matrix during the mining process. As for tensile damage, it is related to linear relationship: where D t p is the tensile damage. Besides, the orientations of these fractures for rock matrixes are also calculated in Figure 22.
Most of the simulated fractures (shear and tensile) are oriented approximately 30° and 150° with respect to the horizontal line of roof, as can be observed in physical model ( Figure 18A).
The modeling findings in this section show good accordance with the large-scale physical model. It indicates that the proposed constitutive models are able to reproduce the actual phenomenon for the 2D complex case in reality. The next section, therefore, moves on to discuss further 3D application of pillar failure, which is difficult to test via experiments.

| Modeling setup
Generally, rock pillar has 3D geometries. Considering the effect of existing geological structure and in situ stress on the pillar behavior, a simply symmetric model ( Figure 23) of a typical block caving pillar in an underground mine is analyzed. The mechanical properties are listed in Table 4. In this case, for boundary conditions, the displacements on each side of the model are fixed perpendicular to the direction of boundaries. The initial stresses imposed in the x-direction, ydirection and z-direction are 10 MPa, 10 MPa, and 15 MPa, respectively. The extracted area shown in Figure 23 will be first extracted, and then sufficient calculation steps (10 000) have been run for stress release and pillar failure.

| Results
To quantitatively assess the degree of damage that pillar undergoes, the ratio of fractured zones to total zones in the pillar is computed and illustrated in Figure 24. The curve in Figure 24 depicts that firstly there is no damage after the extracted area is immediately extracted; then the slope increases due to the progressive generation of fractures; after reaching a certain damage, there is a rapid reduction in the slope; lastly, the slope remains more or less constant. The contours of the maximum principal stress corresponding six points marked in Figure 24 are drawn in Figure 25. It depicts that as calculation step continues, distressed region (corresponding to fractured region in Figure 26) appears at the corner of the pillar surface, then extend to the entire surface and eventually propagate toward the core. Figure 27 shows the progressive outlines of the fractures at different cross sections. We can see from Figure 27 that fractures firstly are generated at the corner of the pillar corresponding to the slabbing in situ and similar results are also obtained by Wagner, 44 then fractures extend to the center of the pillar corresponding to continued spalling in situ; lastly the pillar forms a "hourglass" failure shape ( Figure 26). These findings above are in well keeping with the numerical observations of Sainoki and Mitri 45 and the field observations of Pritchard and Hedley. 46 These results indicate the models proposed are capable of dealing with 3D engineering cases.

| CONCLUDING DISCUSSION
The impact of discontinuities such as fractures, joints and fault on mechanical behaviors has been long known from the experimental and numerical information. These results show discontinuities could strongly affect rock properties including strength, deformability and fracture process, etc. Traditionally, plastic damage are often employed to describe these effects of discontinuities in plastic constitutive models. Generally, this approach used in plastic models is capable of providing good explanations in many topics related to rock mechanics. However, since in reality the contact state of fractures constantly changes from contact to noncontact state (or from noncontact to contact state), these plastic models are not completely reasonable. This study therefore proposes a quasi-continuum model by considering the mechanical behaviors of fractures for both contact and noncontact states. Using the finite-difference code Flac3D as the numerical tool, the present study efficiently simulates the influence of discontinuities on fracture process in a rock. Then, this model is adopted to a 2D complex long-wall mining case and a 3D simple case of a typical block caving pillar. The obtained simulation results are in general agreement with the physical model and field monitoring. The specific conclusions are as follows: 1. Fracture processes including fracture initiation and propagation in three-point bending and uniaxial compression tests have been successfully reproduced with the proposed fracture constitutive model. It indicates this model can be utilized to simulate fracture processes in the tensile, shear, and mixed modes. 2. During the long-wall mining, shear fractures dominate over tensile fractures (mainly generated in form of bedding separations) in terms of rock matrixes. The most majority of fractures are orientated at 30° and 150° with respect to the horizontal roof line. This might be particularly useful in explaining the permeability anisotropy and predicting methane migration in situ but not discussed in detail here.

F I G U R E 2 5
The contour of the maximum compression stress at the points marked on damage curve in Figure 24 F I G U R E 2 6 3D failure shape comparison at point d in Figure 24: (A) numerical fracture shape; (B) zones without fracture; (C) field observation 46 3. The pillar experiences three stages (slow-growth damage; accelerated-growth damage; constant-growth damage) before failure. The obtained results of 3D case indicate the proposed model is capable of simulating the fracture initiation and propagation on 3D scale.
Overall, the proposed model could provide a new way to simulate the fracture initiation and propagation in many fields where fracture behaviors need be taken into account. Further work is required to clearly identify the limits of the proposed model and investigate some more complex engineering applications in particular in fully 3D scale.