Elastoplastic constitutive behavior and strain localization of a low-porosity sandstone during 6 brittle fracturing 7

We investigated the elastoplastic behavior and strain localization of the Zigong sandstone (porosity: 6.5%) during brittle fracturing based on two series of axisymmetric compression experiments. The experiments were conducted under various confining pressures (σ3 = 0 ~ 80 MPa). For each confining pressure, the sandstone specimens were deformed under constant axial and circumferential strain rates, respectively. When σ3 < 60 MPa, the sandstone first undergoes stable deformation in the post-peak stage and then loses its stability. Before the emergence of instability, the mechanical behavior is hardly affected by the controlling method. When the confining is larger, the sandstone manifests a stable failure process during the whole loading stage. The observed elastoplastic behavior was described by a two-yield surface constitutive model established in the framework of generalized plastic mechanics. The proposed constitutive model incorporates two quadratic yield functions, as well as two linearly independent plastic potential functions, to honor the shear yield and volumetric dilatancy, respectively. Via the return mapping algorithm, the proposed constitutive model was verified by comparing the numerical results with experimental results. In addition, the two-yield surface constitutive model, which is equivalent to the model proposed by Rudnicki and Rice,1 was applied to localization analysis. Assuming that the onset of localization occurs at peak stress, frictional coefficient μ and dilatancy factor β were determined from experimental data. The variations of both plastic parameters predict the transition of localization mode from pure dilation bands under uniaxial compression to pure shear bands at high confining pressures, which is consistent with the experimental observations.


Introduction 38
The mechanical behavior of sandstone is intimately related to a variety of geologic processes and 39 engineering applications. As with most of rocks, sandstone exhibits highly-nonlinear pressure 40 dependence of its deformation, such as volumetric dilatancy and strain hardening/softening. 41 Phenomenologically, the deformation process of sandstone subject to compression is typically 42 accompanied by the development of planar deformation bands, also referred to as strain 43 localization. 1 With the increase of pressure, strain localization varies the manifestations from 44 dilation or shear bands in the brittle regime 2 to shear-enhanced compaction or homogenous 45 cataclastic flow when entering the ductile regime. 3

,4 46
In order to understand these abundant mechanical properties, micromechanics and continuum 47 Class II behavior under a low confining pressure and purely Class I behavior under a high confining 115

pressure. 116
Based on these stress-strain curves, we can determine several characteristic stress thresholds, 117 including crack initiation stress σci, crack damage stress σcd, and peak stress σp, to characterize the 118 deformation process. [16][17][18] Among them, the determination of σci is essential to the construction of 119 the constitutive model since it suggests the onset of plastic deformation. However, its determination 120 has usually been subjective and of high uncertainty. 16,17,[19][20][21][22][23] To circumvent this, we utilize an 121 objective method called Lateral Strain Response (LSR) method 24 in this study. In addition, σcd and 122 σp correspond to the point where volumetric strain curve reverses and the highest point of the stress-123 strain curve, respectively, which can be readily identified. Figure 2 shows the determined values of 124 σci, σcd and σp and their positive dependence on σ3. As alluded to above, the consistency between 125 both series for each stress threshold indicates that the two control methods hardly affect the pre-126 peak behavior of Zigong sandstone. In addition, the peak strength of Zigong sandstone can be well 127 captured by Hoek-Brown failure criterion and its extended 3D version. 15,25 128

Elastoplastic Constitutive Properties of Zigong Sandstone 129
Instead of using a single-yield surface model, we adopt the generalized plastic theory 13 to establish 130 a two-yield surface constitutive model for Zigong sandstone. In what follows, a brief introduction 131 of the generalized plastic theory is first provided in Section 3.1. Based on the theoretical framework, 132 we then explore the effect of confining pressure on the elastic shear and bulk moduli in Section 3.2. 133 We propose a shear yield function and a volumetric yield function in Section 3.3 and Section 3.4, 134 respectively, to include shear-and dilatancy-related plastic flow mechanisms. The two-yield 135 constitutive model adopts non-associated flow rules and isotropic hardening rules. 136

Fundamentals of the Generalized Plastic Mechanics 137
Classic elastoplastic theory assumes that the total strain increment consists of elastic and plastic 138 strain increments: 139 in which the elastic strain increment e ij d can be calculated by Hooke's law while the plastic strain increment  p ij d is closely related to the plastic potential function Q by: 142 in which dλ is a non-negative scalar and σij is the stress component. As pointed out by Zheng et al. 13 , 144 Eqn. 2 implies that the direction of the plastic strain increment depends only on stress states rather 145 than stress increment, the latter of which has been experimentally observed in numerous 146 geomaterials. In addition, the plastic potential function Q for geomaterials is usually difficult to 147 determine according to the yield surface, i.e., following non-associated flow rule. 148 In the generalized plastic mechanics, three linearly-independent yield surfaces and corresponding 149 plastic potential functions are assumed to be coexisting at any point in the stress space. In other 150 words, the total plastic strain increment is the sum of the plastic strain increments of all of the yield 151 surfaces. Consequently, the flow rule in the single-yield surface model can be generalized to: 152 in which the plastic potential functions Qk can be any three linearly-independent functions, which 154 jointly contribute to the total plastic strain increment. This advantage avails simple selection of the 155 plastic potential functions and parameter fitting. 156 The magnitude of the plastic strain increment can be determined by the yield functions along with 157 corresponding plastic strains. As alluded to in Section 2.2, the values of σci of all tests depict the 158 initial yield function F(σij) = 0. Following this, the size, center, and shape of the yield surface will 159 change with the ongoing loading, known as subsequent yield surface. To describe how an initial 160 yield surface evolves to the subsequent yield surface, the hardening parameter Hα is imperative, 161 which can be a function of plastic strains changing with deformation. Accordingly  Considering the consistency condition, the stress state should be always located on the yield surfaces 177 to ensure continuous plastic deformation. By differentiating Eqn. 6, we have:

Effect of Confining Pressure on the Elastic Properties
According to the notations in Eqn. 5 to Eqn. 8, we specially transform the principal stresses to von 184 Mises equivalent stress q and mean stress p. For strains, we use the generalized shear strain γ and 185 volumetric strain v  , which are defined as: 186 where ε1, ε2, ε3 are the three principal strains. 189 Consequently, the stress-strain curves in Figure 1 can be plotted as q-γ and p-εv curves, which enable 190 the determination of elastic shear (G) and bulk (K) moduli, respectively, based on the determined 191 values of σci. In Figure 3, the effect of confining pressure on G and K is shown. It is found that, for 192 both elastic moduli, the results of the CSC tests are consistently larger than that obtained from the 193 ASC tests. The mean values of both G and K show an increase-then-decrease trend, which can be 194 fit by the following parabolic functions (black curves in Figure 3 in which the material parameters are given as: aG = -9.568×10 -4 , bG = 0.1165, cG = 8.956; aK = -198 2.135×10 -3 , bK = 0.2178, cK = 12.02. 199 With the determined elastic moduli, we are further able to calculate the plastic shear ( p q  ) and 200 volumetric ( V p  ) strains, which are plotted in Figure 4 in relation to q. In the following two sections, 201 a shear yield function and volumetric yield function will be proposed based on the curves in Figure  202 4. 203

Shear yield surface and hardening function 204
The plastic shear strain p q  is selected as an internal variable for the shear yield mechanism. 205 Specifically, a series of plastic shear strain values are selected, which give the contours of the plastic 206 shear strain in the p-q plane. Note that, the data corresponding to the unstable failure process is not 207 included since the testing system failed to record reliably. As shown in Figure 5a, these contours 208 are essentially the shear yield surfaces. As the plastic shear strain increases, the shear yield surface 209 first increases then decreases (i.e., strain softening, which is not explicitly exemplified in Figure  210 5a). To quantitatively describe these shear yield surfaces, an expression similar to the form of Hoek-211 Brown failure criterion is employed: 212 where ms is the material parameter. By fitting to the test data, it is found that ms can be approximately 214 selected as 10 for Zigong sandstone. 215 For each confining pressure, the average values of the ASC and CSC tests are utilized to obtain the 216 relationship between Hq and p q  . As shown in Figure 5b, such relationship features a transition from 217 strain hardening to softening as p q  increases for all confining pressures. Accordingly, the following 218 function is proposed for hardening rule: 219 where aq, bq, cq, and dq are constants. Based on the curves in Figure 5b, each of these constants is 221 given as a function of σ3: 222

Dilatant Volumetric yield surface and hardening function 227
In the previous section, the hardening effect is ascribed to the plastic shear strain. Since deformation 228 in low-porosity rocks during brittle fracturing is also characterized by volumetric dilatancy, 19 a 229 series of volumetric yield surfaces incorporating dilatant plastic volumetric change is proposed. As 230 implied by Figure 4, Zigong sandstone does not show any volumetric compaction when confining 231 pressure increases up to 80 MPa. Therefore, the volumetric compaction mechanism, specific to the 232 'cap' model for porous rocks 3,26 is not considered in this study. 233 To describe the volumetric yield surface, the plastic volumetric strain V p  is utilized as the internal 234 variable. In Figure 6a, the obtained volumetric yield surfaces are similar to the shear yield surfaces 235 in Figure 5a. Hence, Eqn. 13 can also be used but with different hardening parameter Hv: 236 in which the approximation of ms = 10 is also found to be appropriate. The consistency condition (or Eqn. 8) implies the following loading/unloading conditions: which are often called Kuhn-Tucker complementary conditions. For each yield surface, the first 261 term indicates that the stress state must lie on or within the yield surface while the second term 262 suggests that the plastic strain increment is always non-negative. The last term assures that the stress 263 state remains on the current yield surface during plastic loading. Then the consistency condition can 264 be stated as: When the model is in the elastic regime, D ep is simply equal to D e . 281 The above expressions complete the general mathematical formulation of the proposed two-yield 282 surface constitutive model for Zigong sandstone. Note that, Eqn. 30 can be flexibly reduced to the 283 single surface formulation via leaving out the terms associated with the inactivated yield surface. 284

285
In this section, the proposed constitutive relations are validated by comparing with the experimental 286 data. Unlike the linear elasticity, the elastoplastic constitutive relations are nonlinear and not 287 analytically integrable. Therefore, numerical integration is required to define the resulting 288 mechanical response. Among many techniques for numerical integration, 27 the return mapping 289 algorithm is particularly adopted in this study. 290 Implementation of the return mapping algorithm is essential to find the numerical solution of a 291 nonlinear problem about stress/strain update. Specifically, for a prescribed loading path, we need to 292 find the respective increments of stress σ , strain ε , and hardening parameter Hα (the subscript α = v 293 or q) at each loading step according to the elastoplastic constitutive relations. The problem can be 294 further formulated in discrete forms by using the backward Euler method as: 295     in which the three variables { σ , ε p , Hα}n are known at the current step n and +1 ε p n is the unknown 297 to be determined. In addition, the Kuhn-Tucker conditions are given by: 298 299 To address this nonlinear optimization problem, the return mapping algorithm generally involves 300 two steps: (1) trial elastic prediction; and (2) plastic correction returning the trial state to the current 301 yield surface. During the integration process, the strain increment ε is fixed for each loading step. 302 The trial elastic prediction can be obtained by taking:  represented by Path III. Furthermore, the application of the return mapping algorithm to the shear 317 yield function is shown in Figure 7b as an example. 318 Applying the return mapping algorithm to the proposed elastoplastic constitutive relationship, we 319 are able to perform the numerical simulations. As shown in Figure 8, the numerical results show a 320 good agreement with the experimental results under different confining pressures. The proposed 321 elastoplastic constitutive model is able to capture the mechanical properties (i.e., strain hardening, 322 strain softening, dilatancy, and confining pressure effect) of Zigong sandstone. However, since the 323 unstable post-peak stage and frictional sliding are not included in the constitutive model, the 324 complete stress-strain curves are not shown. 325

Analysis of Strain Localization in Zigong Sandstone 326
In this section, strain localization and its pressure dependence are investigated in the framework of 327 the proposed constitutive model for Zigong sandstone. The proposed two-yield surface model 328 explicitly considers shear yield and volumetric dilatancy during brittle fracturing, which facilitates 329 the description of stress-strain relations. By contrast, Rudnicki and Rice 1 (hereafter referred to as 330

RR model) considers only shear yield in the bifurcation theory while the volumetric dilatancy is 331
included by introducing the dilatancy factor (β). In the following, we first demonstrate the proposed 332 two-yield surface model is equivalent to the RR model in the brittle regime. Then we examine the 333 observed strain localization in Zigong sandstone under different confining pressures. 334

Equivalence of the Proposed Constitutive Model to RR Model 335
In order to compare with the classical RR model, we rewrite Eqn. 20 as: also noteworthy that volumetric compressive mechanism can also be incorporated to model the yield 347 'cap', which, however, needs more complicated algebraic manipulation to apply the bifurcation 348 theory 29 . 349

Effect of Confining Pressure on Strain Localization Mode 350
As indicated by the RR model, the three plastic parameters (μ, β, H) are relevant when applying the 351 bifurcation theory to strain localization. Since the plastic hardening modulus H is a function of the 352 orientation of the potential band and decreases monotonically with the ongoing plastic deformation, 353 When the left inequality is violated, compaction band would emerge; while dilation band would 361 occur when the right inequality is not satisfied.
With these conditions, we are able to investigate the effect of confining pressure on the strain 363 localization mode of Zigong sandstone by evaluating parameters μ and β. In Figure 9 is predicted that the sandstone still has the potential to develop dilation bands when σ3 = 10 and 20 387 MPa. When confining pressure further increases, pure shear bands are predicted, in good accordance 388 with experimental observations. 15 Finally, the sandstone is far from ductile regime under the experimental conditions. Compared with other porous sandstones, Zigong sandstone features a wide 390 brittle range due to its low porosity. 15 391

Conclusions 392
We conducted two series of axisymmetric compression experiments in the low-porosity Zigong 393 sandstone under various confining pressures (0 ~ 80 MPa with a step of 10 MPa). For each confining 394 pressure, sandstone specimens were deformed under constant axial and circumferential strain rates, 395 respectively. It is found that the Zigong sandstone features a combination of Class I (stable) and 396 Class II (unstable) behavior in the post-peak stage under a low confining pressure (< 60 MPa) and 397 purely Class I behavior for a high confining pressure. 398 Based on the theory of generalized plastic mechanics, a two-yield surface elastoplastic constitutive 399 model was proposed to describe the deformation characteristics of the sandstone. The proposed 400 elastoplastic constitutive model employs two quadratic yield functions, along with two linearly 401 independent plastic potential functions, to honor shear yield and volumetric dilatancy, respectively. 402 Numerical integration of the constitutive relations was carried out using the return mapping 403 algorithm. It is found that the resulting stress-strain relations are in good agreement with the 404 experimental results. It can be concluded that the proposed model adequately captures the 405 elastoplastic behavior (i.e., strain hardening, strain softening, dilatancy, and confining pressure 406 effect) of Zigong sandstone in the brittle regime. 407 In the context of brittle fracturing, it is demonstrated that the proposed two-yield surface model is 408 essentially equivalent to the single-yield surface model proposed by Rudnicki and Rice 1 for strain 409 localization analysis. In addition, formulations of three relevant plastic parameters (μ, β, H) were 410 derived according to the proposed constitutive equations. To analyze strain localization, the effects 411 of plastic deformation and confining pressure on parameters μ and β were investigated. As plastic 412 deformation accumulates, β is relatively constant, while μ increases rapidly and reaches a plateau 413 subsequently. With increasing confining pressure, both μ and β decrease, and μ is always larger than 414 β. The theoretical predictions indicate that the localization mode in Zigong sandstone undergoes a 415 transition from pure dilation bands under uniaxial compression to pure shear bands at high confining 416 pressures.