Modeling discrete fractures in continuum analysis and insights for fracture propagation and mechanical behavior of fractured rock

Results in Engineering 4 (2019) 100070

Contents lists available at ScienceDirect

Results in Engineering journal homepage: www.editorialmanager.com/rineng/Default.aspx

Modeling discrete fractures in continuum analysis and insights for fracture propagation and mechanical behavior of fractured rock Goodluck I. Ofoegbu a, *, Kevin J. Smart b a b

GNO Modeling Research LLC., PO Box 769427, San Antonio, TX, 78245-9998, USA Space Science and Engineering Division, Southwest Research Institute®, 6220 Culebra Road, San Antonio, TX, 78238-5166, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Brittle deformation Continuum analysis Discrete fracture Modeling Shear fracture Tensile fracture

This paper presents a procedure for modeling discrete fracture initiation, deformation, and propagation in continuum analysis of rock, and provides numerical simulations to evaluate performance of the procedure against analytical solutions and laboratory data. The procedure is based on aggregating the mechanical behavior of evolving discrete fractures and unfractured rock and provides specific information about fracturing; such as when and where fractures may occur and the type (tensile or shear), geometrical attributes, and aperture of each fracture. By storing fracture geometry as continuum properties that are used internally to generate geometrical parameters needed for analysis, the procedure models the geometry and mechanical behavior of discrete fractures explicitly in continuum analysis without including the fractures in the model grid. Numerical simulations are described for typical laboratory loading conditions such as triaxial compression and extension, and Brazilian testing. The results are compared with laboratory data, analytical solutions, and established concepts to make a case that the procedure performs satisfactorily and can be relied upon for understanding and predicting rock mechanical behavior. Furthermore, the paper discusses insights based on the simulation results for understanding rock fracturing; mechanical behavior of fractured rock; and large-scale behavior such as localization of rock deformation, stress control of fracture mechanisms in failure zones, and conditions for occurrence of multiple parallel fractures.

Introduction Brittle deformation of rocks or rock masses is accommodated through the initiation and growth of new fractures and slip or opening of existing fractures [1–11]. A rock mass consists of both intact rock (with small-scale and discontinuous fractures or cracks) and large-scale discontinuities (e.g., bedding planes, fractures, faults). Brittle response of a rock mass can take the form of intact rock fracturing (initiation and growth of new fractures), or slip or opening of pre-existing small- or large-scale fractures. The deformation mechanisms can occur in combination and affect each other. Recent studies in laboratory rock failure have documented transitions between tensile and shear failure modes (hybrid failure) [12,13], showing that hybrid failure includes interconnected tensile and shear segments. Outcrop-based investigations have recently documented similar hybrid failure in naturally deformed limestone and chalk [14,15]. Modeling rock deformation is important in understanding and controlling processes in civil, mining, petroleum, and environmental

engineering and geology [16–20]. However, available models either typically do not provide sufficient information about fracturing or the ones that can provide such information are too complicated for a typical use. Therefore, a new modeling approach is needed that is reasonably easy to use, can capture the spectrum of failure modes in rock, and provides information regarding typically measurable characteristics of deformed rock (e.g., fracture size, orientation, and density). Two basic categories of numerical models are available: continuum and discontinuum. Continuum models are based on continuum mechanics principles [10] and can be implemented in a variety of numerical techniques (e.g., finite element, finite difference, and boundary element). Typically, these models are used to calculate overall deformation and stress change (based upon specified constitutive relationships) but no specific information about fractures [16,17]. Continuum models allow understanding of patterns of generalized deformation and stress change, and are generally adequate when specific information about fractures is not needed from the calculation [21]. However, they are normally not appropriate for understanding fracture attributes because fractures are

* Corresponding author. E-mail addresses: [email protected] (G.I. Ofoegbu), [email protected] (K.J. Smart). https://doi.org/10.1016/j.rineng.2019.100070 Received 16 August 2019; Received in revised form 15 November 2019; Accepted 20 November 2019 2590-1230/© 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/).

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

not physically manifest in the model. In contrast, discontinuum models provide for representation of the geometry and mechanical behavior of discontinuities explicitly and range from models that consist of an assemblage of intact rock blocks and pre-specified discontinuities to models that represent the rock mass as a combination of bonded particle assemblies and discrete fracture networks [16,17,22–27]. Such models account for slip or opening of existing fractures but typically do not account for fracture initiation and growth. Because discontinuum models with existing fractures typically require a change in model grid for every change in fracture geometry, such models are not well suited for performing multiple analyses as would be needed to evaluate uncertainties in fracture geometry. Furthermore, discontinuum models that include bonded particles generally require extensive effort to calibrate the micromechanical properties, which can be challenging or at times impossible. This paper presents a procedure for modeling the initiation, deformation, and propagation of discrete fractures explicitly in continuum analysis; provides numerical simulations to evaluate performance of the procedure against analytical solutions and laboratory test data from literature; and discusses insights based on the simulation results for understanding rock fracturing and mechanical behavior of fractured rock. Modeling fractures explicitly in continuum analysis Fig. 1. Illustration of brick-shaped element model showing four hypothetical fractures: blue represents unfractured rock; orange, tensile fracture; green, shear fracture; and grey, pre-existing fracture. The maximum, intermediate, and minimum principal compressive stresses, σ 1 , σ 2 , and σ 3 , respectively, are indicated.

The procedure, referred to as CanFrac® (Continuum Analysis of Fracturing) [28], is based on numerical piecewise tensorial aggregation of the governing mechanical relationships for discrete fractures and unfractured rock. Thus, the procedure uses stress-strain relationships that evolve as fractures initiate and deform. In the procedure, fracture geometry is stored numerically as continuum properties used internally to generate geometrical parameters needed for analysis. Thus, while fracture geometry is not included in the model grid, the geometry is represented explicitly in modeling mechanical interactions. Representing fractures this way offers an advantage because a given model grid can be adapted to represent different statistical realizations of fracture geometry and attributes by changing continuum properties of the domain. Each fracture is modeled as segments, each with a specific location and geometrical and mechanical attributes. An elemental volume of the domain can include up to four fracture segments (referred to simply as fractures): two pre-existing fractures and two new (i.e., analysis generated) fractures. The new fractures are limited to one formed in tension and one formed in shear. The pre-existing fractures are unclassified in terms of origin but could be mechanically active or inactive. An inactive fracture does not contribute mechanically until activated, either in shear or tension. Having been formed or activated, a fracture can deform by slip or opening as dictated by the stress state. Every fracture in an element passes through the centroid, cuts through to the boundaries, and has constant orientation and aperture within the element (Fig. 1). However, fracture orientations and aperture can vary between elements. Therefore, fracture tip phenomena, propagation, and coalescence are not based on any predefined relationships but could result from interactions among sub domains as dictated by the evolving mechanical conditions. Results of numerical simulations are provided to make a case that the calculated fracture morphology (e.g., aperture decreasing toward the leading edge) and geometry of failure surfaces due to fracture coalescence are consistent with established concepts.


 q2 ¼ð1=2Þ ðσ 11  σ 22 Þ2 þðσ 22  σ 33 Þ2 þðσ 33  σ 11 Þ2 þðσ 12 Þ2 þðσ 23 Þ2 þ ðσ 31 Þ2 where σ ij is the stress tensor (tension positive convention). The yield function can be specified directly in terms of predefined q  p parameters or indirectly in terms of the Mohr-Coulomb strength parameters: unconfined compressive strength qu and friction angle φ. The latter approach (i.e., specifying the yield function in terms of Mohr-Coulomb parameters) requires an additional parameter (referred to as curvature) to define the shape of the yield surface in q  p space. An example is provided in a subsequent section to illustrate the effect of the curvature parameter on the yield function. The tensile failure criterion is defined in terms of a tensile strength Ts . The strength of a pre-existing fracture is described in terms of a cohesion cs and friction angle, for activation in shear; and a tensile strength for tensile activation. Once a fracture is generated or activated, its contribution to the overall stress-strain response depends on its orientation. Fracture orientation is defined in terms of dip and dip direction and calculated at the time of initiation for new fractures or specified with input data for preexisting fractures. The stress and strain tensors are transformed to the fracture coordinate frame to perform fracture mechanical calculations and transformed back to the general domain coordinate frame for domain-level analyses. The fracture frame calculations depend on the stress state relative to the shear strength (defined in terms of a cohesion and friction angle) and tensile strength. Fracture shear and tensile strength parameters are established based on a parameter referred to as rkDuctility that defines how the strength of intact rock degrades after fracturing (for new fractures) or strength of an inactive fracture degrades after activation (for a pre-existing fracture). The value of rkDuctility can be as large as 1.0 for a rock that retains a large fraction of its strength after fracturing or as small as 106 for a rock that loses most of its strength after fracturing. Meaningful values are in the range of 0.1–1.0 but the model procedure admits smaller values to increase user flexibility.

Fracture initiation and deformation Fracture initiation is governed by shear and tensile strength criteria for the domain. The shear failure criteria are described in terms of a generalized q ¼ qðpÞ yield function: where p and q are the confining and distortional stress intensities, defined as follows. p ¼  ð1 = 3Þðσ 11 þ σ 22 þ σ 33 Þ 2

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

dimensions from Ref. [30]. The absolute specimen dimension would be important if the numerical model included an initial distribution of non-homogeneity. However, the material was considered initially homogenous in the simulation. The specimen for the generic competent rock test has dimensions of 1  1  2:5 (m) and pre-existing fractures were specified throughout the upper two-thirds. For the Coconino sandstone model, the specimen was subjected to an initial state comprised of (1) stress equal to the minimum confining pressure in every direction, (2) boundary pressure on x -normal and y -normal surfaces equal to the confining pressure, and (3) zero z velocity at the top and bottom surfaces. Thereafter, the x -normal boundary pressure was increased in steps to the desired σ 2 value with equilibrium established after each increment. Then, a steady z velocity of 106 m/s was applied at the top surface while monitoring the specimen response. For the generic competent rock model, the planes x ¼ 0, y ¼ 0, z ¼ 0, and z ¼ 2:5 m were held at zero normal velocity; pressure of 10 and 20 MPa was applied instantaneously on the planes y ¼ 1 m and x ¼ 1 m, respectively; and the model was brought to equilibrium. Thereafter, a steady z velocity of 106 m/s was applied at the top surface while monitoring the specimen response. Material properties for the Coconino sandstone model were defined based on [30] in terms of Young's modulus E ¼ 28 GPa, Poisson's ratio ν ¼ 0:25, and tensile strength Ts ¼ 5 MPa (estimated based on the unconfined compressive strength of ~56 MPa). Reference [30] provides the shear strength envelope for Coconino Sandstone in terms of the best-fit relationship for the test data as follows.

Validation of the model procedure This paper provides calculations to make a case for satisfactory performance of CanFrac® as a tool for understanding and predicting the mechanical behavior of rock and rock-like materials. The case is based on comparing results from numerical simulations against analytical solutions, laboratory test data, and behavior from established concepts. The simulations and results are described in section Triaxial compression test of brick shaped specimen for triaxial compression of previously unfractured rock and rock with a regular pattern of pre-existing planes of weakness, in section Triaxial extension test of dog bone shaped specimen for triaxial extension test, and in section Diametral compression of rock disc – The Brazilian test for diametral compression or Brazilian test. The tests were selected to represent failure conditions driven by shear, tensile, or combinations of shear and tensile fracturing. Section Discussion discusses the simulation results as insights for understanding large-scale behavior such as localization of rock deformation, stress control of fracture mechanisms in failure zones, and conditions for occurrence of multiple parallel fractures. All analyses were performed using a CanFrac® implementation that runs as a proprietary plugin with FLAC3D [29]. Triaxial compression test of brick shaped specimen Simulations of triaxial compression tests on brick shaped rock specimens were performed and the calculated behavior is compared with laboratory data in one case and an analytical solution in another. The results from one set of simulations were compared with laboratory data from Ref. [30] for Coconino sandstone. The second set of simulations considered a generic competent rock with a set of pre-existing fractures. The results for this set were compared with an applicable analytical solution ([10], p. 106–107).

τoct ¼ 13:93 þ 0:75σ oct  0:0072σ 2oct where τoct and σ oct represent the octahedral shear and normal stresses, respectively. The shear strength was defined in terms of q  p relationships (Fig. 3). First, the best-fit shear strength envelope was translated pffiffiffi into a q  p relationship using q ¼ ð3 = 2 Þτoct and p ¼ σ oct (referred to as Ma-Haimson best fit in Fig. 3). Second, test results for σ 2 ¼ σ 3 cases from Ref. [30] were used to define an upper bound q  p relationship (Fig. 3). Two sets of simulations were performed based on the data. One set (Group 1) simulates the σ 2 ¼ σ 3 test conditions for σ 3 values of 0, 10, 20, 50, 80, 100, 120, and 150 MPa to compare with laboratory data ([30], Table 2, Fig. 8a). The other set (Group 2) simulates σ 2 > σ 3 test conditions for five ðσ 3 ; σ 2 Þ -pair values that lie close to the Ma-Haimson best fit shear strength envelope (Fig. 3). Material properties for the generic competent rock model were

Test Geometry, loading conditions, and material properties The test specimen is brick shaped with coordinate axes parallel and perpendicular to the brick edges (Fig. 2A). The minimum principal compressive stress σ 3 is applied in the y direction, the intermediate principal stress σ 2 in the x direction, and the maximum principal compressive stress σ 1 in the z direction. For the Coconino sandstone models, the specimen dimensions are 1  1  2:03 (m), which simulate the ratio of length to lateral dimension without using the exact

Fig. 2. Illustration of model configurations for simulations of (A) true triaxial compression test, (B) triaxial extension (dog bone) test, and (C) diametral compression (Brazilian test). Specific geometry description for each test is given in the section under “Test Geometry, Loading Conditions, and Material Properties.” 3

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 3. Relationships for compressive strength of Coconino sandstone showing laboratory data from Ref. [30], fitted q  p relationships for the data, data points selected for numerical simulation, and results calculated based on the simulation.

defined in terms of E ¼ 20 GPa, ν ¼ 0:3, unconfined compressive strength qu ¼ 75 MPa, friction angle φ ¼ 30 , and Ts ¼ 5 MPa. As described in Section Fracture initiation and deformation, a curvature parameter is needed to derive a q  p relationship from the MohrCoulomb shear strength specification. Fig. 4 illustrates the effects of the curvature parameter on the shear strength model. Fig. 4 indicates a cross-over p value at which curvature has no effect on strength, for fixed values of qu and φ. For conditions of low confinement (i.e., p values smaller than the cross-over value), the model shear strength increases as curvature increases. Conversely, for conditions of high confinement (i.e., p values greater than the cross-over value), the model shear strength decreases as curvature increases. Simulations with the generic competent rock model were performed using a curvature of 0.015 (Fig. 4). In addition, every zone (or element) within the upper two-thirds of the

Table 2 Results of simulated triaxial extension test compared with ref. [13] Data. Confining pressure (MPa)

10 20 30 40 50 60 80 100 120 150

Numerical simulation

Ref. [13] data

Peak σ max  σ min (MPa)

θ (deg)

Peak σ max  σ min (MPa)

θ (deg)

13.5 22.6 31.6 40.6 49.6 58.8 76.1 93.8 112.6 138.6

0.0

15.0–15.7 26.3 36.8 47.0 57.5 66.2 82.6 99.6 116.1–116.8 142.7

3.1–6.7

0.0 12.0 29–32

6.1 9.8–11.4 13.4

Table 1 Simulated compressive strength for Coconino Sandstone cases compared with ref. [30] data. σ 3 (MPa)

σ 2 (MPa)

Simulation Test Results

Ref. [30] data

Peak σ 1 (MPa)

Peak σ 1 (MPa)

θ (deg)

θ (deg)

σ 2 ¼ σ 3 test series: Upperbound strength envelop used for numerical simulation 0 10 20 50 80 100 120 150

0 10 20 50 80 100 120 150

113.4 174.4 223.8 339.3 428.2 481.3 529.9 596.3

56–62 47–51 46–52 45–51 49–56 47–52 44–50 42–49

56.1(a) 161.8 235.2 347.2 427.7 465.2 506.4 565.2

80 72 68 62 60 56 54 50

Test for (σ 3 ; σ 2 ) pairs: Best-fit strength envelop used for numerical simulation 0 20 50 120 150

0 90 160 270 360

63.3 241.3 374.3 567.1 620.2

41–52 54–57 55–57 52–56 44–52

56.1 247.3 372.8 557.1 621.5

75 73.6 69.5 52.5 50.5

a Although the calculated σ 1 does not match the test data for this series, the corresponding q  p result fits the failure envelop (Fig. 3), which indicates the upperbound strength envelop is not appropriate for modeling compression behavior under small confining pressure.

Fig. 4. Plots of q  p relationships showing the effects of curvature parameter on the model for a generic competent rock with Mohr-Coulomb parameters of qu ¼ 75 MPa and φ ¼ 30 . 4

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 5. (A) Plots of slip strain on induced shear fractures from triaxial compression test simulations on Coconino sandstone, for simulations with σ 2 ¼ σ 3 (stress level shown for each plot). Each plot shows cross section through the center of the specimen parallel to the YZ-plane. The geometric scale is consistent with Fig. 2A and dimensions given in Section Test geometry, loading conditions, and material properties. Slip strain magnitudes range from 0 (blue) to 0.0025 (red). Dip angles for selected zones of high slip strain are also illustrated. (B) Photographs of Coconino sandstone specimens tested under conventional triaxial compression, σ 2 ¼ σ 3 , with values shown under each photograph (from Ref. [30], Fig. 8a).

specimen was assigned properties representing a plane of weakness that was mechanically inactive initially with cohesive strength of 7.5 MPa, tensile strength of 1.67 MPa (approximately one-third of the intact rock values) and φ ¼ 30 (same as the intact rock value). Several analysis cases were performed, each with a different dip for the weakness planes.

state for Group 1 and Group 2 (Table 1 and Fig. 3). Second, the inclinations of macroscopic failure surfaces from Group 1 simulations were evaluated. The q  p results are reasonably close to the applicable strength envelope: Group 1 results lie on the upperbound envelope and Group 2 results on the best-fit strength envelope, thus quite consistent with the input strength criterion. Furthermore, Table 1 shows the σ 1 results are consistent with the laboratory data, except for the simulated unconfined case that used the upperbound strength criterion. Although the q  p result for this case lies on the strength envelope, the σ 1 value over-estimates the unconfined compressive strength. Therefore, the σ 1

Simulation test results from Coconino Sandstone model Results from Group 1 and Group 2 sets of simulations are compared with ref. [30] data, first in terms of σ 1 and q  p pair values at the peak 5

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Simulation test results from generic competent rock model The generic competent rock model was used to calculate the effects of a regular set of pre-existing planes of weakness on compressive strength. The analytical solution for the case is as follows ([10], p. 106–107).

σ peak  σ 3 ¼

2ðcs þ μw σ 3 Þ ð1  μw cot βÞsin 2 β

where σ peak is the compressive strength, σ 3 the confining pressure, μw the friction coefficient of the weakness plane, and β is the angle between the maximum principal compressive stress σ 1 and the normal to the weakness plane. If σ 1 is vertical, then β ¼ θ, where θ is the dip of the weakness plane. The simulation results match the analytical solution (Fig. 7) for weakness planes with θ in the range 40  θ  80 . For this range, rock behavior is dominated by shear reactivation of pre-existing fractures. If the weakness plane dip is outside this range, behavior of the rock is dominated by intact rock failure. The simulation results differ from the analytical solution for dip values outside of the specified range, because the failure mechanism assumed for the analytical solution (i.e., slip on pre-existing fractures) is not applicable for small dips (e.g., θ < 40 ) or large dips (e.g., θ > 80 ). Thus, the analytical solution predicts infinite compressive strength whereas the simulation predicts strength equal to

Fig. 6. Axial stress increment σ 1  σ 3 versus shortening from the σ 2 ¼ σ 3 (shown for each plot in MPa) simulations with the Coconino sandstone model.

result for this case indicates the upperbound strength envelope is not appropriate for simulating low confinement conditions. However, simulations based on the best fit strength envelope (Group 2 simulations) give results consistent with the data in terms of q  p and σ 1 . The calculated failure surface inclinations from the simulations (Table 1) do not match the experimental data exactly, but there is general agreement. Plots of failure surfaces based on distributions of slip strain on induced shear fractures (Fig. 5A) are similar to failure surface inclinations from laboratory data (Fig. 5B). The simulation results are consistent with the experimental results that the failure angle decreases as the confining pressure increases. The calculated failure surface for the cases with σ 3 values of 80, 100, 120, and 150 MPa are flatter than Fig. 5A suggests, because a through-going failure surface developed first near the top right corner of the specimen but numerically stable post-peak deformation continued thereafter and produced subsequent and steeper “failure surfaces.” Furthermore, the loaddeformation plots from the simulations (Fig. 6) show decreasing slope of load versus deformation, i.e., increasing ductility, as confining pressure increases. The increasing ductility is consistent with the decreasing steepness of failure surfaces in Fig. 5B. Therefore, the calculated failure surface inclinations are generally consistent with laboratory data, more than the numerical comparison in Table 1 indicates.

Fig. 8. (A) Plots of differential stress σmax  σmin versus axial extension strain from triaxial extension test simulations using Berea sandstone model at different confining pressures (shown in the plots in MPa). (B) Stress-strain curves for Berea Sandstone based on laboratory triaxial extension tests (from Ref. [13], Fig. 5). The confining pressure for each test (MPa) is shown on the plot.

Fig. 7. Compressive strength versus dip of weakness plane from simulations with models of the generic competent rock with a regular set of pre-existing planes of weakness and confining pressure of 10 MPa. 6

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

stresses σ max and σ min were calculated as average values for a 5-mm thick disc axially centered in the specimen, representing the neck and simulating a calculation of similar quantities in Ref. [13]. The simulation results (Fig. 8A) compare well with corresponding data (Fig. 8B). Peak values of σ max  σ min from the simulation agree well with the data (Table 2). Also, the post-peak behavior shows greater decrease in stress as confining pressure increases, for both the simulation and laboratory data. However, the calculated stress-strain curves differ from the laboratory data with respect to linearity of the pre-peak behavior. The laboratory pre-peak stress-strain curves (Fig. 8B) show an initial stiffening, interpreted as due to microcracks closing as pressure increased. In contrast, the simulation pre-peak stress-strain curves (Fig. 8A) are linear, which reflects the constant elastic stiffness used in CanFrac® for intact rock with an absence of pre-existing microcracks.

the compressive strength of unfractured rock for small and large dips. Triaxial extension test of dog-bone shaped specimen Triaxial extension tests on dog-bone shaped rock specimens were simulated to compare with independent laboratory test data [12,13]. The laboratory tests show failure surfaces with orientation and morphology that range from characteristics associated with tensile fracturing to those associated with shear fracturing. Ref. [12] interpreted the results as indicating “hybrid fractures” that display characteristics of both tensile and shear fractures but closer to shear fractures as compressive stress increased. CanFrac® does not include any specific “hybrid fracturing” mechanism. Instead, the procedure provides for tensile and shear fractures occurring individually or in combination as determined by the stress state and material properties. The individual fractures coalesce to produce through-going failure surfaces or composite fractures with orientation and morphology that depend on the overall stress and domain boundary conditions. Calculated stress conditions at failure are shown to compare well with the laboratory data. Also, failure stress conditions and interpreted failure surfaces from the simulation are shown to indicate a change in fracture mechanism from tensile fracture domination at low confining pressure, through varying combinations of tensile and shear fractures as confining pressure increases, and shear fracture domination at high confining pressure. The progression of fracturing mechanism with confining pressure indicated by the simulation is consistent with the laboratory data.

Strength envelope The predicted tensile and shear strength based on the simulation are compared with data through plots of peak σ max and σ min (Fig. 9). Both the simulation and laboratory data indicate a tensile strength of e5 MPa based on value of σ min for σ max ¼ 0. For greater values of σ max , the lab data shows tensile resistance increasing to a maximum value and then decreasing as the strength curve transitions to the shear strength criterion. The increase in tensile resistance with increasing σ max in the laboratory data could be due to closure of microcracks, which is consistent with the stiffening of the stress-strain curve as stated earlier. In contrast, the simulation results don't show any increase in tensile resistance. Instead, the results show σ min decreasing asymptotically toward zero as σ max approaches the unconfined compressive strength. Also, the simulation results show that values of σ max greater than the unconfined compressive strength follow a relationship with σ min consistent with the shear strength criterion. That is, the calculated rock strength values are consistent with data for shear behavior but do not match the data for tensile behavior because of the effect of microcracks on tensile resistance. The simulation could conceivably match the tensile behavior better if pre-existing microcracks were represented in the model. A pre-existing fracture distribution could be generated using a statistical model of the microcracks and included in the simulation as part of the input data. The authors expect CanFrac® would perform well with such data and provide results much closer to the laboratory data.

Test Geometry, loading conditions, and material properties The simulations were performed using a dog-bone shaped specimen (Fig. 2B) and were based on laboratory experiments on Berea sandstone conducted at Texas A&M University [13]. The specimen is approximately 10 cm high with a large (shoulder) diameter of approximately 5 cm and a small (neck) diameter of approximately 3 cm. As described in Ref. [13], the radius of curvature of the sample neck is approximately 9 cm. In the laboratory experiment, the specimen was first subjected to equal pressure on all surfaces, with the base restrained axially. Thereafter, the pressure at the specimen top surface was reduced gradually until the specimen failed. For the numerical simulation, static equilibrium was first established with the specimen base held at zero axial velocity and equal pressure (i.e., confining pressure) applied on all other surfaces. Thereafter, the pressure at the top surface was reduced in small increments and static equilibrium re-established with each reduction. When equilibrium could not be re-established – implying that the specified loading was approximately at the failure state – the calculation was continued without changing any applied pressure and the stress-deformation configuration was monitored to a chosen end condition. The procedure was performed for several values of confining pressure ranging from 10 to 170 MPa. Mechanical property values for the simulation were qu ¼ 98:7 MPa and φ ¼ 47 ([13], p. 20); and E ¼ 19 GPa, ν ¼ 0:3, and Ts ¼ 5 MPa, based on general literature regarding Berea sandstone [31–36]. The curvature parameter, needed to derive a q  p relationship from the Mohr-Coulomb strength specification, was assigned a value of 0.015 (e.g., Fig. 4).

Geometry of failure surfaces The geometry of failure surfaces from the simulations (Fig. 10) is compared with similar information based on lab data. The failure surface for the simulation case with 10 MPa confining pressure consists of a

Simulation test results Load versus deformation The calculated load-deformation response is described in terms of plots of differential stress (σ max  σ min ) versus axial extension strain (Δl= l0 ), where σ max and σ min are the maximum and minimum principal compressive stresses, l0 is the initial (zero strain state) length of the specimen, and Δl is the change in specimen length relative to the state of initial static equilibrium under the confining pressure. The principal

Fig. 9. Values of σ max and σ min representing the peak state of the loaddeformation curve from the triaxial extension tests, showing the results of model simulations and Bobich (2005) [13] data. 7

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 10. Distributions of induced fractures showing number and type of fractures at calculation points: from triaxial extension test simulations at confining pressures of 10, 80, 120, 150, and 170 MPa (shown at the top of each plot). Blue indicates unfractured rock; orange, rock with a single tensile fracture; green, rock with a single shear fracture; and red, rock with two fractures (one tensile and one shear). Each plot represents a vertical section through Fig. 2b. Dark lines labeled with degree value indicate authors' interpretation to estimate the inclination of the coalesced failure surface.

t (usually t ¼ D=2). A good review of the test is provided in Ref. [37]. The theoretical basis for the test (e.g. Ref. [10], p. 169) is that diametral compression causes tensile stress of magnitude 2P=π Dt normal to the loaded diametral plane, such that the tensile strength of the rock can be calculated using this relationship if P is the failure load. If the rock is homogenous and isotropic, and the load is applied over a short arc (with arc angle 2α  15 ), then failure should occur through a tensile fracture parallel to the applied load and coplanar with the loaded diametral plane. The fracture initiates at the disc center and propagates toward the loading surface. If the specimen is non-homogenous or non-isotropic or the load application arc is too wide, then failure could initiate at the circumference and propagate toward the interior and the failure load will differ from the assumed relationship with tensile strength. The Brazilian test is well suited for evaluating the CanFrac® procedure because the model output includes whether a fracture has occurred, the type of fracture (tensile or shear) and its location. Thus, the results of a Brazilian test simulation can be evaluated with respect to the type of fracture predicted, how the fracture initiates and propagates, and the failure load relative to the input tensile strength. Two sets of numerical simulation of the Brazilian test were performed: one set on a homogenous and isotropic rock model and the other on a model that includes a regular set of pre-existing planes of weakness. The simulations and evaluation of results are described in this section.

horizontal plane formed from coalescence of coplanar tensile fractures. Failure was also driven by tensile fracturing in the 80 MPa case. However, the tensile fractures for this case are not all coplanar. Although there is a through-going plane of tensile fractures at mid height of the specimen, a series of short non-coplanar but parallel tensile fractures also formed near the specimen outer surface above and below the mid-height fracture. The short fractures could contribute to the eventual failure surface. Thus, the overall failure surface in the 80 MPa case could consist of a composite formed from coalescence of non-coplanar tensile fractures and could be inclined near the specimen outer surface and horizontal in the specimen interior. For the 120 MPa case, both shear and tensile fractures developed near the specimen outer surface but only tensile fractures in the specimen interior. Thus, the failure surface could consist of a composite of shear fractures near the specimen outer surface and tensile fractures in the specimen interior. The composite failure surface would be inclined near the outer surface and horizontal in the interior. For the 150 and 170 MPa cases, failure was driven by shear fractures that propagated around an unfractured interior. Tensile fractures also occurred around the edges of the unfractured interior, aiding separation of the shear-fractured zone from the unfractured interior. The geometry of the failure surface is influenced by the size of the unfractured interior. The overall failure surfaces in the high confining pressure cases are steeper than the failure surfaces in the lower confining pressure cases. However, the failure surface is steeper in the 150 MPa case than the 170 MPa case, because the unfractured interior is smaller in the 170 MPa case than in the 150 MPa case. Generally, failure surfaces from the numerical simulation are steeper as confining pressure increases and failure transitions from tensilefracture driven to shear-fracture driven. Therefore, the effects of confining pressure on the failure surface geometry is generally consistent with data in Refs. [12,13]. However, the numerical simulation predicts another change in trend at high confining pressure because the failure surface deflects around a core of unfractured rock that decreases in size, thus reducing the failure surface inclination, as confining pressure increases.

Test Geometry, loading conditions, and material properties The simulations were performed on a model disc with D ¼ 50 mm, t ¼ 25 mm, and 2α ¼ 10:1 oriented as shown in Fig. 2C. The two ends of the disc (y ¼ 0 and y ¼ 0:025 m) were held at zero y -velocity. Load was applied vertically downward (i.e., z direction) over a narrow strip at the top surface, bisected by the z axial plane such that the arc angle is 10.1 . Vertical support was applied as zero z -velocity over a strip of equal width centered at the opposite end of the z axial plane (i.e., the specimen base surface). The test sequence started with static equilibrium under initial all-round compressive stress of 0.001 MPa and equal normal stress on the specimen cylindrical surface. Thereafter, vertical load on the top strip was increased in steps and static equilibrium re-established after each increment. When equilibrium could not be re-established – implying that the applied loading was approximately at the failure state – the calculation was continued to a chosen end condition without changing the applied load. The stress-deformation configuration was monitored

Diametral compression of rock disc – The Brazilian test The tensile strength of rock can be determined in the laboratory through diametral compression of a rock disc of diameter D and thickness 8

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 11. Plot of vertical load versus downward displacement from Brazilian test simulations for intact rock and rock with regular set of pre-existing weakness planes. The weakness planes have 75% intact rock strength in some analyses and 50% in others. One plot shows annotations to identify mechanical states described in Fig. 12.

were assigned strength equal to 50% of the intact rock strength (i.e., cs ¼ 3:31 MPa, φ ¼ 23:9 , and Ts ¼ 1:25 MPa).

through the simulation. The monitored configuration includes vertical load based on the summation of vertical reaction at the support base. For the second set of simulations, the model included specifications for a regular set of pre-existing planes of weakness everywhere. The weakness planes for each model case have a constant dip in the þ x direction. Simulations were performed for several dip magnitudes: 0, 10, 20, 30, 40, 60, 70, 80, and 90 . A dip magnitude of 0 represents planes of weakness normal to the loading direction, whereas a dip magnitude of 90 represents planes of weakness parallel to the loading direction (consistent with laboratory specimen orientations in Ref. [38]). Material properties for the homogenous and isotropic model represent a generic sandstone with E ¼ 17 GPa, ν ¼ 0:25, qu ¼ 29:5 MPa, φ ¼ 42 , and Ts ¼ 2:5 MPa. The curvature parameter was set to 0.015. Two simulation subsets were performed with the model that includes planes of weakness. In one subset, the planes of weakness were assigned strength equal to 75% of the intact rock strength (i.e., cs ¼ 4:97 MPa, φ ¼ 33:6 , and Ts ¼ 1:875 MPa). In the other subset, the planes of weakness

Simulation test results for homogenous and isotropic model The calculated vertical load (based on support reaction) versus deformation gives a failure load of 5 kN (Fig. 11, case of “No weakness plane”), which implies a tensile strength of 2.55 MPa that compares well with the input tensile strength of 2.5 MPa. The waviness of the pre-peak load-deformation curve is because the plot includes the “march to equilibrium” after each load increment. A straight line would be obtained by plotting only the equilibrium states (highest point of each curve segment). Results show rock failure was driven by tensile fracturing. The tensile fracture initiated on the loaded diametral plane, midway between the load and support surfaces (Fig. 12A). Fracture initiation and propagation in the simulated test can be described as follows:

Fig. 12. Fracture initiation and propagation in simulated Brazilian test for model case with no planes of weakness, showing zones of fractured rock. The exterior outline of the specimen disc is shown, but unfractured rock is excluded from the plots (except for a rectangular plate at each disc end) to facilitate visualization of the fracture zones. Blue indicates unfractured rock; orange, rock with a single tensile fracture; green, rock with a single shear fracture; and red, rock with two fractures (one tensile and one shear). The mechanical state for each plot is annotated on the load-deformation curve (Fig. 11). Plot A represents fracture initiation, while plots B, C, and D represent progressive stages in the fracture propagation. The main fracture is contained in the disc interior in plot A, visible at the disc ends but not at the cylindrical surface in plots B and C, and has propagated to the cylindrical surface in plot D. Auxiliary fractures that initiated at the cylindrical surface are also shown in plots C and D. 9

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

fail before their development. However, the numerical simulation is able to show extended post-peak behavior, thus indicating how rock fracturing might progress if the rock volume was restrained enough to sustain fractures without collapse. Fig. 12D shows several fractures that initiated near the boundary of the disturbed zone and propagated toward the disturbance. In this example, the disturbed zone is defined by stress change around the loaded diametral plane. Several auxiliary fractures initiated near the boundary of the disturbed zone (cylindrical surface) and propagated toward the disturbed zone (diametral plane). Fracture aperture The fracture zones plotted in Fig. 12 represent model sub-domains that contain fractures but the fracture itself is much thinner than the sub-domain. CanFrac® assigns an arbitrary aperture (0.1 mm in the example) to a fracture segment when the fracture is generated. Thereafter, the fracture aperture may increase or decrease, limited to a minimum equal to the initial value, as determined by the mechanical condition. The aperture of the main tensile fracture at the end of the simulation had increased to a maximum of ~0.2 mm at about mid height of the fracture (Fig. 13). The fracture aperture varies laterally and vertically with minimum values at the top and base.

Fig. 13. Fracture aperture distribution for the main fracture shown in Fig. 12D. Aperture varies from  0:1 mm (blue) to  0:2 mm (red).

 Tensile fractures initiate at two locations near the disc ends on the loaded diametral plane, midway between the load and support boundaries. The fractures propagate upwards and downwards, and horizontally toward the interior with an approximately elliptical front (Fig. 12A).  The two fractures coalesce into a single tensile fracture zone that extends horizontally across the specimen and propagates upward and downward with approximately rectangular fronts (Fig. 12B).  The main zone of tensile fracturing reaches a maximum height, but still a finite distance below the loading surface and above the support. Secondary tensile fractures initiate on the cylindrical surface of the specimen, and include two short and isolated fractures near the upper loading surface and two horizontally through-going but vertically short fractures near the lower support surface (Fig. 12C).  Shear fractures form above the main tensile fracture zone near the upper loading surface and, to a lesser extent, below secondary tensile fractures near the lower support. The shear fractures near the loading surface coalesce with the tensile fractures to form a through-going composite failure surface. Several secondary tensile fractures initiate on the cylindrical surface and propagate parallel to the main fracture initially, but then turn to propagate normal to and toward the main fracture zone (Fig. 12D). The occurrence of shear fractures that connect the primary tensile fracturing to the loading surface is consistent with the high compressive stress within close proximity of the loading surface. The simulation predicts that overall specimen failure is driven by the primary tensile fracture that is coplanar with the loaded diametral plane, aided by shear fracturing near the load application boundary. Thus, the failure surface is a composite feature that formed of tensile fractures over most of the surface and shear fractures over parts of the surface in high compressive stress conditions. The tensile fracture initiated midway between the load and support boundaries and propagated upward and downward toward the boundaries. Several of the secondary tensile fractures shown in Fig. 12D likely would not materialize in a laboratory test because the specimen would

Fig. 14. (A) Failure load versus dip of weakness planes based on simulated Brazilian tests. The failure load is normalized relative to the value for a weakness plane dip of 0 (5 kN for weakness planes with 75% intact rock strength and 2.8 kN for weakness planes with 50% intact rock strength). (B) Variation of average failure load, relative to the failure load for loading perpendicular to weak planes, showing cases with decreasing strength as weakness plane dip increases, based on laboratory Brazilian testing of Freiberger Gneiss, Mosel Slate and Herbeumont Slate (from Ref. [38], Fig. 2d). 10

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 15. Fracture initiation and propagation in simulated Brazilian test for model case with weakness planes of 90 dip and 75% intact rock strength, showing zones of fractured rock. The exterior outline of the specimen disc is shown, but unfractured rock is excluded from the plots (except for a rectangular plate at each disc end) to facilitate visualization of the fracture zones. Numbered colors represent fracture type and number of fractures as follows: 0 ¼ no fracture, 1 ¼ new shear fracture only, 2 ¼ new tensile fracture only, 3 ¼ two new fractures (one shear and one tensile), 4 ¼ shear activation only, 5 ¼ tensile activation only, 6 ¼ shear activation plus new fracture, and 7 ¼ tensile activation plus new fracture.

Fracturing in case of weakness plane with 75% intact strength and dip of 90 In the model scenario with 90 dipping weakness planes with 75% intact rock strength, failure initiated with shear activation of weakness planes near the support base and propagated through tensile activation of pre-existing fractures coplanar with the loaded diametral plane (Fig. 15). Evolution of fracturing in this case can be described as follows:

The calculated fracture aperture is largest at about mid height where the fracture initiated and smallest at the top and base representing the fracture tip (or front). This shape of the fracture aperture distribution can be shown to be conceptually consistent with the analytical solution for a penny-shaped crack pressurized internally. As stated in Section Modeling fractures explicitly in continuum analysis, fracture tip phenomenon is not implemented explicitly in CanFrac® but results from mechanical interactions among neighboring sub-domains.

 Two zones of pre-existing fractures activate in shear at the base of the specimen near the support boundary, one on either side of (but not coplanar with) the loaded diametral plane (Fig. 15A).  New tensile fractures form at the upper edge and on the sides of the zones of shear-activated fractures. Additionally, pre-existing fractures activate in tension between and above the shear-activated zones. The tensile-activated zone is approximately coplanar with the loaded diametral plane, initiated near one end of the disc, and propagates upward and laterally toward the other disc end with an approximately elliptical front (Fig. 15B).  The shear-activated zones appear to have terminated but the tensileactivation zone continues to propagate upward and laterally and is coplanar with the loaded diametral plane. Coincident new fractures have appeared over part of the surface of the tensile-activation zone (Fig. 15C).  The zone of tensile activation appears to have attained maximum height and width. Most of the tensile-activated fracture zone is covered with coincident new fractures, but the leading (upper) edge is formed of activated tensile fractures, suggesting that tensile activation of weakness planes is the driving mechanism (Fig. 15D).  Two zones of shear activation develop above and on either side of the zone of tensile activation. The shear activation zones at the specimen

Simulation test results for model with pre-existing planes of weakness The calculated vertical load versus deformation plots (Fig. 11) give a maximum failure load of 5 kN for cases of weakness planes with 75% intact rock strength and 2.8 kN for cases of weakness planes with 50% intact rock strength. The maximum failure load in each subset was obtained for the case of weakness planes normal to the direction of loading (i.e., dip of 0 ). The failure load of 5 kN for the case of weakness planes with 75% intact rock strength and 0 dip is the same as the failure load for the homogenous isotropic rock model. For each model subset (75% or 50% intact rock strength), the failure load decreases as the dip of weakness planes increases (Fig. 14A). The rate of strength decrease is greater for dip in the range of 0  θ  θb and smaller for dip in the range θb  θ  90 , where θb is ~10 for cases with 50% intact rock strength and ~20 for cases with 75% intact rock strength. This trend in the failure load compares well with laboratory data (Fig. 14B). The simulation results indicate the dominant failure mechanism varies with the weakness plane dip. Additional results are provided to describe fracture initiation and propagation for the case of weakness planes with 75% intact rock strength and dips of 90 (Fig. 15) and 60 (Fig. 16). 11

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 16. Fracture initiation and propagation in simulated Brazilian test for model case with weakness planes of 60 dip and 75% intact rock strength, showing zones of fractured rock. The exterior outline of the specimen disc is shown, but unfractured rock is excluded from the plots (except for a rectangular plate at each disc end) to facilitate visualization of the fracture zones. Color legend is the same as explained for Fig. 15.

Fig. 17. Failure surfaces from simulated Brazilian test for model cases including planes of weakness with 75% intact rock strength. The weakness plane dip is shown at the top of the plot. Each plot shows the center cross section normal to the disc axis (XZ cross section in terms of Fig. 2C) and shows the zones of fractured rock using the color legend explained in Fig. 15.

12

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 18. Failure surfaces from simulated Brazilian test for model cases including planes of weakness with 50% intact rock strength. The weakness plane dip is shown at the top of the plot. Each plot shows the center cross section normal to the disc axis (XZ cross section in terms of Fig. 2C) and shows the zones of fractured rock using the color legend explained in Fig. 15.

tively relative to Ref. [38] laboratory data. The comparison is only qualitative because the input for the simulations could not be related to Ref. [38] data. As discussed in Section Fracture initiation and deformation, CanFrac® needs the strength of a pre-existing fracture be specified in terms of tensile strength Ts , cohesion cs , and friction angle φ. In contrast, Ref. [38] defines the strength of weakness planes (or laminations) in terms of an anisotropy ratio (specimen strength parallel to lamination divided by strength normal to lamination). However, based on the relationship of normalized failure load versus dip of weakness planes (Fig. 14A) and comparing the relationship with laboratory data (Fig. 14B), the simulation material appears closest to the behavior of Freiberger Gneiss [38]. Therefore, failure surfaces from the simulation are compared with reference [38] laboratory results for Freiberger Gneiss. The simulation failure surfaces (Fig. 17 for weakness planes with 75% intact strength and Fig. 18 for weakness planes with 50% intact strength) indicate the following: (1) the failure surfaces are parallel to loading, irrespective of weakness plane dip and strength; and (2) the location of a failure surface relative to the loaded diametral plane is influenced by the weakness plane dip and strength. The simulation results are consistent with the laboratory data (Fig. 19) in regard to item (1), including the fact that the direction of the failure surface changes near the specimen top and base. Also, the simulation results can be shown to be consistent with the laboratory data in regard to item (2). However, the comparison regarding item (2) is more subjective. The failure surfaces shown in Figs. 17 and 18 for weakness plane dips of 0 , 30 , 60 , and 90 indicate features similar to corresponding failure surfaces in Fig. 19. For the 0 case as an example, the model for weakness plane with 75% intact rock strength (Fig. 17) shows a vertical diametral failure surface driven by new tensile fracturing, superimposed tensile activation of horizontal weakness planes, and several auxiliary tensile fractures. The horizontal auxiliary fractures were driven by tensile activation with superimposed new fractures, whereas the vertical auxiliary fractures were driven by new tensile fracturing. For the same 0 case, the model for weakness plane with 50% intact rock strength (Fig. 18) shows a generally vertical fracture offset from the diametral plane, with a slight bow in the middle. The fracture morphology shown

base do not grow further but several new tensile fractures form around them (Fig. 15E).  A shear fracture zone develops to connect the main diametral-plane fracture (zone of tensile activation) to the loading surface. The shear zone consists of a combination of shear-activated weakness planes and new shear fractures. Two other tensile fracture zones develop below the shear zone on either side of the diametral-plane fracture zone, parallel to but non-coplanar with the diametral plane, and propagate downward with a rectangular leading edge of tensile-activation (Fig. 15F). The results (Fig. 15) show that fracturing in this case was driven by tensile activation of pre-existing weakness planes along the loaded diametral plane. The failure surface is composite, formed mostly of tensile mechanisms but of shear mechanisms in areas of high compressive stress. Failure initiated in shear near the specimen base, propagated by tension through most of the specimen height, and completed by shear near the specimen top. Fracturing in case of weakness plane with 75% intact strength and dip of 60 In the Brazilian test simulation for this case, fracturing initiated at the specimen top with shear activation of weakness planes near the loading surface (Fig. 16A) and propagated downward through the specimen as a zone of new tensile fractures (Fig. 16A through 16F). The zone of new tensile fractures is coplanar with the loaded diametral plane and parallel to the applied loading, therefore, cuts across the pre-existing weakness planes. The failure surface is composite: formed by shear fracturing in the zone of high compressive stress near the load application boundary, propagated by tension through most of the specimen, and by combination of shear and tension near the support base. The failure surface consists of multiple branches near the specimen top and base (Fig. 16F). Evaluation of Brazilian test simulation results relative to ref. [38] Lab test results Results from the Brazilian test simulations can be evaluated qualita13

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

from the laboratory test for the 0 case (Fig. 19a) indicates features similar to shown in Figs. 17 and 18; such as the middle bow, horizontal auxiliary fracture, and generally vertical failure surface.

understanding large-scale behavior such as localization of rock deformation, stress control of fracture mechanisms in failure zones, and conditions for occurrence of multiple parallel fractures.

Discussion Localization of rock deformation Modeling the mechanical behavior of rock by aggregating the behavior of evolving discrete fractures and unfractured rock is important because fracture deformation and the contribution of fractures to mechanical response are essential to understanding rock behavior. Such modeling offers the capability to account for the behavior of individual fractures and calculate their location and geometrical attributes as part of the output information. The modeling would be almost impractical if the fracture geometry had to be discretized with the domain grid. However, the procedure described in this paper accomplishes the modeling in continuum analysis, thus accounting for fractures without including them explicitly in the domain discretization. The resulting procedure offers promise in expanding understanding of the mechanical behavior of rocks and capability to obtain information about fracturing such as when and where fractures may occur; what type of fractures (new tensile or shear, or activation of pre-existing fracture in tension or shear); and fracture geometrical attributes (dip, dip direction, and aperture) that could be a great benefit in assessing their effects on hydrological characteristics. Additionally, the procedure shows promise in predicting and

The calculated load-deformation responses from triaxial compression and extension and Brazilian tests show elastic deformation up to a peak stress state followed by deformation at decreasing load-bearing capacity (Figs. 6, 8 and 11). The decrease in load-bearing capacity is a consequence of fracturing, because a part of the rock strength is lost at the surface of each fracture when the fracture is generated. Furthermore, the rock becomes more deformable parallel and normal to the fracture at the fracture site. Therefore, unfractured rock adjoining the fracture is stressed more and is more likely to fail. Failure of the adjoining rock follows a pattern determined by the deformation mode of the existing fracture. For example, adjoining rock could fail in tension or shear due to slip on the existing fracture. Thus, new fractures form and join with existing fractures to form an expanding failure zone. As several fractures coalesce to form a failure surface, the rock strength and stiffness decrease along the failure surface. Thus, the individual fractures combine to produce the over-all decrease in load-bearing capacity. The rock weakens and softens in preferential directions (normal and parallel to the individual fractures) at the fractured locations but retains full strength and stiffness at other locations.

Fig. 19. Pictures of Freiberger Gneiss specimens after laboratory Brazilian testing. The specimens before testing had weakness plane angles of (a) 0 , (b) 30 , (c) 60 , and (d) 90 (from Ref. [38], Fig. 7). 14

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 20. Distributions of induced fractures showing number and type of fractures at calculation points: for triaxial compression test simulations using the generic competent rock model. Confining pressures are: for plot A, 1 kPa in Y and 2 kPa in X (essentially unconfined); and for plot B, 10 MPa in Y and 20 MPa in X. Blue indicates unfractured rock; orange, rock with a single tensile fracture; green, rock with a single shear fracture; and red, rock with two fractures (one tensile and one shear).

Similarly, distributions of tensile fracture aperture from two Coconino sandstone model simulations show more tensile fractures and greater fracture aperture for the simulation under zero confining pressure (Fig. 21A) than for the simulation under 80 MPa confining pressure (Fig. 21B). Failure surfaces in brittle rock generally form by tensile and shear fracture mechanisms. Tensile fractures are dominant for failure under tensile loading or low confinement compression loading. For failure under compressive loading, shear fractures are increasingly dominant as confining pressure increases.

Localized weakening and softening in preferential directions result in localized and directional failure. For example, Fig. 5 illustrates the localization as discrete zones of relatively high inelastic deformation accommodated through fracture slip. Such discrete zones could manifest in a rock body as parallel fractures or failure surfaces. Stress control of fracture mechanisms Rock can deform inelastically through tensile or shear fracturing depending on the stress state. The initial induced fracture is tensile if the stress state is dominated by tension (e.g., if the maximum principal stress is tensile) or shear if the stress state is dominated by compression. Subsequent fractures within proximity of the initial fracture could be tensile or shear depending on the fracture deformation mode and prevailing stress. For example, fracture slip could result in adjoining tensile fracturing if the confining stress is small or shear fracturing if the confining stress is high. Generally, fracture slip results in both tensile and shear fractures in adjoining rock, the proportion of shear fractures increasing as the confining pressure increases. Therefore, a failure surface formed under combined tensile and compressive stress conditions generally consists of coalesced tensile and shear fractures. The fractures are dominantly tensile if the stress state is tension-dominated (maximum principal stress is tensile). For a compression-dominant stress state (maximum principal stress is compressive), the proportion of shear fractures increases as confining pressure increases. For example, compare the distribution of induced fractures from simulated triaxial compression of a generic competent rock model under conditions of low and high confinement (Fig. 20). The failure surfaces are driven by shear fractures in both cases. However, whereas failure surfaces in the low confinement test are each surrounded by tensile fracture zones (Fig. 20A), failure surfaces in the high confinement test are sheardominated with sporadic occurrence of tensile fractures (Fig. 20B).

Conditions for occurrence of multiple parallel fractures The simulation results described in this paper suggest an understanding of conditions for development of multiple parallel fracture (or failure) zones in a rock body subjected to external loading. The simulated triaxial compression resulted in multiple parallel shear zones (e.g., Fig. 5), triaxial extension resulted in multiple parallel tensile fractures for simulated tests under intermediate confinement (e.g., Fig. 10 cases for 80 and 120 MPa), and the simulated Brazilian tests resulted in multiple parallel fractures in all cases (e.g., Figs. 12, 15 and 16). In the triaxial extension and Brazilian test cases, the main tensile fracture could not daylight at the model boundary because of an intervening zone of compressive stress that formed a barrier against propagation of the tensile fracture. In the Brazilian test cases, the barrier is a zone of high compressive stress around the loading and support boundaries (Figs. 12, 15 and 16). In the triaxial extension cases (Fig. 10), tensile fractures in the model interior could not propagate to the exterior surface because of a zone of high compressive stress around the model “neck.” In the triaxial compression tests, a shear zone can traverse the entire model but stress relief due to the shear fractures was incomplete because of a small retention of cohesion on shear fractures. The retention was implemented using the parameter rkDuctility to improve numerical stability. Because the applied loading was, thus, not completely relieved 15

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070

Fig. 21. Aperture of tensile fractures from simulated triaxial compression using Coconino sandstone model under σ 2 ¼ σ 3 stress states and with (A) σ 3 ¼ 0 and (B) σ 3 ¼ 80 MPa. Each plot shows cross section through the center of the specimen parallel to the YZ-plane. The geometric scale is consistent with Fig. 2A and dimensions given in Section Test geometry, loading conditions, and material properties. Aperture magnitudes range from  0:1 mm (blue) to  1 mm (red).

Conclusions

through shear failure, multiple parallel failure surfaces occurred instead of a single failure surface as often observed in laboratory specimens (e.g., Fig. 5B). The model simulation differs from the laboratory test in this regard, but does predict a behavior that could explain field observations. Fractures in situ often occur in sets of approximately parallel members, similar to the numerical simulation results. The simulation test results suggest the following conditions as necessary for development of multiple parallel fractures or fracture zones.

Results of numerical simulations performed using CanFrac® are described and evaluated in comparison with laboratory data, analytical solutions, and established concepts. The evaluations lead to the following conclusions. 1. Values of compressive strength based on the simulations are consistent with laboratory data in terms of either q  p relationships or maximum principal compressive stress at the failure state. 2. The simulations predict results consistent with lab data regarding the effects of confinement on failure surface inclinations from triaxial compression and extension loading conditions, though specific inclinations from the simulation may differ from the equivalent laboratory test. 3. Simulation results match the analytical solution for the compressive strength of rock with a regular set of pre-existing weakness planes, except where the driving failure mechanism assumed for the analytical solution is not applicable. 4. The simulation results for peak strength under triaxial extension are consistent with laboratory data for stress conditions where shear fracturing is dominant. However, for tension-dominant conditions, peak strength from the simulations differs from laboratory data because of the effects of microcracks on tensile resistance. The laboratory test data apparently was affected by increase in elastic stiffness and tensile resistance due to microcracks closing as pressure increases under elastic states. CanFrac® could conceivably account for similar effects if pre-existing microcracks are included in the model. 5. The values of failure load from simulated Brazilian test are consistent with the input tensile strength for the case of isotropic and homogenous rock. Also, for rock with a regular set of pre-existing weakness

1. Persistent external loading. 2. The fracture or fracture zone resulting from the driving failure mechanism cannot daylight from the rock body. Thus, the failure zone does not function as a perfect stress relief mechanism. A structure is a perfect stress relief mechanism if it deforms without transmitting stress to its neighborhood, such as a slip surface would under residual strength conditions. 3. The impediment to daylighting could consist of stress or material conditions that do not satisfy the failure criterion, such as a zone of high compressive stress that impedes propagation of a tensile fracture. 4. If daylighting of the failure zone is not impeded, but material behavior within the failure zone permits stress transmission to the neighborhood as the zone deforms, then additional parallel failure zones or fractures are likely. Multiple parallel fractures or fracture zones will likely develop if conditions 1–3 or 1, 2, and 4 are satisfied.

16

G.I. Ofoegbu, K.J. Smart

6.

7.

8.

9.

Results in Engineering 4 (2019) 100070

planes, the variation of failure load with weakness plane dip based on the simulations is consistent with laboratory data. The fracture aperture distribution from the simulated Brazilian test on isotropic homogenous rock is conceptually consistent with the analytical solution for a penny-shaped crack pressurized internally. The failure surface based on the Brazilian test simulations is consistent with laboratory data (e.g., Figs. 17–19) in being generally parallel to the applied loading irrespective of the occurrence and orientation of pre-existing planes of weakness. For rock with pre-existing weakness planes, the Brazilian test simulations indicate that failure initiates by shear activation of the preexisting weakness near the load application or support boundary and propagates by tensile fracturing through most of the specimen. For the case of isotropic and homogenous rock, the Brazilian test simulations indicate that failure initiates in tension on the loaded diametral plane, midway between the load application and support boundaries, and propagates by tensile fracturing through most of the specimen. The authors assessed the effects of mesh density in choosing domain grids for analysis (Fig. 2) and state that a mesh density two times smaller than shown in Fig. 2 for each model will provide similar results as described in this paper without significant loss of accuracy. In contrast, a mesh density four times smaller will provide reasonably similar results but noticeable loss of accuracy.

[5] N.G.W. Cook, The failure of rock, Int. J. Rock Mech. Min. Sci. 2 (1965) 389–403. [6] E. Hoek, Z.T. Bieniawski, Brittle fracture propagation in rock compression, Int. J. Fract. Mech. 1 (1965) 137–155. [7] Z.T. Bieniawski, Mechanism of brittle fracture of rock: part I – theory of the fracture process, Int. J. Rock Mech. Min. Sci. 4 (1967a) 395–406. [8] Z.T. Bieniawski, Mechanism of brittle fracture of rock: part II – experimental studies, Int. J. Rock Mech. Min. Sci. 4 (1967b) 407–442. [9] Z.T. Bieniawski, Mechanism of brittle fracture of rock: part III – fracture in tension and under long-term loading, Int. J. Rock Mech. Min. Sci. 4 (1967c) 425–430. [10] J.C. Jaeger, N.G.W. Cook, Fundamentals of Rock Mechanics, third ed., Chapman and Hall, London, 1979, p. 593. [11] G.I. Ofoegbu, J.H. Curran, Deformability of intact rock, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 29 (1) (1992) 35–48. [12] J.M. Ramsey, F.M. Chester, Hybrid fracture and the transition from extension fracture to shear fracture, Nature 428 (2004) 63–66. [13] J.K. Bobich, Experimental Analysis of the Extension to Shear Fracture Transition in Berea Sandstone, M.S. thesis, Texas A&M University, College Station, 2005, p. 52. [14] D.A. Ferrill, R.N. McGinnis, A.P. Morris, K.J. Smart, Hybrid failure: field evidence and influence on fault refraction, J. Struct. Geol. 42 (2012) 140–150. [15] D.A. Ferrill, R.N. McGinnis, A.P. Morris, K.J. Smart, Z.T. Sickmann, M. Bentz, D. Lehrmann, M.A. Evans, Control of mechanical stratigraphy on bed-restricted jointing and normal faulting: Eagle Ford Formation, south-central Texas, U.S.A, AAPG (Am. Assoc. Pet. Geol.) Bull. 98 (2014) 2477–2506. [16] L. Jing, A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering, Int. J. Rock Mech. Min. Sci. 40 (2003) 283–353, https://doi.org/10.1016/S1365-1609(03)00013-3. [17] S.C. Yuan, J.P. Harrison, A review of the state of the art in modelling progressive mechanical breakdown and associated fluid flow in intact heterogeneous rocks, Int. J. Rock Mech. Min. Sci. 43 (2006) 1001–1022, https://doi.org/10.1016/ j.ijrmms.2006.03.004. [18] H.Y. Liu, J.C. Small, J.P. Carter, Full 3D modelling for effects of tunnelling on existing support systems in the Sydney region, Tunn. Undergr. Space Technol. 23 (2008) 399–420, https://doi.org/10.1016/j.tust.2007.06.009. [19] D. Elmo, D. Stead, E. Eberhardt, A. Vyazmensky, Applications of finite/discrete element modeling to rock engineering problems, Int. J. Geomech. 13 (5) (2013) 565–580, https://doi.org/10.1061/(ASCE)GM.1943-5622.0000238. [20] J. Rutqvist, A.P. Rinaldi, F. Cappa, G.J. Moridis, Modeling of fault activation and seismicity by injection directly into a fault zone associated with hydraulic fracturing of shale-gas reservoirs, J. Pet. Sci. Eng. 127 (2015) 377–386, https://doi.org/ 10.1016/j.petrol.2015.01.019. [21] K.J. Smart, G.I. Ofoegbu, A.P. Morris, R.N. McGinnis, D.A. Ferrill, Geomechanical modeling of hydraulic fracturing: Why mechanical stratigraphy, stress state, and pre-existing structure matter, AAPG (Am. Assoc. Pet. Geol.) Bull. 98 (2014) 2237–2261. [22] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65. [23] P.A. Cundall, R.D. Hart, Numerical modelling of discontinua, in: J.A. Hudson (Ed.), Comprehensive Rock Engineering—Analysis and Design Methods. Volume 2: Rock Testing and Site Characterization, Pergamon Press, Oxford, England, 1993, pp. 231–243. [24] P.A. Cundall, A discontinuous future for numerical modeling in geomechanics? Proc. Inst. Civil Eng. – Geotech. Eng. 149 (1) (2001) 41–47. [25] P.A. Cundall, An Approach to Rock Mass Modelling: SHIRMS Workshop, from Rock Mass to Rock Model Workshop, 09.15.2008, Itasca Consulting Group, 2008. https://www.itascainternational.com/sites/itascacg.com/files/Education/Fairhur stFiles/Cundall%202008%20RockMass%20Workshop.pdf. (Accessed 28 November 2018). [26] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock Mech. Min. Sci. 41 (2004) 1329–1364. [27] K.J. Smart, D.Y. Wyrick, D.A. Ferrill, Discrete element modeling of Martian pit crater formation in response to extensional fracturing and dilational normal faulting, J. Geophys. Res. 116 (2011) E04005, https://doi.org/10.1029/ 2010JE003742. [28] GNO Modeling Research, CanFrac - Continuum Analysis of Discrete Fracture Initiation and Propagation in Rock, 2017. https://gno-canfrac.com/2017/05/21/ welcome-to-canfrac/. (Accessed 15 November 2019). [29] Itasca Consulting Group, Fast Lagrangian Analysis of Continua, Version 8.0, User's Guide, Itasca Consulting Group Inc., Minneapolis, Minnesota, 2016. [30] X. Ma, B.C. Haimson, Failure characteristics of two porous sandstones subjected to true triaxial stresses, J. Geophys. Res.: Solid Earth 121 (2016) 6477–6498, https:// doi.org/10.1002/2016JB012979. [31] B. Wilhelm, W.H. Somerton, Simultaneous measurement of pore and elastic properties of rocks under triaxial conditions, J. Soc. Petrol. Eng. (1967) 283–294. Paper no. 1705. [32] M.J. Aldrich, Pore pressure effects on Berea sandstone subjected to experimental deformation, Geol. Soc. Am. Bull. 80 (8) (1969) 1577–1586. [33] R.E. Goodman, Introduction to Rock Mechanics, John Wiley and Sons, Inc., New York, 1980, p. 478. [34] D.J. Hart, H.F. Wang, Laboratory measurements of a complete set of poroelastic moduli for Berea sandstone and Indiana limestone, J. Geophys. Res. 100 (B9) (1995) 17,741–17,751. [35] W. Zhu, T.F. Wong, The transition from brittle faulting to cataclastic flow: permeability evolution, J. Geophys. Res. 101 (B2) (1997) 3027–3041.

The evaluation of CanFrac® simulations in this paper indicate the procedure can be relied upon for predicting and understanding rock mechanical behavior. The procedure offers promise in expanding understanding of the mechanical behavior of rocks through the capability to obtain information about fracturing: viz., occurrence and location; type (tensile or shear); and geometrical attributes such as dip, dip direction, and aperture of fractures. Furthermore, the procedure shows promise in predicting and understanding large-scale behavior such as localization of rock deformation, stress control of fracture mechanisms in failure zones, and conditions for occurrence of multiple parallel fractures. Authors individual contributions Goodluck I Ofoegbu: Conceptualization; Data curation; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Resources; Software; Supervision; Validation; Visualization; Roles/Writing – original draft; Writing – review & editing. Kevin J. Smart: Conceptualization; Data curation; Funding acquisition; Investigation; Methodology; Resources; Validation; Visualization; Writing – review & editing. Declaration of competing interest None. Acknowledgements This work was supported in part by GNO Modeling Research LLC and in part by Southwest Research Institute internal research and development project R8760. References [1] J. Handin, R.V. Hager, Experimental deformation of sedimentary rocks under confining pressure: tests at room temperature on dry samples, AAPG (Am. Assoc. Pet. Geol.) Bull. 41 (1) (1957) 1–50. [2] W.F. Brace, E.G. Bombolakis, A note on brittle crack growth in compression, J. Geophys. Res. 68 (12) (1963) 3709–3713. [3] J. Handin, R.V. Hager, M. Friedman, J.N. Feather, Experimental deformation of sedimentary rocks under confining pressure: pore pressure tests, AAPG (Am. Assoc. Pet. Geol.) Bull. 47 (5) (1963) 717–755. [4] W.F. Brace, Brittle fracture of rocks, in: W.R. Judd (Ed.), State of Stress in the Earth's Crust, American Elsevier Publishing Company, Inc., New York, 1964, pp. 111–178. 17

G.I. Ofoegbu, K.J. Smart

Results in Engineering 4 (2019) 100070 [38] A. Vervoort, K.-B. Min, H. Konietzky, J.-W. Cho, B. Debecker, Q.-D. Dinh, T. Fruhwirt, A. Tavallali, Failure of transversely isotropic rock under Brazilian test conditions, Int. J. Rock Mech. Min. Sci. 70 (2014) 343–352.

[36] T.F. Wong, P. Baud, Mechanical compaction of porous sandstone, Oil Gas Sci. Technol. Rev. IFP 54 (6) (1999) 715–727. [37] D. Li, L.N.Y. Wong, The Brazilian disc test for rock mechanics applications: review and new insights, Rock Mech. Rock Eng. 46 (2013) 269–287, https://doi.org/ 10.1007/s00603-012-0257-7.

18