Anisotropic subloading surface Cam‐clay plasticity model with rotational hardening: Deformation gradient‐based formulation for finite strain

This study is aimed at developing an anisotropic elastoplastic constitutive model for geomaterials at finite strain and its stress calculation algorithm based on the fully implicit return‐mapping scheme. The Cam‐clay plasticity model is adopted as a specific prototype model of geomaterials. As a pertinent representation of deformation‐induced anisotropy in geomaterials, nonlinear rotational hardening is incorporated into the model in a theoretically reasonable manner by introducing the dual multiplicative decompositions of the deformation gradient tensor. In addition to the usual decomposition into elastic and plastic parts, the plastic part is decomposed further into a part contributing to the rotational hardening and a remainder part. The former part leads to a back stress ratio tensor related to the rotational hardening via a hyperelastic‐type hardening rule. The constitutive theory is thereby formulated on proper intermediate configurations entirely in terms of deformation‐like tensorial variables possessing invariance property, without resort to any objective rates of stress or stress‐like variables. Combining the Cam‐clay plasticity with the concept of subloading surface, a class of unconventional plasticity, enables the model to be capable of reproducing complex hardening/softening accompanied by volumetric contractive/dilative responses. Basic characteristics and predictive capability of the proposed model, as well as the accuracy of the developed numerical scheme, are verified through several numerical examples including monotonic and cyclic loadings.


INTRODUCTION
The general framework of computational plasticity has been established by extensive research efforts since the pioneering research works from the mid 1980s to the early 1990s by Simo, 1-3 Simo and Taylor, 4 Simo and Ortiz 5 , Ortiz and Simo. 6 We refer to the comprehensive monographs by Simo, 7 Simo and Hughes 8 for thorough description. Much progress in this field has been made to develop theoretically sound formulation for constitutive modeling as well as to devise efficient implementation into numerical methods for engineering analysis, especially in finite strain elastoplasticity. [9][10][11] In this context, it would be no exaggeration to say that the theory of finite strain elastoplasticity based on the multiplicative decomposition of deformation gradient tensor with the use of hyperelastic stress-strain relation is the achievement of utmost importance in this research area, although the notion of multiplicative decomposition itself goes back to early works by Kröner, 12 Lee, 13 and Mandel. 14 In view of the historical progress, prior to the emergence of the multiplicative framework, the hypoelastic-based plasticity framework 15,16 were widely employed, in which an additive decomposition of the rate of deformation tensor into elastic and plastic parts is assumed, and the hypoelastic equation in rate form with some proper objective or corotational stress rate tensor is adopted for elastic stress-strain relation. It has been pointed out, however, that the hypoelastic-based framework has several serious drawbacks in theoretical and numerical aspects, and thereby associated controversial issues have been arisen, namely, the arbitrariness or non-uniqueness regarding the choice of an objective rate tensor, 17 pathological model responses such as the oscillatory or unstable behavior of shear stress and/or back stress related to kinematic hardening under large simple shear deformation, [18][19][20][21][22][23] energy dissipation and its accumulation even within purely elastic range during cyclic loadings, 24,25 and the lack of "incremental objectivity" 23 when using improper time-integration schemes. Due to the last mentioned property, a careful treatment is required in the numerical implementation of hypoelastic-based plasticity models to ensure the incremental objectivity. [26][27][28][29][30][31] In contrast, the hyperelastic-based plasticity relying upon the multiplicative decomposition completely bypasses the above mentioned problems encountered in the hypoelastic-based approach. The fundamental requirement of material frame-indifference and the energy conserving property in the elastic range are fulfilled automatically by the intrinsic natures of the multiplicative framework and the hyperelasticity. Nowadays, the hyperelastic-based approach has won the widespread acceptance, and has been adopted for modeling of metallic materials; see for example, the coupling with continuum damage. 32,33 With regard to incorporation of plastic anisotropy into a model, there are mainly two different approaches to formulate the kinematic hardening at finite strain: One approach is a direct extension of the evolution law of the back stress tensor with the use of an objective rate tensor; see Tsakmakis, 34 and Svendsen et al. 35 The other approach is based on a further decomposition of the plastic part of deformation gradient tensor; see for example, Lion, 36 Wallin et al., 37 Wallin and Ristinmaa, 38 Shutov and Kreißig, 39 Henann and Anand, 40 and Vladimirov et al., 41,42 in which the kinematic hardening is described in terms of tensorial internal variables of strain type, and the back stress tensor is derived via a hyperelastic-type relation. In this way, one obtains a thermodynamically consistent formulation in a physically wellbased manner as a multidimensional extension of one-dimensional rheological model. Dettmer and Reese 43 presented a numerical comparison of the above-mentioned different approaches for nonlinear kinematic hardening at finite strain.
In the context described above, Brepols et al. 25 performed a numerical comparison of the inherently different formulations of finite strain elastoplasticity based on the hypoelastic-and hyperelastic-based approaches. In their article, they also presented a detailed review and numerical assessment of the objective integration algorithm for the hypoelastic-based plasticity, as well as of the return-mapping algorithm for the hyperelastic-based plasticity. We also refer to comparative studies on the different finite strain frameworks by Idesman, 44 and Shutov and Ihlemann. 45 As reviewed above, extensive research works on this subject have been done in both theoretical and numerical aspects, and thereby the constitutive framework for finite strain elastoplasticity seems to have already been well established, especially for metallic materials. As for geomaterials, indeed the multiplicative framework has been adopted by many researchers, for example, Simo and Meschke, 46 Borja and Tamagnini, 47 Callari et al., 48 Meschke and Liu, 49 Borja, 50 Ortiz and Pandolfi, 51 Andrade and Borja, 52 Borja and Andrade, 53 and Yamakawa et al., 54,55 but most of their studies are limited to isotropic elastoplasticity, due to the fact that their formulations are mainly based on principal logarithmic strains or two invariants of logarithmic volumetric and deviatoric strains. As far as the authors' knowledge, there are still not many studies that have addressed the constitutive formulation for geomaterials at finite strain within the multiplicative framework amenable to elastic and/or plastic anisotropy involving arbitrary large deformations. There are very early studies on the relevant topic by Jeremić et al., 56 and Rouainia and Muir Wood. 57 Their works, however, dealt with only isotropic models despite that their formulations are applicable to anisotropy. Thereafter, Mosler and Bruhns 58 presented a variational framework for non-associative Drucker-Prager plasticity combined with Hill-type anisotropy, together with volumetric-deviatoric split of elastic and plastic deformations, in which they identified the Mandel stress tensor as a pertinent stress measure for thermodynamically consistent formulation. Gajo 59 presented a finite strain model of saturated porous media with compressible solid and fluid constituents in a Lagrangian description for a volume element moving with the solid skeleton. Recently, Bennett et al. 60 developed a thermodynamically consistent framework to formulate the modified Cam-clay and Drucker-Prager cap plasticity models. During the course of the model development, they also discussed the choice of suitable stress measure in the constitutive modeling of materials exhibiting large plastic volume changes. We also refer to a more recent study by Bennett et al. 61 addressing the extension to include both elastic and plastic anisotropies. They dealt with only anisotropy represented in terms of structure (fabric) tensors, but not plastically induced anisotropy, although they stated in their paper that the provided framework includes the possibility of developing material specific constitutive equations inclusive of induced anisotropy.
Constitutive modeling with anisotropic plasticity for geomaterials has been and is still a very active research area. We refer to several recent research works: anisotropic viscoplastic Cam-clay model for shale by Borja et al., 62 adaptive constitutive modeling with rotational hardening by Dejaloud and Rezania, 63 anisotropic model considering fabric evolution with the aid of DEM numerical tests by Wang et al., 64 and anisotropic model applicable to change in principal stresses and their axes rotation by Xue et al. 65 These studies are; however, limited to the small strain framework.
Motivated by the status of past related researches mentioned above, the objective of this study is aimed at developing a novel Cam-clay plasticity model combined with an induced anisotropy at finite strain within the framework of multiplicative elastoplasticity. As a pertinent expression of the plastically induced anisotropy in geomaterials, 66 the rotational hardening [67][68][69][70][71][72][73][74] is incorporated into the model. To this end, according to the proposal by Lion, 36 we adopt the dual multiplicative decompositions of the deformation gradient tensor; the rotational hardening is thereby described in terms of a strain-type tensorial internal variable, and then the back stress ratio tensor is given via a hyperelastic-type hardening rule defined on a proper intermediate configuration. In this way, the rotational hardening can be incorporated into the model in a well-grounded manner ensuring the frame-indifference property without resort to any objective rates of the back stress ratio tensor.
In order to further improve the capability of the model, particularly for cyclic loadings, we introduce the concept of subloading surface, a class of unconventional plasticity, 75 originally proposed by Hashiguchi and Ueno, 76 and subsequently improved and extended by Hashiguchi. 77 We refer to the comprehensive monographs by Hashiguchi. [78][79][80] The essence of the subloading surface concept is the hypothesis that plastic deformations can occur even at a state of stress below the yield stress (in other words, at a stress point inside the yield surface). Combining the Cam-clay plasticity with the subloading surface concept enables the model to be applicable to a broad class of geomaterials, such as normally and over-consolidated clays, as well as loose and dense sands. A notable feature is that the Cam-clay plasticity equipped with the subloading surface is capable of simulating a smooth transition from initial hardening to softening accompanied by volumetric contractive/dilative responses, typically observed in over-consolidated clays and dense sands; see the references 81,82 for the assessment of predictive capability of the subloading surface model. The organization of this paper is as follows: Section 2 describes two simple small strain examples. One is a onedimensional elastoplastic model with the nonlinear kinematic hardening based on a rheological model. The other is a formulation of the rotational hardening under axisymmetric condition. Analogy of evolution equations appearing in these two examples motivates the development in the subsequent section. Section 3 delineates the entire formulation of general three-dimensional constitutive model. Section 4 presents the implicit stress calculation algorithm using return-mapping for the proposed model. Section 5 presents numerical examples to demonstrate the basic characteristics of the proposed model and to verify the performance of the developed stress calculation algorithm. Finally, Section 6 summarizes the achievement and remaining issues in this study, and mentions potential future developments.
We list notations and symbols used throughout the paper. Bold italic letters denote vectors or second-order tensors. Bold sans-serif letters denote fourth-order tensors. The summation convention is applied over repeated indices. The dot symbol "⋅" denotes the scalar product of two vectors, for example, ⋅ = . A single contraction of adjacent indices of any two second-order tensors is denoted by for example, ( ) = . The colon symbol "∶" denotes the scalar product of two second-order tensors, for example, ∶ = , or the double contraction of two adjacent indices of second-or higherorder tensors, for example, ( ∶ ) = . The symbol "⊗" denotes the dyadic product, for example, ( ⊗ ) = for any two vectors, and ( ⊗ ) = for any two second-order tensors. The symbols "⊗" and "⊗" denote the tensor products of two second-order tensors, for example, ( ⊗ ) = and ( ⊗ ) = . The second-order identity tensor is denoted by , whose Cartesian components are given by the Kronecker delta, that is, ( ) = . The fourth-order identity tensor , the transposing tensor , and the tracing identity tensor are defined by ∶= ⊗ , ∶= ⊗ , and ∶= ⊗ , respectively, denoting the operations ∶ = , ∶ = T , and ∶ = (tr ) for any second-order tensor , F I G U R E 1 Rheological model for Armstrong-Frederick kinematic hardening in metal plasticity where the superscript "T" means the transpose, for example, ( T ) = , and tr = ∶ = is the trace of a secondorder tensor . The symbol ‖ • ‖ denotes the magnitude of a vector or a tensor, for example, ‖ ‖ = ( ⋅ ) 1∕2 for a vector , and ‖ ‖ = ( ∶ ) 1∕2 for a second-order tensor . The superposed dot designates the material time derivative. The subscript "dev" means the deviatoric part of a second-order tensor, for example, dev = − (tr ) ∕3. The superscripts "e" and "p" denote an elastic part and a plastic part, respectively. For simple notation we simply write the inverse and/or transpose of a tensor with the superscript "e" or "p" as for example, ( e ) −1 = e−1 , ( e ) T = eT , and ( where the parentheses are omitted. Remark 1. In the definition of stress and deformation variables, tension is positive and compression is negative as is customary in the continuum mechanics, although the subject of this study is constitutive modeling for geomaterials. In this study on a one-phase medium of the assemblage of soil particles, we shall hereinafter consider the effective (intergranular) stress acting on a soil skeleton.

Motivation: One-dimensional rheological model of nonlinear kinematic hardening
We begin with a one-dimensional constitutive model for small strains with the Armstrong-Frederick type nonlinear kinematic hardening 83,84 based on the rheological model shown in Figure 1. This type of kinematic hardening, together with its interpretation based on the rheological model presented herein, is widely accepted as a suitable plasticity model to capture the characteristic behavior of anisotropic hardening in metallic materials, namely, the Bauschinger effect. 39,[41][42][43] Our discussion here is limited to the rate-independent plasticity, although the rate-dependent viscoplasticity with nonlinear kinematic hardening can also be discussed based on the similar rheological model. We refer to Lion 36 for the development of one-dimensional thermo-viscoplasticity based on a rheological model. The rheological model is constituted by mechanical devices arranged as shown in Figure 1. The whole system of the model has unit length and unit cross-sectional area. The system is a series connection of an elastic spring with stiffness and an inelastic device. The inelastic device consists of a frictional slider with maximum resistance y (≥ 0), and components for plastic hardening, which are connected in parallel. The subscript "y" in y implies the yield stress. In the components for plastic hardening, a dashpot with pseudo-viscosity ∕(̇) and a spring with stiffness are connected serially. This spring is related to kinematic hardening, and so it is called "hardening spring." is a dimensionless constant, anḋ(≥ 0) denotes the magnitude of a slip rate, usually termed "plastic multiplier" in the plasticity theory. Let a stress be applied on the device, and let be a change in the length of device (i.e., strain). The commonly accepted kinematic assumption of the additive decomposition of the strain into elastic and plastic parts is adopted.
where , e , and p are the total strain and its elastic and plastic parts, respectively. The plastic strain is further decomposed additively as where p s is a plastic energy storage part contributing to the kinematic hardening, and p d is a dissipative part. The quantity p d plays a role of internal variable in the model. Differentiating Equations (1) and (2) with respect to time gives the additive decomposition of strain rateṡ=̇e +̇p,̇p =̇p s +̇p d .
The approach based on the rheological model with the idea of the decomposition of plastic strain gives a clear mechanical interpretation of the kinematic hardening, and enables the incorporation of nonlinear kinematic hardening into an elastoplastic model in a physically well-grounded manner without direct use of a rate evolution law for the back-stress as an internal variable. As observed in Figure 1, the stress in the elastic spring and the stress in the hardening spring are given by the following two constitutive equations , .
. F I G U R E 2 Nonlinear kinematic hardening in one-dimensional problem Differentiating Equation (5) with respect to time, and then substituting Equations (7) and (8) lead tȯ where the upper sign corresponds to a positive (tensile) plastic strain rate (̇p > 0), and the lower sign to a negative (compressive) one (̇p < 0). Equation (9) can be rewritten aṡ which is none other than the one-dimensional expression of the well-known evolution equation for the back-stress in the Armstrong-Frederick type nonlinear kinematic hardening model. For further understanding of the hardening characteristic of this model, we consider a monotonic loading process within a time interval [ 0 , ]. Integrating the rate equation (9) over the time interval with an initial condition p = p 0 and = 0 at = 0 leads to with its derivative where we should note that double-sign corresponds in the same order. From Equation (11), the saturation values of the back-stress are The behavior of nonlinear kinematic hardening presented above is illustrated in Figure 2. It should also be noted that, when = 0, the model reduces to Prager's linear kinematic hardening. The notion of the decomposition of inelastic strains facilitates the general three-dimensional constitutive formulation for finite strains combined with kinematic hardening and/or other inelastic mechanisms; it is extended to the multiplicative decomposition of inelastic parts of the deformation gradient tensor. Such approaches have widely been adopted for finite strain formulations of thermo-viscoplastic model, 36 rate-dependent/independent elastoplastic models with nonlinear kinematic hardening, 25,39-43 as well as elastoplastic models coupled to continuum damage. 32

Rotational hardening under axisymmetric condition: Small strain case
The rotational hardening has widely been accepted as a pertinent plastic hardening model to describe induced anisotropy observed in geomaterials. [66][67][68][69][70][71][72][73][74] Figure 3 illustrates a conceptual diagram of the rotational hardening, in which the modified stress tensor with the back stress ratio tensor represented by a deviatoric tensor (i.e., tr = 0), is typically used. We should note here that the stress tensor , as well as its mean normal stress , is defined as positive in tension. Hashiguchi 79 proposed the evolution rule of in the form applicable to finite strains: where • is some proper objective rate of , p dev is the deviatoric part of the plastic deformation rate tensor, r is the evolution coefficient, and r is the constant specifying the stress ratio at a rotation limit. Equation (15) has a form different from the one originally proposed in the earlier studies by him and his co-workers, 68,69,72,73 but it is an improved version consistent with the development presented herein. The small strain counterpart of Equation (15) can be expressed aṡ To simplify the following discussion, we shall consider an axisymmetric condition, in which the matrix representations of stress and strain tensors in terms of Cartesian components are given, respectively, by where the subscripts ' ' and ' ' mean the axial and lateral components, respectively. Here, we assume the axial compression or extension in 1 -direction under axisymmetric stress state with equal lateral confining stresses in 2 -and 3directions. In this case, the mean normal stress and the deviatoric stress invariant are given by The matrix representation of the deviatoric stress tensor dev and its axial and lateral components are Similar relations to Equations (18) and (19) hold for the modified stress tensor mod defined in Equation (14). The volumetric strain v and the deviatoric strain invariant s are given by   .
We should also note that the last equation in Equation (26) defines the evolution rule of The key equations of the nonlinear kinematic hardening and the rotational hardening described above are summarized in Table 1, from which one can observe a similarity in the mathematical structure between these two hardening models. This similarity was discussed also in Hashiguchi 69 based on the rate equations (10) and (15). Although the kinematic hardening model is applied to metallic materials in general, this formal similarity motivates to develop a re-formulation of the rotational hardening in an analogous manner.

CAM-CLAY PLASTICITY COMBINED WITH ROTATIONAL HARDENING AND CONCEPT OF SUBLOADING SURFACE
This section presents a general three-dimensional constitutive formulation of the Cam-clay plasticity model at finite strain combined with the rotational hardening and the subloading surface concept. The present formulation is entirely based on the multiplicative framework. We employ a representation of the plasticity model in terms of the three stress invariants in order to take account of the Lode angle dependencies of both yielding and plastic flow, that is, the effect of intermediate principal stress, commonly observed in geomaterials.

Basic kinematics and deformation and stress tensors in finite strain elastoplasticity
This section briefly summarizes the basic framework of the finite strain elastoplasticity as a preliminary for the constitutive formulation in the subsequent sections. The basic kinematic setting and the definitions of deformation and stress tensors used in the formulation of the present constitutive model are described in this section. A comprehensive description of the deformation and stress variables used in the finite strain elastoplasticity is referred to Hashiguchi and Yamakawa. 85 We employ the widely-accepted kinematic assumption for the finite strain elastoplasticity based on the multiplicative elastic-plastic decomposition of the deformation gradient tensor, 12,13 namely, Following the proposal by Lion, 36 we introduce an additional kinematic assumption that the plastic part p may further be decomposed multiplicatively as where p s denotes a part contributing to the rotational hardening and p d denotes a remainder part. According to the decomposition of Equation (31), a local intermediate configuration at each material point is obtained by releasing the rotational hardening; namely, p d = p s −1 p . Since e = p−1 and p s = p p d −1 , the changes in the stress and a stress ratio variable for rotational hardening in a deformation-driven problem are governed by evolution equations of p and p d . A notable feature of the above-described theory for finite strains based on the multiplicative factorization is that, p , and hence, its factors p s and p d , are invariant under a change in reference frame, as expatiated by Henann and Anand. 40 Therefore, the fundamental requirement of frame-indifference property of the evolution equations for p and p d , and also the constitutive relation for the rotational hardening variable formulated in terms of p s , is automatically satisfied. There has been a long-debated issue concerning the indeterminacy in choosing the intermediate configuration in the multiplicative decomposition (see the references [86][87][88]  Remark 2. The nomenclature "rotational hardening" is not attributed to any physical origin of deformation-induced anisotropy but to a constitutive modeling of anisotropy in frictional materials in such a way that the yield surface rotates around a specific point in the stress space. This is in contrast to the kinematic hardening for metallic materials, in which translation of the yield surface is considered. In cohesionless (purely frictional) materials, the rotation center is usually taken as the origin of the stress space.

Deformation and velocity gradient tensors
In the constitutive formulation in this study, tensorial variables defined on the configurations  0 ,, and  are denoted, respectively, by uppercase letters, overtilded uppercase letters, overlined uppercase letters, and lowercase letters, unless otherwise stated.
The right Cauchy-Green deformation tensor and its plastic and elastic counterparts are given, respectively, by where and p are defined on the reference configuration  0 , whilēe on the intermediate configuration. In a similar manner, the right Cauchy-Green-type deformation tensors related to p d and p s are defined, respectively, by where p d and̃p s are defined on the configurations  0 and, respectively. Taking time-derivative of the both sides of Equation (30), and then multiplying by −1 from the right, we obtain the additive decomposition of the spatial velocity gradient tensor into two terms, namely, in which the velocity gradient tensor and its elastic part e defined on  are given by while we define the plastic velocity gradient tensor on as The second term in the right-hand side of Equation (34) signifies the push-forward operation of̄p from toward .
In a similar way, the plastic velocity gradient tensor on related to p d is defined bỹ

Stress tensors and their invariants
The pull-back operation for the spatial Kirchhoff stress tensor defined on  using e derives the stress tensors defined on, namely,̄∶ = e−1 e−T ,̄∶= eT e−T =̄ē, wherēis the second Piola-Kirchhoff-type stress tensor, and̄is the Mandel stress tensor; they are both defined on. The Mandel stress tensor̄is work-conjugate to the plastic velocity gradient tensor̄p with respect to the plastic work rate; hence, in this study these measures of stress and plastic strain rate are used in the formulation of plasticity based on the description on the intermediate configuration 1 . Note that̄is symmetric by definition, whereas̄is not symmetric in general except the case of elastic isotropy. In addition to the standard stress tensors defined above, we originally introduce the second Piola-Kirchhoff-type (dimensionless) back stress ratio tensor̃r ot defined on related to the rotational hardening. Similar to Equation (38) 2 , the Mandel-type back stress ratio tensor is given bỹr ot ∶=̃p s̃r ot .
This dimensionless tensor̃r ot is a finite strain counterpart of the back stress ratio tensor appearing in Equation (14) in the small strain theory, and will be used in Section 3.1.3 to define the modified stress tensor for rotational hardening via the push-forward operation by Equation (40) 2 . A specific function relating̃p s and̃r ot , as well as̃r ot , which characterizes the rotational hardening, will be presented in Section 3.4. Note that̃r ot is intrinsically symmetric by definition, whereas̃r ot is not symmetric except the limited case when the two symmetric tensors̃p s and̃r ot commute. To be more specific, this is the case wheñr ot is given via an isotropic function of̃p s . In this case,̃r ot and̃p s are coaxial, and thus commute. This is the case in the rotational hardening rule described in Section 3.4.
The push-forward operations with the use of p s converting̃r ot and̃r ot from to derive the back stress ratio tensors defined on, namely, rot ∶= p s̃r ot p s T ,̄r ot ∶= p s −T̃rot p s T =̄p s̄r ot =̄r ot , (40) in which we have used the relation 85̄p s ∶= p s −T̃p s p s −1 = .
The pull-back operations for the above-defined stress tensors and the back stress ratio tensors derive the ones defined on  0 , namely, and where is the second Piola-Kirchhoff stress tensor, and is the (generally non-symmetric) Mandel-type material stress tensor; rot and rot are the back stress ratio tensors in the material description.
Let I , II , and III denote the first, second, and third principal invariants of an arbitrary (generally non-symmetric) second-order tensor , respectively, and let them be defined as where tr , cof , and det denote the trace, cofactor, and determinant of , respectively. The principal invariants of the deviatoric part of , defined by dev ∶= − (1∕3)(tr ) , are thus given by Based on the above definitions, we introduce the three stress invariants of the Mandel stress tensor̄defined on: for IĪd ev ≠ 0, with̄d ev ∶=̄−̄denoting the deviatoric part. In Equations (49)-(51),̄is the mean normal stress,̄is the deviatoric invariant, andΘ is the Lode angle 2 with respect tō. In Equation (51) we take a principal branch of the arccosine function, so thatΘ takes values in the range of 0 ≤Θ ≤ ∕3. The values ofΘ, 0 and ∕3, correspond to the stress states of axisymmetric (triaxial) extension and compression, respectively. In this study, the Lode angle dependencies in plastic yielding and flow property are introduced to describe their differences in the extension and compression sides typically observed in geomaterials. The first and second derivatives of̄,̄, andΘ with respect tōare required in the following constitutive formulation as well as the numerical algorithm based on the return mapping scheme. The explicit expression of these derivatives is given in Appendix A.
Remark 3. From the standpoint of plasticity modeling at finite strain, it is worth noting that the three stress invariants̄, , andΘ are identical to those of , as well as to those of . Let , , and denote the three stress invariants of , and be defined by the form similar to Equations (49)-(51) as follows: with dev = − denoting the deviatoric part of . In a similar way, we define the three stress invariants of as with dev = − denoting the deviatoric part of . Noting the identity in regard to the traces, tr = tr̄= tr ( = 1, 2, 3) according to Equations (38) 2 and (41) 2 , we can readily find the identities among the principal invariants, I = Ī= I , II = IĪ= II , and III = IIĪ= III , and thereby the identities of the three stress invariants of ,̄, and hold, that is, It should be noted, however, that, unlike Equation (54), the identity of the stress invariants does not hold among ,̄, and .

Modified stress tensor for rotational hardening
The schematic diagram of the rotational hardening in the small strain case was shown in Figure 3. In order to introduce the concept of rotational hardening into the finite strain plasticity, we define the Mandel-type modified stress tensor for rotational hardening based on the intermediate configuration as where the deviatoric part of the (dimensionless) back stress ratio tensor,̄r ot dev , defined in Equation (40) 2 , is used. t is the plastic tensile limit hydrostatic stress (0 ≤ t ≤ i ); i is the elastic tensile limit hydrostatic stress, which will be used in 2 The square root of (−3 IĪd ev ) contained in Equation (51) takes a real value, since, as will be seen in Remark 3, (−3 IĪd ev ) = (−3 II dev ) = (3∕2) tr( dev ) 2 = (3∕2) ‖ dev ‖ 2 ≥ 0, in which the equal sign applies only when dev =̄d ev = , that is, hydrostatic stress states. F I G U R E 5 Normal-yield and subloading surfaces, as well as the elastic limit boundary surface, of the Cam-clay plasticity with rotational hardening on the -stress plane the formulation of hyperelasticity in Section 3.3. Note that both i and t are assumed to be constant in this study, whereas t may depend on the plastic deformation inducing degradation or development of bonding between granular particles. 72 The second term in the right-hand side of Equation 55 consists of the two factors, namely, the pressure relative to t and the deviatoric part of the back stress ratio tensor̄r ot dev . This term plays a role like the back stress tensor in the kinematic hardening for metal plasticity, while in the rotational hardening it is proportional to the relative pressure. This thus means that the center of rotation in the stress space is at the plastic tensile limit hydrostatic stress̄= t (≥ 0), and the effect of rotational hardening vanishes at this stress state. When t = 0, the center of rotation is specified to be at the origin of the stress space, similar to the rotational hardening models in past studies. 67,68,70,[72][73][74] The pull-back operation for̄m od leads to the modified stress tensor based on the reference configuration  0 , namely, where rot based on  0 is the counterpart of̄r ot in Equation (55), and it is given by the constitutive equation (107). The three stress invariants of̄m od , denoted bȳm od ,̄m od , andΘ mod , as well as those of mod , denoted by mod , mod , and Θ mod , are defined in similar formats to Equations (49)- (51) and Equation (53), respectively. Similar to Equation (54), the identities of the three stress invariants of̄m od and mod hold, that is, where we note the vanishing trace of the deviatoric tensor, tr̄r ot dev = tr rot dev = 0, in the first equation.

Concept of subloading surface
We introduce the subloading surface 76,78-80 inside the yield surface; see Figure 5. Hereafter, the term "yield surface" commonly used in the conventional plasticity is rephrased as the "normal-yield surface" to clarify the distinction from the subloading surface. The following fundamental properties of the subloading surface are postulated: 1. The subloading surface has a shape and orientation similar to the normal-yield surface. In the present model, the center of their similarity is fixed to the stress point = t on the hydrostatic axis 3 .
2. The ratio of the size of the subloading surface to the size of the normal-yield surface is represented by a scalar variable denoted by (0 ≤ ≤ 1). We call it the "normal-yield ratio." 3. A current stress point always lies on the subloading surface, even in elastic unloading processes. In other words, the subloading surface always passes through the current stress point. This is the crucial difference from the conventional plasticity model. Being endowed with this property, the subloading surface model is capable of describing the evolution of plastic strain even when a current stress point is inside the normal-yield surface. 4. The state 0 ≤ ≤ e is the "quasi-elastic state," where the model behaves elastically with no evolution of plastic strain.
The material constant e (0 ≤ e < 1) defines the value of at the elastic limit. When = 0, the subloading surface degenerates to the stress point of similarity center. The state e < < 1 is referred to as the "sub-yield state," where the subloading surface resides within the normal-yield surface. The state = 1 is called the "normal-yield state," where the subloading surface coincides with the normal-yield surface, and thus the model behaves the same as the conventional plasticity model during plastic loading. 5. The value of increases in the sub-yield state as the plastic strain evolves, approaching the unity corresponding to the normal-yield state. During this process, the subloading surface expands toward the normal-yield surface. The degree of increase in during plastic loading obeys an evolution rule, which will be described in Section 3.7. Once the value of reaches the unity, the evolution rule of becomes inactive, and thereby ceases to evolve, keeping its unity value during plastic loading. 6. The value of is decreased by elastic unloading so that the subloading surface passes through the stress point that is determined by the elastic law. In this case the subloading surface shrinks, while the normal-yield surface maintains its size and shape, since elastic unloading does not induce any change in the variables associated to plasticity.
Remark 4. The state > 1, which could be called a "super-yield state," is the one in which the subloading surface is larger than the normal-yield surface, and a current stress point is outside the normal-yield surface. This is an inadmissible state in the concept of subloading surface for rate-independent plasticity, except when viscoplasticity involving overstress state is considered.

Yield function, subloading function, and plastic potential function
We adopt the subloading function and the plastic potential function expressed in terms of the three invariants of the modified stress tensor̄m od , involving anisotropy by the rotational hardening, as and where c denotes the compression yield stress defined by the isotropic hardening rule, which will be described in Section 3.5. The ellipsoidal subloading surface and plastic potential surface are represented by = 0 and = 0, respectively. The normal-yield function yield is given by putting = 1 in Equation (58), and the elastic limit boundary surface function elast is given by putting = e , namely, yield (̄,̄r ot , c ) ∶= (̄,̄r ot , c , 1), elast (̄,̄r ot , c ) ∶= (̄,̄r ot , c , e ).
The normal-yield surface in the stress space is represented by yield = 0, and the elastic limit boundary surface by elast = 0. (Θ mod ) and (Θ mod ) control the aspect ratio of the ellipsoids, while their mechanical meanings are the critical stress ratio and its counterpart in the plastic potential, respectively. strain. It plays a role as the "elastic-core," at which the model exhibits the most elastic response when a stress point lies on it. The extended subloading surface model has been applied to the constitutive modeling of geomaterials with rotational hardening. 68 are taken. When TC = TC and = , the functions and are identical, and thus the plastic flow rule is associative, whereas the non-associative plastic flow rule is often adopted for geomaterials 4 .
It should also be noted, in the denominator of the first term in Equations (58) and (59), that the square of (Θ mod ) and (Θ mod ) is subtracted by (̄r ot ) 2 , wherēr ot is a measure of the degree of rotational hardening defined by the deviatoric invariant of̄r ot , that is,̄r and (0 ≤ ≤ 1) is the coefficient. This subtraction induces distortion, or flattening, of the ellipsoids accompanying their rotation. Figure 6 illustrates the rotation of the ellipsoidal normal-yield and subloading surfaces accompanied by the distortion (flattening). The slope of the rotation axis is represented bȳr ot . When̄r ot = 0 and/or = 0, the shape of the subloading surface is identical to the standard modified Cam-clay ellipsoid without rotation or distortion. When the rotation axis coincides with the critical state line, the ellipsoid flattens completely, and the subloading surface eventually degenerates to a line segment having the same slope as the critical state line. This form of yield function, involving anisotropy pertinent to geomaterials, was originally proposed by Dafalias. 67 He derived the anisotropic yield function in a theoretically natural manner based on the rate equation of plastic work dissipation, which is consistent with the methodology in the critical state soil mechanics. The resulting yield function was a rotated and distorted ellipse representing anisotropy induced by deviatoric plastic deformations 5 .
Remark 5. Another expression of Equations (58) and (59) can be obtained by replacinḡm od ,̄m od , andΘ mod directly with mod , mod , and Θ mod , respectively. Owing to Equation (57), this alternative expression in terms of the invariants of mod based on  0 is physically equivalent to Equations (58) and (59). 4 The non-associativity of the plastic flow rule is sometimes classified into two types. 52 In the present model, the difference between TC and TC concerns the volumetric non-associativity, while the difference between and the deviatoric non-associativity. 5 Dafalias et al. 70 subsequently extended this form of function to a non-associative case. Then, Dafalias and Taiebat 74 generalized this model to the multiaxial stress space and applied to the simulation of experiments in various conditions. It is worth mentioning that they adopted the non-associative form with TC > TC , that is, volumetric non-associativity, to obtain a better fit to experimental data. We should also refer to Zhang et al., 71 in which they independently developed a similar functional form, and adopted it for their successful simulation of cyclic mobility. Remark 6. Equations (58) and (59) take the form of quadratic function of the modified stress tensor via its invariants to represent the ellipsoidal shape in the stress space. They can be rewritten in a physically equivalent form having the order of stress. The subloading function (58) is rewritten as The plastic potential function (59) is also rewritten in the same way. Other alternative expressions having mathematically different but physically equivalent form can of course be available. However, the authors' numerical experience has shown that the above form provides a better convergent property in the iterative solution process for the return-mapping equation.

Hyperelastic model with pressure-dependent bulk and shear moduli
In this study, we propose a new hyperelastic model for geomaterials with variable elastic bulk and shear moduli in the description on the intermediate configuration. Hyperelastic models expressed in terms of variables defined on the intermediate configuration satisfy the principle of material frame-indifference. The pioneering study on hyperelasticity for geomaterials is the one by Houlsby. 89 He pointed out that simplistic elasticity models with variable tangent moduli dependent on stress can result in non-conservative elastic response, and proposed an energy-conserving nonlinear elasticity model for clays with pressure-dependent bulk and shear moduli based on a two-invariant elastic potential function within the small strain framework. Borja et al. 90 assessed the predictive capability of the two-invariant nonlinear elasticity model by comparing with measured data in clays. Borja and Tamagnini, 47 and Callari et al., 48 independently extended Houlsby's small strain nonlinear elasticity model to include finite strains using the principal logarithmic strains. Extensions to make this class of elasticity model to be applicable in tensile stress range were proposed by Tamagnini et al. 91 within the small strain theory, and by Yamakawa et al. 54,55 within the finite strain theory. Niemunis and Cudny 92 performed a comparative study on several elasticity models for clays in view of stress and energy recoverable properties during cyclic deformations, as well as homogeneity degrees of the constitutive functions. Later, Houlsby et al. 93 proposed an improved version of his original hyperelastic model 89 also within the small strain framework, in which an elastic stress-strain relation, together with variable tangent moduli as power functions of the mean normal stress, is derived from an elastic potential function in terms of two invariants of the elastic strain tensor, or alternatively, from a complementary energy function in terms of two invariants of the stress tensor.
The hyperelastic model by Houlsby et al. 93 underlies the theoretical basis of the model proposed herein. We extend their small strain nonlinear elasticity model to include finite strains. Although hyperelastic models taking account of anisotropy have been developed by several recent research works, 94-97 we shall limit the present formulation to isotropic hyperelasticity. Limiting to isotropic elasticity, the model formulation is realized by a two-invariant expression in terms of the elastic logarithmic strain tensor and the appropriate conjugate stress tensor, and thus it is entirely based on the description on the intermediate configuration. Therefore, the proposed hyperelastic model can readily be combined with plasticity models for finite strains based on the multiplicative framework. In addition, a further extension is introduced for the model to be applicable in tensile stress range.
We use the logarithmic elastic strain tensor defined on: where the elastic stretch tensor̄e is given bȳe = (̄e) 1∕2 . It can also be obtained via the right polar decomposition of e , namely, e = ēe , where e is the elastic rotation tensor. The volumetric and deviatoric invariants of̄e log are defined by e v ∶= tr̄e log = log(det̄e) = log e , where e ∶= det̄e = det e is the elastic volume change, and̄e dis ∶= ( e ) −1∕3̄e is the isochoric (or purely distortional) part of̄e, which is a unimodular tensor having unit determinant, det̄e dis = 1.
We define the following form of hyperelastic potential function expressed in terms of the two strain invariants: in which we put The material constants in Equation (66) are the elastic compressibility index 6 * (> 0), the reference elastic shear modulus e 0 (> 0), the elastic tensile limit hydrostatic stress i (≥ t ≥ 0), and the parameter (0 ≤ ≤ 1) controlling the pressure dependency of elastic shear modulus. The reference mean normal stress at a reference elastic volumetric strain e v0 is denoted by 0 (< i ).
The second equation in Equation (66) for = 1 is derived by considering the limit (2− )∕(2−2 ) → exp as → 1 − . Note that = 1, = 0, and thus e ( e v0 , 0) = 0 hold at the reference state ( e v , e s ) = ( e v0 , 0). When 0 ≤ < 1, ≥ 0 holds for any values of e v ∈ (−∞, +∞) and e s ∈ [0, +∞); → +∞ as e v → ±∞ and/or e s → +∞, since 0 < i . When = 1, → −∞ as e v → +∞ and e s ↛ +∞; → +∞ as e v → −∞ and/or e s → +∞. The first derivatives of Equation (66) with respect to e v and e s respectively give the two stress invariants, and , defined in Equation (52) 1,2 , namely, In Equations (68) and (69), one can confirm = 0 and = 0 at the reference state ( e v , e s ) = ( e v0 , 0). The second derivatives of Equation (66) give the hyperelastic tangent moduli in the two-invariant expression. For 0 ≤ < 1, we 6 In the small strain theory, the elastic compressibility index (also commonly called the swelling index) usually means the slope of the unloading/reloading line, or the swelling/recompression line, on the semi-logarithmic -log(− ) plane or the -log(− ) plane, where and denote the void ratio and the specific volume, respectively ( = 1 + ). Note, in contrast, that the elastic compressibility index * in the present formulation for finite strains means the slope of the linear double-logarithmic log -log(− ) relation. [98][99][100][101][102] It is thus different from the one in the small strain theory. Note further that in the present formulation denotes the mean normal Kirchhoff stress, not the mean normal Cauchy stress used in the small strain theory. This issue was mentioned in the references. [47][48][49] have from which we have the two stress invariants One can observe in Equations (74) and (75) that, when = 0, is given only by e v , and only by e s , leading to vanishing offdiagonal components of the hyperelastic tangent moduli matrix in Equation (72). This means no volumetric-deviatoric coupling. Further, the stress-strain relation is linear with constant bulk and shear moduli. Thus, → ±∞ for infinite elastic volumetric expansion/contraction e v → ±∞ (double-sign corresponds), irrespective of the value of i . Next, we consider the case of = 1. In this case, the volumetric and deviatoric responses are fully coupled. From Equations (68) and (69)  Finally, we consider purely volumetric elastic deformations with e s = 0. In this case, no deviatoric stress is induced ( = 0), while the mean normal stress is given by The hyperelastic tangent moduli for 0 ≤ < 1 are from which, taking a limit → 1 − , we have for = 1: One can observe in Equations (77) and (78) that, during purely volumetric processes, the hyperelastic bulk and shear moduli are both proportional to the -th power of the normalized pressure, with their reference values e vv = −( 0 − i )∕ * and e ss = 3 e 0 at = 0 . As the stress approaches the elastic tensile limit, namely when → − i , the moduli approach null values, i.e. e vv → 0 + and e ss → 0 + . It is also worth noting that Equation (76) is consistent with Equation (124) representing the unloading/reloading response in purely volumetric processes. This point will be discussed again later in Section 3.5.1.
According to the standard expression, the hyperelastic constitutive equation in terms of the stress and elastic deformation tensors defined on is given bȳ where we have used Equations (68) and (69), as well as the identity of the stress invariants argued in Equation (54). Also used in the above are the first derivatives of e v and e s with respect tōe given by while the derivation is omitted here. Since the logarithm of a second-order tensor is an isotropic function of the argument tensor, (loḡe) dev and̄e −1 are coaxial and both symmetric, and thus they commute. Hence, e s ∕̄e is symmetric. Using Equation (38) 2 , we can convert Equation (79) to the expression in terms of the Mandel stress tensor defined on: One can note in Equation (82) that̄is symmetric, since the present hyperelastic model is isotropic. Under purely volumetric elastic processes ( e s = 0 and̄= 0), the deviatoric part represented by the second term in the right-hand side of Equation (82) vanishes, and the stress tensor is composed only of the spherical part.
Taking material-time derivative of Equation (79) The second derivatives of e v and e s with respect tōe appearing in the above expression are given by in which we have used the relation tr(loḡe) = log(det̄e) in the derivation. The fourth-order tensor̄e possesses the major and minor symmetry properties, namely,̄ē̄̄̄=̄ē̄̄̄and̄ē̄̄̄=̄ē̄̄̄=̄ē̄̄̄. Detailed description on the numerical computation of a tensor logarithm and its derivative can be referred to the references. [103][104][105][106] For e s = 0, an approximate computation using numerical differentiation technique is utilized in this study, instead of the analytical formulae (81) and (86).
The pull-back operation for Equation (82) using Equation (41) 2 , with the help of Equation (54), leads to the expression on  0 : where we should note that and p are not necessarily coaxial and thus do not commute. Therefore, is not symmetric in general, even if̄in Equation (82) where the derivative of tensor logarithm appearing in the last term is calculated as The two strain invariants e v and e s , as well as their derivatives, can be expressed in terms of and p−1 as follows: and

Rotational hardening rule
Extending the basic idea described in Section 2.2, we define the constitutive equation for rotational hardening as a function of̃p s that gives the back stress ratio tensor̃r ot . The formulation presented herein is made up on the configuration in a hyperelasticity-like format, while the essential assumption is that only plastic deviatoric deformations contribute to the rotational hardening. In other words, we assume that the rotational hardening is not induced by plastic volumetric deformations nor any elastic deformations. We use the logarithmic strain measure defined on: where the stretch tensor related to rotational hardening,̃p s , is given bỹp s = (̃p s ) 1∕2 , or the right polar decomposition p s = p s̃ps . Similar to Equations (64) and (65) where p s ∶= det̃p s = det p s is the volume change related to p s , and̃p s dis ∶= ( p s ) −1∕3̃p s is the isochoric part of̃p s .
As mentioned above, however, the volumetric invariant p s v does not contribute to the rotational hardening. It is thus excluded in the following formulation of rotational hardening, even if p s contains volumetric component.
We define the potential function for the rotational hardening in terms of p s s : where r (≥ 0) is the dimensionless coefficient of rotational hardening, and the parameter (0 ≤ ≤ 1) controls the nonlinearity of rotational hardening. The second equation in Equation (97) At the reference state p s s = 0, Equations (98) and (99) take the values rot = 0 and rot ss = 3 r , respectively, irrespective of the value of . When = 0, Equation (98) in which we have used the identity of the invariants, rot =̃r ot , as well as the following derivative having a similar form to Equation (81): The tensorial factors in Equation (102), (log̃p s ) dev and̃p s −1 , are coaxial and both symmetric, and thus they commute. Hence, p s s ∕̃p s is symmetric, and thereby one can confirm that̃r ot in Equation (101) One can observe that̃r ot is symmetric in the present isotropic model. Further,̃r ot is deviatoric; this is consistent with the fundamental assumption for the rotational hardening stated at the beginning of this section. Remember that Equation (101) as well as Equation (103) is the extended counterpart of Equation (29) in the axisymmetric description within the small strain theory as the multidimensional rotational hardening rule including finite strains. Taking material-time derivative of Equation (101) The pull-back operation for Equation (103) using Equation (42) 2 leads to the expression on  0 : in which the identity of invariants,̃r ot = rot , has been used. Note that p and p d are not necessarily coaxial, and thus do not necessarily commute. Hence, rot is not symmetric in general, even if̃r ot in Equation (103) where the derivative of tensor logarithm appearing in the last term in each of the above two equations are calculated as and The two strain invariants p s v and p s s can be expressed in terms of p and p d as follows: Their derivatives with respect to p and p d are given by p s s These derivatives will be used in the Jacobi matrix of the return-mapping equation expressed in terms of p and p d , which will be presented in Section 4.5.

Isotropic hardening rule based on plastic volumetric strain
The isotropic hardening rule in terms of the plastic volumetric strain and the pre-consolidation pressure is one of the essential ingredients of the Cam-clay plasticity model. This isotropic hardening rule includes softening, that is, the decrease in pre-consolidation pressure, induced by expansive plastic volumetric strain. Another feature to be mentioned is that volumetric expansion/compression due to isotropic unloading/reloading processes is modeled as a nonlinear elasticity with pressure-dependent bulk modulus. Is this study, the isotropic hardening rule is formulated as the linear doublelogarithmic log -log(− ) relation. [98][99][100][101][102] In accordance with the multiplicative decomposition in Equation (30), the volume change ∶= det is factorized to where the elastic and plastic parts are defined by e ∶= det e = (det̄e) 1∕2 and p ∶= det p = (det p ) 1∕2 , respectively. Taking logarithm of both sides of Equation (118) leads to the additive decomposition of the logarithmic volumetric strain with its rate form: where the logarithmic volumetric strain and its plastic part is defined by v ∶= log and p v ∶= log p , respectively, while the elastic part has already been defined in Equation (64) and used in the formulation of hyperelasticity. Their rates are given bẏv =̇∕ ,̇e v =̇e∕ e , anḋp v =̇p∕ p . In the following, we develop the double-logarithmic relation between the specific volume and the mean normal stress . Referring to Figure 7 and noting = ∕ 0 , e = e ∕ 0 , e 0 = e 0 ∕ 0 , and p = ∕ e , we have We define the following form of isotropic hardening rule expressing the compression yield stress c (tension is positive) as a function of p v : where * (> 0) is the plastic compressibility index, and c0 is the reference value of c at p v = 0 ( p = 1). The limit of c at infinite plastic volumetric expansion p v → +∞ ( p → +∞) is designated by the constant c+ ( c0 < c+ ≤ 0), while the limit in compression side is c → −∞ as p v → −∞ ( p → 0 + ). Therefore, c takes values in the range −∞ < c < c+ . Taking time-derivative of Equation (121) leads to the rate form of the isotropic hardening rule: in which −( c − c+ )∕ * (> 0) can be regarded as the instantaneous modulus of isotropic hardening. Using Equation (120) 3 , we can rewrite Equation (121) from which we note that this isotropic hardening rule represents a linear relation between p versus ( c − c+ ) on the double-logarithmic scale. Next, we consider elastic response under purely volumetric deformation. Recalling Equation (76), we have its inverse relation in which ( e v − e v0 ) has been rewritten in terms of the specific volume with the use of Equation (120) 2 . We then note that the second equation for = 1 in Equation (124)  Taking the sum of Equations (123) and (124) gives the total logarithmic volumetric strain as the function of and c , namely, in which the first and second terms in each of the right-hand sides correspond to the elastic and plastic parts, respectively. Let us consider the elastoplastic compression process under normal compression condition. In this case, we can put e 0 = c0 , = c , = c , and 0 = c0 in Equation (125). For simplicity of the discussion, we shall consider a specific case of = 1. The left-hand side of Equation (125) can be rewritten in terms of the specific volume, and we then have for = 1: where the first and second terms in both sides correspond to the elastic and plastic parts, respectively. The above equation implies that the normal compression curve as the c versus c relation on the double-logarithmic scale is nonlinear in general even when = 1, while it is linear only when i = c+ = 0 (Note that i ≥ 0 and c+ ≤ 0 by their definitions). In this special case, the plastic part of Equation (126) where we have defined the elastoplastic compressibility index * ∶= * + * > 0 which signifies the slope of the normal compression line (see Figure 7), satisfying * > * , since * > 0. The inverse relation of Equation (127)

Extended isotropic hardening rule dependent on deviatoric plastic strain and stress ratio
The isotropic hardening rule described in the previous section is dependent only on the plastic volumetric strain p v . In this section, we extend the isotropic hardening rule by incorporating the influence of the plastic deviatoric strain p s in order to improve the capability of the model. Here, we call it "deviatoric hardening." Let Equation (122) be extended tȯ where the second term in the right-hand side has been added so as to include the effect of the plastic deviatoric strain p s . The dimensionless coefficient ℎ s denotes the deviatoric hardening modulus, whose definition is where ℎ s1 (≥ 0), (≥ 0), (> 1), (−∞ < < +∞), and (≥ 0) are material constants. Equation (130)  .
The above definition is applicable to the range −∞ < < t including tensile stresses; takes a value within the range 0 ≤ < +∞. d (Θ) denotes the boundary stress ratio of the deviatoric hardening and softening; the Lode angle dependency is taken into account by d (Θ) ∶= d,TC ∕ d (Θ) with the boundary stress ratio in the axisymmetric (triaxial) compression side, d,TC , and the function d (Θ) ∶= (Θ; d ) describing the Lode angle dependency with a parameter d . The specific form of the function (Θ; d ) will be described later.
Here, we assume that the model exhibits deviatoric hardening when the stress ratio ∕(− + t ) is higher than d (Θ), that is, > 1, while deviatoric softening when the stress ratio ∕(− + t ) is lower than d (Θ), that is, < 1. No deviatoric hardening or softening occurs on the boundary, namely, ℎ s = 0, when ∕(− + t ) = d (Θ), that is, = 1. This idea is illustrated on the -stress plane in Figure 8 (a). In the principal stress space, the boundary of the deviatoric hardening/softening takes a form of conical surface having its apex at = t on the hydrostatic axis, while its cross-sectional shape on the deviatoric stress plane is a distorted circle due to the Lode angle dependency.
The deviatoric hardening/softening property dependent on the stress ratio assumed above is modeled by the factor containing in Equation (130). Here, a hyperbolic form is employed to model nonlinear hardening/softening property. Figure 8 (B) shows several plots of the deviatoric hardening modulus ℎ s as a function of for different values of the index . When ℎ s1 > 0, the sign of ℎ s varies depending on the value of as follows: In addition, the influence of and p s are considered in Equation (130). When ℎ s1 > 0, > 0, and > 0, the limits of ℎ s are as follows: ℎ s → +∞ as → −∞; ℎ s → 0 as → − t and/or p s → +∞. It should be mentioned here that, unlike Equation (122) leading to Equation (121), the extended version of the isotropic hardening rule is given by the rate form (129), and it is not analytically integrable except when ℎ s is constant. Hence, a numerical time-integration scheme is required to compute an updated value of c . This will be described in Section 4.2.

Plastic evolution rule
The plastic evolution rule in the present model is composed of the plastic flow rule and the evolution rule for rotational hardening variable. This section describes these two evolution rules. The former specifies the evolution of p , which directly determines the plastic deformation. The latter specifies the evolution of p d ; then, p s is determined by p s = p p d −1 , which is used to calculate the back stress ratio tensor̃r ot via the rotational hardening rule described in Section 3.4.

3.6.1
Plastic flow rule According to Equation (36), the evolution equation for the plastic deformation gradient tensor p is given bẏ in which̄p prescribes the evolution of p . The plastic part of the velocity gradient tensor̄p based on the intermediate configuration can be decomposed additively as follows:̄p =̄p +̄p,̄p ∶= sym̄p,̄p ∶= skew̄p, wherēp denotes the symmetric part of̄p, called the plastic rate of deformation tensor on, and̄p the skew-symmetric part, named the plastic material spin tensor on after Dafalias. 88 We define the plastic flow rule for the symmetric part̄p in the non-associative form using the function , namely, wherė≥ 0 is the plastic multiplier. The direction of plastic flow is given bȳp, denoting the symmetric part of the outward normal to the plastic potential surface at a current stress point. The symmetric plastic flow direction tensor p is normalized so that it has a unit magnitude satisfying ‖̄p‖ = 1; hence, the magnitude of plastic flow is given by ‖̄p‖ =̇. Figure 9 illustrates the flow direction at typical stress points on the rotated plastic potential surface. As for the determination of the intermediate configuration, we opt for the simplest choice of a spinless intermediate configuration by specifying the skew-symmetric part as̄p = . (136) The pull-back operation for Equation (135), with the help of the formulȧp = 2 pT̄p p , leads to the expression on the reference configuration  0 : This form will be used in the stress calculation algorithm described in Section 4.
Remark 7. In the case of isotropic hyperelasticity,̄and̄e are coaxial, and thus̄is symmetric. It follows that ∕ī n Equation (135) is symmetric. In this case, the following equality holds (the proof is omitted but straightforward using a pull-back operation): which implies that both ( ∕ ) T p and p ( ∕ ), appearing in Equation (137) as the argument of symmetrization operator, are symmetric, although ∕ is not symmetric or coaxial with p , and therefore they do not commute in general.

Evolution rule for rotational hardening variable
According to Equation (37), the evolution equation foṙp d is given bẏ The velocity gradient tensor̃p d on related to the rotational hardening can be decomposed additively into the symmetric and skew-symmetric parts, namely, where the latter equation specifies that the configuration is spinless, similar to Equation (136) for the intermediate configuration.
The back stress ratio tensor̃r ot dev on is given by the hyperelastic-type constitutive equation (103), which intrinsically provides a deviatoric tensor. The limit of rotational hardening is specified by the stress ratio at a rotation limit, defined by r (Θ) ∶= r,TC ∕ r (Θ) = r,TC ∕ (Θ; r ) with its value at the axisymmetric (triaxial) compression state, r,TC , and the parameter representing the Lode angle dependency, r . Note that the identity of the Lode angle, like Equation (54) 3 , also applies toΘ on.  . (142) Remark 8. As is the case mentioned in Remark 7, if the constitutive equation for rotational hardening described in Section 3.4 is isotropic,̃r ot dev is symmetric (this is the case in the specific form of Equation (103)). In this case, the product rot dev p d appearing in Equation (142) is symmetric (it can readily be proved via the pull-back of̃r ot dev , while omitted here). This symmetricity holds even though rot dev is not symmetric or coaxial with p d , and thus they do not commute in general.

Evolution rule of normal-yield ratio for subloading surface
The normal-yield ratio plays a crucial role in the subloading surface concept; hence, the evolution rule of is of the utmost importance in the constitutive formulation incorporating the subloading surface. Based on the fundamental assumption that the normal-yield ratio evolves along with the progress of plastic strain, we define the evolution rule of in the following form [78][79][80]85 :̇=̇( where the plastic multiplieṙrepresents the magnitude of plastic strain rate, and ( ) is the function that controls the degree of evolution satisfying ( ) > 0, ′ ( ) < 0 ( e < < 1: sub-yield state); ( ) → 0 + ( → 1 − : approaching to normal-yield state); and (1) = 0 ( = 1: normal-yield state).
The material constant e (0 ≤ e < 1) denotes the value of normal-yield ratio at an upper limit of quasi-elastic state. Figure 10 illustrates a typical behavior of the function ( ). When the material is in a state within the range 0 ≤ ≤ e , the right-hand side of Equation (143) must take an indeterminate form with ( ) → +∞ anḋ= 0 so thaṫin the left-hand side takes a finite value. This means that the material exhibits purely elastic response when 0 ≤ ≤ e . When → 1 − , we have ( ) → 0 + implying the convergence of , which means that the material approaches a normal-yield state taking = 1. We here employ the following specific form of ( ) which satisfies the above-stated function property: with two material constants (> 0) and e . The evolution coefficient regulates the degree of evolution of with plastic strain rate, and thereby regulates the smoothness of the transition from elastic to plastic state. The angle brackets ⟨•⟩ denote the ramp function of a scalar argument , namely, which is also often denoted by the Macaulay brackets { }. According to Hashiguchi and Mase, 73 and Zhang et al., 71 we extend the evolution coefficient in Equation (145) by introducing the functional dependence of on the mean normal stress , the Lode angle Θ mod , and the accumulated logarithmic plastic deviatoric strain p s,accum , in order to improve the model capability to describe cyclic behavior. This extension leads to 7 in which we define 0 (> 0) is the reference evolution coefficient; ∞ (0 ≤ ∞ < 0 ) regulates the lower limit of ; and , , and are nonnegative material constants denoting the effects of , Θ mod , and which mean that the value of function becomes smaller as the pressure becomes smaller and/or the deviatoric plastic strain accumulates, leading to the stagnation of evolution of . In addition, the dependence of on the Lode angle Θ mod describes the difference in the evolution of in the axisymmetric (triaxial) compression and extension sides. More specifically, the critical stress ratio (Θ mod ) is generally larger in the compression side than in the extension side. Therefore, when > 0, the function takes smaller values in the compression side than in the extension side, leading to slower evolution of and inducing larger plastic strains due to plastic loading in the compression side. It should also be noted in Equation (148)   .
The numerical integration scheme for the evolution rule (150) to compute an updated value of will be described in Section 4.3.

Incremental form of plastic evolution rules
The incremental form of the plastic evolution rules is derived in this section by integrating the rate equations, presented in Section 3.6, over a time interval [ , +1 ] corresponding to an incremental step in the stress calculation. This provides the update formulas for deformation variables associated to plasticity. The plastic evolution rules (137) and (142) can be expressed symbolically aṡ where ( ) denotes a symmetric tensor symbolizing p in Equation (137), and p d in Equation (142), while ( ) denotes a generally non-symmetric tensor representing their evolution coefficient tensors. When is a constant tensor, that is, ( ) ≡ , Equation (151) is analytically integrable; integration over the time interval [ , +1 ] exactly gives in which we have put ∶= ( ), +1 ∶= ( +1 ), and Δ ∶= | +1 − |. The tensor exponential function, denoted by , plays a role as the integration operator. Numerical computation of the exponential function of a tensor, as well as its derivative, can be referred to the references. [103][104][105]107 In general, ( ) may vary depending on stress states. We thereupon employ an approximate integration like the backward difference scheme with the use of +1 ∶= ( +1 ), which leads to the implicit update formula 8 Based on the above formula, we obtain for p and p d : where we have put the incremental plastic multiplier Δ ∶=̇Δ .

Integration of isotropic hardening rule
As was previously mentioned in Section 3.5.2, the extended version of the isotropic hardening rule, given by the rate form (129), is not analytically integrable except when ℎ s is constant. Hence, an approximate time-integration scheme is required 8 Alternative schemes, such as the forward difference and the generalized mid-point rule, can of course be employed. If the forward difference scheme is adopted, we obtain the explicit update formula.
to compute an updated value of the compression yield stress c at the time +1 . Integrating Equation (129) from which we obtain the update formula for c : Since, in fact, the deviatoric hardening modulus ℎ s ( , , p s ) varies within the time interval, the above result is not an exact integration but an implicit approximation using ℎ s, +1 ∶= ℎ s ( +1 , +1 , p s, +1 ).

Integration of evolution rule of normal-yield ratio
When the evolution coefficient in Equation (145) is constant, the evolution rule of given by (143) is analytically integrable in the following way. Consider again the time interval [ , +1 ] on which an incremental deformation is given. Suppose that the model exhibits elastic response within a subinterval [ , + ] (0 ≤ ≤ 1), and plastic loading in the remainder [ + , +1 ]; that is, the model transits from elastic to plastic state at + . In this case, the inequality ≤ + = e ≤ +1 holds, and we can put Δ ∶= ∫ +1 +̇d . Equation (143), with the use of Equation (145), can be converted to a form in which the two variables and are separated as follows: which is integrated analytically by leading to the update formula for the normal-yield ratio : Note again that + = e at the transition from elastic to plastic state, and thereby the argument of ramp function ⟨•⟩ in Equation (161) vanishes. Otherwise, when the model exhibits plastic loading over the whole interval [ , +1 ], the inequality e ≤ < +1 holds, and thus we can put = 0, + = , and Δ ∶= ∫ +1̇d in Equation (161). The extended version (150) with , defined by Equation (147) taking account of the dependence on stress and plastic strain, is not integrable in analytical ways. In this case, we can integrate Equation (159) in a semi-analytical manner 9 by replacing the constant coefficient in Equation (161) with +1 ∶= ( +1 , Θ mod +1 , p s,accum, +1 ). The updated value +1 is computed directly by the update formula at every iteration during the solution process for the return-mapping equation.
The accumulated logarithmic plastic deviatoric strain  9 Another way to integrate Equation (150) would be the use of forward or backward Euler difference method. However, the authors' numerical experience has shown that the semi-analytical integration presented herein provides a better accuracy.

Elastic predictor and loading/unloading judgment
Let the stress, strain, and all internal variables associated to plasticity be already computed and completely known at the time . Consider a time interval [ , +1 ] where an incremental deformation is imposed onto the known state at by the relative deformation gradient tensor [ , +1] . The deformation gradient tensor is thereby updated by where [ , +1] corresponding to the incremental displacement Δ at a spatial point is given by The right Cauchy-Green deformation tensor +1 is then computed by The elastic trial state is determined under the assumption that this incremental process is purely elastic without any progress of plastic deformations. This can be done by freezing plastic flow, so that we put Once the elastic trial state is determined, we can make a judgment as to whether the incremental process is elastic or plastic according to the loading/unloading criterion 10 : tr +1 > 0 and tr elast, +1 > 0 ⇝ Plastic step. (167) Here we use the values of the subloading function (58) and the elastic limit boundary surface function (60)   .
If the incremental step is judged to be elastic, the variables at the elastic trial state is the actual solution for the updated state, while the value of needs to be updated so that the subloading surface passes through the updated stress point. This can be accomplished by solving the equation +1 = 0 for +1 . Otherwise, if the step is plastic, the return-mapping described in the subsequent section is performed to compute the variables at the updated state.
Remark 9. It is worth noting that, in contrast to the conventional plasticity, the loading/unloading criterion (167) for the subloading surface plasticity utilizes not only the subloading function , as the counterpart of the (normal-)yield function 10 Stress increments pointing inward of the subloading surface do not always lead to purely elastic unloading. In other words, the subloading surface model may exhibit a reverse plastic loading after an elastic unloading when a large stress increment in such direction is imposed. In such cases, the loading/unloading criterion (167), merely based on the elastic trial state, may make an incorrect judgment. To overcome this drawback, improved methods of loading/unloading judgment have been developed, 108,109 while not yet implemented to the analysis code used in this study. in the conventional plasticity, but also the elastic limit boundary surface function elast . This is due to the essential nature of the subloading surface plasticity stated in Section 3.2.1 (particularly, in the items 3, 4, and 6, therein).

Return-mapping
The return-mapping equation is comprised of the incremental form of the plastic evolution rules (154) and the subloading surface function (58), or its variant form (62), namely, which is complemented by Equations (87), (107), (158), and (161). Given that p and p d are symmetric, the set of simultaneous equations (169a)-(169c) has totally 13 unknowns, that is, six independent components of each of the two tensors p +1 and p d +1 , and the scalar Δ . This system of nonlinear equations can be solved by means of the Newton-Raphson iterative method. Let and denote symbolically the unknown vector and the residual vector, respectively. We then have a symbolic expression of the return-mapping equation: Further, let denote the corrector of in each iteration. Linearization of ( ) with respect to , and then equating its result to zero vector lead to where an entity affixed with the superscript "( )" signifies the one at -th iteration. ∕ is the 13 × 13 Jacobi matrix consisting of the derivatives of Equations (169a)-(169c) with respect to the unknown variables, p +1 , p d +1 , and Δ . Its analytical derivation involves cumbersome manipulation of differentiation, or, alternatively, numerical differentiation using a perturbation technique based on a forward difference approximation can be employed. The procedure to update the unknowns in an iteration is as follows: where the convergence criterion is with "TOL" being the tolerance. In the numerical examples presented in Section 5, we set TOL = 1.0 × 10 −12 .

NUMERICAL EXAMPLES
This section presents several numerical examples to demonstrate the verification of the constitutive model and the stress calculation algorithm developed in this study. The material constants and initial values required in this model are listed in Table 2; their values used in the analysis are given in the succeeding texts.

Normal-yield, subloading, and plastic potential surfaces:
Critical stress ratio, TC Aspect ratio of ellipsoidal plastic potential surface, TC Plastic tensile limit hydrostatic stress, t Flattening coefficient of ellipsoidal yield surface,

Basic characteristics of rotational hardening: Effect of parameters
, and This numerical example demonstrates the basic characteristics of the rotational hardening. Another feature of the present model, namely, the subloading surface, was examined in our past study, 54,85,110 as well as in relevant studies by other researchers, 71 so that we here focus on the rotational hardening. To this end, a series of analyses was performed for various values of the two parameters, r,TC and r , related to the rotational hardening. We analyzed axisymmetric (triaxial) compression, followed by extension, under a constant lateral confining stress = −100 kPa in 2 -and 3 -directions. The axial stress in 1 -direction was computed by controlling (1,1)-component of the deformation gradient tensor, 11 , to impose axial compression and then extension. The analysis was started from an initially isotropic stress state of −100 kPa.
The values of material constants and initial values used in the analysis were as follows: * = 0.0065, e 0 = 20 MPa, = 0.60, i = 0.10 kPa, = 0.60, * = 0.0185, c+ = 0 kPa, TC = TC = 1.2, t = 0.0010 kPa, e = 0.010, 0 = 10, ∞ = 0, = = = 0, 0 = −100 kPa, e v0 = 0, c0 = −100 kPa, and 0 = 1.0 (normally consolidated). The deviatoric hardening/softening was not used, that is, ℎ s1 = 0. The Lode angle dependencies of the subloading and plastic potential functions, as well as of the rotational limit stress ratio r (Θ), were taken into account using Equation (B7). Initial isotropy without rotational hardening (̄r ot = 0) was assumed. Figure 11 demonstrates the effect of r,TC . The figures (A) and (B) compare the results obtained for = 1 and 0 to examine the effect of distortion (flattening) of the subloading surface associated with the rotational hardening. Each of the figures (A) and (B) shows three graphs of the deviatoric stress (the difference between axial and lateral stresses), the volume change, and the degree of rotational hardening versus the axial compressive deformation, in which the analysis results for r,TC = 0.2, 0.6, and 1.0, with a fixed value of r = 10, are plotted. As a reference for comparison, the case of r,TC = 0.0 excluding the rotational hardening is also plotted. In (A), a gentle peak of the deviatoric stress is observed in r,TC = 1.0 and 0.6, which can be associated with slight dilative deformation around and after the peak stress. The larger the specified value of r,TC is, the higher the peak stress is observed, while the difference in the peak stress is not very significant, and the stresses after reaching the critical state are almost identical regardless of the value of r,TC . In (B), in contrast, significant difference in the stress at the critical state depending on r,TC is observed. These results are consistent with the interpretation of the model property based on the schematic diagram in Figure 6. With regard to the plots of̄r ot in Figure 11 (A) and (B), one can observe, in the compression process, the evolution of rot toward the specified limit r,TC and its eventual saturation. In the extension process,̄r ot saturates at a value smaller than r,TC due to the Lode angle dependency, namely, r,TC ≥ r,TE = r,TC ∕ r (0), where r (0) = 9∕7 in this case using Equation (B7). Figure 12 shows the analysis results for r = 0, 5, 10, and 20, with a fixed value of r,TC = 1.0. The case of r = 0 is the analysis without rotational hardening. One can observe commonly in the figures (A) and (B) that the stress increases the more rapidly the larger the value of r is, while the stresses after reaching the critical state are almost identical except the case of r = 0 in (B). These material responses are attributed to the fact that, as can be seen in the plots of̄r ot , the rotational hardening evolves the more rapidly the larger the value of r is, and it eventually saturates at the specified rotation limit regardless of the value of r .

Accuracy assessment: Deformation-controlled analysis with large incremental steps
A series of deformation-controlled analyzes was conducted to assess the accuracy of the developed stress calculation algorithm in various deformation conditions. We imposed volumetric-shearing mixed deformation to the model by specifying the deformation gradient tensor, 54 whose matrix representation in terms of Cartesian components is given by where 0 defines an initial state, and * prescribes a deformation imposed on the initial state. Two control parameters, v and s , specify the volumetric and shear deformations, respectively. * can be divided into step equal increments The accuracy of a computed stress was assessed by the relative error with respect to the Kirchhoff stress tensor: where denotes a solution computed by a prescribed number of increments, step , and ref is a reference solution. We used the value of stress computed by sufficiently small increments specified by step = 1000 as a substitute for the reference solution, since exact solution for the present complex constitutive model cannot be available in analytical way. on the -stress plane, and the error maps obtained in Cases A and B, respectively. The graphs on the left and center compare the stresses computed by step = 10 (plotted as markers) with the corresponding reference solution (solid lines). While slight deviations of the computed stress points from the reference solution are observed in some cases with very large strain increments, the error is not significant as a whole and the analysis results demonstrate enough accuracy and robustness for practical use. It is also worthy to note here that, in Case B for OCR = 5.0, a smooth transition from elastic to plastic responses around the peak stress is well simulated by the subloading surface model. The error maps on the right represent the relative error of the stresses computed at the final incremental step (i.e., 10th step) corresponding to v and s , which also convince us of high accuracy in the present stress calculation even with large strain increments.

Simulation of cyclic mobility observed in experiments
Cyclic mobility is typically observed in sands during undrained cyclic loading. Its accurate prediction by constitutive description necessitates highly sophisticated constitutive models with careful parameter determination. 71,73,[111][112][113][114][115][116][117] This phenomenon leading toward eventual flow liquefaction is essentially concerned with instability problem. [118][119][120][121][122] Setting aside the underlying theoretical basis on this aspect, we shall here focus on constitutive simulation of cyclic mobility observed in experiment in a phenomenological level to calibrate the proposed model. The experimental data simulated in the present analysis were obtained by two different testing methods, namely, cyclic torsional shear test on hollow cylindrical specimen and cyclic triaxial test on circular column specimen of a sand sample under undrained condition (CU test). The testing conditions in total three tests are listed in Table 3. The sand sample in all the tests was Gifu silica sand No. 8 with physical properties: the soil particle density s = 2.623 g∕cm 3 , the minimum dry density d,min = 1.125 g∕cm 3 , and the maximum dry density d,max = 1.564 g∕cm 3 (the corresponding maximum and minimum void ratios were max = 1.332 and min = 0.677). The cyclic shear loads applied to the specimens were sinusoidal wave with a frequency of 0.1 Hz.  The material subroutine developed in this study was implemented into a nonlinear finite element code, FINAL-GEO R , 123,124 developed by Obayashi Corporation, and it was used in the present analysis. The loading conditions corresponding to the two different testing methods were simulated in a single finite element model under an assumption of uniform deformation. Prior to performing the analysis, we determined basic parameters in the Cam-clay plasticity according to experimental results obtained by monotonic loading tests. The critical stress ratio TC was determined by a series of triaxial compression tests (CU test) on the same sand sample with relative densities around r = 85.3% ∼ 85.5%, in which monotonic axial loads were imposed under several different lateral confining pressures. The values of * and * were determined by a normal compression curve and a swelling/recompression curve obtained by a consolidation test on the same sample. The pressure dependency parameter in hyperelasticity was determined by the nonlinearity of the swelling/recompression curve. The values of material constants and initial values used in the analysis were as follows: * = 0.00189, e 0 = 20.3 MPa, = 0.60, i = 0.0010 kPa, r,TC = 1.05, r = 90, = 1.0 × 10 −6 , * = 0.02671, c+ = 0 kPa, d,TC = 0.60, ℎ s1 = 4.0, = 10, = 10, = 0, = 0, TC = TC = 1.431, t = 0.0010 kPa, = 1, e = 1.0 × 10 −6 , 0 = 10, ∞ = 0, = 0, = 3.0, = 180, 0 = −100 kPa, e v0 = 0, c0 = −600 kPa, and 0 = 0.167. The Lode angle dependency was introduced by using Equation (B7). Initial isotropy without rotational hardening (̄r ot = 0) was assumed.
The simulation results are compared with the experimental data in Figure 14, in which the time histories of the shear strain and the excess pore water pressure ratio, the shear stress versus shear strain curves, and the stress paths are shown. The constitutive simulation successfully reproduces the cyclic mobility observed in the experiments, even though the same values of material constants were used throughout to simulate the experimental data measured in three different testing methods and conditions. This realistic simulation demonstrates the predictive capability of the proposed constitutive model, while some deviations from the experimental data are observed, particularly in Test 1 in that the overestimation of the decrease in effective mean normal stress toward the zero pressure during several initial cycles. This can be associated with the overestimation of the buildup of excess pore water pressure due to excessive volume contraction in initial cycles. One can also observe in Test 3 that the different responses in compression side and extension side are fairly well simulated by the Lode angle dependency, while deviations from the experiment are still seen especially in the strain evolution in the extension side.

Application to finite element analysis
This example verifies the applicability to finite element (FE) analysis. A plane strain compression test for a rectangular specimen was simulated by FE analysis. The material subroutine developed in this study was implemented into a research numerical code by the authors for nonlinear FE analysis based on the updated Lagrangian formulation. The algorithmic (consistent) tangent moduli tensor was computed by numerical differentiation using a perturbation technique based on the forward Euler difference scheme 125 to construct the global tangent stiffness matrix. The efficiency of this numerical technique has been reported also for constitutive model of geomaterials. [126][127][128] The analysis model is shown in Figure 15 (A). The specimen with a size of width × height = 10 cm × 20 cm was divided into FE mesh consisting of 10 × 20 eight-node quadratic elements using 2 × 2 reduced Gaussian quadrature points. Initially isotropic stress state of 11 = 22 = 33 = −300 kPa was assumed. The nodal displacements in 1 -direction on the top and bottom faces were fully fixed to simulate frictional end boundaries. The specimen was compressed by imposing a downward displacement in 2 -direction uniformly on the top face under a constant lateral confining pressure of 300 kPa, and the reaction force on the top face was computed. The gravity effect was ignored. The values of material constants and initial values used in the analysis were as follows: * = 0.0025, e 0 = 46 MPa, = 1.0, i = 0 kPa, r,TC = 0.50, r = 20, = 1.0, * = 0.0225, c+ = 0 kPa, TC = TC = 0.98, t = 0 kPa, = 0, e = 1.0 × 10 −4 , 0 = 5.5, ∞ = 0, = = = 0, 0 = −300 kPa, and e v0 = 0. The deviatoric hardening/softening was not used, that is, ℎ s1 = 0. The Lode angle dependency was taken into account only for the normal-yield and subloading functions using Equation (B8) with = 0.753, = 1.0, r = 1.0, and d = 1.0. The following two cases were analyzed to examine the effect of difference in initial state of material: • Case A: Normally consolidated soil (OCR = c0 ∕ 0 = 1.0, c0 = −300 kPa, 0 = 1∕OCR = 1.0); • Case B: Over-consolidated soil (OCR = c0 ∕ 0 = 10.0, c0 = −3000 kPa, 0 = 1∕OCR = 0.10). Figures 15 (B) and (C) plot, respectively, the macroscopic response of nominal axial stress and volumetric strain versus the nominal axial strain. Figure 16 shows the contour of shear strain and the distribution of plastic loading/elastic unloading region at an macroscopic axial strain ∕ = 0.1. The analysis result of Case A exhibits a gentle increase in the axial stress along with the contractive volumetric deformation, while a slight decrease in the axial stress is observed after a peak stress. In contrast, Case B exhibits a steep increase in the axial stress, followed by a remarkable decrease after a sharp peak of stress accompanied by the expansive deformation, which can be associated with the localized deformation observed in  Figure 16 (B). Note that, in Figure 16, the X-shaped localized shear zone is more remarkable in Case B than in Case A, and quite large shear strain to reach around 0.3-0.5 is observed. Figure 17 (A) and (B) depict the convergence profile of the global Newton-Raphson iteration at typical three loading steps before, at, and after the peak stress (marked on the curves in Figure 15 (B) and (C)) in the analyses for Cases A and B. In both cases, nicely convergent profiles with quadratic or superlinear rate are observed even in the loading steps involving unstable macroscopic responses due to material softening and localized deformations. This result assures the performance of the stress calculation algorithm developed in this study.

CONCLUDING REMARKS
In this article, we have presented a novel formulation of anisotropic Cam-clay plasticity model at finite strains based on the multiplicative decomposition of the deformation gradient tensor. The rotational hardening has been incorporated into the model as a pertinent expression of the anisotropy induced by plastic deformation in geomaterials in a reasonable manner consistent with the theoretically sound structure of finite strain elastoplasticity. Furthermore, being combined with the subloading surface concept, the model is endowed with the capability of reproducing complex mechanical behavior of geomaterials, such as hardening/softening accompanied by contractive/dilative responses. The basic characteristics and the predictive capability of the proposed model have been demonstrated through numerical examples involving large deformations and cyclic loadings. Critical assessments with regard to the numerical aspect have ensured the accuracy of the stress calculation, as well as the applicability to boundary-value problem analysis using FEM. The proposed constitutive model requires not a few material constants in return for its high capability. Therefore, developing a systematic parameter determination technique with a guidance for practical engineering use should be our future task.