#### Abstract

A grain-based rock model was developed and applied to study mechanical characteristics and failure micromechanics in thick-walled cylinder and wellbore stability tests. The rock is represented as an assembly of tetrahedral blocks with bonded contacts. Material heterogeneity is modeled by varying the tensile strength at the block contacts. This grain-based rock model differs from previous disk/sphere particle-based rock models in its ability to represent a zero (or very low) initial porosity condition, as well as highly interlocked irregular block shapes that provide resistance to movement even after contact breakage. As a result, this model can reach higher uniaxial compressive strength to tensile strength ratios and larger friction coefficients than the disk/sphere particle-based rock model. The model captured the rock fragmentation process near the wellbore due to buckling and spalling. Thin fragments of rock similar to onion skins were produced, as observed in laboratory breakout experiments. The results suggest that this approach may be well suited to study the rock disaggregation process and other geomechanical problems in the rock excavation.

#### 1. Introduction

Massive underground caving operations are becoming ever more prevalent. The development of cave infrastructure in rock strata may lead to excavation instability problems. It is increasingly important to better understand the redistribution of induced stresses and the failure process near underground excavations, and for support design, it is necessary to better predict the related deformation characteristics of brittle failing rock. Due to the natural fractures [1], bedding planes [2], and geological evolution, the rock mass contains lots of different scales of intrinsic defects including microcracks, joints, and cracks [3]. These joints or cracks significantly affect the rock mass mechanical properties such as deformability, strength, and failure process in different ways of providing the crack source for further failure of the rock mass, also leading to the brittle failure of the rock mass by the stress concentration at the crack tips [4, 5]. Investigation of the rock failure process is of great importance for rock stability in engineering exploitations and excavations [6].

Breakout and tensile failure are the typical failure mechanisms of wellbores’ stability. In the early time, Leeman [7] reported stress-induced breakouts causing sidewall spalling in South African gold mine openings. Image logs have been used (Barton and Zoback, [8]) to observe and investigate wellbore failures. Numerical modeling has been used for several decades to predict and simulate wellbore failures [9]. In the late 1980s and early 1990s, several linear elastic and poroelastic models for borehole stability were introduced by Bradley [10], Detournay and Cheng [11], Fuh [12], and Yew and Lui [13]. In the 1990s, nonlinear elastic, plastic, and elastoplastic models were developed by McLean and Addis [14], McLellan and Wang [15], and McLellan [16]. Several thermo- and chemoporoelastic models were also developed by Li [17], Wang and Dusseault [18], and others. Although more complex material behaviors were modeled in the later models, all of the aforementioned models were two-dimensional. Orozco [19] introduced a 3D model which investigated wellbores’ stability through salt zones. They pointed out that complex geomechanical characteristics (i.e., tectonic evolution, fracture distribution, physical and chemical properties of rocks, and the contact of different rocks) such as those in the rock adjacent to creeping salt bodies could induce wellbore instability. Wang and Samuel [20] also investigated wellbores’ stability in salt formations. Righetto [21] investigated wellbore (casing) collapse which could be caused by triaxiality of the stresses, which is not just a single mechanism. Zamiran [22] investigated inclined wellbore stability in anisotropic shale in the Bakken formation using FLAC3D and found that the deformation of the wellbore increased with a decrease in drilling inclination angle. Tang [23] and Liang [24] investigated the three-dimensional failure process of rock subjected to various stress loadings using a finite element method at mesoscopic level. Liu and Zhang [25] simulated the behavior of a cracked rock by introducing a crack into a continuum model with the superfine element division. In most of this literature, continuum codes were used to investigate the rock failures (i.e., stress concentrations). Although continuum approaches are effective to determine the magnitudes and directions of in situ stresses and strains, they have limited capabilities to model postfailure mechanisms such as spalling and bulking. Discontinuum schemes such as the rigid block spring method (RBSM) [26, 27], the discontinuous deformation analysis (DDA) model [28], the displacement discontinuity method (DDM) [29, 30], and the distinct element method (DEM) can simulate these rock failure scenarios, in addition to shear fracturing and tensile failure.

Garza-Cruz [31] developed a grain-based rock model with 3DEC for tunnel stability analysis. Different from particle flow code, when assigning parameters, the model is considered as a continuum and can be assigned macroparameters directly, while the joints of different blocks follow the Mohr-Coulomb criterion in 3DEC [32]. In addition, generally, the model has higher computational efficiency. This paper describes the application of a similar grain-based rock model to study wellbore breakout and thick-walled cylinder (TWC) tests used to predict fragmentation of rock masses. The model was first used to match the properties of typical rocks by simulating the direct tension test and uniaxial compression test. The calibrated model was then applied to study wellbore breakout and the thick-walled cylinder (TWC) test.

#### 2. Numerical Schemes

##### 2.1. Model Description

The numerical code used in this study, 3DEC, is a three-dimensional distinct element code which can be used to simulate the behaviors of discontinuous media such as fractured rock. As shown in Figure 1, the 3DEC modeling domain is discretized into blocks, representing a predefined discontinuity system. Each block can be further discretized into one or more tetrahedral zones in which an explicit finite difference scheme is used for the stress-strain calculation. The discontinuities are modeled by applying various joint models for linear and nonlinear force-displacement relations for movements in both the normal and the shear directions. A user-defined Mohr-Coulomb plug-in by built-in programming language was used to define the contact behavior with the benefit that joint properties can be locally assigned to each joint contact by plug-in joint model, without a constraint on the number of materials that can be defined.

In this study, the domain was directly discretized into tetrahedral blocks, with each block representing one zone or one rock grain. The strength (e.g., tensile strength and cohesion) at block contacts was assigned randomly from a strength distribution curve to introduce strength inhomogeneity in the model. This model has similar characteristics to the linear parallel bond model [33] in PFC3D, except that it exhibits a higher UCS to tensile strength ratio and friction angle because of the use of irregular and highly interlocked block shapes rather than spheres.

##### 2.2. Input Parameters

The wellbore breakout and thick-walled cylinder (TWC) tests are simulated by the proposed model. TWC test is often used to investigate the failure mechanism and form of wellbore, which is of great significance for underground engineering to ensure wellbore’ stability. The blocks and block contacts micromechanical properties used are summarized in Table 1. Note that some of the mesh-size dependent properties, e.g., contact stiffness, differ because the wellbore model and the TWC model have different model sizes and block (or grain) sizes.

##### 2.3. Microtensile Strength

In order to introduce material heterogeneity, a tensile strength value was randomly assigned to each block contact based on the cumulative distribution of microtensile strength in Figure 2, and the cohesion at each contact was fixed at 2.5 times the tensile strength. Finally, the tensile strength and cohesion were set to zero at 10% of the contacts, chosen randomly.

##### 2.4. UCS and Direct Tension Tests for Calibration

UCS and direct tension tests were modeled to benchmark the rock strength for wellbore stability simulations. The sample with the size of 1 m × 0.5 m × 2 m shown in Figure 3 was performed by the uniaxial compression and tension test.

A constant velocity of opposite direction was applied at the top and bottom of the sample to conduct compressing of the sample at a constant strain rate until it failed. The stress-strain response was measured at different locations in the sample to ensure that the test was performed in a quasistatic style. Figure 4 shows the stress-strain curve of the UCS sample measured at the top, the bottom, and the middle of the sample shown in Figure 3. The UCS of the sample is about 14.2 MPa (2060 psi). Similarly, a direct tensile test was performed by pulling each sample apart at a constant velocity until it failed. The tensile strength of this sample is 1.2 MPa, which is less than one-tenth of the UCS. The sample with heterogeneous microtensile and contact strengths exhibits higher macro-UCS strength and much lower macrotensile strength than the maximum of microtensile strengths. When a sample consisting of a collection of bonded irregular and highly interlocked blocks is compressed, strong force chains form to make subcontacts capable of subjecting to a higher level of compression stress. While the sample is tensile, the failure is controlled by the weakest links of the system oriented in the most favorable direction for stress concentration.

#### 3. Thick-Walled Cylinder Test

The TWC test has been widely used for sanding analysis [34, 35]. A TWC test is performed by maintaining a constant borehole pressure on the inner cylindrical wall while applying a stepwise-increasing radial loading stress on the outer cylindrical wall. During the loading process, the external volumetric stain is measured byin which is the volumetric strain; is the initial radius of the outer cylinder before loading; and is the current radius of the outer cylinder. Failure of the sample is indicated by a sudden increase in volumetric strain, and the loading stress at the failure point is identified as the TWC strength.

##### 3.1. Model Setup

The TWC specimen is a hollow cylinder with an axially aligned borehole. Figures 5 and 6 show the dimensions of the two TWC samples modeled in this study: one with an intact rock sample and the other with a rock sample containing an initial crack or flaw. The inner and outer diameters of the sample are 12.7 mm and 38.1 mm, and the height is 38.1 mm [36], modeled by over 42000 tetrahedral blocks with an edge length of about 3 mm.

The initial crack (highlighted in blue in Figure 6) has a dip of 45°, a constant tensile strength of 2 MPa, and cohesion of 5 MPa at each block contact on the crack. In this study, the inner borehole pressure was set to zero and the outer stress was increased with 1 MPa increment.

##### 3.2. Results of the Intact Rock Sample

Figure 7 shows the TWC test loading curve (volumetric strain versus radial loading stress) and the failed rock percentage curve (ratio of failed rock volume to the total sample volume) during the loading process for the TWC test with an intact rock sample. The loading curve shows a gradually increasing slope until reaching the failure point, indicating the increasing compliance of the sample due to the loading. Note that there is no initial strain hardening stage in the curve, possibly because the sample in this study is very compact with nearly zero porosity.

Four states, marked as A, B, C, and D on the loading curve, are selected to present the rock fragmentation profiles. Figures 8 and 9 plot the perspective and top views of the rock fragments at these four states. The colored blocks represent failed rock fragments, while the grey line segments in the plots illustrate the outlines of the unfailed blocks.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

At State A, the sample shows some scattered failed blocks which may occur at local low-strength zones. As the confining stress is increased, the sample fails progressively around the borehole. At State B, the failure near the borehole becomes more pronounced and the connected rock fragments in the perspective view indicate the formation of failure bands. State C is identified as the failure point because the volumetric strain increases dramatically and most of the near-borehole area has experienced failure. The TWC strength can be identified as 14.0 MPa, which is slightly smaller than the UCS (15.7 MPa). State D represents the state of the sample after catastrophic failure, when the sample is crushed and most of the blocks have failed. The sudden decrease in volumetric strain (instead of the expected increase) between State C and State D is due to the measurement methodology. The measurement of the radius of the outer cylinder Ro in (1) is based on the average distance of the outside block surface to the borehole centroid. The sudden catastrophic failure of the sample causes some of the outside blocks to be completely detached from the sample, which is captured as an apparently increased radius even though those blocks are no longer a part of the sample.

Figure 7 shows that the percentage of rock volume which has failed increases gradually up to State C with a sharp increase from State C to State D, confirming the sudden catastrophic failure. At State D, almost 100% of the rock sample has failed.

##### 3.3. Results of the Rock Sample with an Initial Crack

Figure 10 shows the TWC test loading curve and the failed rock percentage curve during the loading process for the TWC test on a rock sample containing an initial crack. The data shown in Figure 7 for the intact rock sample are also shown (dashed lines) for reference. The shape of the loading curve for the sample containing an initial crack is similar to that for the intact sample, with a gradually increasing slope until the failure point is reached. Although the crack and other joints have the same stiffness, the slope of the loading curve for the sample with the crack is steeper, because the lower strength of the crack causes the sample to have a lower effective stiffness. In contrast to the intact case, the measured volumetric strain increases dramatically as failure occurs (between points G and H) because the sample bulges prior to catastrophic failure.

The failed rock percentage in Figure 10 also shows a similar trend before the failed point to that for intact rock. However, the increase from State G to State H is not as sharp as the increase from State C to State D in Figure 7, consistent with the more progressive failure associated with the development of a bulged failure shape shown in Figure 11.

**(a)**

**(b)**

**(c)**

**(d)**

Four states, marked as E, F, G, and H in the loading curve, are selected to present the rock fragmentation profiles for the rock sample containing an initial crack. Figures 11 and 12 show the perspective view and top view of the rock fragments, respectively, at these four states. State E is similar to State A for the intact rock, in which the sample shows failure in a few isolated weak zones. State F is similar to State B for the intact rock, in which the failure near the borehole becomes more pronounced and some failure bands form. State G is identified as the failure point because the volumetric strain starts to increase dramatically and most of the near-borehole area has experienced failure. The TWC strength is 9.0 MPa, which is significantly smaller than the TWC strength for the intact rock sample (14.0 MPa) as expected. In contrast to the sudden catastrophic failure of the intact rock, the cracked sample shows a more progressive failure. The diameter of the failed zone near the crack is larger than in the rest of the sample due to the reduced load-bearing capacity of the crack relative to intact rock as shown in Figure 11.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Wellbore Stability Example and Validation

The numerical scheme was also employed to study wellbore failure around a horizontal well. Both breakout and tensile failure, also known as drilling-induced failure, can be modeled directly.

##### 4.1. Model Setup

Figure 13 shows the 1.2 m × 0.18 m × 1.2 m cuboid model used for the wellbore’ stability test. The horizontal wellbore is represented by a predefined cylinder with a radius of 0.12 m along the *y*-axis. The sample consists of over 200000 tetrahedral blocks with an edge size of approximately 5 cm.

The initial maximum and minimum horizontal stresses of 20 MPa and 18 MPa, respectively, were applied in the *X* and *Y* directions, and the vertical stress of 40 MPa was applied in the *Z* direction. These stresses should be viewed as effective stresses, representing conditions at a depth of approximately 3000 m, that is, SHmax = 50 MPa; SHmin = 48 MPa; Sv = 70 MPa; pore pressure = mud weight pressure = 30 MPa.

##### 4.2. Results

Figure 14 shows and in the plane perpendicular to the horizontal well (*XZ* plane), and Figure 15 shows the rock fragments near the wellbore in the same plane, with colored blocks indicating failed rock fragments. The plots show a similar stress diagram and a breakout angle of about 135 ± 5° which are consistent with FLAC3D modeling using the same conditions and properties. The fragmentation process near the wellbore shows the progressive buckling and spalling of the wellbore, and thin fragments of rock similar to onion skins are produced, as observed in laboratory breakout experiments [37].

**(a)**

**(b)**

#### 5. Conclusions

This paper presents a grain-based model to study the rock failure process. Compared with continuum models, the grain-based model is able to reveal the evolution of the fractures. For particle flow model, the model is easier to calibrate the parameters and more computationally efficient. Numerical simulations of wellbore breakout and thick-walled cylinder (TWC) tests by the proposed model were conducted to predict fragmentation of rock masses. Based on findings from the study, the following general conclusions are drawn:(1)The proposed grain-based rock model can be used for small-scale modeling of low-porosity rocks. The highly interlocked irregular block shapes provide resistance to movement even after contact breakage, enabling simulation of materials with higher uniaxial compressive strength to tensile strength ratios and larger friction coefficients than disk or sphere particle-based rock models.(2)The model has been calibrated using UCS and direct tension tests of samples. The UCS of the sample is about 14.2 MPa, while the tensile strength is 1.2 MPa, which is less than one-tenth of the UCS. The interlocked irregular block shapes that provide resistance to block rotation (moments) after contact breakage contribute the high uniaxial compressive strength to tensile strength ratio.(3)Numerical simulations by the prosed model for thick-walled cylinder and wellbore stability tests show that the effects of both heterogeneous strength properties and larger-scale features such as fractures have been demonstrated. The results of the TWC test with an initial crack indicate the importance of including the effects of natural fracture systems in wellbore stability calculations. The model captured the fragmentation process near the wellbore due to buckling and spalling, and thin fragments of rock similar to onion skins were produced, as observed in laboratory breakout experiments.(4)The effects of more heterogeneous material properties and more extensive natural fracture systems on wellbore stability in naturally fractured media will be presented in future publications. The results suggest that this approach may be well suited to study the rock disaggregation process in other geomechanical problems in the rock excavation.

#### Abbreviations

RBSM: | Rigid block spring method |

DT: | Direct tensile test |

DDA: | Discontinuous deformation analysis |

TWC: | Thick-walled cylinder |

DDM: | Displacement discontinuity method |

3DEC: | Three-dimensional distinct element code |

DEM: | Distinct element method |

PFC3D: | Three-dimensional particles flow code |

UCS: | Uniaxial compression strength |

FLAC3D: | Three-dimensional fast Lagrangian analysis of continua. |

#### Data Availability

This is an open-access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors thank Lee B.T., Mack M., Garza-Cruz T.V., and Pierce M. for the constructive suggestions and useful discussion related to applications of this grain-based rock model in the paper.