Elaborated subloading surface model for accurate description of cyclic mobility in granular materials

The description of the cyclic mobility observed prior to the liquefaction in geomaterials requires the sophisticated constitutive formulation to describe the plastic deformation induced during the cyclic loading with the small stress amplitude inside the yield surface. This requirement is realized in the subloading surface model, in which the surface enclosing a purely elastic domain is not assumed, while a purely elastic domain is assumed in other elastoplasticity models. The subloading surface model has been applied widely to the monotonic/cyclic loading behaviors of metals, soils, rocks, concrete, etc., and the sufficient predictions have been attained to some extent. The subloading surface model will be elaborated so as to predict also the cyclic mobility accurately in this article. First, the rigorous translation rule of the similarity center of the normal yield and the subloading surfaces, i.e., elastic core, is formulated. Further, the mixed hardening rule in terms of volumetric and deviatoric plastic strain rates and the rotational hardening rule are formulated to describe the induced anisotropy of granular materials. In addition, the material functions for the elastic modulus, the yield function and the isotropic hardening/softening will be modified for the accurate description of the cyclic mobility. Then, the validity of the present formulation will be verified through comparisons with various test data of cyclic mobility.


a, b
Material constants for deviatoric hardening b r Material constant for rotational hardening c Elastic core (similarity center of normalyield and subloading surfaces) c e Material constant for evolution of elastic core Normalized outward-normal tensor to elastic-core surface 0 Second-order zero tensor p Pressure R ð 1Þ Normal-yield ratio < c ð vÞ Elastic-core yield ratio u; u; u c Material functions and constants regulating normal-yield ratio U Function for evolution of normal-yield ratio w Continuum spin tensor a Conjugate point in subloading surface to a in normal-yield surface a ð¼ 0Þ Reference point in normal-yield surface, which is fixed to the origin of the stress space b Rotational hardening tensor e Accumulated strain e p Accumulated plastic strain / c Internal friction angle in triaxial compression state / d Border angle of deviatoric hardening and softening / r Limitation angle of rotational hardening n ð\0:5Þ Shifting ratio of yield surface to negative pressurẽ j Inclination of swelling line in both logarithmic plane of volume and pressurẽ k Inclination of normal-consolidation line in both logarithmic plane of volume and pressure _ k Plastic multiplier in terms of stress rate _ K Plastic multiplier in terms of strain rate m Poisson's ratio r Cauchy stress tensor r v Conjugate stress on limit elastic-core surface r y Conjugate stress on normal-yield surface h r Lode's angle # ð [ nÞ Ratio of negative pressure to size of yield surface at which volume becomes infinite v ð\1Þ Maximum value of < c

Introduction
It should be noticed that soils exhibit complex deformation behaviors, which are not observed in metals, e.g., the pressure dependence of the elastic moduli and the plastic deformation characteristics, the plastic compressibility, the plastic volumetric expansion induced by the deviatoric stress, i.e., the dilatancy and the rotational anisotropic hardening instead of the kinematic hardening. Moreover, granular materials such as sands exhibit not only volumetric but also deviatoric isotropic hardening, whereas metals and clays exhibit only the deviatoric isotropic hardening and volumetric isotropic hardening, respectively, usually. Above all, one of the most peculiar phenomena among various elastoplastic deformation behaviors in solids would be the cyclic mobility observed prior to the liquefaction induced in sandy ground during earthquakes. Although the permeability of sands is high in general, the shaking during earthquakes occurs at high frequencies so that the deformation occurs under approximately undrained condition. The butterfly-shaped stress loops in the mean effective and deviatoric stress plane, and the S-shaped deviatoric stress vs. axial strain loops with increasing strain amplitude are engendered in sands subject to cyclic loading of deviatoric stress under undrained condition. The cyclic mobility received the focus of attention after the Chile earthquake in May,1960, the Alaska earthquake in March, 1964, and the Niigata earthquake in June, 1964, while the cyclic mobility is termed as cyclic strain softening in the definition by ASCE [6]. The description of cyclic mobility requires one of the most sophisticated formulations in the constitutive modeling. Numerous research reports have been published on the cyclic mobility, e.g., Zienkiewicz et al. [81], Mroz et al. [51], Prevost and Keane [59], Desai et al. [8], Elgamal et al. [10], Fang et al. [11], Hashiguchi and Chen [34], Oka et al. [60], Zhang et al. [72], Hashiguchi [30], etc., within the framework of the elastoplasticity and Zienkiewicz et al. [82], Iai et al. [44], Akiyoshi et al. [1], Gerolymos and Gazetas [17], Zhang and Wang [73] by the empirical approaches. These extensive research efforts have made significant contribution toward elucidating the fundamental mechanism and establishing accurate prediction methods for cyclic mobility.
The description of the cyclic mobility requires the sophisticated constitutive formulation to describe the plastic deformation during the cyclic loading with a small stress amplitude inside the yield surface. The formulation to fulfill this requirement would be able to be attained in the subloading surface model [19,20,28,30], in which the surface enclosing a purely elastic domain is not assumed, while a purely elastic domain is assumed in the other elastoplasticity models. The subloading surface, which passes through the current stress and has a similar shape and direction to the yield surface (renamed the normalyield surface), is assumed inside the normal-yield surface, and then it is postulated that the plastic strain rate develops as the stress approaches the yield surface, i.e., as the subloading surface expands. Therefore, the smooth transition from the elastic to the plastic state, i.e., the smooth elastic-plastic transition leading to the continuous variation of the tangent stiffness modulus tensor, is always described in this model. The subloading surface model has been applied to the descriptions of the elastoplastic deformation behaviors of various solids, e.g., metals [2-4, 12, 39, 41-43, 45-47, 52], soils [15, 16, 34, 35, 38, 53-56, 58, 66, 70-77, 79], rocks and sediments [13,14,78,80] and friction phenomena [36,37,40]. Here, it should be noted that only the elastic deformation is induced in the unloading process and thus closed hysteresis loop cannot be described in the unloading-reloading process if the similarity center of the normal-yield and the subloading surfaces is fixed to the origin of the stress space. The similarity center is physically regarded as the elasticcore because the purely elastic deformation is induced when the stress coincides with the similarity center inside the yield surface. Therefore, the elastic core must be translated with the plastic deformation in order to describe the cyclic loading behavior. The subloading surface model in which the elastic-core is fixed is called the initial subloading surface model and further the model in which the elastic-core is translated is called the extended subloading surface model. The initial subloading surface model has been applied to the description of the cyclic mobility (e.g. [72]). However, the cyclic mobility can be described by the extended subloading surface model by incorporating the rigorous evolution rule of the elastic-core.
The extended subloading surface model will be elaborated so as to describe the cyclic mobility accurately in this article. First, the rigorous translation rule of the elastic-core is incorporated and the pertinent rotational hardening rule is formulated to describe the induced anisotropy of granular materials. In addition, the material functions for the elastic moduli, the yield function, the isotropic hardening/softening, etc., will be extended for the accurate description of the cyclic mobility. Then, the validity of the present formulation will be verified through comparisons with various test data of cyclic mobility for three kinds of sands under various stress amplitudes and several loading cycles up to eighty cycles.
Throughout this paper, the signs of stress and strain are chosen positive for tension and extension, and all stresses signify the so-called effective stress excluding pore pressure. In addition, the direct notation A : B for A rs B rs , AB for A ir B rj , C : A for C ijrs A rs and A : C for A rs C rsij are used for arbitrary second-order tensors A and B, fourth-order tensors C, where Einstein's summation convention is applied for tensor components with repeated indices taking 1, 2, 3. Further, 0 stands for the second-order zero tensor, I for the second-order identity tensor possessing the Kronecker delta components for the magnitude, and symbol A T for the transposed tensor. The symmetric and the antisymmetric parts of A are denoted by sym½A ðA þ A T Þ=2 and ant½A ðA À A T Þ=2, respectively. ð Á Þ and ð Þ stand for the material-time derivative and the proper objective time-derivative, respectively. The symbol hi signifies the Macaulay's bracket defined by hsi¼ ðsþjsjÞ=2, specifically s\0 : hsi¼ 0 and s ! 0 : hsi¼s for an arbitrary scalar variable s. ðÞ 0 denotes initial values of variables.

Strain rate
Denoting the current position vector and the velocity vector of a material particle by x and v, respectively, the velocity gradient tensor is defined as l ¼ ov=ox. The strain rate tensor and the continuum spin tensor are defined as d sym½l and w ant½l, respectively. Limiting the elastic deformation to be infinitesimal, the strain rate d can be additively decomposed into an elastic strain rate d e and a plastic strain rate d p , i.e.
as verified exactly based on the multiplicative decomposition of the deformation gradient tensor (Hashiguchi [32,33]).
First, let d e be given by the hypoelastic relation (Truesdell [65]): where E is the fourth-order tensor describing the elastic stiffness modulus and r is the Cauchy stress, ð Þ denoting the proper objective time-derivative [30,42]. Let E be given in the Hooke's form as where K and G are the pressure-dependent elastic bulk modulus and shear modulus, respectively, leading to the nonlinear elasticity as will be formulated in Sect. 4.1.

Formulation of constitutive equation of granular materials
The elastoplastic constitutive equation of granular materials for describing cyclic loading behavior is formulated below by modifying and elaborating the past formulations based on the subloading surface model [19,20,26,28,30,34]. The incorporation of the rigorous evolution rule of the elastic-core, in which the most elastic response appears since the subloading surface shrinks to a point in the stress space, is of crucial importance for the description of the cyclic loading behavior. In addition, the incorporation of the rigorous evolution rule of the rotational hardening is required for the description of the induced anisotropy. The basic constitutive equation of granular materials taking account of these phenomena will be formulated in this section.

Normal-yield and subloading surfaces
For granular materials, the conventional yield surface, renamed the normal-yield surface, is given by where F is the stress-valued isotropic hardening function of strain-valued isotropic hardening variable H. The rotational hardening variable of the yield surface incorporated first by Sekiguchi and Ohta [63] is denoted by b, while the yield surface rotates around the origin of the stress space. Consequently, the normal-yield surface expands or contracts and rotates around the fixed reference point a which is fixed to the origin of the stress space, i.e., a ¼ 0 in geomaterials. f is assumed to be a function of r in homogeneous degree-one fulfilling f ðjsjr; bÞ ¼ jsjf ðr; bÞ for an arbitrary scalar variable s. As described in Introduction, the mutual slips of solid particles in materials are not induced simultaneously but induced gradually exhibiting a smooth transition from the elastic to the plastic state, which leads to the continuous variation of the elastoplastic tangent stiffness modulus tensor. Then, it is postulated that the plastic strain rate develops gradually as the stress approaches the yield surface. To describe the approaching degree of stress to the yield surface, the subloading surface which passes through the current stress and has a shape and a direction similar to the yield surface, renamed the normal-yield surface, is incorporated. Then, the approaching degree to the normalyield surface is represented by a scalar variable which is defined by the ratio of the size of the subloading surface to that of the normal-yield surface, while the ratio is called the normal-yield ratio and designated by the symbol R ð0 R 1Þ. Further, to describe the cyclic loading behavior, the similarity center of the normal-yield and the subloading surfaces, denoted by the symbol c, is assumed to translate with the plastic deformation. Here, it is called the elastic-core because the most elastic response is induced when the stress lies on it, i.e., r ¼ c leading to R ¼ 0.
Consequently, the subloading surface for the normalyield surface in Eq. (4) is represented by (see Fig. 1 where the following relation holds by virtue of the similarity of the subloading surface to the normal-yield surface. where It follows from Eq. (6) that a ð1 À RÞc ða À c ¼ Rða À cÞ; a ¼ 0Þ ð 8Þ leading to a is the conjugate (similar) point in the subloading surface to the reference point a ð¼ 0Þ in the normal-yield surface as shown in the ðp; qÞ plane in Fig. 1, where p Àðtr rÞ=3 and q r l À r a , r l and r a representing the lateral and the axial stresses, respectively. r y is the conjugate point on the normal-yield surface to the current stress point r on the subloading surface.
noting Eq. (6), which can be rewritten as

Evolution rule of normal-yield ratio
The rate of the normal-yield ratio R is given by the following equation based on the fundamental concept of the subloading surface described in Sect. 3.1.
where U is a monotonically decreasing function of R, satisfying the following conditions (see Fig. 2).
for R [ 1 : over normal-yield stateÞ Equation (14) with Eq. (15) is incorporated into a consistency condition of the subloading surface so that the stress is controlled to be attracted automatically to the normal-yield surface during the plastic loading process, even if it goes out from the normal-yield surface in numerical calculations, because of _ R\0 for R [ 1 by Eq. (14) with Eq. (15) 4 . Let function U satisfying Eq. (15) be given simply by where u is a material function in general, and its specific form will be presented in Sects. 3.6 and 4.5. The smaller u is, the gentler is the transition from the elastic to the plastic state. Equation (14) with Eq. (16) can be analytically integrated as follows: p Normal -yield surface Compression side: Compression side: Fig. 1 Rotated normal-yield, subloading, elastic-core and limit elastic-core surfaces shown in the (p, q) plane

Evolution rule of elastic-core
Let the rigorous evolution rule of the elastic-core be formulated in this section based on the following facts.
(1) In the physical view point, a smooth elastic-plastic transition is not described if the elastic-core lies on the normal-yield surface at which a remarkable plastic deformation is induced, (2) In the mathematical view point, the subloading surface is not determined uniquely if the elasticcore, i.e., the similarity center lies on the normalyield surface, noting the fact: If the elastic-core lies on the normal-yield surface and the stress coincides with the elastic-core, R is indeterminate as known from the relation 0 ¼ R0 which is induced by substituting r ¼ c ¼ r y into r À c ¼ Rðr y À cÞ based on the similarity of the subloading surface to the normal-yield surface.
Consequently, the elastic-core is not allowed to approach the normal-yield surface unlimitedly.
To avoid the unlimited approach of the elastic-core to the normal-yield surface, first let the following surface, called the elastic-core surface, be introduced as shown in Fig. 1, which passes through the elastic-core c and possesses a similar shape and orientation to the normal-yield surface with respect to the null stress ðr ¼ a ¼ 0Þ.
where the variable < c is the ratio of the size of the elasticcore surface to that of the normal-yield surface, called the elastic-core yield ratio. It plays the role of a measure for the approaching degree of the elastic-core to the normalyield surface. Since the elastic-core must lie inside the normal-yield surface as described above, the elastic-core yield ratio has to be less than unity. Then, the inequality Equations (19) and (21) (rate form) are called the enclosing condition of elastic-core. Let the following relation be adopted.
where c e is the material constant controlling the translational rate of the elastic-core and r v is the stress on the limit elastic-core surface conjugate to the current stress r on the subloading surface (see Fig. 1), i.e.
and r y is the conjugate point on the normal-yield surface to the current stress on the subloading surface (see Fig. 1). The inequality in Eq. (21) noting that of ðc; bÞ=oc is the outward-normal of the elasticcore surface at the current elastic-core c and makes an obtuse angle with r v À c when c lies on the limit elasticcore surface f ðc; bÞ ¼ vFðHÞ as far as it is the convex surface, while r v lies on the limit elastic-core surface. The irrational evolution rule adopting r y instead of r v by which the enclosing condition cannot be satisfied in general was proposed by Hashiguchi [26][27][28] and used by Hashiguchi et al. [41] and Fincato and Tsutsumi [12]. Further, the evolution rule was used by Hashiguchi [30], Hashiguchi and Ueno [39], Iguchi et al. [45][46][47], etc., wherê n c of ðc; bÞ oc which is the normalized outward-normal of the elastic-core surface at the current elastic-core c. It cannot be applicable to the generic deformation behavior, since it depends only on the unit outward-normal tensors n andn c independent of the size and the shape of the normal-yield surface as seen in the right-hand side of Eq. (25). Consequently, the following evolution rule of c is given from Eq. (22) with Eq. (23) as follows: The substitution of Eq. (27) in Eq. (9) yields:

Consistency condition
The substitution of Eq. (28) in Eq. (11) leads to the consistency condition: By virtue of the relations (29) is simplified to the equation:

Plastic strain rate
The associated flow rule is adopted for the subloading surface, i.e.
where _ k is the plastic multiplier, i.e., positive proportionality factor. Here, note that the associated flow rule can be applied to the subloading surface even for soils as will be described in Sect. 3.7.
Substitution of Eqs. (14) and (31) in Eq. (30) yields: The explicit functions of F 0 ; h; and b will be formulated later in Eqs. (68), (77) and (81), respectively, in Sect. 4 for granular materials. It follows from Eq. (32) that The strain rate is given from Eq. (35) along with Eqs. (1) and (2) as follows: from which the proportionality factor described in terms of the strain rate in the flow rule (31), denoted by _ K instead of _ k, is derived as follows: Using Eq. (37), the stress rate is given from Eq. (36) as follows: The loading criterion is given as follows (Hashiguchi [21,22]): or where the judgment of yielding whether the stress reaches the yield surface is not required because the plastic strain rate develops continuously as the stress approaches the yield surface.

Modification of unloading-reloading response
The variation of the accumulated plastic strain e p b À e p a ðe p ¼ R t 0 kd p k dtÞ for a certain variation R b À R a of the normal-yield ratio induced during the plastic deformation process from the state a to the state b is identical regardless of loading processes, e.g. initial loading, reloading and inverse loading, proportional and non-proportional loadings as known from Eq. (17), if the parameter u in Eq. (16) is a constant. It leads to the impertinent description that the reloading stress-strain curve after a partial unloading returns to the preceding stress-strain curve too gently. Therefore, it leads to the inadequate prediction of cyclic loading behavior, resulting in an unrealistically large plastic strain accumulation during cyclic loading. The material parameter u is then extended to describe the generalized Masing effect [50] as follows: where C n n c : n ðÀ1 C n 1Þ ð 42Þ u (average value of u) and u c are material constants and u is a continuous function of variables < c and C n . The forms of the function u for particular states are shown in the bracket. C n ¼ 1, 0 and À1 indicate states for which the plastic strain rate is directed outward-normal, tangential and inwardnormal, respectively, to the elastic-core surface. By this modification, a realistic description is given for the phenomenon in which the reloading curve after a partial unloading returns rapidly to the preceding loading curve and instead the curvature of the inverse loading curve decreases leading to the generalized Masing effect.

Basic characteristics of subloading surface model
The subloading surface model possesses the following distinguished features.
(1) The plastic deformation is not induced abruptly but develops gradually. In fact, mutual slips between material particles, e.g., crystal particles in metals and soil particles in sands and clays, are not induced simultaneously but induced gradually from parts in which mutual slips can be induced easily. Therefore, the plastic strain rate develops continuously as the stress approaches the yield surface. Then, the elastoplastic constitutive equation is required to satisfy the following smoothness condition (Hashiguchi [21,22,24] where H and H are the second-order tensor-valued and the scalar-valued internal variables, respectively, and dðÞ stands for an infinitesimal variation. The rate-linear constitutive equation is described as where the fourth-order tensor M ep is the tangent stiffness modulus tensor, which is the function of the stress and internal variables, and can be described generally by The smoothness condition is satisfied in the subloading surface model, while it is violated in the yield point by the other elastoplastic models because they assume the purely elastic domain. (2) The smooth transition from the elastic to the plastic state and the continuous variation of the tangent stiffness modulus tensor are described by the subloading surface model. However, it is violated at the yield point in any other elastoplasticity models because they assume a purely elastic domain inheriting the conventional elastoplasticity model. (3) The yield judgment whether the stress reaches the yield surface is unnecessary since the plastic strain rate develops continuously as the stress approaches the normal-yield surface. In contrast, a yield judgment is required in the other elastoplasticity models because they assume a surface enclosing a pure elastic domain, while the determination of the yield stress is accompanied with an arbitrariness because the tested stress vs. strain curve are usually smooth. (4) The tangent stiffness modulus changes always continuously, while it changes abruptly at the yield point in the other elastoplasticity model. (5) The plastic strain rate can be described for any small stress variation and for cyclic loading under any small stress amplitudes since a pure elastic domain is not assumed. However, it cannot be described during the stress variation inside the yield surface enclosing a purely elastic domain in the other elastoplasticity models. (6) The automatic stress-controlling function is furnished such that the stress is always attracted to the normal-yield surface. In particular, it is noticeable that the stress is automatically pulled back to the normal-yield surface when it goes over the surface in numerical calculation because of _ R\0 for R [ 1 from Eq. (14) with Eq. (15) 4 (see Fig. 3). In contrast, the particular operation to pull back the stress is required in the other elastoplasticity models. (7) The outward-normal n at the peak stress point in the subloading surface is approximately coincides with the outward-normal of the plastic potential surface adopted in the Drucker-Prager model [9] as shown in Fig. 4 in which n DP designates the outwardnormal of the Drucker-Prager yield surface. Therefore, the associated flow rule can be adopted in the subloading surface model. On the other hand, the non-associated flow must be adopted in the Cap model composed of the Drucker-Prager model in the over-consolidated state and the Cam-clay model [61,62] in the normal-consolidated state in order to suppress the excessive volumetric expansion in the over-consolidated state, which leads to the asymmetry of the elastoplastic tangent modulus tensor.

Material functions for granular materials
The material functions included in the subloading surface model described above are now formulated for a wide range of granular materials involving clays and sands, generalizing the past formulations to achieve the accurate description of cyclic loading behavior up to the cyclic mobility.

Elastic moduli
The elastic bulk modulus K and the elastic shear modulus G are given by Hashiguchi [31], modifying the former equations to describe the pressure-dependency (Hashiguchi [23,32], Hashiguchi and Chen [34]) as follows:  leading to e e v ! 1 (infinite volume expansion) for p ! À#F, while # ð0 #\0:1Þ for common clays and sands (Hashiguchi and Chen [34], Hashiguchi and Mase [35]). d e v is the elastic volumetric strain rate d e v tr d e . The power function of the pressure in Eq. (48) for the shear modulus is referred to Tatsuoka et al. [64]. n ð 1Þ is a material constant, while n ffi 0:5 is chosen for most sands. Consequently, the elastic tangent modulus depends naturally on pressure p and isotropic hardening function F. The elastic moduli in Eq. (48) hold for granular materials consistently in the frameworks of the infinitesimal elasticity, the hypoelasticity and multiplicative hyperelasticity as was verified by Hashiguchi [31].

Yield and subloading functions
Let the normal-yield surface be given by the modified Cam-clay model (Burland [7], Roscoe and Burland [61]): where M is the stress ratio kr 0 k=p in the critical state. To take account of the influence of the third deviatoric stress invariant, M is extended as follows (Hashiguchi [25]): where h r is the so-called Lode angle defined by cos 3h r ffiffi ffi Sands possess quite small strength in the negative pressure range. However, it is of crucial importance to shift the yield surface to the negative pressure range in order to describe rigorously the cyclic mobility, in which the stress changes actively around the null stress state. Here, it should be noticed that the shift of the yield surface is of the importance also in the computational aspect.
To shift the yield surface to the negative pressure range, Eq. (50) is modified by nFðp ! p þ nFÞ as follows (Hashiguchi and Mase [35]).
p À ðð1=2Þ À nÞF F=2 where n is the material constant, while it must fulfill n 1=2 since the tensile yield stress is smaller than the compression yield stress, and further the inequality n\# is required since the volume does not become infinite by the elastic deformation inside the yield surface, satisfying p ! À nF [ À #F. The yield surface in Eq. (54) is shown in Fig. 6 in the pressure-deviatoric stress plane. Equation (54) is extended by taking account of the rotation of the yield surface around the origin of the stress space into account and thus by replacing r 0 with r 0 À pb, where b is referred to as the rotational hardening variable represented by a deviatoric tensor, the evolution rule of which will be formulated in Sect. 4.4. This extension leads to p À ðð1=2Þ À nÞF F=2 where Equation (55) is rewritten as follows: Equation (59) is the quadratic equation of the hardening function F. By solving this equation, it can be described in the separated form to the function f ðp; q _ Þ and the hardening function F as follows: wherẽ n 2ð1 À nÞn; n 1 À 2n; p _ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Note that an analytically separated form cannot be derived from the yield conditions other than the modified Cam-clay model, e.g. the original Cam-clay model (Schofield and Wroth [62]) which is shifted to the negative pressure range by p ! p þ nF leading to p exp Â kr 0 k=ðpMÞ

Isotropic hardening by volumetric and deviatoric plastic strain rates
The hardening function in Eqs. (4) or (5) is given as follows (Hashiguchi [32]): wherek andj stand for the slopes of the normal-consolidation curve and the swelling curve, respectively, in the ln v À lnðp þ #FÞ plane (Hashiguchi [23]), which is based on the following isotropic consolidation characteristics.
where e e v and e p v are the elastic and the plastic logarithmic volumetric strains, respectively; v, v 0 and v are the current volume, the initial volume and the volume unloaded state to the initial pressure, respectively; and p, p 0 , p y and p y0 are the current, the initial and the current yield and the initial yield pressure, respectively, as shown in Fig. 6. Replacing the pressures p y and p y0 to the isotropic hardening function F and its initial value F 0 , respectively, in Eq. (69), one has e e v ¼ Àj ln The isotropic hardening function is given from Eq. (70) by choosing the isotropic hardening variable H to be the plastic volumetric contraction, i.e., H ¼ Àe p v as follows: The rate of isotropic hardening/softening variable H was extended by Nova [57] and Wilde [67] to incorporate the influence of the deviatoric plastic strain rate, called the deviatoric (isotropic) hardening, in addition to the volumetric plastic strain rate described above. Here, aiming at establishing a quantitative description of cyclic mobility, let _ H be extended to depend on the plastic volumetric strain rate d p v and the plastic shear rate d p s for a wide range of stress up to the negative pressure as follows (Fig. 6): where l d , / d ð\/ c Þ, a ð ! 1Þ and b ð [ 1Þ are material constants. Here, / d \/ c holds in general since the deviatoric stress can increase over the critical state (tr n ¼ 0 , tr d p ¼ 0) when the stress ratio increases at a low pressure under the undrained condition even in loose sands. M d in Eq. (76) is formulated analogously to M in Eq. (51). Hardening and softening are induced outside and inside, respectively, the conical surface kr 0 k ¼ M d ðp þ #FÞ which is called the deviatoric hardening boundary surface. It is extended to accommodate the negative pressure p ðÀ#F\p 0Þ by replacing the normalized stress ratio ðkr 0 k=pÞ=M d to v d defined in Eq. (75). The deviatoric hardening rate depends nonlinearly on the modified stress ratio v d with the upper limit l d kd p0 k of the plastic shear rate d p s (Fig. 7).

Rotational hardening: Anisotropic hardening by rotation of yield surface
The rotational hardening rule is formulated by noting that the anisotropic hardening is induced only by the deviatoric plastic strain rate as follows (Fig. 8): where M _ À cos 3h b r is the material constant regulating the rate of the rotation. M _ r is formulated analogously to M in Eqs. (51), (57) and (76), where / r is the material constant denoting the limit of the rotational angle of the normal-yield and subloading surfaces. Note here that Eq. (78) takes a form analogous to the evolution equation of back stress tensor in the nonlinear kinematic hardening rule [5]. For the subloading surface in Eq. (63) based on the modified Camclay model, the rotation does not occur, i.e., b ¼ 0 when the stress lies on the central axis of the subloading surface, fulfilling r 0 ¼ pb, for which n 0 ¼ 0 leading to d p0 ¼ 0 holds as shown in Fig. 1.

Extension of material parameter for evolution of normal-yield ratio
The material parameter u in the evolution equation (41) of the normal-yield ratio is extended as follows (Fig. 9): where e p0 Z t 0 kd p0 k dt ð83Þ u 0 , m and u e are material constants. The material parameter u becomes smaller inducing a larger plastic strain rate as the accumulation of the deviatoric plastic strain rate proceeds, while this trend is more notable for larger values of u e . In addition, u is inversely proportional to M _ and therefore a larger plastic strain rate is induced in the compression side than in the tension side, while this trend is more notable for larger value of m.  4. The difference between the tangent stiffness moduli in reloading and reverse loading is larger and thus the difference between the curvatures of the stress-strain curves in them is larger for larger u c values. 5. Plastic deformation begins sooner after unloading for larger c e values for which the closed hysteresis loop is depicted so that the strain accumulation is suppressed.

Summary of material parameters and their physical meanings
On the other hand, the open hysteresis loop is depicted for c e ¼ 0, returning to the initial subloading surface model. 6. The material constant / d , which regulates the boundary of the deviatoric hardening and softening, is larger for a looser material with a wider range of deviatoric softening. 7. The normal-yield and subloading surfaces rotates in a wider range for a larger / r value. They rotate more Table 1 Physical properties of tested sands, test condition and material parameters in simulations The identifications of the material parameters are commented below.
k andj are determined by the curve-fitting to the test data of the isotropic consolidation. Also, they can be known easily byk ¼ k=ð1 þ e 0 Þ;j ¼ j=ð1 þ e 0 Þ, if there is past data of k and j in the e À log p linear relation.
G 0 and n are determined from the stress vs. strain curve under constant pressure inside the yield surface, while we may use n ¼ 0:5 usually.
/ c is determined by the inclination of the critical state line in the triaxial compression.
n is determined from the pressures in the isotropic compression and the extension test. l d ; / d ; a; and b are determined from the undrained triaxial compression test data.
b r and / r are determined from the loading and unloading process in the triaxial compression or extension test. u c ; u; u e ; and m are inferred from the reloading curve in the drained test. c e and v are determined from the loading-unloading curve, while we may put v ¼ 0:7 usually.
F 0 can be determined from the initial normal-consolidation stress.
b 0 may be chosen to be b 0 ¼ 0 in isotropic-consolidated sands but it can be calculated from K 0 -value in a K 0consolidation.
c 0 is determined from the stress at which the most elastic behavior is observed.

Simulations of test data for cyclic mobility
The constitutive equation of granular materials formulated in the previous sections is applied to the simulations of various test data on the cyclic mobility under the undrained condition with the constant deviatoric stress amplitudes.
All the test data adopted for the simulations were obtained by the cyclic triaxial compression/extension tests with symmetric constant deviatoric stress amplitudes from the isotropic stress state under constant total lateral confining pressures denoted by p c . The initial isotropy, i.e., b 0 ¼ 0 is assumed and the typical values n ¼ 0:5 and v ¼ The physical properties of the tested sands, the test condition and the material parameters used in the simulations are listed in Table 1.
The simulation results of the Toyoura sand with the almost same initial void ratio ðe 0 ¼ 0:722 and 0:718Þ are shown for Figs. 10 and 11, while the same values of the material parameters are used in these simulations. Further, the simulations of the Tone river sand and the Edo river sand are shown for Figs. 12 and 13, respectively. The butterfly-shaped stress paths after the gradual decreases of the pressure and the deviatoric stress vs. axial strain curves with the increasing amplitudes of axial strain are simulated in these figures. The test results are simulated well in these calculations, where various stress ratios and numbers of loading cycles are applied. The maximum amplitudes of the axial strain are simulated well even for the high loading cycles up to the eighty-seven cycles. However, the deceases of the pressure are more gradual in the simulations than in the test results. The amplitudes of the axial strain in the simulations are smaller at the beginnings and later become larger compared with the test results in Figs. 10 and 11 for the Toyoura sand. The stress paths in the simulations are warped inversely to the test results in Figs. 12 and 13 for the Tone river sand and the Edo river sand. Further improvements would be required for a more accurate prediction of cyclic mobility.
Since seismic waves cause the cyclic loading at high frequency, the ground deformation occurs under fully or nearly undrained condition during an earthquake. However, when the earthquake ceases, the vertical load due to the weight of the ground and the gravity of the building floating during the earthquake acts to the soil grounds pushing out the pore water, and thus the soil skeleton deforms under the drained conditions. Therefore, it is required to describe both the cyclic mobility under undrained condition and the monotonic loading behavior under drained condition by a unique elastoplastic constitutive equation with a unified set of material parameters. On the other hand, although there are few test results of cyclic mobility and drained tests on specimens of the same material and the same void ratio, the results of drained tests obtained using Toyoura used in the above-mentioned cyclic mobility are provided from the same experimenter to the authors [69]. The simulation result under the drained condition is shown in Fig. 14, where the tested and the calculated stress and volume strain curves are shown. Therein, since the initial void ratio e 0 ¼ 0:748 in the drained test is sufficiently close to e 0 ¼ 0:722 and 0.718 in the above-mentioned cyclic mobility test, the material parameters shown in Table 1 (which are used in Figs. 10 and 11) are used. A fairly good simulation to the test result is seen in this figure.
Incidentally, for reference, the simulation to the drained test result with the initial void ratio e 0 ¼ 0:924 (looser than in Fig. 14) of the Toyoura sand is shown in Fig. 15, choosing the material parameters as follows: F 0 ¼ 170 kPa; c 0 ¼ À60 I

Conclusion
In this article, the subloading surface model is elaborated to describe the cyclic mobility observed prior to the liquefaction induced by earthquakes. The results obtained in this study are summarized as follows: 1. The elastoplastic constitutive equation of geomaterials is elaborated by the formulations of 2. the rigorous evolution rule of the elastic-core, which is of crucial importance for the description of cyclic loading behavior, 3. the evolution rule of the rotational hardening, which is of crucial importance for the description of the anisotropic hardening behavior of granular materials, 4. the material functions for the elastic moduli extended to the pressure dependence, the yield condition extended to the negative pressure range, the isotropic hardening rule based on not only the volumetric but also the deviatoric plastic deformations for the accurate description of the cyclic loading behavior up to the cyclic mobility. data on the cyclic mobility in three kinds of sands for various numbers of cycles under various stress amplitudes. The butterfly-shaped stress path in the pressure vs. deviatoric stress plane and the S-shaped deviatoric stress vs. axial strain loops were reproduced accurately.
6. The present formulation would provide the substantial foundation for the elastoplastic description of cyclic mobility, although further elaboration is desirable for a highly accurate description of the cyclic mobility.   7. The present constitutive model would be applicable to the predictions of the deformation behavior not only during the earthquake but also after it.
Not a few material constants are incorporated in order to describe the cyclic mobility accurately by the elastoplastic constitutive equation because it is to be one of the most complicated mechanical phenomena in the natural world. Further study is required to clarify their definite physical meanings for formulating a more sophisticated constitutive equation of geomaterials.
The formulation falls within framework of the hypoelastic-based plasticity which is limited to the description of the infinitesimal elastic deformation and requires the cumbersome time-integrations of corotational rates of stress and internal variables. Then, it should be extended to the multiplicative hyperelastic-based plasticity (Hashiguchi [29,33], Hashiguchi and Yamakawa [42]) with a rigorous description of hyperelastic equation of granular materials (Hashiguchi [31]) for further accurate description of cyclic mobility up to the finite deformation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.