An anisotropic viscoplasticity model for shale based on layered microstructure homogenization
AAn anisotropic viscoplasticity model for shale based onlayered microstructure homogenization
Jinhyun Choo a, ∗ , Shabnam J. Semnani b , Joshua A. White c a Department of Civil Engineering, The University of Hong Kong, Hong Kong b Department of Structural Engineering, University of California, San Diego, United States c Atmospheric, Earth, and Energy Division, Lawrence Livermore National Laboratory, United States
Abstract
Viscoplastic deformation of shale is frequently observed in many subsurface applications. Many studies have sug-gested that this viscoplastic behavior is anisotropic—specifically, transversely isotropic—and closely linked to thelayered composite structure at the microscale. In this work, we develop a two-scale constitutive model for shalein which anisotropic viscoplastic behavior naturally emerges from semi-analytical homogenization of a bi-layer mi-crostructure. The microstructure is modeled as a composite of soft layers, representing a ductile matrix formed by clayand organics, and hard layers, corresponding to a brittle matrix composed of sti ff minerals. This layered microstructurerenders the macroscopic behavior anisotropic, even when the individual layers are modeled with isotropic constitutivelaws. Using a common correlation between clay and organic content and magnitude of creep, we apply a viscoplasticModified Cam-Clay plasticity model to the soft layers, while treating the hard layers as a linear elastic material tominimize the number of calibration parameters. We then describe the implementation of the proposed model in astandard material update subroutine. The model is validated with laboratory creep data on samples from three gasshale formations. We also demonstrate the computational behavior of the proposed model through simulation oftime-dependent borehole closure in a shale formation with di ff erent bedding plane directions. Keywords:
Viscoplasticity, Anisotropy, Shale, Homogenization, Creep
1. Introduction
Shale often accumulates a significant amount of irrecoverable deformation over time. This viscoplastic defor-mation plays an important role in a variety of subsurface practices, such as in-situ stress estimation [1, 2], boreholedrilling [3], leakage prevention [4, 5], and hydraulic fracturing [6–8]. This importance has motivated a number ofexperimental investigations into creep of shale at multiple scales ( e.g. [9–20]).Experimental results have shown that the viscoplastic behavior of shale is anisotropic—specifically, transverselyisotropic—and closely linked to the shale composition. For example, creep responses of samples from various shaleformations have commonly exhibited a strong dependence on the relative direction between the di ff erential stress andthe bedding plane ( e.g. [1, 2, 11, 20]). Such anisotropy is a natural consequence of the fact that layering is ubiquitousin sedimentary rocks, including shale—see, e.g. , Fig. 1. Sone and Zoback [11] have also found that creep behavioris apparently correlated with the elastic moduli of shale, although the mechanisms underlying viscoplastic and elasticdeformations are quite di ff erent. They explain this correlation by appealing to stress partitioning in a conceptualbi-layer model, in which one layer is composed of soft constituents such as clay and organics (kerogen), while the ∗ Corresponding Author
Email address: [email protected] (Jinhyun Choo) a r X i v : . [ phy s i c s . g e o - ph ] O c t ther layer is composed of hard constituents such as quartz, feldspar, pyrite (QFP), and carbonates. This bi-layerrepresentation is a simple but promising approach. It is well justified by the fact that shale creep is mainly attributedto the clay and organic constituents rather than the hard minerals [9, 11, 13, 16, 17, 19], and it naturally captures thetransverse isotropy resulting from the depositional process. Green River Shale
Ref.
Mehmani et al. (2016), Fuel ft Figure 1: Optical log of core from the Green River Formation [21]. The log samples a 100 ft (30.5 m) vertical interval through the Mahogany Zone,illustrating the strongly bedded nature of shales. Darker sections contain higher volume fractions of organic content than lighter sections. (Figureadapted from Mehmani et al. [21] with permission from Elsevier.)
Mathematically, time-dependent deformation of geomaterials has been described by viscoelastic / plastic contin-uum models ( e.g. [2, 22]) or empirical relations ( e.g. [23]). For clay-rich materials, a number of anisotropic viscoplas-tic constitutive models ( e.g. [24–29]) have been proposed by extending the Modified Cam-Clay (MCC) model origi-nally developed for isotropic materials [30]. Although most of these models have treated anisotropy in a purely macro-scopic manner, without explicit consideration of material microstructure, a few models have incorporated anisotropybased on microstructural approaches. For example, Pietruszczak et al. [31] have proposed a constitutive frameworkthat introduces two microstructure-based quantities, namely, a microstructure evolution parameter capturing creepdeformation associated with the rearrangement of microstructure, and a microstructure tensor describing the directionof transverse isotropy. Very recently, Borja et al. [32] have developed a two-material framework that represents shaleas a mixture of a softer frame and a harder frame, applying an anisotropic critical state plasticity model [33] to theharder frame to simulate anisotropic viscoplastic deformation.Outside the geomechanics community, homogenization techniques have been developed as an elegant way toincorporate microstructure e ff ects on viscoelastic [34, 35] and viscoplastic [36–38] behavior of materials. However,the application of homogenization techniques to viscoplastic materials entails a few major di ffi culties. First, themathematical description of the macroscopic behavior cannot be determined a priori, because the macroscopic viscousbehavior does not necessarily follow the same law as the constituents [39]. For example, a Maxwell-type behaviorfor the microscopic constituents does not necessarily lead to an overall macroscopic behavior following the Maxwelllaw [40]. Second, exact interactions between constituents can only be derived for linear elastic behavior, while newassumptions and strategies are needed for viscoplastic models.A survey of mathematical homogenization techniques for dissipative materials can be found in Charalambakis et al. [41]. The majority of homogenization techniques developed for viscoplastic creep are based on mean-fieldmethods [42], self-consistent schemes [43], and transformation field analysis [44]. These methods were initiallydeveloped for linear elastic composites based on Eshelby’s solution for an ellipsoidal inclusion in an infinite matrix[45], and were later extended to the nonlinear regime by linearization of the local constitutive equations and definitionof a linear comparison composite for a given deformation state [35, 46–52]. In particular, interaction laws for elasto-2iscoplastic inclusions were developed to extend the Mori–Tanaka and self-consistent schemes [43, 53–57]. Specificassumptions such as piecewise uniform [58] or non-uniform [59] distribution of the strain field in each phase wereused to accommodate viscoplasticity within the transformation field analysis framework [60–62].Other classes of homogenization methods include asymptotic expansion and periodic homogenization methods.Periodic homogenization expresses the strain field as the sum of a macroscopic field and a periodic perturbation[40]. This technique was applied to elasto-viscoplastic composites under simplifying assumptions such as uniformmacroscopic stress and strain in the unit cell [63, 64] and point symmetry of internal distributions [65]. Asymptoticexpansion, which forms the basis of the present work, has been applied to viscoelasticity [34, 66] and viscoplasticty[67] by assuming viscoplastic strains are piecewise constant eigenstrains, i.e. strains not caused by external stresses.The above techniques are generally applicable to isotropic composites made of ellipsoidal inclusions within amatrix, while homogenization-based anisotropic viscoplasticity has received little attention in the literature [68–71].Chatzigeorgiou et al. [68] applied asymptotic expansion techniques to composites with generalized standard materiallaws and demonstrated their application to multi-layered media. Mathematical approaches have also been adopted forperfectly bonded isotropic creep materials [70] and fiber-reinforced laminates [69].In this work, we develop a two-scale constitutive model for shale in which anisotropic viscoplastic behavior nat-urally emerges from semi-analytical homogenization of a bi-layer microstructure. The microstructure is a compositeof soft layers, representing a ductile matrix formed by clay and organics, and hard layers, corresponding to a brittlematrix composed of QFP and carbonates, similar to the conceptual model of Sone and Zoback [11]. In particular,we build the model using an inelastic homogenization framework for layered materials recently proposed by Semnaniand White [71], here extending it to viscoplastic materials. One of the main advantages of this homogenization ap-proach is that the macroscopic behavior is inherently anisotropic, even when the individual layers are modeled withisotropic constitutive laws. This allows us to adopt standard isotropic models with physically meaningful parameters.This provides a significant advantage over many single-scale models that require anisotropy parameters that can bechallenging to calibrate and physically justify. Using a common correlation between clay and organic content andmagnitude of creep, we apply a viscoplastic Modified Cam-Clay plasticity model to the soft layers, while treating thehard layers as a linear elastic material. It will be shown that the proposed model can simulate viscoplastic deforma-tion of shale remarkably well with a relatively small number of parameters. Through a borehole closure example, wewill also demonstrate that the model can be e ffi ciently applied to general initial–boundary value problems involvinganisotropic deformation of shale over time.The proposed model is distinguished from existing shale models in several aspects. First, the majority of existingviscoplasticity models for shale assume isotropic behavior [72, 73]. Although an isotropic viscoplasticity modelmay reproduce uniaxial creep deformation, it requires a di ff erent parameterization when the bedding plane directionis changed; see Haghighat et al. [72] for example. Second, it di ff ers from the anisotropic viscoplasticity modelof shale recently proposed by Borja et al. [32]. That work relies on solid–solid mixture theory, while the presentwork employs asymptotic expansion. Both represent shale as a composite of harder and softer materials. In Borja etal. [32], however, the materials are not layered, requiring at least one of them be modeled by a single-scale anisotropicconstitutive law. The model proposed here simplifies calibration, since anisotropic properties naturally emerge fromthe layered microstructure. Third, the proposed model di ff ers from other homogenization-based models in that itprovides a general and flexible framework to accommodate any desired behavior for the constituents ( e.g. brittle vs.ductile, rate-dependent vs. rate-independent). Its algorithmic implementation in a computational mechanics code isalso quite straightforward. Finally, the viscoplastic deformation considered in this work is an intrinsic process of thesolid matrix, rather than induced by pore fluid di ff usion or other physio-chemical processes studied in other works( e.g. [74–77]).The remainder of the paper is organized as follows. Section 2 describes the model formulation. Section 3 presentsalgorithms for implementing the proposed model in a material update subroutine. Section 4 validates the model with3aboratory creep data of samples from three gas shale formations, namely, Haynesville, Eagle Ford, and Barnett.Section 5 demonstrates the computational performance of the model through numerical simulation of time-dependentborehole closure in a shale formation with di ff erent bedding plane directions. Section 6 closes the work.As for sign conventions and notations, stresses and strains are positive in compression following the geomechanicssign convention. Bold-face letters denote tensors and vectors. The symbol “ · ” denotes single contraction, and “:”denotes double contraction.
2. Formulation
This section begins by describing a bi-layer model of shale, which is a particular class of the inelastic homogeniza-tion framework developed by Semnani and White [71]. For brevity, the complete derivation of the homogenizationframework will be omitted, as it is extensively described in the prior work. Specific constitutive models will thenbe introduced to model the stress-strain response of the two layers at the microscale, aimed at capturing viscoplasticdeformation in shale with minimal ingredients. Quasi-static behavior and infinitesimal deformation will be assumedthroughout.
Consider a body B whose microstructure is described by a spatially periodic unit cell U , as illustrated in Fig. 2.The unit cell is composed of two parallel layers, namely, a soft layer and a hard layer. The soft layer represents thecompliant constituents of shale, such as clay and organic matters. The hard layer corresponds to the sti ff constituentsof shale, including quartz, feldspar, pyrite (QFP) and carbonate. x x x h a r d l a y e r ( Q F P & c a r b o n a t e s ) s o ft l a y e r ( c l a y a n d o r g a n i c s ) y y y shale n θ Figure 2: An illustration of a unit cell comprised of a soft layer and a hard layer. The macroscopic coordinate system is x = { x , x , x } , andthe microscopic coordinate system is y = { y , y , y } . Symbols n and θ denote the unit normal vector of the bedding plane and its angle from thevertical, respectively. A few assumptions are introduced to this bi-layer model. First, the characteristic length-scale of the body(macroscale) and the unit cell (microscale) are clearly separated, a key requirement under asymptotic homogenizationtheory. Second, each layer is a homogeneous, standard continuum. Third, the two layers are perfectly bonded to eachother, such that the displacement field is continuous across the layers as well as within the layers. Note, however, thatcertain strain components may be discontinuous at the layer interface.From a practical point of view, we also note that only quasi-periodicity of the microstructure is strictly required.The material properties of the medium may vary slowly over large-length scales, as long as the local properties at shortscales may be well-approximated via a periodic unit cell. In a boundary value problem, for example, the microstructureat each integration point is assumed to be periodic, but the material properties assigned from one integration point tothe next may vary slowly to capture large-scale heterogeneities like those seen in Fig. 1.4o describe the kinematics of the two-scale model, we introduce a macroscopic coordinate system, x , and amicroscopic coordinate system, y . The two coordinate systems are related as y = x /(cid:15) , (1)where (cid:15) (cid:28) y : = y · n , where n is theunit normal vector of the layers. Let us denote by u ( x , y ) the displacement field, with the dual argument emphasizingthat this field is a function of both the macroscopic and microscopic coordinates. Applying two-scale asymptotichomogenization, we approximate u ( x , y ) as a truncated power series of (cid:15) , i.e. u ( x , y ) ≈ u (0) ( x ) + (cid:15) u (1) ( x , y ) . (2)Here, u (0) represents the macroscopic component of the displacement field, which is only a function of the macroscopiccoordinate. In contrast, u (1) represents a microscopic displacement fluctuation, which varies quasi-periodically acrossunit cells. Accordingly, the total strain field can be additively decomposed as ε = E + e , (3)where E and e are the macroscopic and microscopic parts of the strain field, given by E : = sym (cid:32) ∂ u (0) ∂ x (cid:33) , e : = sym (cid:32) ∂ u (1) ∂ y ⊗ n (cid:33) , (4)with sym( · ) denoting the symmetric part of a tensor. The strain field is potentially discontinuous at the layer interfaces,despite the continuity of displacement field throughout the unit cell. In Semnani and White [71] the possibility of inter-layer slip or normal opening is considered, but here we focus on perfectly bonded layers. The strain field must thenbe continuous in the direction tangential to the layers.The governing equations of the boundary value problem are derived as in Semnani and White [71]. This procedurewill be omitted for brevity. To express the equations, let us introduce layer-wise quantities, using subscripts ( · ) s and( · ) h to denote quantities in the soft and hard layers, respectively. We denote by φ s and φ h the volume fractions of thesoft and hard layers, which satisfy φ s + φ h =
1. We also define the microscopic displacement gradient vectors of thetwo layers as v s : = ∂ u (1)s ∂ y , v h : = ∂ u (1)h ∂ y . (5)To satisfy the bedding-normal traction continuity requirement, we note that the microscopic displacement gradientsmust be constant in each layer [71]. The total strain tensors of the layers are then defined as ε s : = E + sym( v s ⊗ n ) , ε h : = E + sym( v h ⊗ n ) . (6)These strains are related to the stresses in the individual layers as˙ σ s = C s : ˙ ε s , ˙ σ h = C h : ˙ ε h , (7)where σ and C denote the stress tensor and the stress–strain tangent tensor within each phase, respectively. Thegeneral rate form is intended to accommodate potentially nonlinear material behavior. The particular models chosenfor the layers do not impact to the overall homogenization algorithm, so we leave them unspecified for the moment.5e denote by Σ the macroscopic stress tensor, which is defined as the volume average of the stress tensors in eachlayer, i.e. Σ : = φ s σ s + φ h σ h . (8)Introducing the multiscale expansions into the global initial-boundary-value problem, we may arrive at two sep-arate but coupled problems: a microscale problem in the unit cell U and a macroscale problem in the body B . Thegoverning equations of the microscale problem enforce certain stress and displacement field compatibilities: σ s · n − t = in U s (traction continuity in the soft layer) , (9) σ h · n − t = in U h (traction continuity in the hard layer) , (10) φ s v s + φ h v h = in U (periodicity of micro-displacements) . (11)Here, t = Σ · n denotes the traction vector resulting from stresses acting in the bedding normal direction, which mustbe continuous in all layers to satisfy unit cell equilibrium.The governing equation of the macroscopic problem is identical to the linear momentum balance typically em-ployed in single-scale models, div Σ + F = in B . (12)Here, div( · ) is the divergence operator with respect to the macroscopic coordinate system x , and F is the body forcevector. This governing equation is supplemented with appropriate initial conditions in B and boundary conditions on ∂ B to define a complete problem. We may now begin to specialize the general model by defining isotropic constitutive models for the layers. Forthe soft layer composed of clay and organic content, we use a viscoplastic model. This layer is then responsiblefor the time-dependent ductile deformation of the overall shale. The viscoplastic formulation begins by writing thestress–strain relationship as ˙ σ s = C es : ( ˙ ε s − ˙ ε vps ) , (13)where C es is the elastic stress–strain tangent of the soft layer, and ε vps is the viscoplastic part of the strain in the softlayer. For simplicity, we assume that the elastic part of the soft layer deformation is linear elastic. The linear elastictangent can be written as C es = K s ( ⊗ ) + G s (cid:32) I − ⊗ (cid:33) , (14)with K s and G s denoting the bulk and shear moduli of the soft layer, respectively. Also, and I denote the second-orderidentity tensor and the fourth-order symmetric identity tensor, respectively. While the elastic behavior of clay-richmaterials can be highly nonlinear and pressure-dependent [78, 79], the assumption of linear elasticity allows us tominimize the number of model parameters when viscoplastic deformation is of primary interest. We note, however,that the standard Modified Cam-Clay model typically employs a pressure-dependent elastic response rather than alinear model.For the viscoplastic part of deformation, we adopt the modeling approach proposed by Duvaut–Lions [80], whichis one of the most popular approaches to viscoplasticity in the literature ( e.g. [32, 81–84]). The Duvaut–Lions vis-6oplasticity is chosen mainly for its relative simplicity and computational robustness. As shown in Simo et al. [81],the Duvaut–Lions approach can be easily applied to extend any rate-independent plasticity model to viscoplasticity,regardless of the smoothness of yield surface. We note, however, that other well-established approaches to viscoplas-ticity, such as the Perzyna formulation [85], can also be cast into the present homogenization framework. Indeed,Borja et al. [32] have shown that the Duvaut–Lions and Perzyna approaches can produce nearly the same viscoplasticbehavior when their parameters are properly selected.Central to the Duvaut–Lions viscoplasticity is the notion of backbone stress, which is the closest-point projectionof stress onto the yield surface of a rate-independent plasticity model. In a Duvaut–Lions model, the backbone stresstensor is first computed based on a rate-independent model, followed by calculation of the (actual) stress tensor byintegrating the rate equations accommodating viscous e ff ects.For the rate-independent model determining the backbone stress tensor, here we choose the Modified Cam-Claymodel [30], which is widely used for clay-rich materials including shale ( e.g. [33, 86–88]). We denote by p and q thevolumetric and deviatoric stress invariants, respectively, which are defined as p : =
13 tr( σ s ) , q : = (cid:114) (cid:107) σ s − p (cid:107) . (15)Using these two stress invariants, the yield function can be written as f ( p , q , p c ) = q M + p ( p − p c ) ≤ , (16)where M > p c > p c = p c λ ˙ ε vpv , (17)where ε vpv : = tr( ε vp ) is the volumetric part of the viscoplastic strain. As standard in the Modified Cam-Clay model,the plastic flow will be assumed to be associative.The extension of this model to Duvaut–Lions viscoplasticity is carried out through the rate equations [81]˙ ε vps = τ ( C es ) − : ( σ s − ¯ σ s ) , (18)˙ ε vpv = − τ (cid:32) λ p c (cid:33) ( p c − ¯ p c ) , (19)where the bar denotes the inviscid part of a quantity calculated by closest-point projection onto the rate-independentyield surface. Note that while Eq. (18) is general for any plasticity model, Eq. (19) is a particular case of the rateequation for internal variables [81] specialized to the Modified Cam-Clay model. Substituting Eq. (17) into Eq. (19)gives ˙ p c = − τ ( p c − ¯ p c ) . (20)Also, τ is a relaxation time parameter, which may be considered a constant or a function of other quantities asdescribed in Simo et al. [81]. In the context of isotropic Perzyna-type viscoplastic modeling of shale, Haghighat etal. [72] have proposed an exponential relationship between the viscosity coe ffi cient and the volumetric viscoplastic7train ε vpv . Adapting this relationship to Duvaut–Lions viscoplasticity, we consider the following form: τ = τ exp( ζε vpv ) . (21)Here, τ is the initial relaxation time, and ζ is a parameter determining the rate of evolution of the relaxation parameter. The hard layer, which represents a matrix composed of QFP and carbonate minerals, typically deforms in a brittleand rate-independent manner. Because brittle failure of the hard layer is outside the focus of this work, here wesimply model the hard layer as a linear elastic material, reducing the overall number of material parameters. Then themicro-stress tensor of the hard layer, σ h , is related to the micro-strain tensor of the hard layer, ε h , as σ h = C eh : ε h , (22)where C eh = K h ( ⊗ ) + G h (cid:32) I − ⊗ (cid:33) (23)is the elastic stress–strain tangent of the hard layer. Again, K h and G h denote the bulk and shear moduli of the hardlayer, respectively.Note that if the brittle failure of the hard layer is also of interest, the linear elastic model should be replaced by aplasticity (or damage) model for brittle geomaterials. This replacement can be made without any significant changeto the homogenization procedure. Remark . In the absence of viscoplastic deformation, the proposed model becomes equivalent to the bi-layer elasticmodel proposed by Backus [89].
3. Implementation
This section describes how to implement the proposed two-scale model, with a focus on the implementation ofthe material update subroutine. Its implementation in a continuum mechanics code is quite straightforward. Considera load step from time t n to t n + , in which the goal is to determine the macroscopic displacement field u n + for whichboth the macroscopic and microscopic problems are satisfied. When using an implicit finite element method, we adopta two level iteration: a global Newton’s method is used to solve a discretized version of the macroscopic boundaryvalue problem (12), while sub-iterations are used at each material point to solve the microscale balances (9)–(11) todetermine the local stress-strain response.We begin by assuming that an estimate for the discrete macroscale displacement u (0) n + , m is given, where subscript( · ) m is a global iteration counter. The macro-strain E at a quadrature point is therefore known via (4) and used as aninput to a material subroutine. This routine uses the procedure prescribed below to determine the macroscopic stress Σ and tangent sti ff ness C Σ satisfying the local equilibrium and compatibility equations. While the internal computationsdi ff er, the strain-driven interface is the same as for any other material subroutine. Global residual equations may thenbe assembled, and an updated displacement field u (0) n + , m + computed if necessary. The procedure is iterated until allbalance equations are satisfied to a desired tolerance.In the following, we describe a procedure to numerically solve the microscale problem, which is a particular caseof the solution algorithm presented in Semnani and White [71]. We briefly summarize the key steps here, but furtherdetails can be found in the prior work. We also describe the integration procedure for the viscoplastic model of the8oft layer. For brevity, we will denote quantities at t n + without an additional subscript, while distinguishing quantitiesat t n with a subscript ( · ) n . We also use Newton’s method to solve the microscale problem, considering its potential nonlinearity due to vis-coplasticity. Newton’s method for iteratively solving a system of nonlinear equations R ( X ) = proceeds in the twosteps: solving J k ∆ X = − R k , (24)updating X k + = X k + ∆ X . (25)Here, R is the residual vector, X is the unknown vector, ∆ X is the search direction, J is the Jacobian matrix, andthe superscript ( · ) k is a local iteration counter. For this particular problem, the residual vector consists of the threegoverning equations of the microscale problem, namely, Eqs. (9)–(11), and the unknown vector is an array of the threeprimary unknowns, namely, v s , v h , and t . Specifically, R ( X ) : = n · σ s − tn · σ h − t φ s v s + φ h v h → , X : = v s v h t . (26)Note that the micro-stresses in the residual equations are a function of both the microscopic strains and the fixed macroscopic strain via equations (6) and (7). The Jacobian matrix is given by J : = ∂ R ∂ X = n · C s · n − n · C h · n − φ s φ h , (27)where C s and C h are the consistent (algorithmic) tangent operators for the micro-constitutive laws in the soft and hardlayers, respectively. While C h = C eh for the linear elastic hard layer, C s should be distinguished from the continuum(theoretical) tangent of the viscoplasticity model for an optimal convergence rate during the Newton iteration. Aspecific expression for C s will be provided in Eq. (35) later in this section.When the material update is performed within a global iteration algorithm ( e.g. for a stress-driven problem ornonlinear finite element analysis), the macroscopic tangent operator, ∂ Σ /∂ E , is required for convergence of the globaliteration. The macroscopic tangent operator can be written as C Σ : = ∂ Σ ∂ E = φ s ∂ σ s ∂ E + φ h ∂ σ h ∂ E , (28)with ∂ σ s ∂ E = C s − ( C s · n ) · (cid:16) J − · n · C s + J − · n · C h (cid:17) , (29) ∂ σ h ∂ E = C h − ( C h · n ) · (cid:16) J − · n · C s + J − · n · C h (cid:17) . (30)Here, J − · )( · ) is the inverse of the Jacobian matrix (27) with subscripts denoting the appropriate sub-blocks of J − (which has a 3 × C s and C h have emerged from the interaction between the two layers.A complete material update procedure is summarized in Algorithm 1. Note that tol (cid:28) Remark . It should be noted that an additional level of nested iterations is required when using nonlinear materialmodels for the layers (as here for the viscoplastic layer) because sub-iterations are required when calling the up-date routine for the layer stresses. While deep nesting can be expensive, the material subroutine costs remain smallcompared to other components of an implicit finite element solution procedure, and they may be trivially parallelized.
Algorithm 1
Material point update.
Input: E Output: Σ and C Σ Initialize the iteration counter k = X k = . Calculate the total strains ε s = E + sym( v s ⊗ n ) and ε h = E + sym( v h ⊗ n ). Get σ s , σ h , C s and C h by passing ε s and ε h into material subroutines for the microscale constitutive models. Assemble the residual vector R k ( X k ) as Eq. (26). if (cid:107) R k (cid:107) > tol then Solve for the search direction ∆ X = − ( J k ) − R k . Perform a Newton update X k + = X k + ∆ X . Set k ← k + end if Update the macro-stress Σ = φ s σ s + φ h σ h . Compute the macroscopic tangent operator C Σ via Eq. (28). As described above, the solution of the bi-layer problem requires the micro-stress tensor and the consistent tangentoperator of the individual layers. For the hard layer, which is linear elastic, these two quantities can be obtainedanalytically and trivially. For the soft layer, however, they need to be computed by an integration algorithm for aviscoplastic model of Duvaut–Lions type.Time integration of a Duvaut–Lions viscoplastic model proceeds in two steps: (i) integration of the rate-independentplasticity model for the backbone stress tensor, and (ii) integration of the viscoplastic rate equations for the actualstress tensor and the internal variable(s). For brevity, we skip an algorithm for the first step—integration of the rate-independent Modified Cam-Clay model—because it is standard and well described in many references ( e.g. [90]).For the second step, we use the implicit Euler method, except that we evaluate τ in Eq. (21) explicitly using ε vpv atthe previous time step. This semi-implicit approach greatly simplifies the integration, without significant compromisein accuracy because ε vpv does not increase dramatically in a load step. Let us denote the time increment by ∆ t : = t n + − t n . The viscoplastic rate equations (18) and (20) are then integrated as ε vps = ( ε vps ) n + ∆ t τ n ( C es ) − : ( σ s − ¯ σ s ) , (31) p c = ( p c ) n − ∆ t τ n ( p c − ¯ p c ) . (32)10earranging these two equations, we get σ s = σ trs + ( ∆ t /τ n ) ¯ σ s + ∆ t /τ n , (33) p c = ( p c ) n + ( ∆ t /τ n ) ¯ p c + ∆ t /τ n . (34)Here, σ trs : = C es : [ ε s − ( ε vps ) n ] is the trial stress calculated assuming that the strain increment in the load step is fullyelastic. Also, the consistent tangent operator for σ s can be written as [82] C s = C es + ( ∆ t /τ n ) ¯ C s + ∆ t /τ n , (35)where ¯ C s is the consistent tangent operator for the backbone stress ¯ σ s , which can be obtained by an existing algo-rithm for the rate-independent model. The complete procedure for updating the viscoplastic model is summarized inAlgorithm 2. Algorithm 2
Viscoplastic update of the soft layer.
Input: ε s and ∆ t Output: σ s and C s Compute the trial stress σ trs : = C es : [ ε s − ( ε vps ) n ] and set p c = ( p c ) n . Evaluate the yield function f in Eq. (16) with the trial stress. if f < then Elastic step. Set σ s = σ trs and C s = C es . else Viscoplastic step. Compute ¯ σ and ¯ p c using an algorithm for the rate-independent Modified Cam-Clay model. Calculate σ s from Eq. (33), p c from Eq. (34), and C s from Eq. (35). end if4. Validation This section validates the proposed model using creep test data on samples from three gas shale formations,namely, Haynesville, Eagle Ford, and Barnett, which have been obtained by Sone and Zoback [1, 11]. Each loadingstage of the experiments was conducted as follows. First, a cylindrical sample was subjected to an isotropic confiningpressure, which was maintained for ∼ ff erential stress,which was applied incrementally for 60 seconds and then held constant for ∼ ff erential stresses were varied by test, with complete informationavailable in Sone [91]. In what follows, we simulate these tests using the proposed model and algorithms describedabove. To begin, we consider the experimental data of the Haynesville-1V sample in Sone and Zoback [1], which providesthe total (instantaneous plus creep) strain in both the axial and lateral directions after axial loading. The sample wascored along the vertical direction, so it was axially loaded in the direction normal to the bedding plane ( i.e. θ = ◦ ).The confining pressure was 30 MPa, and the di ff erential stress was increased by 29 MPa.11able 1 summarizes the model parameters calibrated to reproduce the experimental data of the Haynesville shalesample. These parameters were obtained through the following procedure. First, the fraction of the soft layer, φ s ,was calculated as the average of the clay and kerogen volume fractions of the Haynesville shale samples in Soneand Zoback [92]. The fraction of the hard layer was then given by φ h = − φ s . The CSL slope, M , was computedusing the sliding friction coe ffi cient reported in Sone and Zoback [11] and its relationship with M under triaxialcompression. For the elasticity parameters of the two layers, we first estimated their potential ranges based on the datain Sone and Zoback [11, 92], and then further calibrated them to fit the instantaneous strains in the axial and lateraldirections. Then, we calibrated the rate-independent plasticity parameters of the soft layer, initial p c and λ , to matchthe final strain of the experimental data, inferring their possible ranges from the estimated in-situ vertical stress inSone and Zoback [92] and Modified Cam-Clay parameters for other shales ( e.g. [32, 33, 88]). It is noted that p c wasalso calibrated due to significant uncertainties in the in-situ horizontal stresses and possible sample disturbance, seeHaghighat et al. [72] for a similar treatment. Lastly, the relaxation time parameters, τ and ζ , of the soft layer werecalibrated to fit the time evolution of creep strain in the experiment.Parameter Symbol Units Soft Layer Hard Layer(Viscoplastic MCC) (Linear elasticity)Volume fraction φ m – 0.455 0.545Bulk modulus K GPa 5.9 41.0Poisson’s ratio ν – 0.28 0.36CSL slope M – 0.9 –Preconsolidation pressure p c MPa 18 –Hardening modulus λ MPa 0.0008 –Initial relaxation time τ s 20 –Relaxation time coe ffi cient ζ – 3400 – Table 1: Parameters for Haynesville shale calibrated using experimental data from Sone and Zoback [1].
Figure 3 compares the simulation results and the experimental data in terms of the axial and lateral strains fromthe beginning of axial loading. The simulation results show an excellent agreement with the experimental data fromthe instantaneous deformation phase (the first 60 seconds) to the creep phase (the remaining period), in both axial andlateral directions. It should be emphasized that although an isotropic viscoplasticity model may also be able to mimicthe axial strain data, it is unable to accurately simulate the lateral strain data.
To further validate the proposed model, we use the experimental data of the Eagle Ford-1 samples in Sone andZoback [1]. Unlike the previous example, the Eagle Ford data were obtained from two di ff erent samples, one coredvertically and the other horizontally. Thus, the vertical sample was axially loaded in the bedding plane normal di-rection ( θ = ◦ ), whereas the horizontal sample was loaded in the bedding plane parallel direction ( θ = ◦ ). Theconfining pressure applied to the two samples were the same (10 MPa, see Sone [91]), but the di ff erential stresseswere slightly di ff erent: 16 MPa for the vertical sample and 17 MPa for the horizontal sample. The model parametersfor the Eagle Ford shale data were obtained through the same procedure as that for the Haynesville shale example andare presented in Table 2.In Fig. 4, the simulation results are compared with the experimental data of the vertical and horizontal Eagle Fordshale samples. The comparison shows that the proposed model captures the di ff erent deformation responses of the12 − S t r a i n ( % ) Figure 3: Haynesville shale: experimental data from Sone and Zoback [11] (open circles) and simulation results using the proposed model (solidlines).
Parameter Symbol Units Soft Layer Hard Layer(Viscoplastic MCC) (Linear elasticity)Volume fraction φ m – 0.260 0.740Bulk modulus K GPa 5.1 49.0Poisson’s ratio ν – 0.39 0.29CSL slope M – 0.87 –Preconsolidation pressure p c MPa 18 –Hardening modulus λ MPa 0.01 –Initial relaxation time τ s 30 –Relaxation time coe ffi cient ζ – 5000 – Table 2: Parameters for Eagle Ford shale calibrated using experimental data from Sone and Zoback [1]. ( θ = ○ ) horizontal ( θ = ○ ) Time (s) S t r a i n ( % ) Figure 4: Eagle Ford shale: experimental data from Sone and Zoback [1] (open circles) and simulation results using the proposed model (solidlines). Note that the di ff erential stresses for the vertical and horizontal samples are di ff erent: 16 MPa and 17 MPa, respectively. As our last validation example, we use the creep test results of the Barnett-1H shale sample in Sone and Zoback [11].The sample was cored horizontally, and it was subjected to an isotropic stress of 20 MPa followed by an increase indi ff erential stress of 48 MPa. Because Sone and Zoback [11] reported creep strains only ( i.e. without instantaneousstrains), the model parameters are calibrated to the creep strains in the axial and lateral directions. The lateral strainin this case is defined as the average of radial strains in the bedding plane parallel and normal directions, in the sameway as in the experimental data of Sone and Zoback [11]. The calibrated parameters for the Barnett shale data arepresented in Table 3.Parameter Symbol Units Soft Layer Hard Layer(Viscoplastic MCC) (Linear elasticity)Volume fraction φ m – 0.470 0.530Bulk modulus K GPa 5.9 43.0Poisson’s ratio ν – 0.12 0.25CSL slope M – 1.2 –Preconsolidation pressure p c MPa 25 –Hardening modulus λ MPa 0.005 –Initial relaxation time τ s 200 –Relaxation time coe ffi cient ζ – 5500 – Table 3: Parameters for Barnett shale calibrated using experimental data from Sone and Zoback [11]. et al. [32] has simulated the same experimental data very well. However, the bi-material model required sevenmore parameters than the present bi-layer model, mainly because the anisotropic elasto-plastic model used in the bi-material model requires six parameters to represent anisotropy. Conversely, the present bi-layer model does not requireany parameters for anisotropy other than the phase volume fractions, because it accommodates anisotropy through thelayered homogenization procedure. The current procedure may be more computationally expensive, however, due tothe need for nested iterations. − ⋅ − axiallateralTime (s) C r ee p s t r a i n ( % ) Figure 5: Barnett shale: experimental data from Sone and Zoback [11] (open circles) and simulation results using the proposed model (solid lines).
5. Borehole Study
This section briefly demonstrates the performance of the proposed model for general initial–boundary value prob-lems. For this purpose, we simulate the the time-dependent closure of a borehole in shale, which is relevant towellbore stability [3] and leakage prevention [4, 5]. We consider a borehole with a radius of 0.1 m, represented in 2Das a circle inside a square domain in plane strain conditions. The shale is modeled using the parameters calibrated forthe Haynesville shale sample in the previous section. To focus on the e ff ect of material anisotropy, we introduce afew simplifying assumptions: the shale is normally consolidated under an isotropic stress field, the borehole is drilledinstantaneously, and the shale deforms under fully drained conditions so that no fluid response must be modeled. Allexternal boundaries are fixed, while the borehole surface is a traction-free boundary. Using the standard displacement-based finite element method in conjunction with Newton’s method, we solve this problem with 1,944 quadrilateralbilinear elements and a constant time increment of 0.1 minute. The deal.II library [93, 94] is used for the finiteelement solution. To examine the e ff ect of shale anisotropy on the borehole closure behavior, we repeat the sameproblem with three di ff erent bedding plane directions: θ = ◦ , θ = ◦ , and θ = ◦ .Figure 6 presents the simulated time evolution of the borehole geometry and the equivalent plastic strain field( √ / (cid:107) ε vp (cid:107) ). It can be seen that in all cases the borehole progressively closes over time, but the closure patterns di ff eraccording to the bedding plane direction. The boreholes squeeze more in the bedding plane normal direction than the15arallel directions to become oval-shaped, even though the initial geometry and stress condition had no significantdirectional dependence. Consequently, the bedding plane direction controls where viscoplastic zones emerge. Theclosure process in this problem nears steady-state after 30 minutes, but the real process in the subsurface wouldbe delayed significantly due to the slow dissipation of pore pressure in low-permeability shale. The fundamentaldeformation patterns would, however, remain similar. θ = ○ θ = ○ θ = ○ equivalent plastic strain (%) Figure 6: Borehole closure simulation: time evolution of the borehole geometry and equivalent plastic strain field with di ff erent bedding planedirections. Displacements are exaggerated by a factor of 200. To demonstrate the computational e ffi ciency and robustness of the proposed model, Fig. 7 presents the conver-gence profiles during the Newton iterations in a few earlier steps where viscoplastic deformations rapidly develop.As can be seen, the relative residual norms are quickly reduced to small values after only four Newton iterations.While not presented, as the viscoplastic process nears completion, the Newton iterations converge even faster. Theconvergence behavior is near-optimal in every step, confirming the correctness of the macroscopic and microscopictangent operators. 16 a) θ = ◦ − − − − − Iteration R e l a t i v e r e s i d u a l n o r m Step 1Step 10Step 30 (b) θ = ◦ − − − − − Iteration R e l a t i v e r e s i d u a l n o r m Step 1Step 10Step 30 (c) θ = ◦ − − − − − Iteration R e l a t i v e r e s i d u a l n o r m Step 1Step 10Step 30
Figure 7: Borehole closure simulation: convergence profiles of Newton iterations.
The results of this example highlight that material anisotropy alone has a profound impact on borehole deforma-tions, and that the proposed model can be easily incorporated into a standard numerical code for initial–boundaryvalue problems. It is also noted that the model could be incorporated into more sophisticated numerical methods forcoupled poromechanical problems ( e.g. [95–98]) to account for fluid flow e ff ects on viscoplastic deformation. Futurework will apply the model to investigate the combined control of the bedding plane direction, stress anisotropy, andfluid flow on the time-dependent behavior of boreholes in shale.
6. Closure
This work has proposed an anisotropic viscoplasticity model for shale, extending a recently developed homoge-nization framework for layered materials [71]. The proposed model conceptualizes the microstructure of shale as acomposite of a viscoplastic layer and an elastic layer. The homogenization approach allows the individual layers to bedescribed by familiar isotropic constitutive models, with anisotropy emerging as an inherent property of the layeredmicrostructure. This combination of layer-wise isotropic models also helps minimize the number of free materialparameters, all of which have clear physical meanings and can be determined without significant di ffi culty. The pro-posed model was validated using experimental creep data from three gas shale formations, showing good agreementfor all strain components. We have also demonstrated that the model can be e ffi ciently incorporated within standardnumerical workflows, with simulation of a borehole closure process as a representative example. The remarkable suc-cess of the proposed model suggests that homogenization is a powerful approach to the modeling of physical behaviorin layered geomaterials, which remains a di ffi cult challenge to accurately predicting the response of geomechanicalsystems. Acknowledgments
JC acknowledges financial support from the Research Grants Council of Hong Kong under Project 27205918. SJSand JAW were supported by Total S.A. through the FC-MAELSTROM project. Portions of this work were performedunder the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under ContractDE-AC52-07-NA27344. The authors wish to thank Yashar Mehmani and Alan Burnham for providing Figure 1and for helpful discussions regarding shale anisotropy. The authors are also grateful to Hiroki Sone for sharing hisexperimental data and for valuable discussions regarding the experimental procedure.17 eferences [1] H. Sone, M. D. Zoback, Time-dependent deformation of shale gas reservoir rocks and its long-term e ff ect on the in situ state of stress,International Journal of Rock Mechanics and Mining Sciences 69 (2014) 120–132.[2] H. Sone, M. D. Zoback, Viscous relaxation model for predicting least principal stress magnitudes in sedimentary rocks, Journal of PetroleumScience and Engineering 124 (2014) 416–431.[3] P. Horsrud, R. M. Holt, E. F. Sonstebo, G. Svano, B. Bostrom, Time dependent borehole stability: laboratory studies and numerical simulationof di ff erent mechanisms in shale, in: Rock Mechanics in Petroleum Engineering, Society of Petroleum Engineers, 1994, pp. SPE–28060–MS.[4] P. Cerasi, E. Lund, M. L. Kleiven, A. Stroisz, S. Pradhan, C. Kjøller, P. Frykman, E. Fjær, Shale creep as leakage healing mechanism in CO sequestration, Energy Procedia 114 (2017) 3096–3112.[5] X. Xie, E. Fjær, E. Detournay, Time-dependent closure of a borehole in a viscoplastic rock, Geomechanics for Energy and the Environment19 (2019) 100115.[6] H. Asala, M. Ahmadi, A. Taleghani, Why re-fracturing works and under what conditions, in: SPE Annual Technical Conference and Exhibi-tion, Society of Petroleum Engineers, 2016, pp. SPE–181516–MS.[7] Y. Yang, M. Zoback, Viscoplastic deformation of the Bakken and adjacent formations and its relation to hydraulic fracture growth, RockMechanics and Rock Engineering 49 (2016) 689–698.[8] V. Chau, C. Li, S. Rahimi-Aghdam, Z. Bažant, The enigma of large-scale permeability of gas shale: Pre-existing or frac-induced?, Journal ofApplied Mechanics 84 (6) (2017).[9] C. Chang, M. D. Zoback, Viscous creep in room-dried unconsolidated Gulf of Mexico shale (I): Experimental results, Journal of PetroleumScience and Engineering 69 (3–4) (2009) 239–246.[10] Y. Li, A. Ghassemi, Creep behavior of Barnett, Haynesville, and Marcellus shale, in: 46th US Rock Mechanics / Geomechanics Symposium,American Rock Mechanics Association, 2012.[11] H. Sone, M. D. Zoback, Mechanical properties of shale-gas reservoir rocks – Part 2: Ductile creep, brittle strength, and their relation to theelastic modulus, Geophysics 78 (5) (2013) D393–D402.[12] K. C. Bennett, L. A. Berla, W. D. Nix, R. I. Borja, Instrumented nanoindentation and 3D mechanistic modeling of a shale at multiple scales,Acta Geotechnica 10 (1) (2015) 1–14.[13] E. Rybacki, J. Herrmann, R. Wirth, G. Dresen, Creep of Posidonia shale at elevated pressure and temperature, Rock Mechanics and RockEngineering 50 (2017) 3121–3140.[14] Z. Geng, A. Bonnelye, M. Chen, Y. Jin, P. Dick, C. David, X. Fang, A. Schubnel, Time and temperature dependent creep in Tournemire shale,Journal of Geophysical Research: Solid Earth 123 (11) (2018) 9658–9675.[15] F. S. Rassouli, M. D. Zoback, Comparison of short-term and long-term creep experiments in shales and carbonates from unconventional gasreservoirs, Rock Mechanics and Rock Engineering 51 (7) (2018) 1995–2014.[16] M. Trzeciak, H. Sone, M. Dabrowski, Long-term creep tests and viscoelastic constitutive modeling of lower Paleozoic shales from the BalticBasin, N Poland, International Journal of Rock Mechanics and Mining Sciences 112 (2018) 139–157.[17] M. Slim, S. Abedi, L. T. Bryndzia, F.-J. Ulm, Role of organic matter on nanoscale and microscale creep properties of source rocks, Journalof Engineering Mechanics 145 (1) (2019) 04018121.[18] S. Mighani, Y. Bernabé, A. Boulenouar, U. Mok, B. Evans, Creep deformation in Vaca Muerta shale from nanoindentation to triaxial experi-ments, Journal of Geophysical Research: Solid Earth 124 (8) (2019) 7842–7868.[19] J. Herrmann, E. Rybacki, H. Sone, G. Dresen, Deformation experiments on Bowland and Posidonia shale—Part II: Creep behavior at in istu p c – T conditions, Rock Mechanics and Rock Engineering 53 (2019) 755–779.[20] C. Li, J. Wang, H. Xie, Anisotropic creep characteristics and mechanism of shale under elevated deviatoric stress, Journal of PetroleumScience and Engineering 185 (2020) 106670.[21] Y. Mehmani, A. K. Burnham, M. D. V. Berg, F. Gelin, H. Tchelepi, Quantification of kerogen content in organic-rich shales from opticalphotographs, Fuel 177 (2016) 63–75.[22] M. Trzeciak, H. Sone, M. Dabrowski, Long-term creep tests and viscoelastic constitutive modeling of lower Paleozoic shales from the BalticBasin, N Poland, International Journal of Rock Mechanics and Mining Sciences 112 (2018) 139–157.[23] M. B. Dusseault, C. J. Fordham, Time-dependent behavior of rocks, in: Rock Testing and Site Characterization, Elsevier, 1993, pp. 119–149.[24] C. Zhou, J.-H. Yin, J.-G. Zhu, C.-M. Cheng, Elastic anisotropic viscoplastic modeling of the strain-rate-dependent stress–strain behavior of K -consolidated natural marine clays in triaxial shear tests, International Journal of Geomechanics 5 (3) (2005) 218–232.[25] N. Sivasithamparam, M. Karstunen, P. Bonnier, Modelling creep behaviour of anisotropic soft soils, Computers and Geotechnics 69 (2015)46–57.[26] J. Jiang, H. I. Ling, V. N. Kaliakin, X. Zeng, C. Hung, Evaluation of an anisotropic elastoplastic–viscoplastic bounding surface model forclays, Acta Geotechnica 12 (2) (2017) 335–348.[27] M. Leoni, M. Karstunen, P. A. Vermeer, Anisotropic creep model for soft soils, Géotechnique 58 (3) (2008) 215–226.[28] Z. Y. Yin, C. S. Chang, M. Karstunen, P. Y. Hicher, An anisotropic elastic-viscoplastic model for soft clays, International Journal of Solidsand Structures 47 (5) (2010) 665–677.
29] M. R. Karim, C. T. Gnanendran, Review of constitutive models for describing the time dependent behaviour of soft clays, Geomechanics andGeoengineering 9 (1) (2014) 36–51.[30] K. Roscoe, J. Burland, On the generalized stress-strain behaviour of ‘wet’ clay, in: J. Heyman, F. Leckie (Eds.), Engineering Plasticity,Cambridge University Press, Cambridge, 1968, pp. 535–609.[31] S. Pietruszczak, D. Lydzba, J.-F. Shao, Description of creep in inherently anisotropic frictional materials, Journal of Engineering Mechanics130 (6) (2004) 681–690.[32] R. I. Borja, Q. Yin, Y. Zhao, Cam-Clay plasticity. Part IX: On the anisotropy, heterogeneity, and viscoplasticity of shale, Computer Methodsin Applied Mechanics and Engineering 360 (2020) 112695.[33] S. J. Semnani, J. A. White, R. I. Borja, Thermoplasticity and strain localization in transversely isotropic materials based on anisotropic criticalstate plasticity, International Journal for Numerical and Analytical Methods in Geomechanics 40 (2016) 2423–2449.[34] Q. Yu, J. Fish, Multiscale asymptotic homogenization for multiphysics problems with multiple spatial and temporal scales: a coupled thermo-viscoelastic example problem, International journal of solids and structures 39 (26) (2002) 6429–6452.[35] N. Lahellec, P. Suquet, E ff ective behavior of linear viscoelastic composites: A time-integration approach, International Journal of Solids andStructures 44 (2) (2007) 507–529.[36] P. A. Fotiu, S. Nemat-Nasser, Overall properties of elastic-viscoplastic periodic composites, International Journal of Plasticity 12 (2) (1996)163–190.[37] S. Nesenenko, Homogenization in viscoplasticity, SIAM Journal on Mathematical Analysis 39 (1) (2007) 236–262.[38] S. Nesenenko, Homogenization of rate-dependent inelastic models of monotone type, Asymptotic Analysis 81 (1) (2013) 1–29.[39] S. Mercier, A. Molinari, Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori–Tanaka schemes, Inter-national Journal of Plasticity 25 (6) (2009) 1024–1048.[40] P. M. Suquet, Elements of Homogenization for Inelastic Solid Mechanics, in: E. Sanchez-Palencia, A. Zaoui (Eds.), Homogenization tech-niques for composite media, Springer-Verlag, 1987, pp. 193–278.[41] N. Charalambakis, G. Chatzigeorgiou, Y. Chemisky, F. Meraghni, Mathematical homogenization of inelastic dissipative materials: a surveyand recent progress, Continuum Mechanics and Thermodynamics 30 (1) (2018) 1–51.[42] S. Berbenni, L. Capolungo, A Mori-Tanaka homogenization scheme for non-linear elasto-viscoplastic heterogeneous materials based ontranslated fields: An a ffi ne extension, Comptes Rendus Mécanique 343 (2) (2015) 95–106.[43] A. Paquin, H. Sabar, M. Berveiller, Integral formulation and self-consistent modelling of elastoviscoplastic behavior of heterogeneous mate-rials, Archive of Applied Mechanics 69 (1) (1999) 14–35.[44] G. J. Dvorak, Transformation field analysis of inelastic composite materials, Proc. R. Soc. Lond. A 437 (1992) 311–327.[45] J. D. Eshelby, Elastic inclusions and inhomogeneities, in: I. Sneddon, R. Hill (Eds.), Progress in Solid Mechanics, North-Holland, Amsterdam,1961, pp. 89–140.[46] A. A.-C. Guéry, F. Cormery, K. Su, J.-F. Shao, D. Kondo, A micromechanical model for the elasto-viscoplastic and damage behavior of acohesive geomaterial, Physics and Chemistry of the Earth, Parts A / B / C 33 (2008) S416–S421.[47] O. Pierard, I. Doghri, An enhanced a ffi ne formulation and the corresponding numerical algorithms for the mean-field homogenization ofelasto-viscoplastic composites, International Journal of Plasticity 22 (1) (2006) 131–157.[48] I. Doghri, L. Adam, N. Bilger, Mean-field homogenization of elasto-viscoplastic composites based on a general incrementally a ffi ne lin-earization method, International Journal of Plasticity 26 (2) (2010) 219–238.[49] N. Lahellec, P. Suquet, On the e ff ective behavior of nonlinear inelastic composites: I. Incremental variational principles, Journal of theMechanics and Physics of Solids 55 (9) (2007) 1932–1963.[50] L. Brassart, L. Stainier, I. Doghri, L. Delannay, Homogenization of elasto-(visco) plastic composites based on an incremental variationalprinciple, International Journal of Plasticity 36 (2012) 86–112.[51] J. Boudet, F. Auslender, M. Bornert, Y. Lapusta, An incremental variational formulation for the prediction of the e ff ective work-hardeningbehavior and field statistics of elasto-(visco)plastic composites, International Journal of Solids and Structures 83 (2016) 90–113.[52] L. Brassart, Homogenization of elasto-(visco)plastic composites: history-dependent incremental and variational approaches, Ph.D. thesis,Université catholique de Louvaino (2011).[53] A. Molinari, Averaging models for heterogeneous viscoplastic and elastic viscoplastic materials, Journal of Engineering Materials and Tech-nology, Transactions of the ASME 124 (1) (2002) 62–70.[54] S. Mercier, A. Molinari, Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes, Interna-tional Journal of Plasticity 25 (6) (2009) 1024–1048.[55] C. Mareau, V. Favier, M. Berveiller, Micromechanical modeling coupling time-independent and time-dependent behaviors for heterogeneousmaterials, International Journal of Solids and Structures 46 (2) (2009) 223–237.[56] C. Mareau, S. Berbenni, An a ffi ne formulation for the self-consistent modeling of elasto-viscoplastic heterogeneous materials based on thetranslated field method, International Journal of Plasticity 64 (2015) 134–150.[57] H. Sabar, M. Berveiller, V. Favier, S. Berbenni, A new class of micro–macro models for elastic–viscoplastic heterogeneous materials, Inter-national Journal of Solids and Structures 39 (12) (2002) 3257–3276.[58] S. Marfia, E. Sacco, Multiscale technique for nonlinear analysis of elastoplastic and viscoplastic composites, Composites Part B: Engineering136 (2018) 241–253.
59] F. Covezzi, S. de Miranda, S. Marfia, E. Sacco, Homogenization of elastic–viscoplastic composites by the Mixed TFA, Computer Methods inApplied Mechanics and Engineering 318 (2017) 701–723.[60] S. Roussette, J. C. Michel, P. Suquet, Nonuniform transformation field analysis of elastic-viscoplastic composites, Composites Science andTechnology 69 (1) (2009) 22–27.[61] S. Kruch, J.-L. Chaboche, Multi-scale analysis in elasto-viscoplasticity coupled with damage, International Journal of Plasticity 27 (12)(2011) 2026–2039.[62] M. Barral, G. Chatzigeorgiou, F. Meraghni, R. Léon, Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: Theory and numerical validation, International Journal of Plasticity 127 (2020).[63] N. Ohno, X. Wu, T. Matsuda, Homogenized properties of elastic-viscoplastic composites with periodic internal structures, InternationalJournal of Mechanical Sciences 42 (8) (2000) 1519–1536.[64] K. Nakata, T. Matsuda, M. Kawai, Multi-scale creep analysis of plain-woven laminates using time-dependent homogenization theory: E ff ectsof laminate configuration, International Journal of Modern Physics B 22 (31n32) (2008) 6173–6178.[65] N. Ohno, T. Matsuda, X. Wu, A homogenization theory for elastic-viscoplastic composites with point symmetry of internal distributions,International Journal of Solids and Structures 3 38 (2001) 2867–2878.[66] P. W. Chung, K. K. Tamma, R. R. Namburu, A micro / macro homogenization approach for viscoelastic creep analysis with dissipativecorrectors for heterogeneous woven-fabric layered media, Composites science and technology 60 (12-13) (2000) 2233–2253.[67] J. Fish, K. Shek, Computational damage mechanics for composite materials based on mathematical homogenization, International Journalfor Numerical Methods in Engineering 45 (1998) 1657–1679.[68] G. Chatzigeorgiou, N. Charalambakis, Y. Chemisky, F. Meraghni, Periodic homogenization for fully coupled thermomechanical modeling ofdissipative generalized standard materials, International Journal of Plasticity 81 (2016) 18–39.[69] T. Matsuda, N. Ohno, H. Tanaka, T. Shimizu, Homogenized in-plane elastic-viscoplastic behavior of long fiber-reinforced laminates, JSMEInternational Journal 45 (4) (2002) 538–544.[70] A. S. Shamaev, V. V. Shumilova, Homogenization of the equations of state for a heterogeneous layered medium consisting of two creepmaterials, Proceedings of the Steklov Institute of Mathematics 295 (1) (2016) 213–224.[71] S. J. Semnani, J. A. White, An inelastic homogenization framework for layered materials with planes of weakness, Computer Methods inApplied Mechanics and Engineering 370 (2020) 113221.[72] E. Haghighat, F. S. Rassouli, M. D. Zoback, R. Juanes, A viscoplastic model of creep in shale, Geophysics 85 (3) (2020) MR155–MR166.[73] S. Xu, F. S. Rassouli, M. D. Zoback, Utilizing a viscoplastic stress relaxation model to study vertical hydraulic fracture propagation in PermianBasin, in: SPE / AAPG / SEG Unconventional Resources Technology Conference 2017, Unconventional Resources Technology Conference(URTEC), 2017.[74] V. Navarro, E. Alonso, Secondary compression of clays as a local dehydration process, Géotechnique 51 (10) (2001) 859–869.[75] P. Cosenza, D. Korošak, Secondary consolidation of clay as an anomalous di ff usion process, International Journal for Numerical and Analyt-ical Methods in Geomechanics 38 (12) (2014) 1231–1246.[76] J. Choo, J. A. White, R. I. Borja, Hydromechanical modeling of unsaturated flow in double porosity media, International Journal of Geome-chanics 16 (6) (2016) D4016002.[77] R. I. Borja, J. Choo, Cam-Clay plasticity, Part VIII: A constitutive framework for porous materials with evolving internal structure, ComputerMethods in Applied Mechanics and Engineering 309 (2016) 653–679.[78] J. Choo, Y.-H. Jung, C.-K. Chung, E ff ect of directional stress history on anisotropy of initial sti ff ness of cohesive soils measured by benderelement tests, Soils and foundations 51 (4) (2011) 737–747.[79] J. Choo, Y.-H. Jung, W. Cho, C.-K. Chung, E ff ect of pre-shear stress path on nonlinear shear sti ff ness degradation of cohesive soils, Geotech-nical Testing Journal 36 (2) (2013) 198–205.[80] G. Duvaut, J. Lions, Inequalities in Mechanics and Physics, Springer, 1972.[81] J. C. Simo, J. G. Kennedy, S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. loading / unloading conditions and numericalalgorithms, International Journal for Numerical Methods in Engineering 26 (10) (1988) 2161–2185.[82] J. W. Ju, Consistent tangent moduli for a class of viscoplasticity, Journal of Engineering Mechanics 116 (8) (1990) 1764–1779.[83] W. Wang, L. Sluys, R. De Borst, Viscoplasticity for instabilities due to strain softening and strain-rate softening, International Journal forNumerical Methods in Engineering 40 (20) (1997) 3839–3864.[84] M. Lazari, L. Sanavia, B. Schrefler, Local and non-local elasto-viscoplasticity in strain localization analysis of multiphase geomaterials,international Journal for Numerical and Analytical methods in Geomechanics 39 (14) (2015) 1570–1592.[85] P. Perzyna, Fundamental problems in viscoplasticity, in: Advances in applied mechanics, Vol. 9, Elsevier, 1966, pp. 243–377.[86] C. Chang, M. D. Zoback, Viscous creep in room-dried unconsolidated Gulf of Mexico shale (II): Development of a viscoplasticity model,Journal of Petroleum Science and Engineering 72 (1–2) (2010) 50–55.[87] J. A. White, A. K. Burnham, D. W. Camp, A thermoplasticity model for oil shale, Rock Mechanics and Rock Engineering 50 (3) (2017)677–688.[88] Y. Zhao, S. J. Semnani, Q. Yin, R. I. Borja, On the strength of transversely isotropic rocks, International Journal for Numerical and AnalyticalMethods in Geomechanics 42 (16) (2018) 1917–1934.[89] G. E. Backus, Long-wave elastic anisotropy produced by horizontal layering, Journal of Geophysical Research 67 (11) (1962) 4427–4440.
90] R. I. Borja, Plasticity Modeling and Computation, Springer, 2013.[91] H. Sone, Mechanical properties of shale gas reservoir rocks, and its relation to the in-situ stress variation observed in shale gas reservoirs,Ph.D. thesis, Stanford university (2012).[92] H. Sone, M. D. Zoback, Mechanical properties of shale-gas reservoir rocks – Part 1: Static and dynamic elastic properties and anisotropy,Geophysics 78 (5) (2013) D381–D392.[93] W. Bangerth, R. Hartmann, G. Kanschat, deal.II – a general purpose object oriented finite element library, ACM Trans. Math. Softw. 33 (4)(2007) 24 / / deal.II library, version 9.1, Journal of Numerical Mathematics 27 (4) (2019)203–213.[95] J. Choo, R. I. Borja, Stabilized mixed finite elements for deformable porous media with double porosity, Computer Methods in AppliedMechanics and Engineering 293 (2015) 131–154.[96] J. Choo, Stabilized mixed continuous / enriched Galerkin formulations for locally mass conservative poromechanics, Computer Methods inApplied Mechanics and Engineering 357 (2019) 112568.[97] Q. Zhang, J. Choo, R. I. Borja, On the preferential flow patterns induced by transverse isotropy and non-Darcy flow in double porosity media,Computer Methods in Applied Mechanics and Engineering 353 (2019) 570–592.[98] J. T. Camargo, J. A. White, R. I. Borja, A macroelement stabilization for mixed finite element / finite volume discretizations of multiphaseporomechanics, Computational Geosciences (2020) 1–18 doi:10.1007/s10596-020-09964-3 ..