Multiplicative Hyperelastic-Based Plasticity for Finite Elastoplastic Deformation/Sliding: A Comprehensive Review

Hyperelastic-based plastic constitutive equation based on the multiplicative decomposition of the deformation gradient tensor is reviewed comprehensively and its exact formulation is given for the description of the finite deformation and rotation in this article. Further, it is extended to describe the general loading behavior including the monotonic, the cyclic and the non-proportional loading behaviors by incorporating the rigorous plastic flow rules and the subloading surface model. In addition, it is extended also to the rate-dependency based on the overstress model, and the exact hyperelastic-based plastic constitutive equation of friction is formulated by incorporating the subloading-friction model. They are the exact constitutive equations describing the monotonic and the cyclic loading behavior up to the finite deformation/rotation and the friction behavior under the finite sliding/rotation with the rate-dependency, which have remained to be unsolved for a long time, although they have been required in the history of elastoplasticity theory.


Introduction
The elastic deformation and the plastic deformation are physically different to each other such that the former is induced by the deformation of material particles themselves but the latter is induced by the mutual slips between the material particles. Therefore, the elastoplasticity is based on the premise that the deformation is decomposed into the elastic and the plastic deformations, so that the irreversible change of substructure is described by the isotropic and the anisotropic hardenings which evolve only by the plastic deformation, while the elastic deformation is irrelevant to the irreversible change of substructure. Here, it should be noticed that the deformation and its rate can be described exactly by the deformation gradient tensor which transforms the reference infinitesimal line element to the current one. Therefore, the exact elastoplastic constitutive equation must be formulated by incorporating the definite decomposition of the deformation gradient tensor into the elastic and the plastic parts which is realized by the multiplicative decomposition of the deformation gradient tensor [52,55,56,[58][59][60]. However, it now passed already a half century after the proposition of the multiplicative decomposition of deformation gradient tensor. In the meantime, unfortunately the hypoelastic-based plasticity has been studied enthusiastically by numerous workers represented by Rodney Hill and James R. Rice after the proposition of the hypoelasticity by Truesdell [83], which is not based on the multiplicative decomposition definitely.
The multiplicative hyperelastic-based plasticity has been studied centrally by Simo and his colleagues (e.g. [68,[74][75][76][77][79][80][81]) in the last century, in which the logarithmic strain has been used mainly leading to the coaxiality of stress and strain rate so that it has been limited to the isotropy. It has been developed actively from this century on by Menzel and Steinmann [63], Wallin et al. [87], Dettmer and Reese [14], Menzel et al. [64], Wallin and Ristinmaa [86], Gurtin and Anand [18], Sansour et al. [73], Vladimirov et al. [84,85], Henann and Anand [47], Hashiguchi and Yamakawa [46], Brepols et al. [6], Hashiguchi [27,31], etc. in which constitutive relations are formulated in the intermediate configuration imagined fictitiously by the hyperelastic unloading to the stress-free state. However, the plastic flow rule with the generality for the elastically anisotropy remains unsolved and only the conventional plasticity model (named by Drucker [15]) with the yield surface enclosing the elastic domain have been incorporated so that only the monotonic loading The precise description of the plastic strain rate induced by the rate of stress inside the yield surface is inevitable for the prediction of cyclic loading behavior which is crucial for the accurate mechanical design of solids and structures in engineering. A lot of works have been executed and various unconventional plastic constitutive (cyclic plasticity) models (named by Drucker [15]) have been proposed aiming at describing the plastic strain rate caused by the rate of stress inside the yield surface after 1960s when the demands of mechanical designs for the vibration proof and the earthquake proof have been highly raised responding to the high development of machine industries and the frequent occurrences of earthquakes in Niigata (Japan) and Chili, etc. Among various unconventional models the multi surface model [51,66], the two surface model [12,53,91] and the superposed-kinematic hardening model [10,67] are well-known. However, they inherit a small yield surface enclosing purely-elastic domain from the conventional plasticity model and are based on the premise that the plastic strain rate develops with the translation of the small yield surface so that they are called the cyclic kinematic hardening model. Therefore, they possess various defects, e.g. (1) the abrupt transition from the elastic to the plastic state violating the continuity and the smoothness conditions [21,22,24], (2) the incorporation of the offset value of the plastic strain at yield, (3) the incapability of cyclic loading behavior for the stress amplitude less than the small yield surface assumed inside the conventional yield surface, (4) the incapability of the non-proportional loading behavior, (5) the incapability of extension to the rate-dependency at high rate of deformation up to the impact loading behavior, (6) the limitation to the description of deformation behavior in metals, (7) the necessity of the additional cumbersome operation to pull-back the stress to the yield surface. On the other hand, only the subloading surface model [19,20,33] does not incorporate the yield surface enclosing a purely-elastic domain so that it possesses the high ability for describing not only the monotonic but also the cyclic and the non-proportional loading behaviors, resolving all the above-mentioned defects in the other elastoplasticity models. In addition, it possesses the automatic controlling function to pull-back the stress to the yield surface when the stress goes out from the yield surface in numerical calculations [19,20,33]. It is capable of describing the elastoplastic deformation behavior in not only metals but also soils [38,42,89] and further the friction behavior rigorously [24,26,40,41]. It would be regarded to be the governing law of the irreversible mechanical phenomena of solids. The subloading surface model has been incorporated to the commercial software Marc 2017. 1

version in Marc
Software Corporation as the standard installation by the name ''Hashiguchi model'', which can be used by all Marc contractors (users). Therefore, it is explained in the Marc user manual [61] in brief. Needless to say, however, it is limited to the formulation in the ordinary current configuration, whereas there is no commercial software installed the multiplicative hyperelastic-based plastic constitutive equation up to date.
First, the basic frameworks of the elastoplastic constitutive equations will be reviewed briefly in order to confirm the necessity of the multiplicative decomposition of the deformation gradient tensor into the elastic and the plastic parts. Therein, the fundamentals in the multiplicative decomposition will be delineated concisely where the interpretation of the isoclinic concept is given thoroughly. Then, the exact hyperelastic-based plastic constitutive equation will be formulated within the framework of the multiplicative decomposition of the deformation gradient tensor, incorporating the rigorous plastic flow rules and the subloading surface model. It will be extended to the ratedependency based on the overstress model by revising the former formulations of the overstress model so as to be applicable to the general rate ranging from the quasi-static to the impact loading behaviors [31,33]. Further, the exact hyperelastic-based plastic constitutive equation of friction is formulated rigorously, in which not only the rotation but also the deformation of the contact surface can be incorporated, although the hypoelastic-based plastic constitutive equation has been formulated in the past [31,33,40,41]. They are the exact constitutive equations describing not only the monotonic and the cyclic loading behavior under the finite deformation/rotation and the cyclic friction behavior under the finite sliding/rotation. Then, the comprehensive review will be given for the multiplicative hyperelastic-based plasticity for finite deformation/rotation and finite friction-sliding/rotation, incorporating the various novel, rational formulations which have not been shown in existing literatures.
The direct notations u Á v for u r v r , u v for u i v j , Av for A ir v r , A : B for A rs B rs , AB for A ir B rj , C : A for C ijrs A rs , A : C for A rs C rsij , A B for A ij B kl and C : N for C ijrs N rskl are used for arbitrary vectors u and v, second-order tensors A and B, fourth-order tensors C and N, where the Einstein's summation convention is applied for letters with repeated indices taking 1, 2, 3. Further, 0 and O stand for the zero vector and tensor, respectively, I and I stand for the second-order and the fourth-order identity tensors possessing the components of the Kronecker's delta d ij ðd ij ¼ 1 for i ¼ j; d ij ¼ 0 for i 6 ¼ jÞ and I ijkl ¼ d ik d jl , respectively, trA ¼ A ij d ij for the trace, A 0 ¼ A À (trAÞI=3 for the deviatoric tensor, A À1 for the inverse tensor satisfying AA À1 ¼ I and kAk ¼ ffiffiffiffiffiffiffiffiffiffi ffi A ij A ij p for the magnitude and the symbol ðÞ T for the transposed tensor. The notations sym½A ðA þ A T Þ=2 and ant½A ðA À A T Þ=2 stand for the symmetric and the antisymmetric parts, respectively. ( Á ) stands for the material-time derivative. The symbol h i designates the Macaulay's bracket defined by hsi ¼ ðs þ jsj=2Þ, i.e. s\0 : hsi ¼ 0 and s ! 0 : hsi ¼ s (s: arbitrary scalar variable).

Classification of Elastoplastic Constitutive Equations
The advantages of the multiplicative hyperelastic-based plasticity should be clarified prior to its formulation which is to be the purpose of this article. The basic frameworks of elastoplastic constitutive equations are classified and their features are described concisely in the following.

Infinitesimal Elastoplasticity
1. The infinitesimal strain is additively decomposed into the elastic and the plastic parts. 2. The reference and the current configurations are not distinguished. 3. The Cauchy stress is used as the stress measure. 4. The stress is related to the elastic strain by the hyperelastic relation based on the strain energy potential function. 5. The plastic strain rate is calculated and its accumulation leads to the plastic strain. 6. The elastic strain is calculated by subtracting the plastic strain from the strain. 7. The stress is calculated by substituting the elastic strain into the hyperelastic equation. 8. The infinitesimal elastic and plastic deformation without a rotation is described.

Hypoelastic-Based Plasticity
1. The symmetric and the anti-symmetric parts of the velocity gradient are defined as the strain rate and the spin, respectively. 2. The strain rate and the spin are additively decomposed into the elastic and the plastic parts. 3. The formulation is performed in the current configuration which is influenced by the material rotation. 4. The elastic strain rate is related to the stress rate by the hypoelastic relation [83], so that the elastic deformation cannot be described exactly. 5. The corotational stress rate is used in order to exclude the influence of material rotation from the materialtime derivative of stress, while the responses by various corotational stress rate have been examined by Dafalias [11], Zbib and Aifantis [92], Gambirasio et al. [17], etc. 6. The cumbersome operation is required for the timeintegration of the corotational stress rate in order to calculate the stress (cf. [13,33,78], etc.). 7. The finite plastic deformation with the finite rotation are described under the restriction of the infinitesimal elastic deformation.
Most of past literatures ( [3,4,13,33,78], etc.) are concerned mainly with the explanation within the framework of the hypoelastic-based plasticity. tion is adopted as the stress measure, which is workconjugate to the strain rate in the intermediate configuration. 5. The second Piola-Kirchhoff stress is related to the right Cauchy-Green deformation tensor in the intermediate configuration by the hyperelastic relation, and the Mandel stress is calculated by multiplying the elastic right Cauchy-Green deformation tensor to the second Piola-Kirchhoff stress. 6. The plastic velocity gradient is calculated and its accumulation leads to the plastic deformation gradient. 7. The elastic deformation gradient is calculated by excluding the plastic deformation gradient from the deformation gradient, from which the elastic right Cauchy-Green deformation tensor is calculated. 8. The second Piola-Kirchhoff stress is calculated by substituting the elastic right Cauchy-Green deformation tensor into the hyperelastic relation. 9. The finite elastic and plastic deformations with the finite rotation can be described exactly.

Multiplicative Hyperelastic-Based Plasticity
Consequently, the exact elastoplastic constitutive equation based on the multiplicative decomposition will be formulated in the subsequent sections.

Multiplicative Decomposition of Deformation Gradient Tensor Based on Isoclinic Concept
The following physical facts should be noticed for the rigorous formulation of elastoplastic constitutive equation.
(1) Decomposition of deformation/rotation into elastic and plastic parts: The physical mechanisms of the elastic deformation and the plastic deformation are substantially different to each other such that the former is induced by the deformation of material particles themselves and the latter is induced by the mutual slips between material particles. Therefore, the deformation must be decomposed into the purely elastic deformation with the potential energy function and the plastic deformation exactly. Then, the irreversible variation of mechanical response is described rigorously by incorporating the internal variables which evolve only by the plastic deformation. If not, the irrational result is caused such that the internal variables evolve even by the elastic deformation. (2) Deformation/rotation (rate) measures based on deformation gradient tensor: The deformation/rotation (rate) of materials can be described exactly by the deformation gradient tensor F. Then, the elastic and the plastic parts must be described rigorously by the definite decomposition of the deformation gradient tensor into the elastic and the plastic parts.
Based on the above-mentioned basic requirements, the multiplicative decomposition of the deformation gradient tensor was advocated by Kroner [54], Lee and Liu [56], Lee [55], Mandel [58][59][60] and Kratochvil [52]. Then, the deformation gradient F ¼ ox=oX, where x and X are the position vectors in the current and the reference configurations, respectively, is multiplicatively decomposed into the elastic deformation gradient F e and the plastic deformation gradient F p .
Here, based on the requirement of the exact decomposition of deformation into the purely elastic and the plastic parts, F e is calculated from the current stress by the hyperelastic relation reflecting the real hyperelastic property of material and then F p is calculated from F and F e . In other words, the plastic deformation gradient F p is obtained by the unloading to the stress-free state along the hyperelastic relation. The unloaded configuration to the stress-free state is called the intermediate configuration.
Equation (1) is regarded as the three-dimensional extension of the multiplicative decomposition of the stretch k ¼ l=l 0 into k ¼ k e k p with k e ¼ l=l p ; k p ¼ l p =l 0 where l is the current length, l 0 is the initial length, and l p is the length in the unloaded state to the stress-free state.
Here, note that solids possess the heterogeneous substructures which are statically-indeterminate (non-static stability) in general. Therefore, the purely-elastic deformation is induced only in the initiation of unloading process and the plastic deformation is slightly induced (removed) in the actual unloading process to the stress-free state. In order to let all material points be released to the real stress-free state, we must give different amounts of destressing to individual material points by cutting a material up into pieces. Therefore, the vector d X of the infinitesimal line-element in the intermediate configuration is merely a fictitious vector differing from an actual position vector because it is merely calculated to fulfill the following equation from the elastic deformation gradient tensor F e or the plastic deformation gradient F p , while dX and dx are the actual infinitesimal line-elements.
In addition, the following equations do not possess any primal physical roles.
The elastoplastic deformation process base on this notion is illustrated in Fig. 1 where the initial, the intermediate and the current configurations are specified by the symbols K 0 , K and K, respectively. The tensors in the current configuration are designated by the lowercase letter as t, the ones in the reference configuration by the uppercase letters as T and the ones in the intermediate configuration by the uppercase letters with the upper bar as T. Note here that Fig. 1 Multiplicative decomposition of deformation gradient (1) A rigid-body rotation is involved in the deformation gradient F in general. (2) The multiplicative decomposition in Eq. (1) is where R is an arbitrary rigid-body rotation added to the intermediate configuration and Therefore, the intermediate configuration is not determined uniquely depending on the rigid-body rotation added to the intermediate configuration.
Then, the extensive debates as to which an elastic deformation gradient or a plastic deformation gradient must involve the rigid-body rotation have been repeated for a long time after the proposition of the multiplicative decomposition. The inclusion of the rigid-body rotation in the plastic deformation gradient has been insisted in many literatures. This situation would be caused by worrying the fact that the elastic distortion is known from the current stress but the rigid-body rotation is unknown and thus it is possible to exclude only the elastic distortion but it is impossible to exclude both of the rigid-body rotation and the elastic distortion from the current configuration in order to get to the intermediate configuration. On the other hand, the inclusion of the rigid-body rotation in the elastic deformation gradient has been insisted by the others.
One has only to assume the intermediate configuration so as to lead to the rigorous formulation of constitutive equation, since it is not actual but merely fictitious. The multiplicative hyperelastic-based plastic constitutive equation can be formulated rigorously in the intermediate configuration, since the strain rate and the spin tensors in the intermediate configuration can be decomposed into the purely elastic and plastic parts exactly as will be shown in Eq. (13) with Eq. (10). Then, the plastic and the elastic deformation gradient tensors are formulated as follows: (1) The plastic deformation gradient F p is formulated by the plastic constitutive equation which, needless to say, is independent of the rigid-body rotation. Here, note that the plastic deformation is independent of the rotation of substructure since the material particles move in parallel along the substructure. (2) The elastic deformation gradient F e is formulated by excluding F p from the deformation gradient F. (3) The elastic orthogonal tensor R e in the polar decomposition F e ¼ R e U e stands for the rotation of substructure.
Consequently, the substructure does not rotate in the intermediate configuration, where, needless to say, the rigid-body rotation is involved in the elastic deformation gradient tensor. It has been called the isoclinic concept by Mandel [58][59][60], while ''isoclinic'' means ''equal inclination'', which is schematically shown in Fig. 2.
Here, it should be noted that the actual configuration is the current configuration but the intermediate configuration is the fictitious one. Therefore, it is natural to formulate first the constitutive equation in the current configuration and after that we transform it to the intermediate configuration, resulting in the multiplicative hyperelastic-based plastic constitutive equation.
Further, F p is decomposed into the plastic storage part F p ks causing the kinematic hardening and the plastic dissipative part F p kd multiplicatively [57] as follows (see Fig. 3): Based on the right Cauchy-Green deformation tensor the following tensors of the storage parts C e and C _ p ks and the dissipative parts C p and C _ p  Fig. 3.
The velocity gradient l in the current configuration K is additively decomposed into the elastic and the plastic parts: where l _ FF À1 ; Further, the velocity gradient L defined as the contravariant-covariant pull-back of the velocity gradient tensor l in the current configuration to the intermediate configuration K can be additively decomposed into the purely elastic and the purely plastic parts exactly as follows: where Multiplicative Hyperelastic-Based Plasticity for Finite Elastoplastic Deformation… Therefore, L and L e ; L p can be adopted in the exact formulation of elastoplastic constitutive equation. Let them be decomposed additively into the symmetric and the antisymmetric parts, i.e.  Fig. 3 Multiplicative decompositions of deformation gradient tensor for material with kinematic hardening D ¼ sym½ The material-time derivative of C e is given from Eqs. (5) 1 and (9) as Further, the plastic velocity gradient L p is additively decomposed into the storage and the dissipative parts for the kinematic hardening as follows: noting The material-time derivative of C _ p ks in Eq. (5) is given by

Stress Measures
Introduce the second Piola-Kirchhoff stress tensor S in the intermediate configuration, which is the contravariant pullback of the Kirchhoff stress tensor s, i.e.
and the Mandel stress Here, note that the work-conjugate stress measure with the strain rare L in the intermediate configuration is the Mandel stress M as known from Further, the contravariant push-forward of the 2nd Piola-Kirchhoff stress-like variable S _ k in the kinematic hardening variable from K _ to K is given by Further, the Mandel-like variable M k for the kinematic hardening variable is given by noting Eq. (6). Note here that the Mandel stress M is not symmetric tensor in general but the Mandel-like kinematic variable M k is the symmetric tensor. The material-time derivative of the kinematic hardening variable S k in the intermediate configuration is given by from Eqs. (16) 1 and (23), noting Further, the material-time derivative of M k is given from Eq. (25) with Eqs. (15) and (24) by

Hyperelastic Constitutive Equations
The 2nd Piola-Kirchhoff stress push-forwarded to the intermediate configuration, S, is given by incorporating the strain energy function wð C e Þ as follows: and the Mandel stress is given as The tensor M satisfies the symmetry M ¼ M T for the particular case that w e is the function of invariants of C e , leading to the elastic isotropy. The material-time derivative of the Mandel stress is given noting Eq. (14) as where L e is the fourth-order hyperelastic tangent modulus tensor given by with Further, let S _ k be formulated incorporating the strain energy function w k ðC _ p ks Þ as noting Eqs. (23) and (24).
The material-time derivative of S _ k is given from Eq. (32) with Eq. (19) as Substituting Eq. (34) into Eq. (26), _ M k is given as follows:

Multiplicative Hyperelastic-Based Plastic Equation for Conventional Model
The multiplicative hyperelastic-based plastic constitutive equation [34] for the conventional model on the premise that the inside of the yield surface is a purely elastic domain, which is applicable only to the description of monotonic loading behavior, will be given based on the equations formulated in the preceding sections.

Flow Rules for Plastic Strain Rate and Plastic Spin
The yield surface with the isotropic and the kinematichardenings is described by Here, let the function f ð MÞ be chosen to be homogeneously degree-one of M. The yield surface in Eq. (37) is shown in Fig. 4. The yield condition in Eq. (37) is described in the current configuration as r is the Cauchy stress, i.e.
The evolution rule of kinematic hardening variable (back stress) a in the infinitesimal elastoplasticity is given by the nonlinear kinematic hardening rule [1,33] as follows: _ e p is the infinitesimal plastic strain rate in the current configuration, noting that the anisotropic hardening is induced only by the deviatoric strain rate, and c k and b k are the material parameters. Note that the kinematic hardening saturates as a ! b kn . Let the plastic strain rate be given by the symmetrized associated flow rule, noting that the plastic strain rate is to be the symmetric tensor, as follows [29,30]: where _ k is the positive plastic multiplier and N is the normalized and symmetrized outward-normal tensor of the yield surface, i.e.
M also holds. The expression in Eq. (44) for the asymmetric tensor M is first adopted by Hashiguchi [29].
The dissipative part in the plastic velocity gradient for the kinematic hardening variable in Eq. (16) is assumed for the nonlinear-kinematic rule in Eq. (42) as follows: Let the plastic spin W p in Eq. (13) and the kinematichardening spin W _ p kd in Eq. (16), which are induced by the plastic strain rate and the dissipative part of kinematic hardening rate, respectively, be given extending the equation in the hypoelastic-based plasticity by Zbib and Aifantis [92] as follows: where g p and g k p are the material parameters, while the flow rules in Eqs. (44) and (46) (44), (46) and (47) hold for the general material with an anisotropy. On the other and, the past formulations of flow rules (e.g. [14,18,46,73,[84][85][86][87] have been concerned only to materials with the symmetry of where N is the normalized outwardnormal to the yield surface. Needless to say, it cannot be allowed to assume easily the flow rule by which the plastic spin is inevitably ignored for the elastically even for the plastically anisotropic material. The velocity gradients are given by substituting Eqs. (44), (46) and (47) into Eqs. (11) and (16) as follows: The flow rule is irrelevant in general. The substitutions of Eq. (48) into Eqs. (29) and (36) yield:

Plastic Strain Rate
The time-differentiation of Eq. (37) leads to the consistency condition as follows: Here, it holds from Eq. (37) that by the Euler's theorem for the homogenous function f ð MÞ of M in degree-one, and then it follows that The substitution of Eq. (53) into Eq. (51) leads tô where F 0 : dF/dH and noting Eq. (44) and the homogeneity of holds for the equivalent plastic strain hardening.
The substitutions of Eqs. (49), (50) and (57) into Eq. (56) from which it follows that where M p N : The substitution of Eq. (49) into Eq. (58) leads to the consistency condition N : using the symbol _ K for the plastic multiplier in terms of the strain rate instead of _ k in terms of the stress rate. The plastic multiplier is expressed from Eq. (61) as follows: The loading criterion is given by which can be given actually as noting the positivity of the denominator in the plastic multiplier in Eq. (62), while the verification of these loading criterion was given by Hashiguchi [24,33] for the current configuration.

Multiplicative Hyperelastic-Based Plastic Equation Incorporating Initial Subloading Surface Model
The conventional elastoplasticity model has been extended to describe the plastic strain rate caused by the rate of stress inside the yield surface, while the extended models are called the unconventional elastoplasticity models by Drucker [15]. The subloading surface model [19,20] as the unconventional elastoplasticity model is based on a quite natural postulate that the plastic strain rate is induced progressively as the stress approaches the yield surface. In this model, the surface, termed the subloading surface, be incorporated, which always passes through the current stress point and maintains the similar shape and same orientation to the yield surface, renamed as the normal-yield surface. Then, let the similarity ratio, i.e. the ratio of the size of subloading surface to that of normal-yield surface, termed the normal-yield ratio and designated by the symbol R ð0 R 1Þ, be introduced as the general measure for approaching degree of stress to the normal-yield surface.

Subloading surface and evolution rule of normal-yield ratio
The subloading surface is represented by the following equation (see Fig. 5).
Assume the following evolution rule of normal-yield ratio [20].
where UðRÞ is the monotonically-decreasing function of the normal-yield ratio fulfilling the conditions (see Fig. 6) ð67Þ R e ð\1Þ is the material constant denoting the value of R below which only an elastic deformation is induced practically. Here, we assume the explicit form of UðRÞ as follows: where u is the material parameter. Equation (66) can be time-integrated analytically as follows:

Plastic Strain Rate
The time-differentiation of the subloading surface in Eq. (65) reads: which can be described aŝ N : The following relation based on the Euler's homogenous function is used for the derivation of Eq. (71) from Eq. (70), noting Eq. (65).
The flow rules for the plastic strain rate, the dissipative part of the kinematic hardening variable and the plastic spin and the spin of the dissipative part of the kinematic hardening are given by Eqs. (44), (46) and (47) themselves. The plastic modulus is given by the following equation instead of Eq. (60).

Fig. 5 Normal-yield and subloading surfaces in intermediate configuration
The loading criterion is given by which can be given actually as in which the fulfillment of the yield condition is not required, although it is required in Eq. (63) or (64) for the conventional model. The subloading surface model possesses the following distinguished abilities.

Smooth transition from elastic to plastic state is
described, which is observed in real material behavior. Therefore, we don't need to suffer from the determination of an offset value (plastic strain value at yield point) influenced by an arbitrariness. In contrast, the abrupt transition from the elastic to the plastic state is depicted and the determination of offset value is required in all of the other models i.e. the conventional model and the cyclic kinematic hardening models (the multi surface model: Mroz [66] and Iwan [51], the two surface model: Dafalias and Popov [12] and Krieg [53] and the superposed-kinematic hardening model: Chaboche et al. [10]) since they assume the small yield surface enclosing the purely-elastic domain. The influence of the material parameter u on the stressstrain curve is depicted in Fig. 7. The larger the material parameter u, the more rapidly the normalyield ratio R increases causing the more rapid increase of stress, i.e. approaching the behavior of the conventional elastoplasticity. 2. The continuity and the smoothness conditions [21,22,24] are satisfied, while these conditions are violated in all of the other models assuming the yield surface enclosing the purely-elastic domain. 3. Plastic strain rate can be described even for low stress level and for cyclic loading process under small stress amplitudes since a purely-elastic domain is not assumed. 4. The yield-judgment whether or not 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, the yield judgment is required in all of the other elastoplastic models since they assume a surface enclosing a purely-elastic domain. 5. The stress is automatically pulled-back to the normalyield surface when it goes out from the surface in numerical calculation because of _ R\0 for R [ 1 from Eq. (66) with Eq. (67) 4 as seen in Fig. 8. In contrast, the particular operation to pull-back the stress to the yield surface is required in all of the other models because they assume a surface enclosing a purelyelastic domain.
For the concise illustration of the above-mentioned feature of the subloading surface model, let the stress Influence of material parameter u on stress-strain curve Stress is automatically controlled to be attracted to yield surface in subloading surface model versus strain relation in the uniaxial loading for the isotropic Mises material with the yield condition ffiffiffiffiffiffiffi ffi 3=2 p r 0 k k ¼ F be shown in the current configuration. The relations of the axial Cauchy stress r a and the normal-yield ratio R versus the axial strain e a : $ d a dt are depicted in Fig. 9. The responses adopting the linear isotropic hardening Fig. 9a and those for the nonlinear isotropic hardening rule FðHÞ ¼ F 0 f1 þ h 1 ½1 À expðÀh 2 e ep Þg are shown in Fig. 9b. The two levels of axial strain increment de a = 0.0006 and 0.0055 are input in the numerical calculations. Here, any special algorithm for pulling back the stress to the yield surface is not introduced. The material parameters are chosen as follows: Material constants:  loading process, while, needless to say, the amplitude is smaller for a smaller input increment of strain. Eventually, the subloading surface model posseses the distinguished high ability for numerical calculation as verified also quantitatively in these concrete examples, which cannot be attained in any other elastoplastic constitutive models including the conventional model and the cyclic kinematic hardening models.

Multiplicative Hyperelastic-Based Plastic Equation Incorporating Extended Subloading Surface Model
The initial subloading surface model is furnished with the distinguished ability improving the conventional elastoplastic model as described in the last section. However, it is incapable of describing the cyclic loading behavior, predicting open hysteresis loops because only elastic deformation is induced in the unloading process. Then, it has been extended such that the similarity-center of the normalyield and the subloading surfaces moves with the plastic deformation, while the similarity-center is regarded as the elastic-core because the most elastic deformation behavior is induced when the stress lies on it, the normal-yield ratio reducing to zero, i.e. the subloading surface reducing to a point.

Basic Feature of Extended Subloading Surface Model in Current Configuration
The uniaxial loading behavior is depicted in Fig. 10 in the current configuration for the Mises material without a hardening (F = F 0 and a ¼ O) for simplicity. Here, a is the conjugate point in the subloading surface to the kinematic hardening variable (back stress) a in the normalyield surface. The elastic-core c goes up following the stress by the plastic strain rate in the initial loading process as seen in Fig. 10a, b. The subloading surface shrinks and thus only elastic strain rate is induced until the stress goes down to the elastic-core in the unloading process as seen in Fig. 10c. After that the subloading surface begins to expand and thus the plastic strain rate in the compression is induced in the unloading-inverse loading process whilst the similarity-center goes down following the stress by the plastic strain rate as seen in Fig. 10d. Again only the elastic strain rate is induced until the stress goes up to the similarity-center in the reloading process from the complete unloading as seen in Fig. 10e. After that the subloading surface begins to expand and thus the plastic strain rate is induced whilst the similarity-center goes up following the stress by the plastic strain rate as seen in Fig. 10f. The expanded figure of Fig. 10f is shown in the lowest part of Fig. 10. Consequently, the closed hysteresis loop is depicted realistically as shown in this figure.
The extended subloading surface model would describe the cyclic loading behavior realistically as illustratively shown in Fig. 11. It does not contain any drawbacks in the cyclic plasticity models based on the kinematic hardening concept, while the continuity and the smoothness conditions [21,22,24] are satisfied only in this model. Then, it has been applied to the descriptions of rate-independent and rate-dependent elastoplastic deformation behavior of not only metals but also soils and further the friction phenomena between solids as will be described in detail in the later sections.
The elastoplastic models other than the subloading surface model, i.e. the conventional model and the cyclic kinematic hardening models are incapable of describing monotonic loading behavior realistically and also cyclic loading behavior appropriately as was described in the foregoing. In addition, it would be incapable of formulating multiplicative hyperelastic-based plastic equation based on them. Only the extended subloading surface model with the translation of the elastic-core is capable of describing the monotonic and cyclic loading behavior rigorously and can be led to the multiplicative hyperelastic-based plasticity as will be described in detail in the subsequent sections.
The normal-yield, subloading and elastic-core surfaces in the current configuration is shown in Fig. 12. The variables in the hypoelastic-based plasticity correspond to the variables in the intermediate configuration for the multiplicative hyperelastic-based plasticity as follows: where c is the material parameters and n is the normalized outward-normal of the subloading surface, i.e. n of ð rÞ or . of ð rÞ a is the conjugate point in the subloading surface to a in the normal-yield surface. The elastic-core surface which passes through the elastic-core and is similar to the normalyield surface with respect to the back-stress a is given as follows (Fig. 12): < c ð0 < c 1Þ is called the elastic-core yield ratio designating the ratio of the size of the elastic-core surface to that of the normal-yield surface, i.e. the approaching-degree of the elastic-core to the normal surface, and vð\1Þ is the material parameter designating the limit value of < c .n c is the normalized outward-normal of the elastic-core surface, i.e. n c of ðĉÞ oc The following inequality holds for Eq. (78).
Therefore, the elastic-core does not go out from the limit elastic-core surface f ðĉÞ ¼ vFðHÞ (see Fig. 12). The validity of the extended subloading surface model in the hypoelastic-based plasticity has been verified widely (cf. e.g. [33,43,45]).

Multiplicative Decomposition of Plastic Deformation Gradient for Elastic-Core
Analogously to the multiplicative decomposition of the plastic deformation gradient for the kinematic hardening in Eq. (3), decompose F p into the plastic storage part F p cs causing the translation of the elastic-core and its plastic dissipative part F p cd multiplicatively as follows [29,30]: The configurations based on these decompositions are illustrated in Fig. 13.
The following tensors of the storage part Ĉ p cs F pT cs F p cs and the dissipative part Ĉ p cd F pT cs F p cd are defined.
where one has C p cs F pÀT cs Ĉ p cs F pÀ1 Tensor variables in the elastic-core intermediate configuration are specified by adding the hat symbol ð^Þ. Further, the following additive decomposition of the velocity gradient holds for the elastic-core analogously to Eqs. (15)- (19) for the kinematic hardening variable. noting The material-time derivative of Ĉ p ks in Eq. (5) is given by Fig. 13 Multiplicative decompositions of deformation gradient tensor for material with translations of kinematic hardening variable and elasticcore Multiplicative Hyperelastic-Based Plasticity for Finite Elastoplastic Deformation… 613

Elastic-Core
The contravariant push-forward of the 2nd Piola-Kirchhoff stress-like variable for the elastic-core Ŝ c from K to K is given by Further, the Mandel-like variable M c for the elastic-core is given by noting Eq. (87). Note here that the Mandel stress M is the asymmetric tensor in general but the Mandel-like elasticcore M c is the symmetric tensor. The material-time derivative of M c is given analogously to Eq. (26) for that of the kinematic hardening as follows:

Hyperelasticity for Elastic-Core
Further, let Ŝ c be formulated incorporating the strain energy function w c ðĈ p cs Þ as from which S c and M c are given by noting Eqs. (87), (93) and (96).
The material-time derivative of Ŝ c is given from Eq. (96) with Eq. (92) as where Substituting Eq. (98) into Eq. (95), _ M c is given as follows: 8.5 Normal-Yield, Subloading and Elastic-Core Surfaces The subloading surface in Eq. (65) for the initial subloading surface model is extended as follows: in the intermediate configuration, which are depicted in Fig. 14.
The subloading surface is given from Eq. (101) with Eq. (77) as follows: from which the normal-yield ratio R is calculated.
leading to The elastic-core surface which passes through the elastic-core M c and is similar to the normal-yield surface with respect to the back-stress M k in the hyperelastic-based plasticity is given noting Eq. (81) as follows:

Plastic Flow Rules
The plastic strain rate is given in the following associated flow rule proposed by Hashiguchi [29,30].
where _ k is the positive plastic multiplier and which is the normalized and symmetrized tensor. If the strain energy functions w e is given by invariants of C e leading to the elastic isotropy, the symmetries of the In order to describe a higher tangent modulus in the reloading process than in the unloading process, e.g. the generalized Masing effect [62], the material parameter u is extended to where u and u c are material constants, < c is given by Eq. (105) and C r is given as follows [27,31]: C r is larger when the direction of the plastic strain rate is nearer to the outward-normal of the elastic-core surface, and thus u is larger in the reloading process and smaller in the reverse loading process. Let the plastic spin W p ; the kinematic-hardening spin W _ p kd and the elastic-core spin Ŵ p cd be given extending the equation in the hypoelastic-based plasticity by Zbib and Aifantis [92] and incorporating Eqs. (106), (108) and (109) as follows: where g c p is the material parameter.

Plastic Strain Rate
The elastic constitutive equation is given by Eqs. (27)- (29). The plastic strain rate will be formulated in this section. The formulations given in this section is not necessary in the numerical calculation by the return-mapping based on the closet-point projection.
The time-differentiation of Eq. (101) leads to the consistency condition of the subloading surface as follows: It holds from Eq. (101) that by the Euler's theorem for the homogenous function f ð MÞ of M in degree-one, and then it follows that N : The substitution of Eq. (120) into Eq. (118) leads to The further substitution of Eq. (104) into Eq. (122) leads to N : _ M À N : Furthermore, substituting the relation Equation (123) is rewritten as N : _ M À N : where based on Eqs. (57) and (66) from which it follows that where M p N : The substitution of Eq. (115) into Eq. (128) leads to the consistency condition: N : using the symbol _ K for the plastic multiplier in terms of the strain rate instead of _ k in terms of the stress rate. The plastic multiplier is given from Eq. (131) as follows: The loading criterion is given by which can be given actually as

Material Functions
Material functions contained in the constitutive equations formulated in the preceding sections are given for metals and soils in this section.

Metals
The hyperelastic equation and the yield functions for metals are shown below in the intermediate configuration.

Hyperelastic Equations
The hyperelastic equations are given for the stress, the kinematic hardening variable and the elastic-core.

Stress
The following strain energy function of elastic deformation in the Modified Neo-Hookean elasticity may be adopted [84,85].
where l and K are the material constants. The substitution of Eq. (135) into Eq. (27) reads: It is followed from Eq. (136) that where g ¼ F eÀT GF eT is the metric tensor in the current configurations.
The hyperelastic tangent modulus tensor C e in Eq. (31) is given for Eq. (136) by where C e is the fourth-order tensor defined as which possesses the minor symmetry C e ijkl ¼ C e jikl ¼ C e ijlk and the major symmetry C e ijkl ¼ C e klij .

Kinematic Hardening Variable
Assume the following strain energy function for kinematic hardening, which possesses the identical form to the shear part in Eq. (135).
where C k is material constant. It is derived from Eq. (32) that Then, the Mandel-like kinematic hardening variable is given from Eqs. (33) and (143) as follows: (35) is given for Eq. (143) as 9.1.1.3 Elastic-Core Assume the following strain energy function for elastic-core analogously to the kinematic hardening variable described in 9.1.2.
where C c is material constant. The following relations hold.

Yield Functions
Assume the von Mises yield function and the plastic equivalent hardening, i.e.
for which we have ffiffi ffi 3 2 from which the normal-yield yield ratio is given by

Soils
The hyperelastic equation and the yield functions for soils are shown below in the intermediate configuration.

Hyperelastic Equation
The hyperelastic constitutive equation of soils was proposed first by Houlsby [48] and it has been extended to various multiplicative hyperelastic equations by Borja and Tamagnini [5], Callari et al. [9], etc. using the Hencky (logarithmic) strain in the current configuration and by Yamakawa et al. [90] using the second Piola-Kirchhoff stress in the intermediate configuration, while they involve various impertinences or incompleteness. The hyperelastic equation within the framework of the multiplicative finite strain theory for soils was shown by Hashiguchi [33] but it involves some impertinences. The exact equation will be shown in the following [37]: The isotropic hyperelastic equation is given by introducing the function of the variables lnJ e and tr C e (tr C e À 3 in detailed expression) which stand for the volumetric strain and deviatoric strain, respectively, as follows: where F e vol is the volumetric part and F e is the so-called unimodular tensor designating the isochoric (constant volume, i.e. deviatoric) part of F e . In addition, tr C e ¼ 3 (F e ¼ F e vol Þ and ln J e ¼ 0 ðF e ¼ F e Þ is required in the purely volumetric deformation and the purely deviatoric deformation, respectively.
The following partial derivatives hold.
The substitution of Eq. (159) into Eq. (157) reads: Let the following strain energy function be assumed. noting where P M0 is the initial value of the pressure defined in terms of the Mandel stress M, i.e. P M À(1/3)tr M.k and j are the material constants standing for the inclinations of the isotropic and the swelling lines in the both logarithmic linear relation of pressure and the volume (cf. [33]). #(\1/2) is material constant, while the volume becomes infinite as the pressure approaches #F, while the hardening function F coincides to the preconsolidation pressure in the isotropic consolidation process. G 0 is the initial value of elastic shear modulus, n standing for the pressure-dependence of the shear modulus.
The following partial derivatives hold for Eq. (161).
noting oJ e Àn=j =o ln J e ¼ Àðn=jÞJ eÀn=j . Equation (160) with Eq. (162) reads: from which the Mandel stress is given as follows: noting J eÀ2=3 C e0 ¼ C e0 by virtue of Eq. (158) 7 . It is follows from Eq. (164) that i.e. lnJ e ¼ Àj ln and which would describe appropriately the basic characteristics in the volumetric and the deviatoric deformations. Equations (164), (166) and (167) are rewritten as with lnJ e ¼ Àj ln and in the current configuration, where s is the Kirchhoff stress tensor, p s0 is the initial value of p s ÀðtrsÞ=3 and b e is the elastic unimodular left Cauchy-Green deformation tensor, i.e.
The above-mentioned elastic equation is reduced to the infinitesimal strain theory by adopting the strain energy function as follows: resulting in and where p is the pressure, i.e. p Àðtrr)/3 and its initial value is denoted by p 0 . e e is the infinitesimal elastic strain tensor and e e v ffi tre e ; e e d ¼ e e0 k k.
The elastic balk modulus K and the elastic shear modulus G adopted in the hypoelasticity for the above-mentioned elastic equations are given as follows: where d e v trl e ; d e sym½l e , while n ffi 0.5 can be chosen in most of soils [82].
The elastic constitutive equations of soils formulated above possess the following physical validities.
1. It is applicable up to the finite deformation/rotation. 2. It is applicable up to the negative pressure range which depends on the preconsolidation pressure. 3. The shear modulus increases depending on the pressure. 4. They are consistent for the multiplicative hyperelasticity, the infinitesimal hyperelasticity and the hypoelasticity.

Yield Functions
The modified Cam-clay yield surface [8,72], which is exhibited as the ellipsoid passing through the origin and possessing the long axis coinciding to the hydrostatic axis in the three-dimensional stress space, was extended by Hashiguchi and Mase [39] as follows (Fig. 15): where F is the isotropic hardening function corresponding to the long axis of the ellipsoidal yield surface, n h ð\#Þ is the material constant describing the translation to the negative pressure range by n h F along the hydrostatic direction and M is the ratio of the short axis to the long axis, which is the function of the Lode angle as follows [25]): where / c is the friction angle in the tri-axial compression sate, while M c 2 ffiffi ffi 6 p sin / c =ð3 À sin / c Þ is the value of M in that state (h M = p/3), and Equation (178) is expressed in the separated form of the function f ð P M ; XÞ of the stress and the hardening function wherẽ The isotropic hardening function is given as [23].
where J p ¼ det F p and F 0 is the initial value of F.

Calculation Procedures
The calculation procedure by the above-mentioned formulations is described in this section.
First, the plastic multiplier _ K is calculated by the input of the velocity gradient L into Eq. (62), while L is calculated from the current velocity gradient l by Eq. (10). Then, substituting it into Eq. (114), the plastic and the dissipative parts L p , L p kd and L p cd are calculated. On the other hand, the plastic multiplier _ K is calculated directly from the plastic flow rule in Eq. (106) and the stress rate is calculated by Eq. (29) under L ¼ O in the plastic corrector process in the return-mapping method. Thereafter, the stress and the tensor-valued internal variables are calculated by the process described below.
The deformation gradient tensor is updated by where with the displacement vector u, designating r x n oð Þ=ox n and noting The rates of the plastic gradient and its dissipative parts are given from Eqs. (10) 3 , (18) and (91) as follows: where L p , L p kd and L p cd are given by Eq. (114). The storage parts F e ; F p ks and F p cs of the deformation gradient are given by substituting the results of the time-integrations of Eq. (188) into The plastic constitutive equation with the plastic modulus in Eq. (130) is not necessary to be used in the numerical calculation by the return-mapping in which the plastic strain rate is calculated by use of only the plastic flow rule in Eq. (106) and then the stress and internal variables are calculated by the procedures described above.
The time-integrations of Eq. (188) for the tensors F p , F p kd and F p cd can be executed in high efficiency by the tensor exponential method [46,65,88] which is delineated below.
The equations in Eq. (188) are collectively described in terms of an arbitrary second-order tensors T and Z as follows: Consider the following candidate as the numerical solution of Eq. (190) for the time-interval Dt ¼ ½t n ; t nþ1 .
provided that Z and T n are constant during the time-interval. The time-differentiation of Eq. (191) leads to

Loading Criterion in Return-Mapping for Subloading Surface Model
The subloading surface model is capable of describing not only the monotonic but also the cyclic and the non-proportional loading behaviors, since the yield surface enclosing the purely-elastic domain is not incorporated in this model. However, the particular attention is required for the return-mapping method in numerical calculations for this model, in which rather large loading increments are input in the elastic trial step. In fact, there exists the possibility of the occurrence of plastic strain increment even when the stress increment is directed inward of the subloading surface. Therefore, a particular loading criterion after the elastic trial step and a particular calculation procedure in the initiation of plastic corrector step are required in the return-mapping projection for the subloading surface model. This fact has not been recognized in the past formulations [2,16,28,33,46,90] and the inexact loading criterions after the elastic trial step have been shown in the past [32,49]. The closet point projection in the return-mapping method for the extended subloading surface model was proposed first by Anjiki et al. [2] and subsequently Iguchi et al. [49] in the multiplicative decomposition and followed the results of Anjiki et al. [2] later by Fincato and Tsutsumi [16]. The exact loading criterion after the elastic trial step and the initial treatment in the plastic corrector step for the subloading surface model will be formulated below extending the former equations in the current configuration [36] to the equations in the intermediate configuration within the framework of the multiplicative elastoplaticity.

Formulation of Loading Criterion
Analogously, the following equation holds at the end of the elastic trial step for the Mises material.
which is rewritten as from which R trial nþ1 is given by The plastic corrector step is executed until the subloading surface equation is satisfied within the prescribed tolerance. The plastic corrector step is started for the difference from the subloading surface f ð M n Þ ¼ R n FðH n Þ in the end of the step n to the subloading surface f ð M trial nþ1 Þ ¼ R trial nþ1 FðH n Þ in the elastic trial step for the case 1) ii). On the other hand, it must be executed for the case 2) ii) b) by the special procedures described in the following.
The subloading surface once shrinks inducing only elastic strain increment and turns to expand inducing the plastic strain increment after it contacts tangentially to the stress increment D M trial nþ1 in the process of the stress variation from M n as shown in Fig. 16. Then, we need to calculate the normal-yield ratio R 0 nþ1 in the state that the subloading surface contacts tangentially to the line-element. Here, the following relations must be satisfied at the contact point, where the stress at the contact point, i.e. the plastic-loading initiation stress is denoted by M 0 nþ1 from which the subloading surface changes from the contraction to the expansion.
s(0 s 1) is the unknown scalar parameter which must be determined so as to satisfy Eq. (207) at the stress M 0 nþ1 . The two unknown variables s and R 0 nþ1 are involved in Eq. (207). In what follows, the explicit calculation procedure of them will be shown for the Mises material.
Equation (207) is explicitly described for the Mises as follows: The upper equation in Eq. (213) is expressed as The substitution of Eq. (214) into Eq. (213) 2 leads to resulting in The solution of R 0 nþ1 in Eq. (215) is given by where A S 2 ca À S ss M The plastic strain increment is induced during the variation of normal-yield ratio from R 0 nþ1 to R nþ1 which is related to the plastic strain increment e p À e p 0 as follows: The plastic strain increment in the initiation of the plastic corrector step is calculated based on the overstress M

Subloading-Overstress Model Based on Multiplicative Decomposition
The deformation gradient F is multiplicatively decomposed into the elastic deformation gradient F e and the plastic deformation gradient F vp instead of the plastic deformation gradient F p in the multiplicative elastoplasticity described in the preceding sections. Then, we first adopt the following equation instead of Eq. (1).
Further, the viscoplastic deformation gradient F vp is multiplicatively decomposed into the viscoplastic storage part F vp ks causing the kinematic hardening and its dissipative part F vp kd and into the viscoplastic storage part F vp cs causing the elastic-core and its dissipative part F vp cd multiplicatively as follows: Then, the following right Cauchy-Green deformation tensors for the viscoplastic deformation are introduced.
The velocity gradient l in the current configuration is additively decomposed into the elastic and the viscoplastic parts: where l _ FF À1 ; Further, the velocity gradient L in the intermediate configuration is additively decomposed into the elastic and the plastic parts as follows: where from which it follows that where The viscoplastic velocity gradient L vp is additively decomposed for the kinematic hardening and the elasticcore into the storage and the dissipative parts as follows: The viscoplastic strain rate is given extending the overstress model [70,71] as follows: where N is given by Eq. (107), l, n and R m ([ 1) are the material parameters, while R m is the maximum value of R, and The viscoplastic spin W p , the kinematic-hardening viscoplastic spin W p kd and the elastic-core viscoplastic spin W p cd are given extending Eq. (113) as follows: where g vp , g vp k and g vp c are the material parameters. The velocity gradients are given by substituting Eqs. (237) or (239), (243) and (245) into Eqs. (227) 3 , (231) 2 and (234) 2 as follows: The velocity gradient is given by Eq. (225) with Eqs. (226) and (246) as follows: which is expressed in the incremental form as follows: The rate of the viscoplastic gradient is given from Eq. (224) 3 as follows: The time-integration for F vp can be performed effectively by the tensor exponential method as described for the plastic deformation gradient F p in the end of Sect. 10. The elastic deformation gradient F e is given by substituting the time-integration F vp in Eq. (188) into Then, C e is calculated by Eq. (5) and further the stresses S and M are calculated by Eqs. (27) and (28)

Hyperelastic-Based Plastic Constitutive Equation for Subloading-Friction Model
The reduction of the friction coefficient from the static to kinetic friction coefficient and the recovery of the friction coefficient under the reduction of sliding velocity are the fundamental characteristics in friction phenomena, which are recognized widely. Difference of the static and kinetic frictions often reaches up to several tens of percent. The formulation of friction phenomenon taken account of these fundamental friction behavior is of importance for accurate analyses of practical problems in engineering. The constitutive equation of friction describing these behavior has been formulated based on the subloading surface concept [19] and thus it is called the subloading-friction model [33,40,41,44]. However, it has been formulated based on the hypoelastic-based plasticity (e.g. [7,33,40,41,69]), in which the elastic sliding velocity is limited to be infinitesimal and the cumbersome time-integration of corotational contract stress rate is required.
The subloading-friction model will be improved based on the hyperelastic-based plasticity in this section. It is capable of describing exactly the finite sliding behavior under the finite rotation of contact surface without the cumbersome time-integration of corotational contact stress rate by formulating to relate the contact stress directly to the elastic sliding displacement by the hyperelastic equation. Further, it is extended to be applicable to the general case that the contact surface undergoes the rotation and the deformation.

Sliding Displacement and Contact Traction Vectors
The sliding displacement vector u, which is defined as the sliding displacement of the counter (slave) body to the main (master) body, is orthogonally decomposed into the normal sliding displacement vector u n and the tangential sliding displacement vector u t to the contact surface as follows: where u n ¼ ð u Á nÞn ¼ ðn nÞ u ¼ À u n n u t ¼ u À u n ¼ ðI À n nÞ u ( ð252Þ n being the unit outward-normal vector of the surface of main body and u n Àn Á u n ¼ Àn Á u ð253Þ Furthermore, u is decomposed into the elastic sliding displacement u e and the plastic sliding displacement u p as follows: where u e n ¼ ð u e Á nÞn ¼ ðn nÞ u e ¼ À u e n n u e t ¼ u e À u e n ¼ ðI À n nÞ u e ( ð256Þ The minus sign is added for u e n to be positive when the counter body approaches the main body. On the other hand, the plastic sliding flow rule will be formulated to fulfill u p n ¼ 0 in Sect. 13.6.
The contact traction vector f acting on the main body is additively decomposed into the normal traction vector f n and the tangential traction vector f t as follows: The minus sign is added for f n to be positive when the compression is applied to the main body by the counter body.
The contact traction vector f, f n and f t are calculated from the Cauchy stress r applied to the contact surface by virtue of the Cauchy's fundamental theorem (cf. [33]) as follows:

Elastic Sliding Displacement
Let the contact traction vector f be given by the hyperelastic relation with the elastic sliding displacement energy function uð u e Þ as follows: The simplest function ouð u e Þ is given by the quadratic form: where the second-order symmetric tensor E is the elastic contact tangent modulus tensor fulfilling the symmetry E ¼ E T . The substitution of Eq. (264) into Eq. (263) leads to The inverse relation of Eq. (265) is given by The elastic contact tangent modulus tensor E is given for the isotropy on the contact surface as follows: E ¼ a n n n þ a t ðI À n nÞ E À1 ¼ 1 a n n n þ 1 a t ðI À n nÞ where a n and a t are the normal and tangential contact elastic moduli, respectively. Their values are quite large usually as 10 2 À 10 5 GPa/mm 3 for metals because the elastic sliding is caused by elastic deformations of the surface asperities. Equations (265) and (266) with Eq. (267) leads to Now, introduce the normalized rectangular coordinate system ðê 1 ;ê 2 ;ê 3 Þ ¼ ðê 1 ;ê 2 ; nÞ fixed to the contact surface, which changes with the rotation of the contact surface. The elastic sliding displacement and the contact traction are described as follows: Hence, Eq. (265) is described in the simple form as follows: The sliding velocity vector _ u is the objective vector, since it is not an absolute velocity vector but the mutual velocity vector between surface points on the master and the counter bodies. Therefore, it is not necessary to use a corotational velocity vector but we only have to use the time derivative for the sliding velocity vector. Further, note that one does not need to adopt a corotational rate but one has only to use the time derivative for the contact traction vector f by the fact: The contact traction f is calculated from the hyperelastic equation with the substitution of the elastic displacement u e which is obtained by subtracting the plastic displacement vector u p from the displacement vector u as will be explained in Sect. 13.7.

Normal Sliding-Yield and Sliding-Subloading Surfaces
Assume the following sliding-yield surface with the isotropic hardening/softening, which describes the slidingyield condition.
f ðfÞ ¼ l ð271Þ l is the isotropic hardening/softening function denoting the variation of the size of the sliding-yield surface. The friction-yield stress function f ðfÞ for the Coulomb friction law is given by for which l specifies the coefficient of friction. Then, in order to introduce the measure of approaching degree to the sliding-yield surface, renamed the normal sliding-yield surface, let the following subloading-sliding surface passing through the current contact stress and maintaining a similarity to the normal sliding-yield surface be introduced, which plays the general measure of approaching degree of the contact stress to the normal sliding-yield surface (see Fig. 17).
f ðfÞ ¼ rl ð273Þ where rð0 r 1) is the ratio of the size of the subloading surface to that of the normal sliding-yield surface and called the normal sliding-yield ratio, playing the role of the measure of the approaching degree of the contact stress to the normal sliding-yield surface.

Evolution Rule of Normal Sliding-Yield Ratio
The evolution rule of the normal sliding-yield ratio is given as follows [33]: where m p Àj l l k À 1 r þ UðrÞl; m c n 1 À l l s r ð ! 0Þ ð286Þ are relevant to the plastic and the creep sliding velocity, respectively. It is obtained from Eqs. (281) and (285) that Substituting the rate form of Eqs. (266) and (287) into the rate form of Eq. (254), the sliding velocity is given by The plastic multiplier in terms of the sliding velocity, denoted by the symbol _ K, is given from Eq. (288) as The inverse relation Eq. (288) is given by substituting the rate form of Eq. (254) with Eq. (289) into the rate form of Eq. (265) as follows: The loading criterion is given as follows: Some numerical experiments for the linear sliding process for subloading-friction model formulated above are shown in the following. The seven material constants of l s , l k , j, n,ũ, a n and a t and the initial value l 0 of the friction coefficient are chosen as follows: l 0 ¼ l s ¼ 0.4, l k ¼ 0.2 j ¼ 10 mm À1 ; n ¼ 0.01/s u ¼ 1000 mm À1 a n ¼ a t ¼ 1000 kN/mm 3 The influence of the sliding velocity on the relation of the traction ratio f t =f n versus the tangential sliding displacement u t are shown in Fig. 18, where _ u t k _ u t k. Smooth transitions from the static friction to the kinetic friction and the decreases of the friction coefficient are depicted. Faster decrease of friction coefficient is shown for higher sliding velocity.
The recovery of the static friction coefficient from the kinetic friction with the elapsed time t after the stop of sliding is shown in Fig. 19. In the calculation, the constant sliding velocity _ u t ¼ 0:1 mm/s is given in the first stage reaching the kinetic friction and then the tangential contact traction is unloaded to zero. After the cessation of sliding for several elapsed times, the same sliding velocity in the first stage is given again. The recovery is larger for a longer stationary time.

Calculation Procedure
The exact calculation of the contact stress can be performed through the hyperelastic relation, while the cumbersome calculation procedure for the time-integration of the corotational contact stress rate in the hypoelastic relation is not required. It will be shown for the sliding processes under the rotation and under the rotation/ deformation of the contact surface in the following.

Sliding Process Under Rotation of Contact Surface
The sliding displacement is calculated by the time-integration of the input value of the sliding rate, i.e. u ¼ R _ udt. The plastic sliding displacement is calculated by the timeintegration of the plastic sliding rate, i.e.
The elastic sliding displacement U e in the reference configuration is calculated by the following equation, from which the elastic sliding displacement ue in the current configurations is known.
The deformation behaviour described by the subloadingoverstress friction model is schematically shown in Fig. 21 in which the normalized overstress divided by the normal contact stress f n is written simply by the term ''overstress''. The dynamic loading friction-yield ratio r coincides with the subloading friction-yield ratio r s , i.e. r ¼ r s in quasistatic sliding and increases above r s with increasing sliding velocity. However, r does not rise above r m , i.e. r r m , with the equality r ¼ r m being realised only for impact sliding.
The decrease of friction resistance with the increase of sliding velocity is observed in the dry friction and called the negative rate-sensitivity. Inversely, the increase of friction resistance with the increase of sliding velocity is observed in the fluid friction and called the positive ratesensitivity.

Concluding Remarks
The formulation of the exact multiplicative hyperelasticbased plastic constitutive equation has been desired earnestly during a half century after the proposition of the multiplicative decomposition of the deformation gradient tensor in 1960s. Now, it would have been attained by incorporating the rigorous plastic flow rules and the subloading surface model in this article.
The main results in this article are summarized in the following.
1. The exact multiplicative hyperelastic-based plastic constitutive equations are formulated for the conventional elastoplastic model with the yield surface enclosing the purely elastic domain, the initial subloading surface model and the extended subloading surface model which is capable of describing not only the monotonic but also the cyclic loading behavior under the finite deformation and rotation. Among various cyclic plasticity models only the subloading surface model can be formulated based on the framework of the multiplicative elastoplasticity. 2. These equations are extended to the description of the rate-dependent deformation behavior by incorporating the notion of the overstress. They are applicable to the deformation behavior at the general rate ranging from the quasi-static to the impact deformation behavior. 3. The hyperelastic-based plastic constitutive equation for the subloading-friction model is formulated by extending the hypoelastic-based plastic one. It is applicable to the sliding process in which the contact surface undergoes the rotation. It is further extended to be applicable to the sliding process in which the contact surface undergoes not only the rotation but also the deformation. 4. These equations possess the distinguished ability in numerical calculations by the forward-Euler method such that the stress is automatically pulled-back to the yield surface when it goes out from the yield surface, which is the inherent ability of the subloading surface model. 5. The general loading criterion is formulated, which is applicable not only to the monotonic loading process but also to the reverse loading process, and further the initial treatment in the plastic-corrector step is formulated for the return-mapping method adopting the subloading surface model.
These formulations possess the high numerical efficiency in addition to the rigorous description of mechanical behaviors up to the finite deformation/sliding and rotation. They have been long-awaited in the history of elastoplasticity theory. They will highly contribute to the steady developments of the continuum mechanics and the mechanical design technique in industries.
Acknowledgements The author would like to express his sincere gratitude to Professor Genki Yagawa (Emeritus Professor, University of Tokyo and Toyo University) for encouraging and inviting the author to contribute to this journal, highly appreciating this monograph on the multiplicative hyperelasticity-base plasticity incorporating the subloading surface model. He is also indebted to Prof. Yuki Yamakawa, Tohoku University for valuable comments and suggestions on the finite strain elastoplasticity.

Compliance with Ethical Standards
Ethical Standards The author declares that this article complies the ethical standard.

Conflict of interest
The author declares that he has no conflict of interest. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.