An Inelastic Homogenization Framework for Layered Materials with Planes of Weakness
AAn Inelastic Homogenization Framework for Layered Materials with Planesof Weakness
Shabnam J. Semnani a, ∗ , Joshua A. White b a Department of Structural Engineering, University of California, San Diego, United States b Atmospheric, Earth and Energy Division, Lawrence Livermore National Laboratory, United States
Abstract
Many geologic materials have a composite structure, in which macroscopic mechanical behavior is determinedby the properties, shape, and heterogeneous distribution of individual constituents. In particular, sedimentaryrocks commonly exhibit a layered microstructure, with distinct bedding planes that can also form planes ofweakness. In this work, we present a homogenization framework for modeling inelastic layered media. Theproposed constitutive model allows for distinct micro-constitutive laws for each layer, explicit representationof layer distributions, as well as incorporation of imperfect bonding at the interface between adjacent layers.No a priori assumptions are needed regarding the specific consitutive models used for the layers and interfaces,providing significant modeling flexibility. The overall framework provides a simple and physically-motivatedway of defining anisotropic material behavior as an emergent property of the layered microstructure. Themodel is calibrated using triaxial and true-triaxial experimental data to demonstrate its ability to describeanisotropic deformation and multiple modes of failure.
Keywords:
Transverse isotropy, Asymptotic homogenization, Composite media, Plasticity, Shale
1. Introduction
Many geologic materials are multi-constituent composites. For example, sedimentary rocks commonlyconsist of distinct material layers that have been deposited over geologic time. An interesting question is towhat degree the composite material behavior can be predicted from knowledge of the individual constituents.For example, elastic properties of layered rocks [1] or soils [2] can be readily derived via homogenization of theproperties within a representative volume element. Layer-to-layer interactions induce a transversely isotropicmaterial response even if the separate layers are isotropic. A more challenging question is how to homogenizebehavior in the presence of plasticity and inelastic deformation. To date, most constitutive models for plasticbehavior of transversely isotropic materials take a single-scale viewpoint, with no explicit representation of thematerial microstructure (e.g. [3, 4]). Rather, the e ff ect of layering is introduced through anisotropy parametersor fabric tensors postulated at the macroscopic level.Another common limitation is that many constitutive models ignore debonding and sliding that may occurat layer interfaces. It is clear that these interfaces can significantly a ff ect the mechanical deformation [3, 5, 6].In particular, two main modes of failure for sedimentary rocks include sliding along weak discontinuities andshear localization within matrix layers. To capture these multi-mode failure mechanisms, a variety of failurecriteria for transversely isotropic rocks have been proposed over the years [5, 7–9]. These models, however,typically focus on describing just the brittle failure envelope, and do not predict the complete stress-strainresponse.In this work, we propose a novel strategy for modeling anisotropic, layered media. We rely on an asymp-toptic homogenization framework, in which the governing momentum and traction balance equations for thelayered composite are expanded in a multi-scale representation. This leads to a system of microscale equationsthat determine the homogenized stress, strain, and tangent sti ff ness for the composite. We present an algorithmto solve these microscale equations that has the same structure as an implicit return-mapping algorithm famil-iar from classical plasticity theory. As a result, the proposed model may be readily incorporated in standardfinite element analysis software as a material point subroutine with no special modification. The framework is ∗ Corresponding author
Email addresses: [email protected] (Shabnam J. Semnani), [email protected] (Joshua A. White)
Preprint submitted to arXiv March 26, 2020 a r X i v : . [ phy s i c s . g e o - ph ] M a r odular is the sense that arbitrary constitutive models may be adopted for the individual layers and interfaces.This provides significant flexibility to describe a wide array of layered composite material behavior.The paper is organized as follows. We first provide a brief overview of the homogenization literature inSection 2 to put the current work into context. The multi-scale conceptual model is described in Section 3,followed by a two-scale asymptotic homogenization of the governing equations in Section 4. The key resultis a numerical material point integration algorithm presented in Section 4.5. The resulting model is calibratedusing triaxial and true-triaxial experimental data for several rock types, and numerical results are providedin Section 5 to demonstrate the capability of the present framework. We end the contribution with a fewconcluding remarks.
2. Prior work
Numerous homogenization techniques have been developed in the past few decades to establish a connec-tion between microscopic and macroscopic behavior, determine e ff ective properties, and predict macroscopicresponse of heterogeneous media. Homogenization methods can generally be divided into analytical, semi-analytical and computational techniques, with comprehensive reviews of these methods found in [10–13].The literature on computational homogenizations methods is particularly rich—e.g. FFT-based techniques[14–16], methods based on transform field analysis [17–20], the Ritz-Galerkin approach [21], and the gradient-enhanced scheme [22]. These methods allow for significant modeling flexibility, though with higher imple-mentation complexity and computational cost. Analytical and semi-analytical homogenization techniques, onthe other hand, are appealing in their simplicity, though the underlying conceptual model must be su ffi cientlysimple to lend itself to analytical analysis. In general, (semi-)analytical methods can be further sub-dividedinto two categories. In one category, the material is assumed to be periodic or quasi-periodic, introducing astrong regularity assumption. The method proposed here falls in this category. The second category addressescases where the microstructure is not explicitly regular. In this case, one can look to at least two options. Thefirst includes approximate methods, e.g. the orientational averaging method [23], the self-consistent schemesdeveloped by Hill [24, 25] and Kr¨oner [26] for elastic domains, or e ff ective medium approximations estab-lished by Eshelby [27]. The second option also includes methods based on variational principles to determineupper and lower bounds for e ff ective properties [28–30].For periodic media, asymptotic homogenization methods [31–34] and mean-field methods [35, 36, 36–41](also referred to as average-field methods) are among the most popular techniques. The asymptotic frameworkis particularly appealing for the layered composite considered here. It is based on introducing a fast (or local)coordinate and posing a series expansion of the solution in terms of a scale parameter (cid:15) . The coe ffi cients ofthe homogenized equations can then be written in terms of the solution on a unit cell [42]. We note that highcontrast among the constituent properties may result in a nonlocal structure in the homogenized constitutiverelations. This particular case has been addressed in the literature by a number of researchers [43, 44].Methods based on asymptotic expansion were initially developed for linear systems with regular or nearlyregular structures, and were referred to as mathematical homogenization [45]. Chung et al. [46] developeda generic recipe for the asymptotic homogenization approach and its implementation within the Finite Ele-ment Method (FEM). Pinho-da-Cruz et al. [47] and Oliveria et al. [48] derived the strong form of asymptotichomogenization for linear elasticity, and discussed its numerical implementation using FEM. The asymp-totic homogenization techniques were later extended to nonlinear problems, e.g. composites with nonlinearconstituents [49], nonlinear elastic laminates with imperfect interface contact conditions [50], interaction ofmicrocracks [51], layered thermoelectric composites [52], quasilinear transport equations [53], damage evo-lution [54, 55], and large deformations [56]. Several other researchers have developed advanced variations ofthe asymptotic homogenization technique [57–61]. Among them, Fish and Belsky [57, 62] developed a variantof this technique that removes the assumption of uniform macroscopic fields within the unit cell, which is ad-vantageous for high-gradient regions of the macroscopic field. Ram´ırez-Torres et al. [58] extended asymptotichomogenization to three-scale linear elasticity. Rezakhani and Cusatis [59] combined asymptotic homogeniza-tion with a lattice discrete particle model for quasi-brittle materials to include rotational degrees of freedom.Their method was later applied to modeling of the anisotropic behavior of shale by Li et al. [60]. Higher-orderor strain gradient terms in asymptotic expansion homogenization method for elastic media were accounted forin [61, 63, 64].The majority of the literature on homogenization pertains to elastic materials or plasticity of polycrystals,while problems involving plasticity or dissipation in composites and geomaterials have received less attention.Geological materials often deform beyond the elastic regime, thus limiting the application of multi-scale and2omogenization methods to these materials. The work by Suquet [65] was one of the first to present a ho-mogenization framework for periodic composites with rigid-plastic and elastic-perfectly plastic constituents.In Suquet’s approach, micro-stress and strain fields are written as a sum of an average field and a fluctuatingfield. In addition, macroscopic internal energy and dissipation (or plastic work) are written as an average oftheir microscopic counterparts. Suquet [65] presented the qualitative structure of the macroscopic constitutivelaw, in which macro-strain and the whole field of plastic strains are considered as the state variables. This leadsto an infinite number of internal variables. Therefore, simplifying assumptions were made to derive approxi-mate models such as piecewise constant plastic strains [65, 66]. Subsequently, a number of researchers builton this general framework [66–70]. Among them, Pruchnicki [66] discretized the unit cell into subregionswith constant plastic microstrain tensors, and used a Fourier series approach to solve the integral equation.Construction of macroscopic plastic constitutive equations for arbitrary loading typically requires simpli-fying assumptions regarding stress or strain fields within the individual phases of composites [71]. The Methodof Cells [72] and Generalized Method of Cells [73] are analytical methods based on a first-order representationof displacement field in each subcell, thus piece-wise uniform stress and strain fields in the periodic unit cell.The transformation field analysis approach is another method used in the literature to reduce the number ofmacroscopic internal variables by assuming a piecewise uniform [17] or a non-uniform [19, 20, 74] distribu-tion of microscopic fields of internal variables. Fourier series approximation of stress and strain fields withineach unit cell is another strategy adopted by a number of researchers [75, 76]. Recently, rigorous mathematicalhomogenization of plasticity equations has also been presented [43, 44, 77–83].A number of researchers have extended previously developed homogenization methods to elasto-plasticityproblems, e.g. asymptotic homogenization [45, 46, 68, 84–86], mean-field methods [36, 41, 87], as well asmethods based on variational principles [28–30, 30, 88–92]. Fish et al. [45] generalized the asymptopticmethod by regarding all inelastic strains as eigenstrains in an otherwise elastic domain and approximatingdisplacements and eigenstrains (i.e., plastic strains) by an asymptotic power series expansion. Chatzigeorgiouet al. [93] applied asymptotic expansion homogenization to thermo-mechanical coupling of multi-layereddissipative standard generalized materials. Variational asymptotic homogenization was developed by Yu andTang [94] to combine the advantages of variational methods and asymptotic methods by asymptotic expansionof the energy functional for a unit cell, and this approach was later applied to elasto-plastic problems [95].A number of researchers have adopted a di ff erent view for homogenization involving plasticity, and focusedon deriving the homogenized or e ff ective yield limit [67, 96–103]. Among them, Sawicki [99] formulated thee ff ective yield surface for two-phase layered materials; however, they did not consider its evolution due toplastic strains. Gl¨uge [97] used stress concentration tensors to obtain e ff ective yield surface and e ff ective plas-tic flow for two-phase laminates, and showed that the e ff ective yield surface evolves with plastic deformations,even if each layer is considered isotropic and elasto-perfectly plastic. Gl¨uge [98] examined the application oforientation averaging to derive homogenized plastic properties of laminates. Ponte Casta˜neda and deBotton[101] focused on obtaining bounds and estimates for e ff ective yield strength of rigid-perfectly plastic two-phase composites using a variational approach. Shen et al. [103] derived closed forms of the macroscopicyield criterion for porous materials with Drucker-Prager type plastic matrix and spherical voids.A particular class of periodic materials are layered microstructures, which is the specific focus of thepresent framework. A number of researchers have made di ff erent simplifying assumption to address plasticityin layered materials [104–106]. Among them, El Omri et al. [104] derived a homogenization approach forrigid-plastic and elastic perfectly plastic layered composites, assuming plastic strains are unweighted averagesof the microscopic ones. This assumption, however, holds if the multi-layered medium has homogeneouselastic properties [67, 105]. He and Feng [105] derived exact closed-form solutions for macroscopic behaviorof elastic perfectly plastic periodic composites with perfectly bonded layers. They derived the formulationin the general anisotropic case, and defined the macroscopic elastic energy and plastic dissipation power asvolume averages of their microscopic counterparts.Layered geometry allows for simplifying the general homogenization plasticity framework of Suquet [65].It has been shown that for multi-layered media with homogeneous layers, the microscopic stresses and strainsare constant in each layer, as a direct result of equilibrium equations, compatibility and constitutive relations[66, 67]. The macroscopic quantities can therefore be written as volume averages of the microscopic quantities,in the absence of localization. Following the approach presented by Suquet [65], a number of researchers haveassumed a decomposition of the strain field into a macroscopic average term and a fluctuating term, and takenadvantage of uniformity of stresses and strains in each layer to develop homogenization procedures for elasto-plastic behavior of layered media [67–70]. Among them, Lourenc¸o [70] derived a matrix formulation forhomogenization of elasto-plastic periodic layered media, assuming a homogeneous state of stress and strain ineach layer, piecewise constant inelastic strains, and no relative displacement in the interface between layers.3 x x y yy Ω Y h h y y y n n n n u a ~ u a1,2 u ~ a h N hh k1 ~~ ~~~~~~ h k-1k+1 h u a1,2 ~ u a ~ k-1,k u a ~ k,k+1 n n Γ Y (a) (b) (c) (a) Y W ¶W Y G Y n y y y y n h m Y m Y m+ m , m+ Let u ( x ,y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x ,y ) = r s ˙ u ( x ,y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x ,y ) = C i ( x ,y ) : ˙ " i ( x ,y ) , (7)where C ( x ,y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x ,y ) @ x ⌧ @ C ( x ,y ) @y and @ D ( x ,y ) @ x ⌧ @ D ( x ,y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x ,y,
0) = u ( x ,y ) x ⌦ , y Y (initial displacement) (10e)6 Y Layer m … … Y M x x x y m (b) Figure 1: (a) Schematic representation of unit cell, Y, with multiple parallel homogeneous layers, which compose Ω , the heterogeneousperiodic macroscopic medium. n is the unit vector normal to layers. The macroscopic and microscopic coordinate systems are denotedby x and y , respectively. (b) Schematic of a unit cell with M layers and local scalar coordinate system y . Layer m is shown as Y m , and S m , m + is the interface between layers m and m + They derived an equivalent return mapping algorithm for the average strains and average stresses, in whichthe Jacobian is calculated numerically. Pruchnicki and Shahrour [67] adopted a thermodynamical approach toderive a macroscopic constitutive law for multi-layered media with elastic perfectly plastic constituents. Theydefined the macroscopic free energy and dual potential as volume averages of their microscopic counterparts,and showed that the macroscopic yield surface undergoes kinematic hardening as a result of microscopicplastic deformations. The work of Pruchnicki and Shahrour [67] was later extended by Pruchnicki [68] aswell as Ensan and Shahrour [69] to introduce imperfect interfaces which allow for slippage at the interface oflayers.In the present work, we describe how the asymptotic homogenization technique can be extended to addressinelasticity and imperfect bonding in layered materials, by developing the formulation in rate form and takingadvantage of the layered geometry of the unit cell. We will show that the same form proposed by Suquet [65]for the stress and strain fields, i.e., sum of a fluctuating term and the average term, emerges from the asymptotichomogenization approach. Based on the layered microstructure of the material, we develop a representationand solution strategy for the microscopic unknowns. Moreover, a numerical procedure is presented that canaccommodate general micro-constitutive laws for layers and interfaces.
3. Conceptual model
Consider a macroscopic body Ω composed of a spatially periodic unit cell Y as shown in Figure 1a. Theunit cell consists of parallel layers Y m with index m ∈ L = { , , ..., M } . The material is Y-periodic in thedirection normal to the layers, n , resulting in a unidirectional microstructure. The surface S nm , identified withdouble index nm ∈ I = { , , ..., M } , separates adjacent layers n and m . Note that from periodicity thefinal surface S M wraps around from the last layer to the first, and the double index proves convenient fortracking such details. The displacement field is assumed to be continuous within each layer Y m but may bediscontinuous at each interfaces S nm .Following asymptotic homogenization theory, two separate spatial coordinates are considered: a macro-scopic coordinate, x , and a microscopic coordinate, y . Separation of spatial scales holds, with a non-dimensionalperiod (cid:15) defined as (cid:15) = lL , (cid:15) (cid:28) , (1)where L and l denote the macroscopic and microscopic characteristic length scales, respectively. The twocoordinate systems are related as y = x /(cid:15) . (2)Due to the uni-directional periodicity of the unit cell, it is convenient to define a scalar micro-coordinate axis, y , as (see Figure 1b) y = y · n = (cid:15) ( x · n ) . (3)4 nm S u u [[ u ] ] Figure 2: Decomposition of the total displacement field u into a continuous part ¯ u and a jump [[ u ]] in the vicinity of an interface S nm . The goal of two-scale asymptotic homogenization is to determine the macroscopic behavior of the systemwhen the macroscopic and microscopic scales are su ffi ciently separate; i.e., for (cid:15) →
0. In this case, it isimplicitly assumed in the formulation that the two scales, x and y are independent. As a result, x is treated asa parameter during integration and di ff erentiation with respect to microscopic coordinate y .All material properties are now assumed to be unidirectional, Y-periodic fields that satisfy f ( x , y ) = f ( x , y + kl ) ∀ k ∈ Z , (4)where Z is the space of integers. While here we describe a strictly periodic assumption, in practice it issu ffi cient to have a quasi- periodic medium. Material properties may vary with the macroscopic coordinate,provided that the variations with respect to x are much slower than their variations with respect to the micro-scopic coordinate y . Thus, macroscopic heterogeneities can still be accommodated in the proposed framework.In practice, the material model developed here would be applied at an integration point within a finite elementsimulation, representing material behavior in a local neighborhood of that point. Other integration pointscould be assigned di ff erent properties to capture a slow, macroscopic variation of material microstructure andorientation. Let the domain Ω be delimited by boundary ∂ Ω . The boundary is split as ∂ Ω = ∂ Ω D ∪ ∂ Ω N , with ∂ Ω D ∩ ∂ Ω N = ∅ , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. Atthe microscale, a periodic set of layers and layer surfaces S is included as shown in Figure 1b. We introduce aY-periodic displacement field u ( x , y , t ) with corresponding velocity field ˙ u ( x , y , t ). Because we are interestedin nonlinear material behavior, it is most convenient to work with the rate form of the governing equations toderive the homogenized model. Later, we will switch back to a displacement-based formulation for the actualcomputational implementation.Within the conceptual geometry, the global initial / boundary value problem is to find ˙ u ( x , y , t ) such that ∇ · ˙ σ + ˙ f = in Ω \ S (balance of linear momentum) (5a)[[ ˙ σ ]] · n = on S (traction continuity) (5b)˙ u = ˙ u D on ∂ Ω D (Dirichlet boundary condition) (5c)˙ σ · n N = ˙ t N on ∂ Ω N (Neumann boundary condition) (5d) u ( t ) = u in Ω (initial displacement) . (5e)The discussion here is limited to quasi-static behavior and infinitesimal deformations. For brevity, in theremainder of the paper will drop the argument t when specifying the dependent variables of a field, and focusonly on the spatial dependencies x and y , as they are of primary interest. It should be clear from the discussionwhen a field is time-dependent. We observe, using Equation (3), that the total derivative of u ( x , y ) with respect to x can be computed asd u d x = ∂ u ∂ x + (cid:15) ∂∂ y ( u ⊗ n ) . (6)As a convenient notation, let ∇ x u = ∂ u ∂ x and ∇ y u = ∂∂ y ( u ⊗ n ) , (7)5escribe macro- and micro-scale gradient operators, respectively. In our conceptual model, the displacementfield within a given layer Y m is continuous, but it is potentially discontinuous at layer interfaces. Consideringa small neighborhood in the vicinity of a surface S nm as in Figure 2, we may express the total displacementfield u in the particular form [107, 108] u = ¯ u + H nm [[ u ]] nm . (8)Here, ¯ u ( x , y ) is a continuous field, H nm ( y ) is a Heaviside function centered at the interface under consideration,and [[ u ]] nm is the displacement jump, [[ u ]] nm = u ( x , y + ) − u ( x , y − ) , (9)where y + and y − denote the micro-coordinates at each side of the discontinuity. We emphasize that the Heav-iside function only depends on the micro-coordinate y , as the jump location is determined by the unit cellgeometry.Noting that the strain field ˆ ε is the symmetric gradient of this displacement field, we may writeˆ ε = ε x + (cid:15) ε y + δ nm (cid:15) ε nm , (10)with δ nm the Dirac delta function centered at the interface, and strain components ε x = ∇ sx ¯ u , ε y = ∂∂ y ( ¯ u ⊗ s n ) , and ε nm = [[ u ]] nm ⊗ s n . (11)Here, ∇ s and ⊗ s denote symmetric gradient and symmetric dyadic product, respectively. Due to the appearanceof the Dirac delta, the last term in Equation (10) is singular, but it is only present at the interface itself. Awayfrom an interface, the strain field is regular and bounded, and can be written as ε = ε x + (cid:15) ε y . (12)Each layer Y m in the unit cell is assumed to be homogeneous. The extension of the model to non-homogeneous layers will be discussed in Section 4.7. We assume the constitutive behavior of individuallayers can be expressed in a generic rate form as˙ σ m = C m : ˙ ε m ∀ m ∈ L , (13)where ˙ σ m denotes the stress rate, ˙ ε m the strain rate, and C m the sti ff ness tangent operator within layer m . Thisgeneric form can represent a wide variety of elasto-plastic-damage models for the constituent layers.The traction rate at an interface is ˙ t = ˙ σ · n . While the strain rate at an interface may be singular, thetraction rate must be bounded in order for the traction equilibrium (5b) to be satisfied. To be consistent withEquations (10)–(13), we assume that the interface constitutive model can be expressed in the form,˙ t nm = (cid:15) D nm · [[ ˙ u ]] nm ∀ nm ∈ I , (14)where D nm is the sti ff ness tangent operator of the interface. We will see that the presence of the 1 /(cid:15) weightinghere is crucial within the asymptoptic analysis to preserve a proper balancing of terms at various scales inthe expansion. A scale-weighted form is consistent with traction expressions proposed by other authors in[68, 109–111]. This generic form also allows for a variety of specific interface models in order to capturedi ff erent types of contact behavior. For example, the model in [112] describes a complete framework fordescribing friction, dilation, and damage on rough interfaces.Equations (13) and (14) let us handle the continuous and discontinuous portions of the displacement fieldseparately, and alleviate the need to directly work with a singular strain field. Therefore, in the remainderof this paper, we will only use the regular part of the strain rate tensor from Equation (10). It is noteworthythat the homogenized model will now be developed based on the generic micro-constitutive behaviors givenin Equations (13) and (14) without any further specification. We will illustrate the incorporation of specificmaterial models through numerical examples in Section 5, allowing for a diverse spectrum of macroscopicbehavior.We again remark that in practical applications all fields can be modeled as quasi-periodic; that is, thesti ff ness tangents C and D are functions of both macroscopic and microscopic scales. Their variations at themacroscopic scale are, however, su ffi ciently slow so that separation is retained. In this case we have ∂ C ( x , y ) ∂ x (cid:28) ∂ C ( x , y ) ∂ y and ∂ D ( x , y ) ∂ x (cid:28) ∂ D ( x , y ) ∂ y , (15)6llowing for homogenization on the unit cell. In this case, all response functions, including displacements,stresses, and strains, are also quasi-periodic.Before concluding this section, we remark that the divergence of a continuous second order tensor σ canbe similiarly expanded as ∇ · σ = ∇ x · σ + (cid:15) ∇ y · σ , (16)with ( ∇ x · σ ) i j = ∂σ i j ∂ x j and (cid:16) ∇ y · σ (cid:17) i j = ∂σ i j ∂ y n j . (17)
4. Two-scale asymptotic homogenization
Our goal now is to transform the governing formulation into a two-scale system of equilibrium equations,from which the homogenized material behavior may be derived.
Following asymptotic homogenization theory, the unknown velocity field can be represented by an expan-sion in a power series of (cid:15) as ˙ u = ∞ (cid:88) i = (cid:15) i ˙ u ( i ) = ˙ u (0) + (cid:15) ˙ u (1) + (cid:15) ˙ u (2) + O (cid:16) (cid:15) (cid:17) . (18)Similarly, the regular portion of the strain rate is expanded as˙ ε = ∞ (cid:88) i = (cid:16) (cid:15) i ˙ ε ( i ) x + (cid:15) i − ˙ ε ( i ) y (cid:17) = ∞ (cid:88) i = − (cid:15) i ˙ ε ( i ) , (19)in which we have introduced the notation˙ ε ( − = ˙ ε (0) y and ˙ ε ( i ) = ˙ ε ( i ) x + ˙ ε ( i + y for i = { , , , ... } . (20)The stress rate then follows as˙ σ = C : ∞ (cid:88) i = − (cid:15) i ˙ ε ( i ) = ∞ (cid:88) i = − (cid:15) i σ ( i ) with ˙ σ ( i ) = C : ˙ ε ( i ) . (21)Similarly, the traction rate is ˙ t = ∞ (cid:88) i = − (cid:15) i ˙ t ( i ) with ˙ t ( i ) = D · [[ ˙ u ( i + ]] . (22)Comparing the stress and traction expansions, we may infer that the relationship ˙ t ( i ) = ˙ σ ( i ) · n also holds. Using the macroscopic spatial divergence operator defined in Equation (16), the balance of linear momen-tum given in (5a) can be expanded as ∇ x · ˙ σ + (cid:15) ∇ y · ˙ σ + ˙ f = . (23)Replacing the expanded expression for the stress rate tensor given in Equation (21) into Equation (23) leads to (cid:15) − ∇ y · ˙ σ ( − + (cid:15) − (cid:16) ∇ x · ˙ σ ( − + ∇ y · ˙ σ (0) (cid:17) + (cid:15) (cid:16) ∇ x · ˙ σ (0) + ∇ y · ˙ σ (1) + ˙ f (cid:17) + O ( (cid:15) ) = . (24)Each term multiplied by a power of (cid:15) in Equations (24) is now separately set to zero to ensure that the asymp-totic series approximation is valid as (cid:15) →
0. The results obtained from setting terms of order (cid:15) − , (cid:15) − , and (cid:15) to zero imply ∇ y · ˙ σ ( − = , (25a) ∇ x · ˙ σ ( − + ∇ y · ˙ σ (0) = , (25b) ∇ x · ˙ σ (0) + ∇ y · ˙ σ (1) + ˙ f = . (25c)7n addition, traction continuity as expressed in Equation (5b) needs to be satisfied at all scales. Therefore, wehave [[ ˙ σ ]] · n = ∞ (cid:88) i = − (cid:15) i [[ ˙ σ ( i ) ]] · n = , (26)from which we infer, [[ ˙ σ ( − ]] · n = , (27a)[[ ˙ σ (0) ]] · n = . (27b)We now examine each of the equilibrium equations in turn. Consider the unit cell Y as shown in Figure 1b. We now multiply both sides of (25a) by ˙ u (0) , integrate overthe unit cell, and use the definition of the microscale divergence operator in equation (17), as follows0 = (cid:88) m ∈ L (cid:90) Y m ˙ u (0) · ∂ ˙ σ ( − ∂ y · n d y = (cid:88) m ∈ L (cid:104) ˙ u (0) · ˙ σ ( − · n (cid:105) Y + m Y − m − (cid:88) m ∈ L (cid:90) Y m ∂∂ y (cid:16) ˙ u (0) ⊗ n (cid:17) : ˙ σ ( − d y = (cid:88) nm ∈ I [[ ˙ u (0) ]] ⊗ ˙ t ( − − (cid:88) m ∈ L (cid:90) Y m ∂∂ y (cid:16) ˙ u (0) ⊗ n (cid:17) : ˙ σ ( − d y = (cid:88) nm ∈ I [[ ˙ u (0) ]] · D · [[ ˙ u (0) ]] − (cid:88) m ∈ L (cid:90) Y m ∂∂ y (cid:16) ˙ u (0) ⊗ s n (cid:17) : C : ∂∂ y (cid:16) ˙ u (0) ⊗ s n (cid:17) d y . (28)Both terms must go to zero to satisfy this relationship. If C is positive definite, we may infer ∂ s ˙ u (0) ∂ y = ⇒ ˙ u (0) = ˙ u (0) ( x ) . (29)That is, ˙ u (0) is constant in y and only depends on x . This leads to ˙ ε (0) y = and consequently, ˙ σ ( − = .Similarly, if D is positive definite, we have [[ ˙ u (0) ]] = . Therefore, the medium is continuous at the macroscale(as desired) and we only need to consider here velocity jumps on the microscale surfaces in S .It should be noted that the assumptions regarding positive definiteness of C and D hold in the absenceof localization, when separation of scales holds. Slip may still occur at the microscopic interfaces, but ina pervasive and periodic manner. Of course, a softening mode is likely to be unstable and may eventuallylocalize to one surface in a non-periodic manner. Such a surface would then need to be explicitly includedwithin the macroscopic boundary-value-problem ([[ ˙ u (0) ]] (cid:44) ) since a separation-of-scales violation has takenplace. Applying the result of equilibrium equation obtained in the previous section to Equation (25b) and usingthe divergence operator defined for scalar axis y in Equation (16), we can write ∇ y · ˙ σ (0) = ∂∂ y (cid:16) ˙ σ (0) · n (cid:17) = . (30)Therefore, ˙ σ (0) · n = ˙ t (0) ( x ) is constant in unit cell Y . As a result, the traction continuity condition in Equa-tion (27b) is trivially satisfied. Consequently, we obtain the following microscale balance equation for eachsubdomain Y m (cid:16) C m : ˙ ε (0) (cid:17) · n = ˙ t (0) ( x ) ∀ m ∈ L . (31)The equilibrium at interface S nm is similarly obtained from Equation (22) as D nm · [[ ˙ u (1) ]] nm = ˙ t (0) ( x ) ∀ nm ∈ I . (32)In Section 4.4 we will use the solution of the microscale balance equations (31) and (32) to determine ˙ σ (0) .8 .2.3. Equilibrium equation (25c) Let us define the averaging operator over the unit cell, (cid:104) • (cid:105) = | Y | (cid:90) Y ( • ) d y , (33)in which | Y | is volume of the cell. Applying the averaging operator to both sides of Equation (25c) gives1 | Y | (cid:34)(cid:90) Y ∇ x · ˙ σ (0) d y + (cid:90) Y ∇ y · ˙ σ (1) d y + (cid:90) Y ˙ f d y (cid:35) = . (34)Here, x is treated as a parameter during integration with respect to microscopic coordinate y . Applying thedivergence theorem to the second term on the left hand side of Equation (34) leads to:1 | Y | (cid:34) ∇ x · (cid:90) Y ˙ σ (0) d y + (cid:90) Γ Y ˙ σ (1) · n d Γ + (cid:90) Y ˙ f d y (cid:35) = . (35)Due to periodicity of the unit cell, the second integral on the left hand side of (35) vanishes, and we obtain themacroscopic balance equation ∇ x · ˙ Σ + ˙ F = in Ω . (36)Here, the macroscopic stress and body force rates are˙ Σ ( x ) = (cid:68) ˙ σ (0) (cid:69) and ˙ F = (cid:68) ˙ f (cid:69) . (37) The previous developments have employed the rate form of the governing equations, which simplifiesthe presentation in the nonlinear material model context. In practice, however, we often prefer a non-rate,displacement-based form for implementation purposes. The two are essentially equivalent, so we may readilyswitch back and forth.As mentioned earlier, we are interested in the solution to the problem as (cid:15) →
0. A two-scale expansion issu ffi cient to provide an approximate solution to the initial / boundary value problem, with an error of O ( (cid:15) ). Letus therefore write the displacement field as, u ( x , y ) ≈ u (0) ( x ) + (cid:15) u (1) ( x , y ) . (38)The field u (0) provides the macroscopic displacement field in Ω , while u (1) provides a microscale “fluctuation”that is periodic in Y . The strain field ε (0) is decomposed as ε (0) = E + e , (39)with macroscale and microscale strain tensors identified, respectively, as E = ∇ sx u (0) and e = ∇ sy u (1) . (40)The two-scale initial / boundary value problem is then: • Macroscale Problem : Find u (0) ( x ) such that the following macroscopic conditions are satisfied: ∇ x · Σ + F = in Ω (balance of linear momentum) (41a) u (0) = u D on ∂ Ω D (Dirichlet B.C.) (41b) Σ · n N = t N on ∂ Ω N (Neumann B.C.) (41c) u (0) ( t ) = u in Ω (initial displacement) (41d)where the macroscopic stress Σ = (cid:68) σ (0) (cid:69) at point x satisfies the microscale problem: • Microscale Problem : Given u (0) ( x ), find u (1) ( x , y ) such that σ (0) · n = t (0) ( x ) in Y m for each layer m ∈ L (42a) t (0) nm = t (0) ( x ) on S nm for each interface nm ∈ I (42b)Periodicity conditions requires that: u (1) is periodic in Y (42c) σ (0) · n is anti-periodic on Γ Y (42d)9 … … y n M-1 M M Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 u u (1) E . n h A H y y A B B A Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 A B Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 u (1) Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 Let u ( x ,y )representthedisplacementfieldintheheterogeneousbody⌦.Theinfinitesimalstrainrate tensorisdefinedas˙ " ( x ,y )= r s ˙ u ( x ,y ) , (5) whereoperator r s denotessymmetricgradient,asdefinedinAppendixA,and˙ ⇤ denotesderivativewith respecttotime. Atthesurface S thefollowingvelocityjumpvectorisdefined: [[˙ u ]]=˙ u + ˙ u (6) EachlayerintheunitcellYisassumedtobehomogeneous.Theextensionofthemodeltonon- homogeneouslayersisdiscussedinSection4.9.Theconstitutivebehavioroflayershasthegeneralform of˙ i ( x ,y )= C i ( x ,y ):˙ " i ( x ,y ) , (7) where C ( x ,y )denotesthesti↵nesstangentoperatorofthelayer.Atsurface S betweentwoadjacent layers,velocityjumpvector[[˙ u ]]isrelatedtorateoftractionvectoractingonthesurface,˙ t =˙ · n ,as ˙ t = D · [[˙ u ]] , (8) where D isthesti↵nesstangentoperatoroftheinterface S .Variousconstitutivemodelscanbedefined todescribedi↵erenttypesofcontactandinterfacebehavior.Forexample,themodelprovidedin[108] accountsforfrictionandevolvingdilationatthejointinterface.Itisnoteworthythattheformulationis developedherebasedonthegeneralforms(7)and(8)forthemicro-constitutivebehaviors,withoutany furtherassumptions.Subsequently,wewilldemonstrateapplicationofspecificmaterialmodelsinthe numericalexamplesinSection5. Alsonotethatthepresentsystemisquasi-periodic;thatis,thesti↵nesstangents C and D arefunctions ofboththemacroscopicandmicroscopicscales.Thevariationsatthemacroscopicscaleare,however, su cientlyslowsothattheseparationofscalesisretained.Inthiscasewehave @ C ( x ,y ) @ x ⌧ @ C ( x ,y ) @y and @ D ( x ,y ) @ x ⌧ @ D ( x ,y ) @y, (9) thushomogenizationcanbeperformedonaunitcell.Itisalsoacceptabletoassumethatallresponse functions,includingdisplacements,stresses,andstrains,arealsoquasi-periodic.Inotherwords,their localvariationswithrespecttomicro-coordinate y isperiodic,whiletheiraveragevaluescanvarywith respecttomacro-coordinate x . Letboundary @ ⌦bepartitionedas @ ⌦= @ ⌦ D [ @ ⌦ N ,with @ ⌦ D \ @ ⌦ N = ; ,wheresubscripts D and N denoteDirichletandNeumannconditions,respectively.Theinitial/boundaryvalueproblemcanbe statedasfollows: r· ˙ +˙ f = in⌦(balanceoflinearmomentum)(10a) ˙ u =˙¯ u ( x )on @ ⌦ D (Dirichletboundaryconditions)(10b) ˙ · ¯ n =˙¯ t ( x )on @ ⌦ N (Neumannboundarycondition)(10c) [[˙ · n ]]= on S (tractioncontinuity)(10d) u ( x ,y, u ( x ,y ) x ⌦ ,y Y (initialdisplacement)(10e) 6 AB u (1) Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 Let u ( x ,y )representthedisplacementfieldintheheterogeneousbody⌦.Theinfinitesimalstrainrate tensorisdefinedas˙ " ( x ,y )= r s ˙ u ( x ,y ) , (5) whereoperator r s denotessymmetricgradient,asdefinedinAppendixA,and˙ ⇤ denotesderivativewith respecttotime. Atthesurface S thefollowingvelocityjumpvectorisdefined: [[˙ u ]]=˙ u + ˙ u (6) EachlayerintheunitcellYisassumedtobehomogeneous.Theextensionofthemodeltonon- homogeneouslayersisdiscussedinSection4.9.Theconstitutivebehavioroflayershasthegeneralform of˙ i ( x ,y )= C i ( x ,y ):˙ " i ( x ,y ) , (7) where C ( x ,y )denotesthesti↵nesstangentoperatorofthelayer.Atsurface S betweentwoadjacent layers,velocityjumpvector[[˙ u ]]isrelatedtorateoftractionvectoractingonthesurface,˙ t =˙ · n ,as ˙ t = D · [[˙ u ]] , (8) where D isthesti↵nesstangentoperatoroftheinterface S .Variousconstitutivemodelscanbedefined todescribedi↵erenttypesofcontactandinterfacebehavior.Forexample,themodelprovidedin[108] accountsforfrictionandevolvingdilationatthejointinterface.Itisnoteworthythattheformulationis developedherebasedonthegeneralforms(7)and(8)forthemicro-constitutivebehaviors,withoutany furtherassumptions.Subsequently,wewilldemonstrateapplicationofspecificmaterialmodelsinthe numericalexamplesinSection5. Alsonotethatthepresentsystemisquasi-periodic;thatis,thesti↵nesstangents C and D arefunctions ofboththemacroscopicandmicroscopicscales.Thevariationsatthemacroscopicscaleare,however, su cientlyslowsothattheseparationofscalesisretained.Inthiscasewehave @ C ( x ,y ) @ x ⌧ @ C ( x ,y ) @y and @ D ( x ,y ) @ x ⌧ @ D ( x ,y ) @y, (9) thushomogenizationcanbeperformedonaunitcell.Itisalsoacceptabletoassumethatallresponse functions,includingdisplacements,stresses,andstrains,arealsoquasi-periodic.Inotherwords,their localvariationswithrespecttomicro-coordinate y isperiodic,whiletheiraveragevaluescanvarywith respecttomacro-coordinate x . Letboundary @ ⌦bepartitionedas @ ⌦= @ ⌦ D [ @ ⌦ N ,with @ ⌦ D \ @ ⌦ N = ; ,wheresubscripts D and N denoteDirichletandNeumannconditions,respectively.Theinitial/boundaryvalueproblemcanbe statedasfollows: r· ˙ +˙ f = in⌦(balanceoflinearmomentum)(10a) ˙ u =˙¯ u ( x )on @ ⌦ D (Dirichletboundaryconditions)(10b) ˙ · ¯ n =˙¯ t ( x )on @ ⌦ N (Neumannboundarycondition)(10c) [[˙ · n ]]= on S (tractioncontinuity)(10d) u ( x ,y, u ( x ,y ) x ⌦ ,y Y (initialdisplacement)(10e) 6 BA y h H y y u u (1) ~ ~ ~ ~ h h u (1) Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 Let u ( x ,y )representthedisplacementfieldintheheterogeneousbody⌦.Theinfinitesimalstrainrate tensorisdefinedas˙ " ( x ,y )= r s ˙ u ( x ,y ) , (5) whereoperator r s denotessymmetricgradient,asdefinedinAppendixA,and˙ ⇤ denotesderivativewith respecttotime. Atthesurface S thefollowingvelocityjumpvectorisdefined: [[˙ u ]]=˙ u + ˙ u (6) EachlayerintheunitcellYisassumedtobehomogeneous.Theextensionofthemodeltonon- homogeneouslayersisdiscussedinSection4.9.Theconstitutivebehavioroflayershasthegeneralform of˙ i ( x ,y )= C i ( x ,y ):˙ " i ( x ,y ) , (7) where C ( x ,y )denotesthesti↵nesstangentoperatorofthelayer.Atsurface S betweentwoadjacent layers,velocityjumpvector[[˙ u ]]isrelatedtorateoftractionvectoractingonthesurface,˙ t =˙ · n ,as ˙ t = D · [[˙ u ]] , (8) where D isthesti↵nesstangentoperatoroftheinterface S .Variousconstitutivemodelscanbedefined todescribedi↵erenttypesofcontactandinterfacebehavior.Forexample,themodelprovidedin[108] accountsforfrictionandevolvingdilationatthejointinterface.Itisnoteworthythattheformulationis developedherebasedonthegeneralforms(7)and(8)forthemicro-constitutivebehaviors,withoutany furtherassumptions.Subsequently,wewilldemonstrateapplicationofspecificmaterialmodelsinthe numericalexamplesinSection5. Alsonotethatthepresentsystemisquasi-periodic;thatis,thesti↵nesstangents C and D arefunctions ofboththemacroscopicandmicroscopicscales.Thevariationsatthemacroscopicscaleare,however, su cientlyslowsothattheseparationofscalesisretained.Inthiscasewehave @ C ( x ,y ) @ x ⌧ @ C ( x ,y ) @y and @ D ( x ,y ) @ x ⌧ @ D ( x ,y ) @y, (9) thushomogenizationcanbeperformedonaunitcell.Itisalsoacceptabletoassumethatallresponse functions,includingdisplacements,stresses,andstrains,arealsoquasi-periodic.Inotherwords,their localvariationswithrespecttomicro-coordinate y isperiodic,whiletheiraveragevaluescanvarywith respecttomacro-coordinate x . Letboundary @ ⌦bepartitionedas @ ⌦= @ ⌦ D [ @ ⌦ N ,with @ ⌦ D \ @ ⌦ N = ; ,wheresubscripts D and N denoteDirichletandNeumannconditions,respectively.Theinitial/boundaryvalueproblemcanbe statedasfollows: r· ˙ +˙ f = in⌦(balanceoflinearmomentum)(10a) ˙ u =˙¯ u ( x )on @ ⌦ D (Dirichletboundaryconditions)(10b) ˙ · ¯ n =˙¯ t ( x )on @ ⌦ N (Neumannboundarycondition)(10c) [[˙ · n ]]= on S (tractioncontinuity)(10d) u ( x ,y, u ( x ,y ) x ⌦ ,y Y (initialdisplacement)(10e) 6 u (1) Let u ( x , y ) represent the displacement field in the heterogeneous body ⌦. The infinitesimal strain ratetensor is defined as ˙ " ( x , y ) = r s ˙ u ( x , y ) , (5)where operator r s denotes symmetric gradient, as defined in Appendix A, and ˙ ⇤ denotes derivative withrespect to time.At the surface S the following velocity jump vector is defined:[[ ˙ u ]] = ˙ u + ˙ u (6) Each layer in the unit cell Y is assumed to be homogeneous. The extension of the model to non-homogeneous layers is discussed in Section 4.9. The constitutive behavior of layers has the general formof ˙ i ( x , y ) = C i ( x , y ) : ˙ " i ( x , y ) , (7)where C ( x , y ) denotes the sti↵ness tangent operator of the layer. At surface S between two adjacentlayers, velocity jump vector [[ ˙ u ]] is related to rate of traction vector acting on the surface, ˙ t = ˙ · n , as˙ t = D · [[ ˙ u ]] , (8)where D is the sti↵ness tangent operator of the interface S . Various constitutive models can be definedto describe di↵erent types of contact and interface behavior. For example, the model provided in [108]accounts for friction and evolving dilation at the joint interface. It is noteworthy that the formulation isdeveloped here based on the general forms (7) and (8) for the micro-constitutive behaviors, without anyfurther assumptions. Subsequently, we will demonstrate application of specific material models in thenumerical examples in Section 5.Also note that the present system is quasi-periodic; that is, the sti↵ness tangents C and D are functionsof both the macroscopic and microscopic scales. The variations at the macroscopic scale are, however,su ciently slow so that the separation of scales is retained. In this case we have @ C ( x , y ) @ x ⌧ @ C ( x , y ) @y and @ D ( x , y ) @ x ⌧ @ D ( x , y ) @y , (9)thus homogenization can be performed on a unit cell. It is also acceptable to assume that all responsefunctions, including displacements, stresses, and strains, are also quasi-periodic. In other words, theirlocal variations with respect to micro-coordinate y is periodic, while their average values can vary withrespect to macro-coordinate x . Let boundary @ ⌦ be partitioned as @ ⌦ = @ ⌦ D [ @ ⌦ N , with @ ⌦ D \ @ ⌦ N = ; , where subscripts D and N denote Dirichlet and Neumann conditions, respectively. The initial/boundary value problem can bestated as follows: r · ˙ + ˙ f = in ⌦ (balance of linear momentum) (10a)˙ u = ˙¯ u ( x ) on @ ⌦ D (Dirichlet boundary conditions) (10b)˙ · ¯ n = ˙¯ t ( x ) on @ ⌦ N (Neumann boundary condition) (10c)[[ ˙ · n ]] = on S (traction continuity) (10d) u ( x , y,
0) = u ( x , y ) x ⌦ , y Y (initial displacement) (10e)6 Let u ( x ,y )representthedisplacementfieldintheheterogeneousbody⌦.Theinfinitesimalstrainrate tensorisdefinedas˙ " ( x ,y )= r s ˙ u ( x ,y ) , (5) whereoperator r s denotessymmetricgradient,asdefinedinAppendixA,and˙ ⇤ denotesderivativewith respecttotime. Atthesurface S thefollowingvelocityjumpvectorisdefined: [[˙ u ]]=˙ u + ˙ u (6) EachlayerintheunitcellYisassumedtobehomogeneous.Theextensionofthemodeltonon- homogeneouslayersisdiscussedinSection4.9.Theconstitutivebehavioroflayershasthegeneralform of˙ i ( x ,y )= C i ( x ,y ):˙ " i ( x ,y ) , (7) where C ( x ,y )denotesthesti↵nesstangentoperatorofthelayer.Atsurface S betweentwoadjacent layers,velocityjumpvector[[˙ u ]]isrelatedtorateoftractionvectoractingonthesurface,˙ t =˙ · n ,as ˙ t = D · [[˙ u ]] , (8) where D isthesti↵nesstangentoperatoroftheinterface S .Variousconstitutivemodelscanbedefined todescribedi↵erenttypesofcontactandinterfacebehavior.Forexample,themodelprovidedin[108] accountsforfrictionandevolvingdilationatthejointinterface.Itisnoteworthythattheformulationis developedherebasedonthegeneralforms(7)and(8)forthemicro-constitutivebehaviors,withoutany furtherassumptions.Subsequently,wewilldemonstrateapplicationofspecificmaterialmodelsinthe numericalexamplesinSection5. Alsonotethatthepresentsystemisquasi-periodic;thatis,thesti↵nesstangents C and D arefunctions ofboththemacroscopicandmicroscopicscales.Thevariationsatthemacroscopicscaleare,however, su cientlyslowsothattheseparationofscalesisretained.Inthiscasewehave @ C ( x ,y ) @ x ⌧ @ C ( x ,y ) @y and @ D ( x ,y ) @ x ⌧ @ D ( x ,y ) @y, (9) thushomogenizationcanbeperformedonaunitcell.Itisalsoacceptabletoassumethatallresponse functions,includingdisplacements,stresses,andstrains,arealsoquasi-periodic.Inotherwords,their localvariationswithrespecttomicro-coordinate y isperiodic,whiletheiraveragevaluescanvarywith respecttomacro-coordinate x . Letboundary @ ⌦bepartitionedas @ ⌦= @ ⌦ D [ @ ⌦ N ,with @ ⌦ D \ @ ⌦ N = ; ,wheresubscripts D and N denoteDirichletandNeumannconditions,respectively.Theinitial/boundaryvalueproblemcanbe statedasfollows: r· ˙ +˙ f = in⌦(balanceoflinearmomentum)(10a) ˙ u =˙¯ u ( x )on @ ⌦ D (Dirichletboundaryconditions)(10b) ˙ · ¯ n =˙¯ t ( x )on @ ⌦ N (Neumannboundarycondition)(10c) [[˙ · n ]]= on S (tractioncontinuity)(10d) u ( x ,y, u ( x ,y ) x ⌦ ,y Y (initialdisplacement)(10e) 6 M1 (a) (b) Figure 3: Schematic of the total displacement field u and the periodic fluctuation u (1) in (a) bi-layer and (b) multi-layer unit cell. The macroscale problem is identical to a typical solid mechanics formulation solved via, e.g., a finite elementdiscretization. The constitutive relationship between the macroscopic stress Σ and macroscopic strain E atan integration point, however, is computed by solving the microscale problem. Given a macroscopic field u (0) with associated strain E , we solve the micro-equilibrium equations to determine σ (0) , from which themacroscopic stress Σ may be computed by averaging over Y . A procedure for doing so is detailed in the nextsection.We remark that for a unit cell located at the macroscopic point x , Equation (39) can be used to obtain theapproximate displacement vector u at the microscale. For this purpose, Equation (39) is integrated in terms of y , which leads to: u ≈ E ( x ) · y + u (1) . (43)Equations (39-43) link the asymptotic homogenization approach to a group of other methods used in theliterature, for example [65, 105], in which the microstrain field is written as sum of an average value and aperiodic fluctuation. Within the material point update, the macroscale strain E is viewed as a known input, from which we seekto derive the macroscale stress Σ and the consistent tangent operator C . The general formulation is derived fora multi-layer unit cell with M layers, as shown in Figure 3. To clarify the formulation further, specific formsare also presented for the simple subcase of a bi-layer unit cell consisting of two parallel layers, A and B , andtwo interfaces, S AB and S BA .The volumetric ratio for each phase, m , is defined as φ m = h m l with (cid:88) m ∈ L φ m = , (44)where l is the total thickness of the unit cell. If layers are individually homogeneous, microscopic strains andstresses are uniform in each layer, as a direct result of equilibrium equations, compatibility, and constitutiverelations [66, 67]. We will indicate by ε m and σ m the restriction of the strain field ε (0) and stress field σ (0) toeach layer m . 10niformity of the strain tensor e in each layer directly indicates that the microscopic periodic displacement u (1) is a piecewise linear function of y . To satisfy the periodic conditions given in (42c), micro-displacementsmust satisfy the following compatibility condition: (cid:88) m ∈ L (cid:90) Y m ∂ u (1) ∂ y d y + (cid:88) nm ∈ I [[ u (1) ]] nm = . (45)See Figure 3 for a schematic visualization. It is logical then to define the following primary unknowns: v m = ∂ u (1) m ∂ y (constant for each layer m ∈ L ) , w nm = [[ u (1) ]] nm (constant for each interface nm ∈ I ) . (46)Within the solution algorithm, we will solve for one displacement gradient vector v m for each layer, and onedisplacement jump vector w nm for each interface. Note, however, that symmetry requires w nm = w mn if twomaterials come into contact multiple times due to periodicity, reducing the number of independent unknowns.This is easiest to see in the bi-layer case in Figure 3a. Surfaces S AB and S BA both experience the same loadingconditions, and therefore w AB = w BA .Replacing (46) into (45) gives the simplified compatibility condition (cid:88) m ∈ L φ m v m + (cid:88) nm ∈ I w nm = . (47)Using these unknowns, the total strain in layer m is ε m = E + v m ⊗ s n , (48)from which the layer stress σ m ( ε m ) can be computed by calling a material model update routine specific to thatlayer. Similarly, the traction t nm ( w nm ) can be computed by calling the desired traction model. The micro-scalebalance equations to be satisfied are σ m · n − t (0) ( x ) = for all layers m ∈ L , (49) t nm − t (0) ( x ) = for all interfaces nm ∈ I . (50)The normal traction t (0) ( x ), which is constant in Y , is also unknown and must be determined by the materialpoint algorithm. For simplicity, in this section we set t = t (0) ( x ). In total, we have M layer balance equations(49), N surface balance equations (50), and one compatibility condition (47). This provides a fully determinedsystem for the M displacement gradients v m , the N displacement jumps w nm , and the overall traction t .In the following, we now switch to Voigt representation of the tensor quantities to clarify how the algorithmis implemented. We use brackets {•} to indicate a symmetric tensor that has been unrolled into 6 × { ε m } = { E } + Bv m , (51) { ε m } = (cid:104) ε ε ε ε ε ε (cid:105) T , (52) v m = (cid:104) v v v (cid:105) T . (53)Here, B is a 6 × n , B = n n
00 0 n n n n n n n . (54)Its transpose can be used to express balance equation (49) as B T { σ m } − t = , (55)11 σ m } = (cid:104) σ σ σ σ σ σ (cid:105) T . (56)For the general multilayer system, we can assemble the complete set of nonlinear residual equations as R ( X ) = B T { σ } − tB T { σ } − t ... B T { σ M } − tt − tt − t ... t M − t (cid:80) m φ m v m + (cid:80) nm w nm = with X = v v ... v M w w ... w M t . (57)Note that any redundant slip unknowns should be condensed out to avoid an ill-posed system. The solutionto these equations can then be found using, e.g., Newton’s method. Let J = ∂ R /∂ X denote the Jacobianmatrix that arises from linearization of the nonlinear residual. Given a current estimate X k for the solution, animproved estimate can be found by: solving J k ∆ X = − R k , (58)updating X k + = X k + α ∆ X . (59)Here, α ∈ (0 ,
1] is a suitably chosen line-search or backtracking parameter than can be used to increase therobustness of the search when far away from the neighborhood of convergence. The sequence is repeated untilthe residual norm drops below a desired tolerance. At convergence, we then obtain the macroscopic stresstensor as the volumetric average: Σ = (cid:88) m ∈ L φ m σ m (60)It is perhaps easiest to understand the algorithm structure by examining the special case of a bi-layer medium—i.e. with two layers, A and B , and two interfaces, AB and BA . We first identify the symmetry w AB = w BA , sothat one slip unknown is made redundant. The unknowns and residual equations are then R ( X ) = B T { σ A } − tB T { σ B } − tt AB − t φ A v A + φ B v A + w AB = with X = v A v B w AB t . (61)Note the appearance of the 2 is the compatibility condition, as we have two surfaces with identical slip w AB .The Jacobian for this system is J = ∂ R ∂ X = B T C A B − B T C B B − D AB − φ A φ B , (62)in which denotes a diagonal identity matrix, and is two times the identity matrix. Here, C A and C B denotethe 6 × consistent (or algorithmic) sti ff ness tangents of layers A and B , which should beavailable from the layer material model subroutines. Similarly, D AB is the tangent from the traction model.Throughout this section, all tangent operators are consistent operators, rather than the continuum versions, inorder to maintain optimal convergence rates.At convergence, we may derive the homogenized consistent tangent operator, C = ∂ { Σ } ∂ { E } = (cid:88) m ∈ L φ m ∂ { σ m } ∂ { E } , (63)in which ∂ { σ m } ∂ { E } = ∂ { σ m } ∂ { E } (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X + ∂ { σ m } ∂ X (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E ∂ X ∂ { E } = C m + C m B ∂ v m ∂ { E } . (64)12o compute the final derivative, we examine the converged residual R ( X ∗ ) = , in which X ∗ denotes thesolution to the local iteration. At this configuration, ∂ R ∂ { E } = ∂ R ∂ { E } (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X + ∂ R ∂ X (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E ∂ X ∂ { E } = . (65)Let G = J − be the inverse of the local tangent evaluated at X ∗ . Solving (65) gives ∂ X ∂ { E } = − G ∂ R ∂ { E } (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X , (66)from which we infer ∂ v m ∂ { E } = − (cid:88) n ∈ L G mn B T C n . (67)Here, G mn denotes the 3 × G whose row and columns indices correspond to the residualequations for layer m and n , respectively. Combining the previous equations, the homogenized tangent operatorin matrix form is finally obtained as, C = (cid:88) m ∈ L φ m C m − (cid:88) m ∈ L (cid:88) n ∈ L φ m C m B G mn B T C n (68)It is interesting to observe that this operator is composed of a volume weighted average term plus an inter-layerinteraction term. This latter e ff ect induces anisotropic behavior, even if the layers are isotropic materials. The final material point update is summarized in Algorithm 1. The state is fully known at time t n − , and thegoal is to compute an updated state at time t n given a new macroscopic strain E n as input. We use κ n − m and κ n − nm to denote, respectively, internal state parameters for layer m and interface nm from the previous timestep. Notethat this routine can be coded in such a way that the micro-constitutive models are called in a polymorphicfashion, so that one homogenization routine can employ a variety of micro-constitutive models on demand. Algorithm 1.
Material point update.
Input: E n Output: Σ n and C n
1. Set local iteration counter k = X k = .2. For each layer m ∈ L : • From X k , compute micro-strain e m . • Call material subroutine for layer m :Inputs: ε m = ( E n + e m ), ε n − m , κ n − m Outputs: σ m , C m , κ m
3. For each interface nm ∈ I : • From X k , compute jump w nm . • Call material subroutine for interface nm :Inputs: w nm , w n − nm , κ n − nm Outputs: t nm , D nm , κ nm
4. Assemble residual R k (cid:16) X k (cid:17)
5. If (cid:107) R k (cid:107) < (cid:15) tol go to Step 8.6. Else, perform Newton update X k + = X k − α (cid:16) J − (cid:17) k R k .7. Set k ← k + Σ n = (cid:80) m φ m σ m
9. Compute the consistent tangent operator C n via equation (68).10. Save the sub-model states in preparation for the next timestep: • For each layer: ε nm ← ε m , κ nm ← κ m • For each interface: w nnm ← w nm , κ nnm ← κ nm W ¶W x x x n Y y y y z z z n Z Y Z Z ¶ Y ¶ Z ¶ Z n Z Figure 4: Schematic representation of a three-scale domain. The macroscopic domain Ω is a periodic domain consisting of a multi-layerunit cell Y with period (cid:15) y . Each layer is itself a periodic domain consisting of multi-layer unit cell Z with period (cid:15) z . Consistency of mechanical work calculated at both scales must be assured [41] leading to a Hill-Mandelcondition [65]. For the present problem, this would require Σ : E = (cid:88) m ∈ L φ m σ m : ε m + (cid:88) nm ∈ I t nm · w nm . (69)We can use Equations (47)–(50) and (60) to expand the right hand side as (cid:88) m ∈ L φ m σ m : ε m + (cid:88) nm ∈ I t nm · w nm = (cid:88) m ∈ L φ m σ m : ( E + v m ⊗ s n ) + t · (cid:88) nm ∈ I w nm = Σ : E + (cid:88) m ∈ L φ m ( σ m · n ) · v m + t · (cid:88) nm ∈ I w nm = Σ : E + t · (cid:88) m ∈ L φ m v m + (cid:88) nm ∈ I w nm = Σ : E , (70)confirming the Hill-Mandel condition is indeed satisfied. The formulation derived in the previous sections is based on the assumption that all layers are individuallyhomogeneous. The constituent model for each layer, however, could itself arise from a homogenization pro-cess. For example, Figure 4 shows a macroscopic domain Ω composed of a multi-layer unit cell Y with period (cid:15) y and layer normal n y . Each layer in unit cell Y , however, is itself a periodic domain unit cell Z with period (cid:15) z and layer normal n z . The two-scale approach presented in the previous sections can then be recursivelyextended to multiple scales, provided that all scales are su ffi ciently separate. That is, (cid:15) z (cid:28) (cid:15) y (cid:28) . (71)In this case, Step 2 in the numerical algorithm described in Algorithm 1 would consist of a recursive call tothe homogenization routine to obtain the response of each Z-periodic layer in Y . We note that a number ofresearchers have extended the asymptotic homogenization approach beyond two scales [33]. See [58] for areview of multi-scale asymptotic homogenization.
5. Numerical examples
We now present several numerical examples using the proposed constitutive modeling framework. The aimis to demonstrate its performance in capturing mechanical response of a variety of rock types, as observed fromexperiments. This includes anisotropy in strength and elastic properties, as well as changing modes of failure.To demonstrate the framework’s flexibility, we will incorporate a number of di ff erent micro-constitutive lawsfor the layers and interfaces. These sub-models are briefly described here, but more extensive discussions oftheir implementation and extensions can be found elsewhere [112, 113]. More sophisticated formulations canbe readily included as needed. 14 .1. Specific micro-constitutive models Each layer is assumed to be homogeneous, isotropic, and elastoplastic. The strain tensor is additivelypartitioned into elastic and plastic components as ε = ε e + ε p . (72)The stress is computed from an isotropic, linear model as σ = C e : ε e , with the elastic moduli C e a function ofthe bulk modulus K and Poisson ratio ν . Evolution of plastic strains in each layer is determined by a convexyield surface and an associative flow rule, f ( σ , κ ) = , ˙ ε p = ˙ λ ∂ f ∂ σ , (73)where κ and ˙ λ ≥ f ( σ m , c ) = q + (tan ϕ ) p − c (74)in which p = tr ( σ ) / s = σ − p is the deviatoric stress, and q = √ / (cid:107) s (cid:107) is the von Misesstress. In addition, ϕ and c denote the friction angle and cohesion, respectively. Hardening (or softening) isintroduced in the model as ˙ c = h ˙ λ with modulus h controlling the cohesion-hardening rate.The second yield criterion considered is the Modified Cam-Clay surface, f ( σ , p c ) = q M + p ( p − p c ) (75)in which p c is the pre-consolidation pressure. The pre-consolidation pressure depends on volumetric plasticstrains as ˙ p c = h tr(˙ ε p ).For each interface, the slip w is additively partitioned into elastic and inelastic components as w = w e + w p . (76)The interface sti ff ness follows a linear model t = D e · w e , with the elastic moduli D e a function of the normalsti ff ness k and shear sti ff ness µ of the interface. Choosing very large sti ff ness coe ffi cients will mimic anincompressible joint with rigid / plastic shear response. Inelastic slip is governed by a Coulomb yield surfaceand non-associative flow rule, F ( t , c ) = t s + (tan ϕ ) t n − c , ˙ w p = ˙ λ ∂ t s ∂ t . (77)Here, t n = t · n is the normal traction magnitude, t s = t − t n n is the shear traction, and t s = (cid:107) t s (cid:107) is the sheartraction magnitude. This model assumes a dilation-free inelastic slip. See [112] for additional details andextensions. Tien et al. [114] used di ff erent weight ratios of cement, kaolinite, and water to produce two cementitiousmaterials with di ff erent strength and sti ff ness values. They created a synthetic rock by layering these twomaterials in an alternating fashion. Subsequently, triaxial tests were conducted at various confining pressures Parameter Symbol Units Material A Material B Interface(Drucker-Prager) (Drucker-Prager) (Coulomb)Phase fraction φ m – 0.5 0.5 –Bulk modulus K MPa 13,395 6,840 –Poisson’s ratio ν – 0.23 0.21 –Normal sti ff ness k MPa – – ∞ Shear sti ff ness µ MPa – – ∞ Friction angle ϕ ◦
35 18 18Cohesion c MPa 90 35 11Hardening modulus h MPa 0 0 0Table 1: Parameters calibrated using data from the Tien et al. [114] experiments on synthetic layered rock. ormal stress on plane of weakness (MPa)0 10 20 30 40 50 60 S hea r s t r e ss on t he p l ane o f w ea k ne ss ( M P a ) Figure 5: Shear strength versus normal stress acting on an interface, for bedding plane orientation θ = ◦ , for the synthetic layered rock. Bedding plane orientation, degrees M a x i m u m s t r eng t h , M P a p = 6 MPap = 14 MPap = 35 MPap = 6 MPa (experiment)p = 14 MPa (experiment)p = 35 MPa (experiment) Figure 6: Variation of maximum strength with respect to bedding plane orientation and confining pressure for the synthetic layered rock. C on f i n i ng p r e ss u r e , M P a Failure along interfaces (simulation) Tien and Kuo criterion Failure along interfaces (experiment) Failure in rock matrix (experiment) Tensile fracture in rock matrix (experiment) Tensile splitting along interfaces (experiment)
Figure 7: Comparison of the failure modes as predicted by Tien and Kuo’s criterion [5, 114], experiments, and simulations using thepresent method. The zone inside the marked regions corresponds to failure along discontinuities. θ , defined as the angle between the plane of the layers and the horizontal axis.These experiments are an appealing test because the material microstructure is highly controlled, and multipleexperiments in di ff erent orientations can be readily performed.To mimic this synthetic microstructure, we use Drucker-Prager plasticity models for the two materials, anda rigid-plastic Coulomb model for the interface. Calibrated material parameters are summarized in Table 1.Phase ratios for materials A and B are φ A = φ B = .
5. Since θ = ◦ is reported as a case where failure occursalong the interface and not through the rock layers, we have used test results at this orientation to calibrate theCoulomb friction and cohesion parameters (Figure 5). The intercept and slope on this plot correspond to thecohesion and tangent of friction angle, respectively.Figure 6 compares peak strength results across a range of bedding orientations and confining pressures. Itcan be seen that the simple model presented here captures the complex strength anisotropy of the synthetic rockand its pressure dependence reasonably well. We note, however, that Tien et al.’s experimental results underzero confining pressure are not included in this figure, as tensile fractures in the rock matrix were observedduring unconfined tests. This particular failure mode is not possible using the Drucker-Prager model, thoughin the future a tensile failure criterion could be added to capture this particular mode if desired.Di ff erent failure modes observed from the experiments are shown in Figure 7. The regions correspondingto failure along an interface, as predicted by the proposed model and Tien and Kuo’s criterion [5, 114] arealso shown. It can be seen that the simulation results are consistent with the experimentally observed failurealong the interfaces. At low confining pressures, the stress acting on the interface becomes critical and leadsto failure for bedding plane orientations of ∼ ◦ to 80 ◦ . At increased confining pressure, the strength of theinterfaces increases due to the higher normal stress, and the range of bedding plane orientations at which theshear stress reaches a critical value decreases. At su ffi ciently high confining pressures, failure in the rockmatrix will dominate. We now explore triaxial test data on Vaca Muerta shale from Ambrose [115]. Vaca Muerta shale has a veryhigh organic content that is well-dispersed throughout the shale matrix. X-ray tomography images obtainedby Ambrose show laminations and weak planes in the samples, while on thin sections one dominant rockcomposition can be observed. Therefore, we propose to model the shale with a single matrix layer (phasefraction φ = .
0) with Drucker-Prager plasticity. Matrix layers are then separated by weak interfaces treatedusing a Coulomb frictional model, with elastic normal and shear compliance.The calibrated model parameters are summarized in Table 2. Figure 8 shows the variation of maximumstrength with respect to bedding plane orientations at confining pressures p = . , . , . , and 137 . ν and ν . Except for ν , for which the experimental dataappear very scattered, the other two elasticity parameters exhibit a smooth ascending trend as bedding planeorientations vary from 0 ◦ to 90 ◦ , which is captured reasonably well by the present model. We note that theselected elastic model has no intrinsic pressure dependence. It is di ffi cult to discern a consistent trend heregiven the experimental scatter in the data, however, so we have ignored this feature in the comparison. Parameter Symbol Units Matrix Interface(Drucker-Prager) (Coulomb)Bulk modulus K MPa 17,390 –Poisson’s ratio ν – 0.27 –Normal sti ff ness k MPa – 70,000Shear sti ff ness µ MPa – 52,500Friction angle ϕ ◦
47 26Cohesion c MPa 70 18Hardening modulus h MPa 0 0Table 2: Parameters calibrated using experimental data from Ambrose [115].
In this example, we use true triaxial test data from Mogi [116] for a schist rock from Chichibu on HonshuIsland, Japan. This metamorphic crystalline rock has a densely foliated structure [117]. For modeling purposes,17
10 20 30 40 50 60 70 80 90Bedding plane orientation, degrees50100150200250300350400450500 M a x i m u m s t r eng t h , M P a p = 6.9 MPap = 17.2 MPap = 34.5 MPap = 137.9 MPap = 6.9 MPa (experiment)p = 17.2 MPa (experiment)p = 34.5 MPa (experiment)p = 137.9 MPa (experiment) Figure 8: Variation of maximum strength with respect to orientation of bedding planes obtained from triaxial tests on Vaca Muerta shale.
Bedding plane orientation, degrees E l a s t i c m odu l u s , M P a Bedding plane orientation, degrees P o i ss on ' s r a t i o , Bedding plane orientation, degrees P o i ss on ' s r a t i o , p = 6.9 MPa p = 17.2 MPa p = 34.5 MPa p = 137.9 MPaSimulation Figure 9: Variation of macroscopic (apparent) elastic modulus in the loading direction (axis 3) and Poisson’s ratios, with respect toorientation of bedding planes obtained from triaxial tests on Vaca Muerta shale. arameter Symbol Units Matrix Interface(Drucker-Prager) (Coulomb)Bulk modulus K MPa 16880 –Poisson’s ratio ν – 0.3 –Normal sti ff ness k MPa – ∞ Shear sti ff ness µ MPa – ∞ Friction angle ϕ ◦ Intermediate stress, MPa M a x i m u m s t r e ss , M P a Mode 1: =60 , =0 (experiment)Mode 2: =60 , =45 (experiment)Mode 3: =60 , =90 (experiment)Mode 4: =0 (experiment)Mode 1Mode 2Mode 3Mode 4 Y Y Y s s s b w Figure 10: Variation of maximum stress, σ , with respect to intermediate principal stress, σ . Minimum principal stress is σ =
50 MPafor all cases.
80 100 120 140 160 180 200 220 240Normal stress, MPa6080100120140 S hea r s t r e ss , M P a Mode 1 (simulation)Mode 2 (simulation)Mode 1 (experiment)Mode 2 (experiment)
Figure 11: Shear stress versus normal stress acting on the plane of weakness at failure, corresponding to Mode 1 and Mode 2.
19 single Drucker-Prager matrix type and rigid-plastic Coulomb interfaces are considered, with parametersprovided in Table 3.The orientation of the bedding planes in true triaxial tests is determined by two angles, β and ω , whichare defined with respect to the principal stress directions σ > σ > σ . Figure 10 shows variations ofmaximum strength in the direction of σ versus intermediate stress σ , while σ =
50 MPa is held constantin all experiments. The results obtained from simulations are compared with experiments in four di ff erentmodes, showing good agreement across a variety of failure mechanisms. Mode 1 ( β = ◦ , ω = ◦ ) and Mode2 ( β = ◦ , ω = ◦ ) lead to failure along the interfaces, while in Mode 4 ( β = ◦ ) failure is only observedin the rock matrix. Mode 3 ( β = ◦ , ω = ◦ ) exhibits a transition from failure along interfaces at smallervalues of σ to failure in the rock matrix at higher values of σ . Note that the data points marked with anarrow in Figure 10 correspond to a transitional mode of failure and are disregarded in the calibration process.Figure 11 shows the variation of shear stress versus normal stress acting on the interfaces at failure for Modes1 and 2. For Mode 1, we have ω =
0. Changes in the intermediate stress σ in the simulations do not a ff ectthe stress conditions on the interface, and we obtain only one point corresponding to this mode (red circle inFigure 11). Experimental data corresponding to Mode 1 shows some small variation in the maximum strengthwith respect to σ , while these changes are much more pronounced for other modes with non-zero angle ω . In this final example, we focus on capturing a transition from brittle to ductile behavior in geologic mate-rials. This example is motivated by oil shale experiments [118] which show a clear transition from brittle toductile behavior with increasing organic content. For this purpose, we present results for a hypothetical ma-terial with two perfectly bonded layers (no weak interface). Material A has a ductile behavior modeled withModified Cam-Clay. Layer B has a brittle behavior modeled with Drucker-Prager. Material parameters areindicated in Table 4. We will vary the phase fraction of each material to control the fraction of brittle versusductile constituents.Triaxial simulations are conducted with a constant confining pressure of p =
10 MPa and bedding planeorientations of θ = ◦ , ◦ , and 90 ◦ . Figure 12 illustrates variation of the macroscopic deviatoric stress versusmacroscopic axial and radial strains, as well as variation of macroscopic volumetric strain versus axial strainfor various ratios of phase A , from φ A = . φ A = .
9. It should be noted that the two components of radialstrains in two di ff erent directions are not identical, and only one is shown here for clarity. As expected, uponincreasing the ratio of the ductile phase, the material shows a more ductile overall macroscopic behavior. Wealso observe a significant reduction in compaction of the material for θ = ◦ and 90 ◦ upon increasing the ratioof the ductile phase. A transition from macroscopic compaction to dilation is observed for θ = ◦ , when theductile phase comprises around 70% of the volume. Parameter Symbol Units Material A Material B(Modified Cam-Clay) (Drucker-Prager)Bulk modulus K MPa 26.7 40Poisson’s ratio ν – 0.25 0.25CSL slope M – 1.5 –Pre-consolidation pressure p c MPa 10 –Friction angle ϕ ◦ – 50Cohesion c MPa – 5Hardening modulus h MPa 5000 -200Table 4: Parameters for the brittle-ductile material example.
6. Concluding remarks
We have presented an inelastic homogenization framework for layered materials with planes of weakness.Key features of the proposed approach include explicit representation of layered microstructures, potentialfor slip along interface discontinuities, and the ability to include arbitrary combinations of micro-constituentmodels within the same numerical framework. This approach was applied to a diverse variety of experimentaldata—e.g. layered synthetic rocks, Vaca Muerta shale, and Chichibu schist. The proposed model performswell in capturing both strength and sti ff ness anisotropy. Moreover, it was shown that the models can describefailure modes in the rock matrix and along planes of weakness in a manner consistent with experimental20
15 -10 -5 0 5 10 15
Axial and radial millistrain D e v i a t o r i c s t r e ss ( M P a ) Axial millistrain -505101520 V o l u m e t r i c m illi s t r a i n A =0.1 A =0.3 A =0.5 A =0.7 A =0.9 = 0 o -15 -10 -5 0 5 10 15 Axial and radial millistrain D e v i a t o r i c s t r e ss ( M P a ) Axial millistrain -505101520 V o l u m e t r i c m illi s t r a i n A =0.1 A =0.3 A =0.5 A =0.7 A =0.9 = 45 o -20 -10 0 10 Axial and radial millistrain D e v i a t o r i c s t r e ss ( M P a ) Axial millistrain V o l u m e t r i c m illi s t r a i n A = 0.1 A = 0.3 A = 0.5 A = 0.7 A = 0.9 = 90 o Figure 12: Variation of deviatoric stress versus axial and radial strains (left) and volumetric strain versus axial strain (right) for phase ratiosof the ductile layer varying between 0 . .
9. Panels from (top) to (bottom) indicate bedding plane orienations of θ = ◦ , ◦ , and 90 ◦ . observations. Numerical examples for a hypothetical bi-layer material also demonstrate the ability to describecomplex phenomenology for materials with mixed brittle and ductile constituents.A number of failure criteria for transversely isotropic rocks have been developed in the literature [5, 7–9].A useful example is Jaeger’s plane of weakness model [7], which employs a Mohr-Coulomb failure criterionfor the rock matrix and sliding on planes of weakness. The present method can exactly mimic Jaeger’s modelwhen a unit cell with one rock layer and one rigid-plastic interface is used. This model, however, will describethe complete stress-strain response of the material, not just the failure envelope.Anisotropic elasto-plastic continuum models have also been developed in the literature to capture anisotropyin the mechanical response of transversely isotropic materials—e.g. by extension of isotropic Drucker-Prager[4] and Cam-Clay type models [3, 119]. These models have also been applied to strength anisotropy of rocks[4, 115, 120]. One advantage of the present framework is that the overall anisotropic behavior emerges nat-urally from the interaction of isotropic layers within a clearly defined microstructure. We therefore view theproposed homogenization methodology as a natural framework for designing anisotropic material models,allowing model parameters to be directly connected to a microstructural model in a physically motivated way. Acknowledgements
Funding for this work was provided by Total S.A. through the FC-MAELSTROM Project. Portions of thiswork were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore NationalLaboratory under Contract DE-AC52-07-NA27344. 21 eferencesReferences [1] M. D. Salamon, Elastic moduli of a stratified rock mass, International Journal of Rock Mechanics and Mining Sciences & Geome-chanics Abstracts 5 (6) (1968) 519–527. doi:10.1016/0148-9062(68)90039-9 .[2] A. Sawicki, On application of e ff ective moduli theory to layered soil, Rozprawy Hydrotechniczne 39 (1978) 3–13.[3] S. J. Semnani, J. A. White, R. I. Borja, Thermoplasticity and strain localization in transversely isotropic materials based onanisotropic critical state plasticity, International Journal for Numerical and Analytical Methods in Geomechanics 40 (18) (2016)2423–2449. doi:10.1002/nag.2536 .[4] W. G. Pariseau, Plasticity Theory For Anisotropic Rocks And Soil, in: The 10th US Symposium on Rock Mechanics (USRMS),1968.[5] Y. M. Tien, M. C. Kuo, A failure criterion for transversely isotropic rocks, International Journal of Rock Mechanics and MiningSciences 38 (3) (2001) 399–412.[6] L. T. Drzal, The role of the fiber-matrix interphase on composite properties, Vacuum 41 (7-9) (1990) 1615–1618. doi:10.1016/0042-207X(90)94034-N .[7] J. C. Jaeger, Shear failure of anistropic rocks, Geological Magazine 97 (1) (1960) 65–72. doi:10.1017/S0016756800061100 .[8] E. Hoek, Estimating Mohr-Coulomb friction and cohesion values from the Hoek-Brown failure criterion, International Journal ofRock Mechanics and Mining Sciences & Geomechanics Abstracts 27 (3) (1990) 227–229.[9] E. Hoek, C. Carranza-Torres, B. Corkum, Hoek-Brown failure criterion-2002 edition, Proceedings of NARMS-Tac 1 (2002) 267–273.[10] P. Kanout´e, D. P. Boso, J. L. Chaboche, B. A. Schrefler, Multiscale methods for composites: A review, Archives of ComputationalMethods in Engineering 16 (1) (2009) 31–75. doi:10.1007/s11831-008-9028-8 .[11] M. G. D. Geers, V. G. Kouznetsova, K. Matouˇs, J. Yvonnet, Homogenization Methods and Multiscale Modeling: NonlinearProblems, in: E. Stein, R. de Borst, T. J. Hughes (Eds.), Encyclopedia of Computational Mechanics Second Edition, John Wiley &Sons, Ltd., 2017. doi:10.1002/9781119176817.ecm2107 .[12] G. A. Pavliotis, A. M. Stuart, Multi-scale methods: Averaging and homogenization, Springer Science & Business Media, NewYork, 2008.[13] J. M. Ortolano Gonz´alez, J. A. Hern´andez Ortega, X. Oliver Olivella, A comparative study on homogenization strategies formulti-scale analysis of materials, Centre Internacional de M`etodes Num`erics en Enginyeria (CIMNE), 2013.[14] J. Vondrejc, FFT-based method for homogenization of periodic media: Theory and applications FFT-based micro-mechanicalsolver View project, Ph.D. thesis, Czech Technical University in Prague (2013). doi:10.13140/RG.2.1.2534.2489 .URL [15] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear mechanical properties of composites, C.R Acad. Sci. Paris 318 (1994) 1417–1423.[16] J. C. Michel, H. Moulinec, P. Suquet, E ff ective properties of composite materials with periodic microstructure: A computational ap-proach, Computer Methods in Applied Mechanics and Engineering 172 (1-4) (1999) 109–143. doi:10.1016/S0045-7825(98)00227-8 .[17] G. J. Dvorak, Transformation field analysis of inelastic composite materials, Proc. R. Soc. Lond. A 437 (1992) 311–327.URL http://rspa.royalsocietypublishing.org/ [18] S. Marfia, E. Sacco, Computational homogenization of composites experiencing plasticity, cracking and debonding phenomena,Computer Methods in Applied Mechanics and Engineering 304 (2016) 319–341. doi:10.1016/j.cma.2016.02.007 .[19] V. Sepe, S. Marfia, E. Sacco, A nonuniform TFA homogenization technique based on piecewise interpolation functions of theinelastic field, International Journal of Solids and Structures 50 (5) (2013) 725–742. doi:10.1016/j.ijsolstr.2012.11.005 .[20] J. C. Michel, P. Suquet, Nonuniform transformation field analysis, International journal of solids and structures 40 (25) (2003)6937–6955. doi:10.1016/S0020-7683(03)00346-9 .[21] S. Wulfingho ff , S. Reese, E ffi cient computational homogenization of simple elasto-plastic microstructures using a modified Ritz-Galerkin approach, in: COMPLAS XIII: proceedings of the XIII International Conference on Computational Plasticity: fundamen-tals and applications, CIMNE, 2015, pp. 956–964.[22] V. Kouznetsova, M. G. Geers, W. A. Brekelmans, Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme, International Journal for Numerical Methods in Engineering 54 (8) (2002) 1235–1260. doi:10.1002/nme.541 .[23] A. Lagzdinˇs, Orientational averaging in mechanics of solids, Vol. 265, Longman Scientific and Technical, 1992.[24] R. Hill, Continuum micro-mechanics of elastoplastic polycrystals, Journal of the Mechanics and Physics of Solids 13 (2) (1965)89–101. doi:10.1016/0022-5096(65)90023-2 .[25] R. R. Hill, A self-consistent mechanics of composite materials, Mech. Phys. Solids 13 (1965) 213–222.[26] E. Kr¨oner, Zur plastischen verformung des vielkristalls, Acta metallurgica 9 (2) (1961) 155–161. doi:10.1016/0001-6160(61)90060-8 .[27] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. Lond. A241 (1226) (1957) 376–396. doi:10.1098/rspa.1957.0133 .[28] P. Ponte Casta˜neda, The e ff ective mechanical properties of nonlinear isotropic composites, Journal of the Mechanics and Physicsof Solids 39 (1) (1991) 45–71. doi:10.1016/0022-5096(91)90030-R .[29] P. Ponte Casta˜neda, Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: I - Theory,Journal of the Mechanics and Physics of Solids 50 (4) (2002) 737–757. doi:10.1016/S0022-5096(01)00099-0 .[30] P. M. Suquet, Overall potentials and extremal surfaces of power law or ideally plastic composites, Journal of the Mechanics andPhysics of Solids 41 (6) (1993) 981–1002. doi:10.1016/0022-5096(93)90051-G .[31] N. Bakhvalov, G. Panasenko, Homogenisation: Averaging Processes in Periodic Media: mathematical problems in the mechanicsof composite materials, Kluwer, Dordrecht, The Netherlands, 1989. doi:10.1007/978-94-009-2247-1 .[32] R. Penta, A. Gerisch, An Introduction to Asymptotic Homogenization, in: A. Gerisch, R. Penta, J. Lang (Eds.), Multiscale Modelsin Mechano and Tumor Biology, Springer, 2018. doi:10.1007/978-3-319-73371-5_1 .[33] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, Vol. 374, American Mathematical Soc.,2011.
34] J. L. Auriault, C. Boutin, C. Geindreau, Homogenization of Coupled Phenomena in Heterogenous Media, John Wiley & Sons,2010. doi:10.1002/9780470612033 .[35] M. Hori, S. Nemat-Nasser, On two micromechanics theories for determining micro-macro relations in heterogeneous solids, Me-chanics of Materials 31 (1999) 667–682. doi:10.1016/S0167-6636(99)00020-4 .[36] I. Doghri, A. Ouaar, Homogenization of two-phase elasto-plastic composite materials and structures Study of tangent operators,cyclic plasticity and numerical algorithms, International Journal of Solids and Structures 40 (7) (2003) 1681–1712. doi:10.1016/S0020-7683(03)00013-1 .[37] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta Metallurgica21 (5) (1973) 571–574. doi:10.1016/0001-6160(73)90064-3 .[38] S. Mercier, A. Molinari, S. Berbenni, M. Berveiller, Comparison of di ff erent homogenization approaches for elastic-viscoplasticmaterials, Modelling and Simulation in Materials Science and Engineering 20 (2) (2012) 024004. doi:10.1088/0965-0393/20/2/024004 .[39] S. Nemat-Nasser, M. Hori, Micromechanics: overall properties of heterogeneous materials, Elsevier, 1993.[40] S. Nemat-Nasser, Averaging theorems in finite deformation plasticity, Mechanics of Materials 31 (8) (1999) 493–523. doi:10.1016/S0167-6636(98)00073-8 .[41] E. S. Perdahco˘glu, H. J. M. Geijselaers, Constitutive modeling of two phase materials using the mean field method for homoge-nization, International Journal of Material Forming 4 (2) (2011) 93–102. doi:10.1007/s12289-010-1007-6 .[42] A. L. Kalamkarov, I. V. Andrianov, V. V. Danishevs’kyy, Asymptotic Homogenization of Composite Materials and Structures,Applied Mechanics Reviews 62 (3) (2009) 030802. doi:10.1115/1.3090830 .[43] G. Allaire, Homogenization and two-scale convergence, SIAM Journal on Mathematical Analysis 23 (6) (1992) 1482–1518. doi:10.1137/0523084 .[44] V. V. Zhikov, On an extension of the method of two-scale convergence and its applications, Sbornik: Mathematics 191 (7) (2000)973. doi:10.1070/SM2000v191n07ABEH000491 .[45] J. Fish, K. Shek, M. Pandheeradi, M. S. Shephard, Computational plasticity for composite structures based on mathematicalhomogenization: Theory and practice, Computer Methods in Applied Mechanics and Engineering 148 (1997) 53–73. doi:10.1016/S0045-7825(97)00030-3 .[46] P. W. Chung, K. K. Tamma, R. R. Namburu, Asymptotic expansion homogenization for heterogeneous media: computationalissues and applications, Composites Part A: Applied Science and Manufacturing 32 (9) (2001) 1291–1301. doi:10.1016/S1359-835X(01)00100-2 .[47] J. Pinho-da Cruz, J. A. Oliveira, F. Teixeira-Dias, Asymptotic homogenisation in linear elasticity. Part I: Mathematical formulationand finite element modelling, Computational Materials Science 45 (4) (2009) 1073–1080. doi:10.1016/j.commatsci.2009.02.025 .[48] J. A. Oliveira, J. Pinho-da Cruz, F. Teixeira-Dias, Asymptotic homogenisation in linear elasticity. Part II: Finite element proceduresand multiscale applications, Computational Materials Science 45 (4) (2009) 1081–1096. doi:10.1016/j.commatsci.2009.01.027 .[49] S. Jansson, Homogenized nonlinear constitutive properties and local stress concentrations for composites with periodic internalstructure, International Journal of Solids and Structures 29 (17) (1992) 2181–2200. doi:10.1016/0020-7683(92)90065-2 .[50] J. C. L´opez-Realpozo, R. Rodr´ıguez-Ramos, R. Guinovart-D´ıaz, J. Bravo-Castillero, L. P. Fern´andez, F. J. Sabina, G. A. Mau-gin, E ff ective properties of non-linear elastic laminated composites with perfect and imperfect contact conditions, Mechanics ofAdvanced Materials and Structures 15 (5) (2008) 375–385. doi:10.1080/15376490801977742 .[51] X. Markensco ff , C. Dascalu, Asymptotic homogenization analysis for damage amplification due to singular interaction of micro-cracks, Journal of the Mechanics and Physics of Solids 60 (8) (2012) 1478–1485. doi:10.1016/j.jmps.2012.04.004 .[52] Y. Yang, F. Y. Ma, C. H. Lei, Y. Y. Liu, J. Y. Li, Nonlinear asymptotic homogenization and the e ff ective behavior of layeredthermoelectric composites, Journal of the Mechanics and Physics of Solids 61 (8) (2013) 1768–1783. doi:10.1016/j.jmps.2013.03.006 .[53] J. J. Telega, S. Tokarzewski, A. Gałka, E ff ective Conductivity of Nonlinear Two-Phase Media: Homogenization and Two-PointPad´e Approximants, Acta Applicandae Mathematicae 61 (1-3) (2000) 295–315.[54] F. Devries, H. Dumontet, G. Duvaut, F. L´en´e, Homogenization and damage for composite structures, International Journal forNumerical Methods in Engineering 27 (2) (1989) 285–298. doi:10.1002/nme.1620270206 .[55] J. Fish, Q. Yu, K. Shek, Computational damage mechanics for composite materials based on mathematical homogenization, Interna-tional Journal for Numerical Methods in Engineering 45 (1999) 1657–1679. doi:10.1002/(SICI)1097-0207(19990820)45:11<1657::AID-NME648>3.0.CO;2-H .[56] J. Fish, R. Fan, Mathematical homogenization of nonperiodic heterogeneous media subjected to large deformation transient load-ing, International Journal for Numerical Methods in Engineering 76 (7) (2008) 1044–1064. doi:10.1002/nme.2355 .[57] J. Fish, V. Belsky, Multi-grid method for periodic heterogeneous media Part 2: Multiscale modeling and quality control inmultidimensional case, Computer Methods in Applied Mechanics and Engineering 126 (1) (1995) 17–38. doi:10.1016/0045-7825(95)00812-F .[58] A. Ram´ırez-Torres, R. Penta, R. Rodr´ıguez-Ramos, J. Merodio, F. J. Sabina, J. Bravo-Castillero, R. Guinovart-D´ıaz, L. Preziosi,A. Grillo, Three scales asymptotic homogenization and its application to layered hierarchical hard tissues, International Journal ofSolids and Structures 130 (2018) 190–198. doi:10.1016/j.ijsolstr.2017.09.035 .[59] R. Rezakhani, G. Cusatis, Asymptotic expansion homogenization of discrete fine-scale models with rotational degrees of freedomfor the simulation of quasi-brittle materials, Journal of the Mechanics and Physics of Solids 88 (2016) 320–345. doi:10.1016/j.jmps.2016.01.001 .[60] W. Li, R. Rezakhani, C. Jin, X. Zhou, G. Cusatis, A multiscale framework for the simulation of the anisotropic mechanicalbehavior of shale, International Journal for Numerical and Analytical Methods in Geomechanics 41 (14) (2017) 1494–1522. doi:10.1002/nag.2684 .[61] V. P. Smyshlyaev, K. D. Cherednichenko, On rigorous derivation of strain gradient e ff ects in the overall behaviour of periodic het-erogeneous media, Journal of the Mechanics and Physics of Solids 48 (6-7) (2000) 1325–1357. doi:10.1016/S0022-5096(99)00090-3 .[62] J. Fish, V. Belsky, Multigrid method for periodic heterogeneous media Part 1: Convergence studies for one-dimensional case,Computer Methods in Applied Mechanics and Engineering 126 (1-2) (1995) 1–16. doi:10.1016/0045-7825(95)00811-E .[63] N. Triantafyllidis, S. Bardenhagen, The influence of scale size on the stability of periodic solids and the role of associated higher rder gradient continuum models, Journal of the Mechanics and Physics of Solids 44 (11) (1996) 1891–1928. doi:10.1016/0022-5096(96)00047-6 .[64] I. V. Andrianov, V. I. Bolshakov, V. V. Danishevs’kyy, D. Weichert, Higher order asymptotic homogenization and wave propagationin periodic composite materials, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 464 (2008)1181–1201. doi:10.1098/rspa.2007.0267 .[65] P. M. Suquet, Elements of Homogenization for Inelastic Solid Mechanics, in: E. Sanchez-Palencia, A. Zaoui (Eds.), Homogeniza-tion techniques for composite media, Springer-Verlag, 1987, pp. 193–278.[66] E. Pruchnicki, Homogenized nonlinear constitutive law using fourier series expansion, International Journal of Solids and Struc-tures 35 (16) (1998) 1895–1913. doi:10.1016/S0020-7683(97)00128-5 .[67] E. Pruchnicki, I. Shahrour, A macroscopic elastoplastic constitutive law for multilayered media: Application to reinforced earthmaterial, International Journal for Numerical and Analytical Methods in Geomechanics 18 (7) (1994) 507–518. doi:10.1002/nag.1610180705 .[68] E. Pruchnicki, Homogenized elastoplastic properties for a partially cohesive composite material, Zeitschrift f¨ur angewandte Math-ematik und Physik ZAMP 49 (4) (1998) 568–589. doi:10.1007/s000000050109 .[69] M. N. Ensan, I. Shahrour, A macroscopic constitutive law for elasto-plastic multilayered materials with imperfect interfaces:Application to reinforced soils, Computers and Geotechnics 30 (2003) 339–345. doi:10.1016/S0266-352X(03)00007-7 .[70] P. B. Lourenc¸o, A matrix formulation for the elastoplastic homogenisation of layered materials, Mechanics of Cohesive–frictionalMaterials: An International Journal on Experiments, Modelling and Computation of Materials and Structures 1 (3) (1996) 273–294. doi:10.1002/(SICI)1099-1484(199607)1:3<273::AID-CFM14>3.0.CO;2-T .[71] J. Aboudi, M. J. Pindera, S. M. Arnold, Higher-order theory for periodic multiphase materials with inelastic phases, InternationalJournal of Plasticity 19 (6) (2003) 805–847. doi:10.1016/S0749-6419(02)00007-4 .[72] J. Aboudi, A continuum theory for fiber-reinforced elastic-viscoplastic composites, International Journal of Engineering Science20 (5) (1982) 605–621. doi:10.1016/0020-7225(82)90115-X .[73] M. Paley, J. Aboudi, Micromechanical analysis of composites by the generalized cells model, Mechanics of materials 14 (2) (1992)127–139. doi:10.1016/0167-6636(92)90010-B .[74] F. Covezzi, S. de Miranda, S. Marfia, E. Sacco, Homogenization of elastic-viscoplastic composites by the Mixed TFA, ComputerMethods in Applied Mechanics and Engineering 318 (2017) 701–723. doi:10.1016/j.cma.2017.02.009 .[75] P. A. Fotiu, S. Nemat-Nasser, Overall properties of elastic-viscoplastic periodic composites, International Journal of Plasticity12 (2) (1996) 163–190. doi:10.1016/S0749-6419(96)00002-2 .[76] K. P. Walker, A. D. Freed, E. H. Jordan, Thermoviscoplastic analysis of fibrous periodic composites by the use of triangularsubvolumes, Composites science and technology 50 (1) (1994) 71–84. doi:10.1016/0266-3538(94)90127-9 .[77] B. Schweizer, M. Veneroni, The needle problem approach to non-periodic homogenization, Networks and Heterogeneous Media6 (4) (2011) 755–781. doi:10.3934/nhm.2011.6.755 .[78] M. Heida, B. Schweizer, Non-periodic homogenization of infinitesimal strain plasticity equations, ZAMM Zeitschrift fur Ange-wandte Mathematik und Mechanik 96 (1) (2016) 5–23. doi:10.1002/zamm.201400112 .[79] A. Visintin, On homogenization of elasto-plasticity, in: Journal of Physics. Conference Series, Vol. 22, 2005, pp. 222–234. doi:10.1088/1742-6596/22/1/015 .[80] B. Schweizer, M. Veneroni, Homogenization of plasticity equations with two-scale convergence methods, Applicable Analysis94 (2) (2015) 375–398. doi:10.1080/00036811.2014.896992 .[81] S. Nesenenko, Homogenization in viscoplasticity, SIAM Journal on Mathematical Analysis 39 (1) (2007) 236–262. doi:10.1137/060655092 .[82] G. Francfort, A. Giacomini, On periodic homogenization in perfect elasto-plasticity, Journal of the European Mathematical Society16 (3) (2014) 409–461. doi:10.4171/JEMS/437 .[83] K. Sab, Homogenization of non-linear random media by a duality method. Application to plasticity, Asymptotic Analysis 9 (4)(1994) 311–336. doi:10.3233/ASY-1994-9402 .[84] N. Ohno, X. Wu, T. Matsuda, Homogenized properties of elastic-viscoplastic composites with periodic internal structures, Interna-tional Journal of Mechanical Sciences 42 (8) (2000) 1519–1536. doi:10.1016/S0020-7403(99)00088-0 .[85] A. Ram´ırez-Torres, S. D. Stefano, A. Grillo, R. Rodr´ıguez-Ramos, J. Merodio, R. Penta, An asymptotic homogenization approachto the microstructural evolution of heterogeneous media, nternational Journal of Non-Linear Mechanics 106 (2018) 245–257. doi:10.1016/j.ijnonlinmec.2018.06.012 .[86] P. W. Chung, K. K. Tamma, R. R. Namburu, A computational approach for multi-scale analysis of heterogeneous elasto-plasticmedia subjected to short duration loads, International Journal for Numerical Methods in Engineering 59 (6) (2004) 825–848. doi:10.1002/nme.880 .[87] I. Doghri, M. I. El Ghezal, L. Adam, Finite strain mean-field homogenization of composite materials with hyperelastic-plasticconstituents, International Journal of Plasticity 81 (2016) 40–62. doi:10.1016/j.ijplas.2016.01.009 .[88] K. Danas, M. I. Idiart, P. Ponte Casta˜neda, A homogenization-based constitutive model for isotropic viscoplastic porous media,International Journal of Solids and Structures 45 (11-12) (2008) 3392–3409. doi:10.1016/j.ijsolstr.2008.02.007 .[89] P. Suquet, E ff ective Properties of Nonlinear Composites, in: P. Suquet (Ed.), Continuum Micromechanics, Springer Vienna, Vienna,1997, pp. 197–264. doi:10.1007/978-3-7091-2662-2_4 .[90] D. R. S. Talbot, J. R. Willis, Variational principles for inhomogeneous non-linear media, IMA Journal of Applied Mathematics35 (1) (1985) 39–54. doi:10.1093/imamat/35.1.39 .[91] P. Ponte Casta˜neda, Exact second-order estimates for the e ff ective mechanical properties of nonlinear composite materials, J. Mech.Phys. Solids 44 (6) (1996) 827–862. doi:10.1016/0022-5096(96)00015-4 .[92] M. Agoras, P. Ponte Casta˜neda, Homogenization estimates for multi-scale nonlinear composites, European Journal of Mechanics,A / Solids 30 (6) (2011) 828–843. doi:10.1016/j.euromechsol.2011.05.007 .[93] G. Chatzigeorgiou, N. Charalambakis, Y. Chemisky, F. Meraghni, Periodic homogenization for fully coupled thermomechanicalmodeling of dissipative generalized standard materials, International Journal of Plasticity 81 (2016) 18–39. doi:10.1016/j.ijplas.2016.01.013 .[94] W. Yu, T. Tang, Variational asymptotic method for unit cell homogenization of periodically heterogeneous materials, InternationalJournal of Solids and Structures 44 (11-12) (2007) 3738–3755. doi:10.1016/j.ijsolstr.2006.10.020 .[95] L. Zhang, W. Yu, Variational asymptotic homogenization of elastoplastic composites, Composite Structures 133 (2015) 947–958. doi:10.1016/j.compstruct.2015.07.117 .
96] G. Bouchitte, P. Suquet, Homogenization, plasticity and yield design, in: Composite media and homogenization theory, Birkh¨auser,Boston, 1991, pp. 107–133. doi:10.1007/978-1-4684-6787-1_7 .[97] R. Gl¨uge, E ff ective plastic properties of laminates made of isotropic elastic plastic materials, Composite Structures 149 (2016)434–443. doi:10.1016/j.compstruct.2016.04.029 .[98] R. Gl¨uge, E ff ective yield limits of microstructured materials, Composite Structures 176 (2017) 496–504. doi:10.1016/j.compstruct.2017.05.051 .[99] A. Sawicki, Yield conditions for layered composites, International Journal of Solids and Structures 17 (10) (1981) 969–979. doi:10.1016/0020-7683(81)90035-4 .[100] P. De Buhan, A. Taliercio, A homogenization approach to the yield strength of composite materials, European Journal of Mechan-ics. A, Solids 10 (2) (1991) 129–154.[101] P. Ponte Casta˜neda, G. deBotton, On the homogenized yield strength of two-phase composites, Proc. R. Soc. Lond. A 438 (1903)(1992) 419–431. doi:10.1098/rspa.1992.0116 .[102] G. deBotton, P. Ponte Casta˜neda, On the ductility of laminated materials, International Journal of Solids and Structures 23 (19)(1992) 2329–2353. doi:10.1016/0020-7683(92)90219-J .[103] W. Q. Shen, J. Zhang, J. F. Shao, D. Kondo, Approximate macroscopic yield criteria for Drucker-Prager type solids with spheroidalvoids, International Journal of Plasticity 99 (2017) 221–247. doi:10.1016/j.ijplas.2017.09.008 .[104] A. El Omri, A. Fennan, F. Sidoro ff , A. Hihi, Elastic-plastic homogenization for layered composites, European Journal of Mechan-ics, A / Solids 19 (4) (2000) 585–601. doi:10.1016/S0997-7538(00)00182-0 .[105] Q. C. He, Z. Q. Feng, Homogenization of layered elastoplastic composites: Theoretical results, International Journal of Non-LinearMechanics 47 (2) (2012) 367–376. doi:10.1016/j.ijnonlinmec.2011.09.018 .[106] K. Poulios, C. F. Niordson, A homogenization method for ductile-brittle composite laminates at large deformations, InternationalJournal for Numerical Methods in Engineering 113 (2018) 814–833. doi:10.1002/nme.5637 .[107] R. I. Borja, R. A. Regueiro, Strain localization in frictional materials exhibiting displacement jumps, Computer Methods in AppliedMechanics and Engineering 190 (20-21) (2001) 2555–2580. doi:10.1016/S0045-7825(00)00253-X .[108] J. C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelasticsolids, Computational Mechanics 12 (5) (1993) 277–296. doi:10.1007/BF00372173 .[109] F. Lene, D. Leguillon, Homogenized constitutive law for a partially cohesive composite material, International Journal of Solidsand Structures 18 (5) (1982) 443–458. doi:10.1016/0020-7683(82)90082-8 .[110] S. Shkoller, A. Maewal, G. A. Hegemier, A dispersive continuum model of jointed media, Quarterly of applied mathematics 52 (3)(1994) 481–498.URL [111] H. Murakami, G. A. Hegemier, Development of a nonlinear continuum model for wave propagation in joined media: theory forsingle joint set, Mechanics of Materials 8 (2-3) (1989) 199–218. doi:10.1016/0167-6636(89)90012-4 .[112] J. A. White, Anisotropic damage of rock joints during cyclic loading: Constitutive framework and numerical integration, Interna-tional Journal for Numerical and Analytical Methods in Geomechanics 38 (2014) 1036–1057. doi:10.1002/nag.2247 .[113] R. I. Borja, Plasticity, Springer, 2013.[114] Y. M. Tien, M. C. Kuo, C. H. Juang, An experimental investigation of the failure mechanism of simulated transversely isotropicrocks, International Journal of Rock Mechanics and Mining Sciences 43 (8) (2006) 1163–1181. doi:10.1016/j.ijrmms.2006.03.011 .[115] J. Ambrose, Failure of Anisotropic Shales under Triaxial Stress Conditions, Ph.D. thesis, Imperial College London (2014).[116] K. Mogi, Flow and fracture of rocks under general triaxial compression, in: 4th ISRM Congress. International Society for RockMechanics and Rock Engineering, 1979.[117] M. Kwa´sniewski, Mechanical behaviour of rocks under true triaxial compression conditions-Volumetric strain and dilatancy,Archives of Mining Sciences 52 (3) (2007) 409–435.[118] 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. doi:10.1007/s00603-016-0947-7 .[119] R. Nova, An extended Cam Clay model for soft anisotropic rocks, Computers and Geotechnics 2 (2) (1986) 69–88. doi:10.1016/0266-352X(86)90005-4 .[120] Y. Zhao, S. J. Semnani, Q. Yin, R. I. Borja, On the strength of transversely isotropic rocks, International Journal for Numerical andAnalytical Methods in Geomechanics 42 (2018) 1917–1934. doi:10.1002/nag.2809 ..