Numerical analysis of anisotropic stiffness and strength for geomaterials

Full Length Article Numerical analysis of anisotropic stiffness and strength for geomaterials Fei Song a , Manuel A. González-Fernández b , Alfonso Rodriguez-Dono a , c , * , Leandro R. Alejano b a Department of Civil and Environmental Engineering, Universitat Politècnica de Catalunya (UPC), Barcelona, 08034, Spain b CINTECX, GESSMin Group, Department of Natural Resources and Environmental Engineering, University of Vigo, Vigo, 36310, Spain c International Centre for Numerical Methods in Engineering (CIMNE), Barcelona, 08034, Spain article info Article history: Received 19 October 2021 Received in revised form 31 January 2022 Accepted 14 April 2022 Available online 27 June 2022 Keywords: Rock mechanics Anisotropy Numerical modelling CODE_BRIGHT Shale Slate abstract In numerical modelling, selection of the constitutive model is a critical factor in predicting the actual response of a geomaterial. The use of oversimplified or inadequate models may not be sufficient to reproduce the actual geomaterial behaviour. That selection is especially relevant in the case of anisotropic rocks, and particularly for shales and slates, whose behaviour may be affected, e.g. well stability in geothermal or oil and gas production operations. In this paper, an alternative anisotropic constitutive model has been implemented in the finite element method software CODE_BRIGHT, which is able to account for the anisotropy of shales and slates in terms of both deformability and strength. For this purpose, a transversely isotropic version of the generalised Hooke’s law is adopted to represent the stiffness anisotropy, while a nonuniform scaling of the stress tensor is introduced in the plastic model to represent the strength anisotropy. Furthermore, a detailed approach has been proposed to determine the model parameters based on the stressestrain results of laboratory tests. Moreover, numerical analyses are performed to model uniaxial and triaxial tests on Vaca Muerta shale, Bossier shale and slate from the northwest of Spain (NW Spain slate). The experimental data have been recovered from the literature in the case of the shale and, in the case of the slate, performed by the authors in terms of stress-strain curves and strengths. A good agreement can be generally observed between numerical and experimental results, hence showing the potential applicability of the approach to actual case studies. Therefore, the presented constitutive model may be a promising approach for analysing the anisotropic behaviour of rocks and its impact on well stability or other relevant geomechanical problems in anisotropic rocks. Ó2023 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/). 1. Introduction Within the framework of geomechanics, excavation or well stability problems can be associated with poor design approaches, often related to improper understanding, perception or use of rock constitutive models (Rodriguez-Dono, 2011;Mánica et al., 2016, 2021;Mánica, 2018;Alonso et al., 2021a, b;Song, 2021). Isotropic mechanical behaviour approaches have been widely studied and analysed, including elastic, viscoelastic, perfectly-plastic, purely brittle or even strain-softening constitutive models (Alejano and Alonso, 2005;Arzúa and Alejano, 2013;Wang et al., 2014;Chen et al., 2019;Song et al., 2020,2021a, b;Song, 2021;Song and Rodriguez-Dono, 2021). In these cases, there are generally valid and widely accepted approaches to estimating in a reasonably accurate manner the material parameters representative of the behaviour of the rock mass at stake. However, many sedimentary and metamorphic rocks tend to exhibit anisotropy due to bedding or foliation, i.e. planes of weakness that control the mechanical response of these materials, thus different responses were observed for different directions of stress application (Alejano et al., 2021). Furthermore, interpreting and modelling the anisotropy of rocks are still insufficiently understood (Hudson, 2008;Alejano et al., 2021;Mánica et al., 2021). If the anisotropic properties are not accounted for, significant errors may be introduced in stress and displacement analyses (Barla, 1972;Kirkgard and Lade, 1993; Abelev and Lade, 2004;Kanfar et al., 2015). Hence, having reliable studies and proper simulations of the mechanical behaviour of this type of rock, especially accounting for anisotropy, remains a critical issue in the field of rock mechanics. *Corresponding author. Department of Civil and Environmental Engineering, Universitat Politècnica de Catalunya (UPC), Barcelona, 08034, Spain. E-mail address: [email protected] (A. Rodriguez-Dono). Peer review under responsibility of Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Contents lists available at ScienceDirect Journal of Rock Mechanics and Geotechnical Engineering journal homepage: www.jrmge.cn Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338 https://doi.org/10.1016/j.jrmge.2022.04.016 1674-7755 Ó2023 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Production and hosting by Elsevier B.V. This is an open access article under theCCBYNC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Laboratory tests show that many types of rocks, such as shale, slate or schist, but also in a less relevant manner, gneiss and sandstone, exhibit a marked anisotropic behaviour (Ambrose, 2014; Ding et al., 2020;Alejano et al., 2021). The mechanical properties, such as stiffness and strength, are significantly influenced by planes of weakness (Mánica, 2018;Ismael and Konietzky, 2019;Ismael et al., 2019;Ding et al., 2020;Alejano et al., 2021;Mánica et al., 2021). The anisotropy of stiffness has been rigorously addressed in the past, thus we are today capable of interpreting and modelling this elastisc part of the behaviour accurately (Barla, 1972; Worotnicki, 1993;Amadei, 1996;Chen et al., 1998;Cho et al., 2012; Nejati et al., 2019). Moreover, many standard numerical codes used in rock mechanics allow simulating anisotropic deformability (Gonzaga et al., 2008). In the most general case of anisotropic stiffness, 21 independent constants are needed to completely define the anisotropy of elastic stiffness through the generalised Hooke’slaw(Barla, 1972;Wittke, 1990). When considering elastic symmetry, the number of constants can be reduced to nine for orthotropic materials or five for the transversely isotropic materials (Barla,1972;Cho et al., 2012). In conclusion, anisotropic deformability, especially transverse or cross-anisotropy, has been well understood and analysed. However, there is still room for further understanding the strength anisotropy. In the last century, anisotropic failure criteria were developed. Jaeger (1960) proposed the Jaeger’s plane of weakness (JPW) strength criterion, which seems to be the most widely used strength criterion for cross-anisotropic rocks (Cho et al., 2012; Ambrose, 2014). The concept behind the JPW strength criterion is associating two potential failure mechanisms to rock strength: one associated with intact rock and the other associated with shear failure along pre-existing planes of weakness having particular orientation. The parameters needed to represent this criterion are the cohesion and friction of the intact rock and those of the planes of weakness. Thus, they have a clear physical meaning and can be obtained based on a limited number of tests. Moreover, Ambrose (2014) proposed a robust approach for estimating the strength constants of the JPW model based on statistical approaches, and thus, facilitating the empirical estimate of the parameters needed for the application of this model. It is easy to find in literature practical applications of this JPW failure criterion approach implemented in numerical simulations to model laboratory triaxial tests and underground stability problems (Chang and Konietzky, 2018;Zimmerman et al., 2018;Xu et al., 2021). Nevertheless, observed experimental results show that for some rocks, and notably for slate, the strength varies continuously through all anisotropy angles ( b ), something well reflected by the JPW criterion. However, two significantly different strength levels were observed for samples cut perpendicular and parallel to the foliation, respectively, which the JPW cannot account for (Alejano et al., 2021). Thus, at least in this point, the JPW model meets one of its limitations for representing the actual anisotropic strengths of rocks. As opposed to the ‘discontinuous’JPW model, Pariseau (1968) proposed a criterion that can predict a smooth, continuous variation of strength with anisotropy angles. However, this approach is not commonly used in practice since the needed parameters might have not a clear physical meaning, and its estimate should be performed in a rigorous statistical approach and usually based on many (often unavailable) data (Ambrose, 2014). Alternativelly, Mánica et al. (2016) proposed a cross-anisotropic strength model (named Mánica’s model or Mánica’s approach in the following), by introducing a nonuniform scaling of the stress tensor. Mánica’s model can simulate a continuous variation of the strength at different anisotropy angles ( b ). Nevertheless, unlike the JPW model, it allows different strength levels between samples with b ¼0  and b ¼90  . Additionally, Mánica’s approach can be easily incorporated into the implemented isotropic constitutive models with minor revisions (Mánica et al., 2016). Based on the above advantages, Mánica’s approach is adopted in this research to represent the anisotropy of strength characteristic of shale and slate samples. Moreover, a relatively simple approach based on laboratory tests is proposed in this study, to calibrate the parameters in the anisotropic constitutive model in a rigorous manner. It is worth noticing that the presented approach can be used for modelling actual cases of well stability or small-scale (laboratory) rock response. Nevertheless, for modelling large-scale excavations, such as tunnels or slopes, the parameters should be recalibrated based on the geotechnical quality of the rock masses formed by these rocks. In summary, the use of inappropriate anisotropic constitutive models is one relevant limitation in numerical modelling to predict the actual response of rocks. In this paper, an alternative numerical approach is presented and programmed in the finite element method software CODE_BRIGHT (Olivella et al., 2021), which can reproduce the anisotropy of stiffness and strength for different types of rocks. Moreover, a simple approach is proposed to calibrate the model parameters. Furthermore, numerical simulations of laboratory uniaxial and triaxial tests are performed, considering a wide range of confining stress levels, load orientations and three types of rocks, i.e. two types of shales with data recovered from the literature (Ambrose, 2014) and a slate tested by the authors and results in terms of stiffness and strength properties reported in the previous study (Alejano et al., 2021). Note that the stress-strain curves of slate samples are presented in this article for the first time. Comparisons are carried out between numerical predictions and experimental data in terms of the stressestrain curves, stiffnesses and strengths, which result in a level of accuracy well over the average in rock mechanics applications. Accordingly, the proposed parameter estimation and numerical approaches can be considered a useful tool for predicting the anisotropic behaviour of the rocks at stake. 2. Theoretical background It has long been recognised that geomaterials may exhibit anisotropic behaviour, i.e. geomaterials may have different properties at different directions (Donath, 1964;Cho et al., 2012; Alsuwaidi et al., 2021). As an example, in slate rocks (see Fig. 1), the mechanical properties heavily depend on the anisotropy angle ( b ), including the axial stressestrain curves (Fig. 1a), the apparent elastic moduli (Fig.1b) and the peak compressive strengths (Fig.1c). Additionally, previous studies found that the compressive strength of anisotropic rocks could differ up to an order of magnitude depending on the directions of application of stress with respect to the orientation of the planes of weakness (Ambrose, 2014). Therefore, proper simulation of anisotropic behaviour, and notably of strength, is crucial to reliably represent the actual behaviour of some types of rocks. In general, the stressestrain behaviour of rocks can be referred to as cross-anisotropic or transversely isotropic. As shown in Fig. 1b, the theoretical solutions of the cross-anisotropic elastic model can reasonably represent apparent elastic moduli from experimental results. Therefore, in this study, the cross-anisotropic constitutive model is adopted to reproduce the anisotropic stiffness of geomaterials. Five independent elastic constants are needed (Barla, 1972;Amadei, 1996), i.e. two Young’s moduli for the directions parallel (E) and perpendicular (E 0 ) to the isotropic plane, two corresponding Poisson’s ratios ( n and n 0 ), and the shear modulus G 0 for F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338324 shear loading in the isotropic plane. The adopted cross-anisotropic elastic model is described in detail in Section 3. Concerning peak compressive strengths, various failure models have been developed to reproduce the anisotropy of strength. Among all these models, JPW model may be the most commonly used one in interpreting experimental results and numerical modelling. Fig. 2a shows the conceptual model of rock samples with planes of weakness. As shown in Fig. 2b, two groups of failure criteria are adopted in the JPW model: one associated with the failure through intact rock and the other associated with the failure along the plane of weakness. The minimum strength value obtained from both criteria plays a critical role in the failure mode. The so-called JPW model was initially considered for jointed rocks, but it has been extended for foliated rocks, such as shales and slates (Cho et al., 2012). However, concerning intact rock strength, the predicted strength from the JPW model is forced to be the same when the major principal stress is parallel ð b ¼90  ) or perpendicular ( b ¼ 0  ) to the foliation, respectively, which is different from the observed results. Different failure mechanisms result in different strength levels for samples with varied anisotropy angles ( b ). As observed in laboratory tests of slate from the northwest of Spain (NW Spain slate) and other rocks (see Fig. 3)(Donath, 1964;Alejano et al., 2021), for the samples cut perpendicular to foliation ð b ¼0  ), a double-cone failure can be observed, as in typical isotropic rocks. However, for samples cut parallel to foliation ( b ¼90  ), failure through vertical planes combined with local shear bands can be Fig. 2. (a) Conceptual representation of a cross-anisotropic geomaterial, (b) Peak compressive strength versus anisotropy angle ( b ) based on the JPW theory, and (c) Fictitious scaled yield surfaces in the principal stress space in the Mánica’s model, based on the works by Jaeger (1960),Mánica (2018) and Alejano et al. (2021). s 1 is the major principal stress, s 2 is the intermediate principal stress, and C N and C S are the nonuniform scaling factors in the Mánica’s model. Fig. 1. Graphs representing (a) stressestrain curves for slate samples with different anisotropy angles ( b ¼30,45 and 90), (b) apparent elastic moduli versus b , and (c) peak compressive strengths ( s 1,max ) versus b (Modified from Alejano et al., 2021). s 3 is the minor principal stress. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338 325 observed. In general, for the rock samples with b ¼0  and 90  ,the strengths are normally different although the failure mechanism of both cases could be referred to as failure through intact rock (Jaeger, 1960). Additionally, for the rock sample with b ¼45  ,thefailure mode was typically clean sliding through a foliation plane (Cho et al., 2012;Alejano et al., 2021). To overcome this limitation of the JPW model, Alejano et al. (2021) proposed the so-called 2 MC-JPW and 2HB-JPW strength models, where two different strength levels can be adapted for samples with b ¼0  and b ¼90  . However, in these improved models, the strength variation along different anisotropy angles is still discontinuous (as in the original JPW model). Moreover, 2 MCJPW and 2HB-JPW models have not been implemented yet in any code, and hence, it still needs further research. As a result of all stated above, strength anisotropy is still particularly challenging from the modelling point of view. On the other hand, some other models have been developed to represent anisotropy of strength. Among all these models, Mánica et al. (2016) and Mánica (2018) proposed a general crossanisotropic plastic model based on a nonuniform scaling of the stress tensor. According to this model, the anisotropic stress space can be obtained by introducing the nonuniform stress scaling factors (C N and C S ) to modify the isotropic yield and, thus, account for the cross-anisotropy of strength, as shown in Fig. 2c(Mánica, 2018). A more detailed explanation of scaling factors (C N and C S ) can be found in Section 3. The advantages of Mánica’s approach are that (1) it is able to represent different levels of strength for samples cut perpendicular and parallel to foliation, (2) it is able to represent a continuous variation of the strength along different anisotropy angles. In addition to the contribution of Mánica et al. (2016) and Mánica (2018), this article (see Section 3) proposes a simple but robust approach to calibrating model parameters and introducing the concept of the dilatancy angle in the constitutive model, which can facilitate its application in rock mechanics. Furthermore, in this study, the smoothed Mohr-Coulomb model is adopted in the numerical implementation, with the aim of improving the numerical reliability and efficiency. Finally, numerical simulations are performed for three types of rocks, i.e. Vaca Muerta shale, Bossier shale and NW Spain slate. Numerical predictions are compared with experimental results, which display a good agreement, validating the proposed constitutive model and approach for calibrating model parameters for the case of the rocks under scrutiny. In addition, the presented model is developed within the framework of the viscoelasticeviscoplastic (VEVP) series constitutive models (Song et al., 2020,2021a, b) in the coupled Thermo-HydroMechanical (THM) CODE_BRIGHT, and therefore, the presented model can be easily extended to simulate coupled THM as well as time-dependent problems. 3. Anisotropic constitutive model for rocks The viscoelastic-viscoplastic (VEVP) series constitutive models have been proposed by Song (2021) and Song et al. (2020, 2021a, b). However, this type of VEVP models can only represent isotropic properties of geomaterials. Hence, this study constitutes a development of the previous VEVP models proposed by Song (2021) and Song et al. (2020, 2021a, b), by introducing anisotropic properties in terms of stiffness and strengths. Note that further improvements may be introduced into these models to account for the anisotropy in viscous dashpots and Kelvin models. 3.1. Constitutive model In order to represent the anisotropic stiffness and strength for geomaterials, an alternative anisotropic constitutive model is presented in this study. As shown in Fig. 4, the adopted anisotropic constitutive model consists of an elastic spring and Perzyna’s viscoplastic model. The generalised Hooke’s law is used to describe the reversible elastic constitutive relationship (Lechnickij and Lekhnit  skiĭ,1963; Wittke,1990), while Perzyna’s viscoplastic model is used to represent the irreversible strain. Global ðx;y;zÞand local ð1;2;3Þcoordinate systems are needed to describe the crossanisotropic behaviour of geomaterials, as shown in Fig. 5. The global coordinate system ðx;y;zÞfirst rotates around the z-axis for the angle of a , and then around the y-axis for the angle of b ,to obtain the local coordinate system ð1;2;3Þ. Therefore, for cases of bedding plane parallel to the 1-o-2 plane as referred in this article, b should be the same as the anisotropy angle of rock samples. Experience showed that the irreversible strains and the associated stress redistribution might heavily depend on time (Wittke, 1990). Thus, the viscoplastic model may be appropriate to represent the behaviour of geomaterials. Moreover, the time-dependent viscoplastic model can be simplified to the time-independent constitutive model by adjusting the viscosity (Fig. 4), and therefore, the viscoplastic model is used. Due to the lack of experimental data representing the time-dependent behaviour of the studied specimens, no time-dependency is considered in this research. However, the proposed approach can be easily extended to Fig. 3. Different failure modes of rock samples (NW Spain slate) with different anisotropy angles: (a) b ¼0, (b) b ¼45, and (c) b ¼90. The confining pressure is 0 MPa. Fig. 4. Adopted anisotropic constitutive model for rocks. h vp is the viscosity in the viscoplastic model. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338326 simulate the time-dependent cross-anisotropic behaviour of geomaterials. The total strain rate tensor of the proposed anisotropic constitutive model, dε=dt, can be decomposed into components describing the reversible (caused by the elastic spring, ε e ) and irreversible (caused by the viscoplastic model, ε vp ) parts, as shown in Eq. (1). Note that a viscoplastic solution which is close to the ‘true’purely plastic solution can be obtained by sufficiently decreasing the viscosity h vp to be 0 (Alonso et al., 2005). dε dt¼dεe dtþdεvp dt¼dεe dtþ1 h vp h F ðFÞivG v s (1) F ðFÞ¼ F ðFÞð F ðFÞ0Þ 0ð F ðFÞ<0Þ(2) where Gis the potential, Fis the failure criterion, and F ðFÞ¼F m with mrepresenting stress power (Perzyna, 1966;Mánica, 2018; Song, 2021). The generalised Hooke’s law could describe the elastic constitutive relationship of cross-anisotropic rocks (Payne and Lekhnitskii, 1964;Wittke, 1990;Alejano et al., 2021), as shown in Eq. (3). For cases where the global coordinate system ðx;y;zÞcoincides with the local one ð1;2;3Þ, the local elastic compliance matrix (C local e ) is the same as the global elastic compliance matrix (C global e ). Therefore, Eq. (3) can be simplified to Eq. (4) for this case. εglobal e¼Cglobal e s global ¼TT DClocal eTD s global (3) εglobal e¼Clocal e s global (4) where Clocal e¼ 2 6 6 6 6 6 6 4 1=E n 0=E0 n =E00 0  n 0=E01=E0 n 0=E000 0  n =E n 0=E01=E00 0 0001 =G000 00002ð1þ n Þ=E0 00000 1 =G0 3 7 7 7 7 7 7 5 where ε global e (or s global ) is the elastic strain (or the current stresses) in the global coordinate system; C global e (or C local e ) is the compliance elastic matrix with respect to the global (or local) coordinate system; T D is the transformation matrix, which depends on the two anisotropic angles ( a and b in Fig. 5). The detailed expression of T D is shown in Appendix A. To represent the strength anisotropy, three stress spaces ( s global , s local and s ani ) need to be introduced (Mánica et al., 2016;Mánica, 2018). The stress space in the global coordinate system ( s global ) can be transferred to the local stress space ( s local ), using the transformation matrix ‘a’(Eq. (5)). Moreover, the anisotropic stress space ( s ani ) can be obtained by introducing the nonuniform stress scaling factors (C N and C S ), as shown in Eq. (6). Note that the anisotropic strength can be simplified to the isotropic strength when C N ¼C S ¼1(Mánica, 2018). s local ¼a s globalaT(5) s ani ¼2 6 6 6 4 s local 11 .CN s local 12 CS s local 13 s local 12 s local 12 .CNCS s local 23 CS s local 13 CS s local 23 CN s local 33 3 7 7 7 5 (6) where a¼2 4 cos a cos b sin a cos a sin b sin a cos b cos a sin a sin b sin b 0 cos b 3 5 A non-associated flow rule is used in this model. The expressions of failure criterion and plastic potential of the adopted MohrCoulomb model are shown in Eqs. (7) and (8), respectively. Fani MC ¼panisin4þffiffiffiffiffiffiffi Jani 2 qcos q ani 1 ffiffiffi 3 psin4sin q aniccos4 (7) Gani MC ¼panisin j þffiffiffiffiffiffiffi Jani 2 qcos q ani 1 ffiffiffi 3 psin j sin q ani(8) where cand 4are the cohesion and friction angle, respectively; j is the dilatancy angle; p ani ,J ani 2 and q ani represent the mean effective stress, the second invariant of the deviatoric stress tensor, and the Lode angle in the anisotropic stress space ( s ani ), respectively. 3.2. Calibration of model parameters In the adopted anisotropic constitutive model, there are five independent elastic constants (E,E 0 , n , n 0 ,G 0 ), and four independent strength parameters (c,4,C N ,C S ). Time-dependency is not considered, and thus, a small enough value of h vp should be adopted (Song, 2021). Based on numerous tests, h vp is adopted as 100 MPa 5 s in this research, with m¼5. No dilatancy is considered. Fig. 5. (a) Global coordinate system ðx;y;zÞ, (b) Local coordinate system ð1;2;3Þ, and (c) Angles in the transformation matrix, based on the work of Wittke (1990),Mánica et al. (2016) and Mánica (2018). a represents the first rotation around z-axis, and b represents the subsequent rotation around the y-axis. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338 327 However, at least in CODE_BRIGHT simulations, null dilatancy ( j ¼0  ) might result in numerical difficulties. Therefore, a very small value of dilatancy ( j ¼0:01  ) has been adopted in these examples. A simple approach for calibration of the model parameters is proposed and described in this section. As shown in Figs. 6 and 7, the stressestrain curves of anisotropic rock samples can be obtained using strain gauges glued in the appropriate directions. The relationships between the five transversely isotropic elastic constants (E,E 0 , n , n 0 ,G 0 ) and D ε= Ds are displayed in Eqs. (9)-(11). D εx Ds z¼sin2ð2 b Þ 41 Eþ1 E01 G0 n 0 E0cos4 b þsin4 b (9) D εz Ds z¼sin4 b Eþcos4 b E0þsin2ð2 b Þ 42 n 0 E0þ1 G0(10) D εy Ds z¼ n Esin2 b  n 0 E0cos2 b (11) Using the same approach described in former publications (Barla, 1972;Worotnicki, 1993;Amadei, 1996;Chen et al., 1998; Talesnick and Bloch-Friedman, 1999;Hakala et al., 2007;Cho et al., 2012;Alsuwaidi et al., 2021), two elastic moduli (E,E 0 ) and two Poisson’s ratios ( n , n 0 ) can be calculated from rock samples of b ¼0  (Fig. 6a) and b ¼90  (Fig. 6c). After that, the stressestrain curves obtained from samples of 0  < b <90  (Fig. 6b) can be used to determine the shear modulus perpendicular to the isotropic plane (G 0 ). Finally, five independent elastic constants (E,E 0 , n , n 0 ,G 0 ) can be determined. A detailed description of the approach used to derived the elastic constants can be found in the literature (Barla, 1972;Worotnicki, 1993; Amadei, 1996;Chen et al., 1998;Talesnick and Bloch-Friedman, 1999;Hakala et al., 2007;Cho et al., 2012;Alsuwaidi et al., 2021). Then, the four independent plastic parameters (c,4,C N ,C S )canbe evaluated using the Generalised Reduced Gradient non-linear algorithm (GRG method) (Maia et al., 2017), based on all available test results for each type of rock under study, as shown in Fig. 8.As commented in Mánica et al. (2016) and Mánica(2018),C N controls the difference between the strengths of cases for b ¼0  and b ¼90  , while C S affects the strength for samples of 0  < b <90  .Noteagain that the anisotropic strength model can be simplified to the isotropic one when C N ¼C S ¼1. As described in Fig. 8, experimental data ( b ; s global ij ) can be obtained for each uniaxial or triaxial laboratory test. Then, groups of data ( b ; s local ij ) and ( b ; s ani ij ) in the local and anisotropic stress spaces can be sequentially determined, respectively. Consequently, the value of F k can be obtained through inputting the k-th group of experimental data in Eq. (7), i.e. F k represents the value of Fusing the k-th group experimental data. Finally, F sum is chosen as the objective function, where F sum is total value of F k using all groups of experimental data, as shown in Eq. (12). Note that F sum varies with input values of the plastic constants (c,4,C N ,C S ). Fsum ¼X n k¼1 Fk(12) In the determining process, an iteration of the GRG method is carried out, which varies the plastic parameters c,4,C N and C S to obtain the optional (final) plastic constants, where F sum is closest to zero. Consequently, the final values of plastic constants (c,4,C N ,C S ) can be output. Note that the output values of plastic parameters may depend on the initial values and boundary conditions set for the GRG method. However, the finally obtained values are visually checked to be reliable, physically acceptable numbers, consequent with observations in test results. 3.3. Numerical implementation The presented anisotropic constitutive model (Fig. 4) has been implemented into finite element method software CODE_BRIGHT (Olivella et al., 2021). CODE_BRIGHT is developed at the Universitat Politècnica de Catalunya (UPC) and works in combination with the pre-/post-processor GID. GID has been developed by the International Centre for Numerical Methods in Engineering (CIMNE). The total strain rate of the presented anisotropic elasticviscoplastic model can be decomposed into elastic (dε global e =dt) and viscoplastic (dε global vp =dt) parts, with respect to the global coordinate system, as shown in Eq. (13). dεglobal dt¼dεglobal e dtþdεglobal vp dt¼dεglobal e dtþ1 h vpDFani MCmEvGani MC v s global (13) The strain rate of the elastic spring can be expressed as Fig. 6. Conceptual model of rock samples with different anisotropy angles: (a) b ¼0, (b) 0< b <90, and (c) b ¼90. Fig. 7. The recommended locations of strain gauges in the laboratory tests. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338328 dεglobal e dt¼Cglobal e d s global dt¼TT DClocal eTD d s global dt(14) The derivative of the strain rate of the viscoplastic model with respect to the global stress tensor ( s global ) can be expressed as v v s global dεglobal vp dt! ¼1 h vp "mFani MCm1vFani MC v s global vGani MC v s global þFani MCmv2Gani MC v s global2#(15) where vFani MC v s global ¼2 6 40 B @ vFani MC vpani vpani v s ani þvFani MC vffiffiffiffiffiffiffi Jani 2 q vffiffiffiffiffiffiffi Jani 2 q v s ani þvFani MC v q ani v q ani v s ani1 C A T v s ani v s global3 7 5 T vGani MC v s global ¼2 6 40 B @ vGani MC vpani vpani v s ani þvGani MC vffiffiffiffiffiffiffi Jani 2 q vffiffiffiffiffiffiffi Jani 2 q v s ani þvGani MC v q ani v q ani v s ani1 C A T v s ani v s global3 7 5 T v v s global vGani MC v s global! ¼v2Gani MC vffiffiffiffiffiffiffi Jani 2 qv q ani v q ani v s global0 @ vffiffiffiffiffiffiffi Jani 2 q v s global1 A T þvGani MC vffiffiffiffiffiffiffi Jani 2 q v2ffiffiffiffiffiffiffi Jani 2 q v s global2þ vGani MC v q ani v2 q ani v s global2þv2Gani MC v q anivffiffiffiffiffiffiffi Jani 2 q vffiffiffiffiffiffiffi Jani 2 q v s global v q ani v s global!T þv2Gani MC v q ani2 v q ani v s global v q ani v s global!T þv2Gani MC vffiffiffiffiffiffiffi Jani 2 q2 vffiffiffiffiffiffiffi Jani 2 q v s global0 @ vffiffiffiffiffiffiffi Jani 2 q v s global1 A T Due to the gradient discontinuities of the yield surface and the plastic potential, numerical calculations may meet numerical inefficiency when using the standard Mohr-Coulomb model in Eqs. (7) and (8) (Song et al., 2020). Therefore, the smoothed MohrCoulomb model (Abbo and Sloan, 1995;Song et al., 2020)is adopted to improve numerical efficiency and reliability, as shown in Eqs. (16) and (17). Fani MC ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jani 2K2 MC þa2 MCsin24 qþpanisin4ccos4(16) KMC ¼8 > > < > > : AMC þBMCsin3 q ani  q ani> q ani T cos q ani 1 ffiffiffi 3 psin 4sin q ani  q ani q ani T(17) Fig. 8. Process of determining the strength parameters. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338 329 where AMC ¼1 3cos q ani T3þtan q ani Ttan3 q ani T þ1 ffiffiffi 3 ph q aniihtan3 q ani T3tan q ani Tisin4 BMC ¼ 1 3cos3 q ani Th q aniisin q ani Tþ1 ffiffiffi 3 psin4cos q ani T h q anii¼8 < : þ1 q ani 0 1 q ani <0 in which q ani T is a specified transition angle in the smoothed theory (Abbo and Sloan, 1995). The typical values of a MC and q ani T in the smoothed Mohr-Coulomb model are 0 a MC 0:25 and 25   q ani T 30  (Abbo and Sloan, 1995). Note that, if not specified, a MC and q ani T are adopted as 0.14  and 28  , respectively, in this article. Since the second derivative of the viscoplastic potential should also be continuous, the C2 smoothed theory (Abbo et al., 2011;Song et al., 2020) is used to smooth the potential of the viscoplastic model. The adopted smoothed potential can be expressed in Eq. (18). The alternative form of K G ð q ani Þin the vicinity of the singularities can be expressed in Eq. (19). Gani MC ¼panisin j þJani 2K2 G(18) KG¼8 > > < > > : DGþEGsin3 q aniþFGsin23 q ani q ani> q ani T cos q ani 1 ffiffiffi 3 psin j sin q ani  q ani q ani T(19) where FG¼cos3 q ani Tcos q ani T1 ffiffiffi 3 psin jq ani 3h q aniisin3 q ani Th q aniisin q ani T þ1 ffiffiffi 3 psin j cos q ani T18cos33 q ani T EG¼h q aniisin6 q ani Tcos q ani T1 ffiffiffi 3 psin j h q aniisin q ani T 6cos6 q ani Th q aniisin q ani Tþ1 ffiffiffi 3 psin j cos q ani T 18cos33 q ani T DG¼1 ffiffiffi 3 psin j h q aniisin q ani TEGh q aniisin3 q ani T FGsin23 q ani Tþcos q ani T Eqs. (16)e(19) are adapted from the smoothed approximation (C1 and C2 smoothed theories) by Abbo and Sloan (1995) and Abbo et al. (2011). 4. Numerical simulations of anisotropic behaviour of rocks In this section, numerical simulations of uniaxial and triaxial tests are performed on three types of rocks using the anisotropic constitutive model presented in Section 3. Comparisons between numerical and experimental results are carried out for Vaca Muerta shale (Section 4.1), Bossier shale (Section 4.2)(Ambrose, 2014) and NW Spain slate (Section 4.3)(Alejano et al., 2021), in terms of stressestrain curves, apparent elastic modulus (E q ), apparent Poisson’s ratios ( n q and n 0 q ) and peak compressive strengths ( s 1;max ). In this section, the sign conventions are defined as negative for tension and positive for compression. For consistency with experimental results and in line with the typical convention adopted in rock mechanics studies, positive is defined as opposite to the direction of the coordinate axis. Concerning sample preparation, it is challenging for shales and slates. During the process of shale sample preparation, shale cores with significant microcracks were treated with low-viscosity epoxy externally, vacuum suctioned to fill the external microcracks, and cured. Slate samples were prepared with the assistance of a saw disk machine (CEDIMA model CTS-265, 400 mm radius disk), a drilling machine (WEKA, model DK22) and a grinding machine. To obtain these slate cores, 40 kg slab-like blocks have to be the first cut to produce a base that provides the desired schistosity orientation, and then these properly positioned blocks were cored. It is relevant to mention that cutting of these oriented samples was not an easy task, and it was time-consuming. Thus, in the process of coring, numerous cores were broken, particularly when the foliation formed angles of 15  and 30  with the sample bases since the drilling process generate shear stresses that produce breaking of the samples. More than twice the typical quantity of rock material (with granite or sandstone) was needed to produce the roughly 90 tested samples. Note again that the experimental data of stressstrain curves of NW Spain slate are presented in this article for the first time. 4.1. Vaca Muerta shale Ambrose (2014) carried out 21 uniaxial and triaxial tests for Vaca Muerta shale. In this section, numerical simulations using CODE_BRIGHT are performed to fit these laboratory results from Ambrose (2014). The numerical model used (Fig. 9a) is a threeFig. 9. Uniaxial and triaxial numerical tests: (a) Basic features and boundary conditions (conceptual model), and (b) Mesh with 16,980 tetrahedral elements. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338330 dimensional (3D) cylindrical model with a diameter of 0.055 m and a height of 0.11 m. Concerning boundary conditions, along the bottom boundary of the model, all displacements are restrained. Along the top boundary, a constant z-displacement rate 1 10 7 m/s is applied, consistent with that typically applied in laboratory tests, while displacements in the xand y-directions are restrained. Note that the obtained peak strength is independent on the loading rate in this case since time-dependency has not been considered in the simulations. Finally, on the lateral boundary, a constant radial confining stress ( s 3 ) is applied, representing the actual confinement in the cylindrical sampler walls, applied through fluid pressure on the sample enclosed in a rubber sleeve. A mesh with 16,980 tetrahedral elements has been considered (Fig. 9b), sufficiently accurate and resulting in an acceptable computation time. Table 1 lists the input parameters. Fig. 10aed shows the comparison of apparent elastic modulus (E q ), apparent Poisson’s ratio ( n q ), apparent Poisson’s ratio ( n 0 q ), and peak compressive strength ( s 1;max ), respectively, for Vaca Muerta shale between numerical predictions (in terms of lines) and experimental data (in terms of points). Note that apparent elastic modulus (E q ) is the observed stiffness response of the sample, which can be computed as Ds z = D ε z ; the apparent Poisson’s ratios ( n q and n 0 q ) can be computed as n q ¼1=ð D ε z = D ε y Þand n 0 q ¼1= ð D ε z = D ε x Þ. A good agreement concerning stiffnesses and strengths between CODE_BRIGHT results and experimental data can be observed. The root mean square error (RMSE) value of numerical predictions to fit experimental results of Vaca Muerta is 18.45 MPa. As a comparison, the RMSE values of the JPW model and the Pariseau’s model to match experimental data of the Vaca Muerta shale are 12.8 MPa and 16.5 MPa, respectively (Ambrose, 2014). The smaller the values of the RMSE is, the better the agreement between experimental data and the values of prediction is. Therefore, based on the obtained RMSE values, it can be concluded that the adopted anisotropic strength model can reasonably represent the stressestrain anisotropic behaviour of Vaca Muerta shale. In addition, Fig. 11 presents the comparisons of stressestrain curves of Vaca Muerta shale between numerical and experimental results for six samples. A good agreement of stressestrain curves between numerical predictions obtained from CODE_BRIGHT simulations and experimental data obtained from laboratory tests can be observed, validating that the proposed numerical approach and the adopted transversely isotropic stiffness constitutive model can reasonably represent the anisotropic deformability of Vaca Muerta shale. 4.2. Bossier shale Ambrose (2014) carried out 36 uniaxial and triaxial tests for Bossier shale. In this section, numerical simulations are performed to fit the laboratory results of Bossier shale from Ambrose (2014). The geometry, conditions and mesh of the numerical model are the same as those described in Section 4.1 (Fig. 9). The input parameters are listed in Table 1. Fig. 12 presents the comparison of Bossier shale in terms of stiffness (apparent elastic modulus and apparent Poisson’sratios) and peak strength between numerical predictions (in terms of Fig. 10. Comparisons of (a) apparent elastic modulus (E q ), (b) apparent Poisson’s ratio ( n q ), (c) apparent Poisson’s ratio ( n 0 q ), and (d) peak strength results between numerical and experimental results. ‘C_B’represents the CODE_BRIGHT results. ‘Exp’represents the experimental results. Experimental data are obtained from Ambrose (2014). F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338 331 Song, F., Rodriguez-Dono, A., Olivella, S., Zhong, Z., 2020. Analysis and modelling of longitudinal deformation profiles of tunnels excavated in strain-softening timedependent rock masses. Comput. Geotech. 125, 103643. Song, F., 2021. Modelling Time-dependent Plastic Behaviour of Geomaterials. PhD Thesis. Universitat Politècnica de Catalunya (UPC), Spain. Song, F., Rodriguez-Dono, A., 2021. Numerical solutions for tunnels excavated in strain-softening rock masses considering a combined support system. Appl. Math. Model. 92, 905e930. Song, F., Rodriguez-Dono, A., Olivella, S., 2021a. Hydro-mechanical modelling and analysis of multi-stage tunnel excavations using a smoothed excavation method. Comput. Geotech. 135, 104150. Song, F., Rodriguez-Dono, A., Olivella, S., Gens, A., 2021b. Coupled solid-fluid response of deep tunnels excavated in saturated rock masses with a timedependent plastic behaviour. Appl. Math. Model. 100, 508e535. Talesnick, M.L., Bloch-Friedman, E.A., 1999. Compatibility of different methodologies for the determination of elastic parameters of intact anisotropic rocks. Int. J. Rock Mech. Min. Sci. 36, 919e940. Wang, Y.H., Lau, Y.M., Gao, Y., 2014. Examining the mechanisms of sand creep using DEM simulations. Granul. Matter 16, 733e750. Wittke, W., 1990. Rock Mechanics, Theory and Applications-With Case Histories. Springer-Verlag. Worotnicki, G., 1993. CSIRO triaxial stress measurement cell. In: Hudson, J.A. (Ed.), Comprehensive Rock Engineering, vol. 3. Pergamon Press Ltd, Oxford, pp. 329e 394. Xu, G., Gutierrez, M., He, C., Wang, S., 2021. Modeling of the effects of weakness planes in rock masses on the stability of tunnels using an equivalent continuum and damage model. Acta. Geotech. 16, 2143e2164. Zimmerman, R.W., Ambrose, J., Setiawan, N.B., 2018. Failure of anisotropic rocks such as shales, and implications for borehole stability. In: Proceedings of the ISRM International Symposium - 10th Asian Rock Mechanics Symposium. ISRM-ARMS10-2018-263. Dr. Fei Song obtained his PhD in Geotechnical Engineering in 2021 from Universitat Politècnica de Catalunya (UPC), Spain. He is currently a postdoctor at the same university. His research interests include constitutive modelling, rock mechanics, and developing numerical and analytical approaches for the design of underground engineering. F. Song et al. / Journal of Rock Mechanics and Geotechnical Engineering 15 (2023) 323e338338