Elastoplastic model of metals with smooth elastic–plastic transition

The subloading surface model is based on the simple and natural postulate that the plastic strain rate develops as the stress approaches the yield surface. It therefore always describes the continuous variation of the tangent modulus. It requires no incorporation of an algorithm for the judgment of yielding, i.e., a judgment of whether or not the stress reaches the yield surface. Furthermore, the calculation is controlled to fulfill the consistency condition. Consequently, the stress is attracted automatically to the normal-yield surface in the plastic loading process even if it goes out from that surface. The model has been adopted widely to the description of deformation behavior of geomaterials and friction behavior. In this article, it is applied to the formulation of the constitutive equation of metals by modifying the past formulation of this model. This modification enables to avoid the indetermination of the subloading surface, to make the reloading curve recover promptly to the preceding loading curve, and to describe the cyclic stagnation of isotropic hardening. The applicability of the present model to the description of actual metal deformation behavior is verified through comparison with various cyclic loading test data.


Introduction
The conventional elastoplasticity model [9], premised upon the idea that the interior of a yield surface is an elastic region, has contributed to the prediction of the elastoplastic deformation behavior of solids such as metals, geomaterials, and concretes. However, it cannot describe cyclic loading behavior, the application being limited to the description of monotonic loading behavior. Various elastoplasticity models aimed at describing cyclic loading behavior have been proposed to date. They require a pertinent description of the plastic strain rate induced by the rate of stress inside the yield surface. They are termed the unconventional plasticity model by Drucker [9] and also designated as the cyclic plasticity model [32]. They are classified into models based on the concepts of kinematic hardening, i.e., the "translation" (parallel displacement) of (sub)yield surface enclosing a purely elastic region and the concept of expansion/contraction, i.e., "size-variation" of the loading surface. The multi surface model [36,43], the two surface model [5,40], and the nonlinear kinematic hardening model [1,3,4] belong to the former category, and only the subloading surface model belongs to the latter category.
The direct notations AB for A ir B r j , A : B for A rs B rs , : A for Γ i jrs A rs and A : for A rs Γ rsi j are used for arbitrary second-order tensors A, B and fourth-order tensor , where Einstein's summation convention is applied for a letter of the repeated suffix taking 1, 2, 3.

Strain rate
Designating the current configuration of the material particle as x and the current velocity as v, the velocity gradient is described as L = ∂v/∂x by which the strain rate and the continuum spin are defined through D ≡ (L + L T )/2 and W ≡ (L − L T )/2, respectively, () T standing for the transpose. Limited to the infinitesimal strain, let the strain rate D be decomposed additively into the elastic strain rate D e and the plastic strain rate D p as where in the context of the infinitesimal elastic deformation D e be given as the hypoelastic equation [48] in the Hooke's type as where E is the fourth-order tensor describing the elastic tangent modulus and σ is the Cauchy stress, with ( • ) denoting the proper objective corotational rate tensor for an arbitrary second-order tensor T, where ω is the spin of the rigid-body rotation of the material, i.e., the spin of the substructure of the material. The continuum spin W can be used as the corotational spin ω for infinitesimal deformation, whereas W is calculated from the velocity gradient. The elastic tangent modulus tensor E in Hooke's form is . (4) E and ν are the Young's modulus and shear modulus, respectively. δ i j is the Kronecker's delta, i.e., δ i j = 1 for i = j and δ i j = 0 for i = j. Here, the hypoelastic-plastic constitutive equation incorporating proper corotational rates for stress and anisotropic hardening variables holds only for infinitesimal deformation but it holds even for the finite rotation.

Refinement of formulation in the subloading surface model
The discussion of this Section is devoted to the formulation of the constitutive equation of metals, modifying the past formulation of the subloading surface model by avoiding the indetermination of the subloading surface, by making the reloading curve recover promptly to the preceding loading curve and by introducing the cyclic stagnation of isotropic hardening. The conventional yield surface, renamed as the normal-yield surface, for metals is given as whereσ ≡ σ − α. (6) F is the isotropic hardening function exhibiting the size of the normal-yield surface, which is the function of the isotropic hardening variable H . α is the kinematic hardening variable, i.e., the back stress. Here, it is assumed that f (σ) is the function ofσ in the homogeneous degree-one fulfilling f (|s|σ) = |s| f (σ) for an arbitrary scalar variable s. Therefore, the normal-yield surface maintains a similar shape.
To describe the plastic strain rate induced by the rate of stress inside the yield surface, it is assumed that the plastic strain rate is induced progressively as the stress approaches the yield surface. First, let the surface, called 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 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, called the normal-yield ratio and designated by the symbol R(0 ≤ R ≤ 1), be introduced as the general measure of approaching to the normal-yield surface. The subloading surface is represented by the following equation (see Fig. 1): whereσ ≡ σ −ᾱ = R(σ y − α) =σ + Rŝ which is rewritten asN noting the following Euler's theorem for the homogeneous function in degree-one ofσ : designates the magnitude, i.e., T = √ T : T. The variable • α in Eq. (10) is described from Eq. (8.3) as The most elastic deformation behavior is induced in the state in which the normal-yield ratio is zero, i.e., R = 0 in which the stress lies on the similarity-center, i.e., σ = s. Consequently, the similarity-center s is interpreted as the most elastic stress state. Here, the most elastic stress state, i.e., the similarity-center s approaches the normal-yield surface, following stress σ as the plastic deformation proceeds. However, from physical perspective, the similarity-center should not approach the normal-yield surface without limitation, although the small yield surface enclosing a purely elastic region is allowed to go up unlimitedly to the fully plastic state in the other cyclic plasticity models based on the concept of kinematic hardening. In addition, from the mathematical perspective, the subloading surface is not determined uniquely if the stress coincides with the similarity-center lying just in the normal-yield surface.
To formulate the evolution rule of the similarity-center to avoid the unlimited approach of the similarity-center to the normal-yield surface, first, let the following surface, called the similarity-center surface, be introduced as presented in Fig. 1, which passes through the similarity-center s and has a similar shape and orientation to the normal-yield surface with respect to the back stress α, where the variable s (0 ≤ s ≤ 1) represents the ratio of the size of the similarity-center surface to that of the normal-yield surface, called the similarity-center yield ratio. It plays the role of the measure for the approaching degree of the similarity-center to the normal-yield surface. Since the similarity-center must lie inside the normal-yield surface as described above, the similarity-center yield ratio must be less than unity. Then, the following inequality must hold: where χ (<1) is the material constant exhibiting the maximum value of s . The time differentiation of Eq. (15) at the limit state that s lies on the outermost surface in which the relation {∂ f (ŝ)/∂ŝ}:ŝ(= f (ŝ)) = χ F for s = χ on account of Euler's homogeneous function f (ŝ) in degree-one ofŝ is taken into account. Let the following equation be assumed, which fulfills the inequality (16): where it holds that c is the material constant affecting the rate of similarity-center. s y is the conjugate point on the normal-yield surface to the point s on the similarity-center surface with respect to α, i.e., The similarity-center approaches the current stress for s = 0 under the non-hardening state • α = 0. The fulfillment of inequality (16) by Eq. (17) is justified by the fact that the vector σ y −s y forms an obtuse angle to the outward-normal vector ∂ f (ŝ)/∂ŝ of the similarity-center surface leading to {∂ f (ŝ)/∂ŝ} : (σ y − s y ) ≤ 0 in the six-dimensional stress space as presented in Fig. 1. It holds from Eqs. (8.4) and (19) that Inserting this relation in Eq. (17), the translation rule of the similarity-center is given by Substituting Eq. (20) into Eq. (13), one obtains Inserting Eq. (21) into Eq. (10), it follows that N : Making use of the relationsσ Furthermore, let the evolution rule of the normal-yield ratio R be given by where U is a monotonically decreasing function of R, satisfying the following conditions: where u is a material function in general and is the Macaulay bracket implying s = (s + |s|)/2, i.e., s = s for s ≥ 0 and s = 0 for s<0 for an arbitrary scalar variable. Equation (24) with Eq. (26) can be integrated analytically as under the initial condition ε p = ε p 0 : R = R 0 , where ε p ≡ D p dt (t: time). Substituting Eq. (24) into Eq. (23), the following consistency condition is obtained: In what follows, we introduce the associated flow rule where λ is the proportionality factor. Inserting Eq. (29), Eq. (28) is rewritten as N: It holds from Eq. (30) that λ =N : The strain rate is given from Eqs. (1), (2), and (33) as from which the proportionality factor described in terms of the strain rate, denoted by Λ instead of λ, in the flow rule (29) is given as follows: Using Eq. (35), the stress rate is given from Eq. (34) as follows: The loading criterion is given by Hashiguchi [21], Hill [34,35]: where the judgment of yielding is not necessary because the current stress always lies on the subloading surface, which plays the role of a loading surface, although a judgment is necessary in the other plasticity models.
For the normal-yield state R = 1(ᾱ = α, U = 0), the plastic strain rate in Eq. (33) is reduced to (38) which is no more than the plastic strain rate in the conventional plasticity model with the isotropic and kinematic hardening. For u → ∞ leading to the abrupt decrease in the function U from U → ∞ for R < 1 to U = 0 for R = 1 in Eq. (26), the plastic modulus M p in Eq. (38) drops abruptly from the infinite value to the value M p in Eq. (38) so that the present model behavior is reduced to the conventional elastoplasticity model behavior, thereby exhibiting an abrupt transition from the elastic to the plastic state. However, as u becomes smaller, a more gentle transition from the elastic to the plastic state is described. Therefore, u plays the role to alleviate the abrupt transition from the elastic to the plastic state. The plastic strain rate formulated above depends only on the stress rate component normal to the yield/subloading surface, called the normal stress rate, but it is independent of the component tangential to the yield/subloading surface, called the tangential stress rate, because it is derived simply based on the consistency condition. However, the inelastic strain rate induced by the tangential stress rate, called the tangential-inelastic strain rate [19], cannot be ignored in the non-proportional loading process, which is necessary also for the analysis of plastic instability problems [23,25,26,29,38,39].
where • σ n ≡ (n ⊗n ) : I andĪ denote the identity tensors in the second-order and the fourth-order, respectively. () in the second-order tensor implies the deviatoric part, i.e., is the elastic shear modulus, and ξ and τ (≥1) are material constants. The tangential-inelastic strain rate has no loading criterion and is generated if and only if the deviatoric tangential stress rate is added, falling within the framework of hypoelasticity for which the complete integrability condition does not hold, so that the time integration depends on the loading path. Here, it is generated gradually with the increase in the normal-yield ratio R as sown in the third terms of the right-hand sides in Eqs. (39) and (40), so that both the continuity and the smoothness conditions are fulfilled always. In contrast, not only the smoothness but also the continuity conditions are violated, resulting in losing the uniqueness of the solution in the other plasticity models if the tangential-inelastic strain rate is incorporated because it is predicted suddenly when the stress reaches the surface enclosing the purely elastic region. Therefore, no other cyclic plasticity model has been extended to describe the tangential-inelastic strain rate.

Enhancement of the description of cyclic loading behavior
The unique relation ε p − ε (24) is the function of only the normal-yield ratio R as shown in Eq. (27) as an example. Therefore, ε p induced during a certain change of R in the monotonic loading process is identical irrespective of loading processes, e.g., the initial loading, the reloading and the inverse loading, and of the proportional and non-proportional loadings. This fact engenders the description that the returning of the reloading stress-strain curve to the previous loading curve is unrealistically gentle. Therefore, it engenders the impertinent prediction of cyclic loading behavior, i.e., the prediction of the unrealistically large plastic strain accumulation during the cyclic loading process. The above-described insufficiency in the subloading surface model for the description of cyclic loading behavior has been criticized intensely by Dafalias [7] stating that the subloading surface model fits the description of soil behavior but is inapplicable to that of metal behavior predicting "quite unrealistic, basically due to the strong undershooting." Here, it is noteworthy that the bounding surface model with a radial mapping proposed later by Dafalias and Herrmann [6] has a structure similar to the subloading surface, as was written "It appears that the first time a radial mapping formulation was proposed, it was in reference to granular materials by Hashiguchi and Ueno [12]" by Dafalias [7]. However, the above-described insufficiency in the description of deformation behavior by the present formulation in the subloading surface model would not originate from the intrinsic nature of this model contrary to the assessment by Dafalias [7]. In what follows, the past formulation of the subloading surface model will be modified to supply this insufficiency.
First, note the following facts: 1. The difference between the curvatures in the reloading and the inverse loading curves becomes larger as the plastic state develops. This fact is known as the Masing rule [42]. 2. The similarity-center to be the most elastic stress state approaches the normal-yield surface as the plastic state develops, and the approaching degree of the similarity-center to the normal-yield surface is expressed by the similarity-center yield ratio s in Eq. (14) as described in Sect. 3. 3. The transition from the elastic to the plastic state is more abrupt, i.e. the curvature of the stress-strain curve is larger for a larger value of the material parameter u in the function U (R) in Eq. (26), as described in Sect. 3. Therefore, u is expected to be larger for the reloading state but smaller for the inverse loading state. 4. Therefore, the difference between the values of u for the reloading and the inverse loading states should be larger for the larger value of s . 5. The degree of reloading can be judged by how the stress is directed outward from the similarity-center surface when it is observed from the similarity-center. It can be expressed by the scalar product of the direction vectorσ r / σ r and the normalized outward-normal vectorn s of the similarity-center surface in the six-dimensional stress space, whileσ r is the modified stress stemming from the similarity-center to the current stress as described below. The scalar product is presented as S σ can be regarded as the measure of the direction of deviatoric stress observed from the similarity-center surface.
Eventually, introducing the variables s and S σ , let the material parameter u in Eq. (26) be extended as follows: Fig. 2 Stress-plastic strain curve predicted by the modified subloading surface model: rapid recovery to preceding monotonic loading curve u is the continuous function of the variables s and S σ , while the forms of u for the particular values of them are designated in the bracket. S σ = 1, 0 and −1 designate the states that the similarity-center lies in the outward-normal, tangential, and inward-normal direction, respectively, of the similarity-center surface. s = 1 and 0 designate the states that the similarity-center lies on the center of the normal-yield surface and on the normal-yield surface, respectively.ū and u s are the material constants, whereas the former denotes the average value of u. Then, u increases in the loading direction leading to the reloading but inversely it decreases in its opposite direction leading to the inverse loading. By this modification, the phenomenon that the stress returns rapidly to the previous loading curve in the reloading process after a partial unloading is described realistically as depicted in Fig. 2, and thus, the plastic strain accumulation for the cyclic loading process is suppressed realistically.

Material functions of metals
The functions of the normal-yield, the subloading, and the similarity-center surface for the von Mises yield condition with the kinematic hardening are adopted for metals as follows: The hardening function H in Eq. (5) is given for metals as follows [14,17]): where h 1 , h 2 are the material constants, and ε E p is referred to as the equivalent plastic strain. The evolution equation of the kinematic hardening variable α is given as follows [17,22]: where a α and r α are material constants. In fact, Eq. (54) reduces to dα a = a α (r α F ∓ α a )dε p a in the uniaxial loading process, where the upper and lower signs signify extension and compression, respectively, ε p a denoting the axial plastic strain. Here, α translates to approach the conjugate point √ 3/2r α FN(= √ 3/2r α Fσ / σ ) on the limit surface σ = r α F of kinematic hardening, i.e., α → √ 3/2r α Fσ / σ (α a → r α F for the uniaxial loading).
Introducing Eq. (52) into Eqs. (14) and (50), one obtains Inserting Eqs. (51), (53), and (54) into Eqs. (20) and (31), the evolution rule of the similarity-center and the plastic modulus read The normal-yield ratio R is calculated by solving the equation of the subloading surface as follows: Substituting Eqs. (8.1) and (51.2) into Eq. (7), the subloading surface is expressed as The normal-yield ratio R is calculable using the following equation that is derived by solving the quadratic equation (58): where only the known variables σ, α, s and F are included in the right-hand side.

Cyclic stagnation of isotropic hardening in metals
It has been observed through experiments for metals that isotropic hardening stagnates despite the development of the equivalent plastic strain for a certain range in the beginning stage of re-yielding after reverse loading. This phenomenon considerably affects the cyclic loading behavior in which the reverse loading is repeated. In particular, the isotropic hardening saturates finally in the cyclic loading under a constant strain amplitude.
To describe this phenomenon, the concept of the non-hardening region, i.e., cyclic stagnation of isotropic hardening, was proposed by Ohno [44], modifying the memory surface of plastic strain proposed by Chaboche et al. [3] (see also [4]). The concept holds that isotropic hardening does not proceed when the plastic strain lies inside a certain region, called the non-hardening region, in the plastic strain space. The non-hardening region expands and translates when the plastic strain lies on the boundary of the region, and the plastic strain rate is induced in the outward direction of the region. It is similar to the notion of the yield surface based on the assumption that the plastic strain rate is not induced when the stress lies inside that surface, while the plastic strain and the isotropic hardening variable for the non-hardening region correspond to the stress and the plastic strain rate, respectively, for the yield surface. Thereafter, the formulation that the isotropic hardening stagnates when the back stress lies inside the certain region of stress space was proposed by Yoshida and Uemori [51,52], where the nonlinear kinematic hardening rule is adopted. However, it cannot describe the stagnation behavior of isotropic hardening of metals without kinematic hardening, while isotropic hardening and the kinematic hardening would be mutually independent phenomena. The following defects are involved in these formulations.
1. Isotropic hardening is induced suddenly when the plastic strain (or the back stress in Yoshida and Uemori's [51,52] formulation) reaches the boundary of the nonhardening region, violating the smoothness condition [15,16,18,21]. Consequently, the smooth stress-strain curve cannot be described except for a particular treatment. 2. A judgment of whether or not the plastic strain (or the back stress) reaches the boundary of the nonhardening region is necessary.

3.
A numerical operation to pull back the plastic strain (or the back stress) to the boundary of the nonhardening region so as not to move outward from the nonhardening region is necessary. 4. The direct time integration of plastic strain rate itself is physically meaningless in the general deformation process with a rigid-body rotation, but it has been used in the formulations of the cyclic stagnation of isotropic hardening [3,44].
In the following discussion, a pertinent formulation without these deficiencies is presented for the cyclic stagnation of isotropic hardening, based on the concept of a subloading surface.
First, define the basic variable for the formulation of the cyclic stagnation of isotropic hardening. The plastic strain rate D p is independent of the rigid-body rotation possessing the objectivity, and thus, it is related to the corotational stress rate • σ as the constitutive relation. Then, consider the state variable whose corotational rate is the plastic strain rate D p and let it be termed the accumulation of plastic strain rate, abbreviated as the plastic strain below for the sake of simplicity, letting it be denoted by ε p which may be regarded as a sort of measure of plastic deformation. By replacing T, Therefore, the plastic strain is calculated by the equation The continuum spin W can be adopted as the corotational spin ω describing the spin of rigid-body rotation of material for the infinitesimal strain but finite rotation as described in Sect. 2. Here, the rationality of Eq. (60) would be captured also by the example of Prager's [46] linear-kinematic hardening rule  (60) is to be one of the most fundamental and direct application of the corotational rate tensor, while it should be noted that this notion is limited to the infinitesimal deformation since the additive split of the strain rate into the elastic and the plastic parts does not hold in the finite deformation.
Equations (60) and (61) are indispensable in order to apply the cyclic stagnation theory of isotropic hardening to practical engineering problems that material rotates as seen often in machine tools such as gears, screws, wheel shafts, and pulleys. The past formulations for the cyclic stagnation of isotropic hardening have been limited to the deformation without rigid-body rotation.
Assuming that the isotropic hardening stagnates when the plastic strain ε p lies inside a certain region in the plastic strain space, let the following surface, called the normal-isotropic hardening surface, be introduced as stated by Chaboche et al. [3]: K and α designate the size and the center, respectively, of the normal-isotropic hardening surface, the evolution rules of which will be formulated later. The function f ( ε p ) in Eq. (62) is explicitly given by for Eq. (51) in the Mises yield condition.

Fig. 3 Normal-isotropic and sub-isotropic hardening surfaces
Furthermore, we introduce the surface, called the sub-isotropic hardening surface, which always passes through the current plastic strain ε p and which has the similar shape and orientation to the normal-isotropic hardening surface (see Fig. 3). It is expressed by the following equation: where R(0 ≤ R ≤ 1) is the ratio of the size of the sub-isotropic hardening surface to that of the normalisotropic hardening surface. It plays the role as the measure for the approaching degree of plastic strain to the normal-isotropic hardening surface. Then, R is referred to as the normal-isotropic hardening ratio. It is calculable from the equation R = f ( ε p )/K in terms of the known values ε p , α and K .
The material-time derivative of Eq. (65) leads to the consistency condition for the sub-isotropic hardening surface: Now, for the formulation of the evolution rule of the size K of the normal-isotropic hardening surface, let it be assumed that 1. The normal-isotropic hardening surface expands if and only if the sub-isotropic hardening surface expands, i.e. the plastic strain rate is directed outward of the sub-isotropic hardening surface. It therefore holds that 2. The expansion rate of the normal-isotropic hardening surface increases as the plastic strain approaches the normal-isotropic hardening surface, i.e. as the normal-isotropic hardening ratio increases. Then, • K is the monotonically increasing function of R. 3. In order that the normal-isotropic hardening surface expands continuously, its rate must be zero when the plastic strain lies just on its center, i.e.
• K = 0 for R = 0 because the rate is zero during the process that the plastic strain returns to the center, i.e.
Based on these assumptions, let the evolution rule for the size of the normal-isotropic hardening surface be given by where C and ζ(≥1) are the material constants. Now, for the formulation of the evolution rule of the center α of the normal-isotropic hardening surface, let the following be assumed.

The center of the normal-isotropic hardening surface translates if and only if the sub-isotropic hardening
surface expands, i.e. if the plastic strain rate is directed outward of the sub-isotropic hardening surface and thus it holds that Here, it is known from the consistency condition in Eq. (66) that N must be used for α The translation rate of the center of the normal-isotropic hardening surface increases as the plastic strain approaches the normal-isotropic hardening surface, i.e. as the normal-isotropic hardening ratio increases.
Then, α • is the monotonically increasing function of R. 3. In order that the normal-isotropic hardening surface translates continuously, the translational rate must be zero when the plastic strain lies just on its center, i.e. α • = 0 for R = 0 because the rate is zero during the process that the plastic strain returns to the center, i.e. α 4. The direction of translation is outward-normal N of the sub-isotropic hardening surface.
Based on these assumptions, let the following evolution rule for the center of the normal-isotropic hardening surface be given by where A is the material constant which will be related to the material constant C below. Then, substituting Eqs. (67) and (69) into Eq. (66), one obtains Taking account of R • = 0 for R = 1 in Eq. (70), the following must hold: Substituting Eq. (71) into Eq. (69), the evolution rule for the center of the normal-isotropic hardening surface is determined finally as The normal-isotropic hardening surface expands without the translation in case of C = 1 but inversely translates without the expansion in case of C = 0. Here, assume that the normal-isotropic hardening surface evolves in between these extreme cases leading to Then, R • is given by in which it holds that noting Eq. (73) (see Fig. 4). The evolution equation (24) of the normal-yield ratio R is first assumed and then it is incorporated into the consistency condition of the subloading surface to formulate the plastic modulus. On the other hand, the evolution equation (74) of the normal-isotropic hardening ratio R is not formulated first but is instead obtained from the consistency condition of the sub-isotropic hardening surface with the substitutions of evolution rules for the size and the translation of the normal-isotropic hardening surface. The calculation is controlled to fulfill the consistency condition (66). Therefore, the plastic strain is attracted automatically to that surface even if it goes out from that surface by virtue of the inclusion of the inequality R • < 0 for R > 1 as shown in Eq. (75). Furthermore, the judgment of whether or not the plastic strain lies on the normal-isotropic hardening surface is not necessary in the present model. In contrast, the incorporations of the computer algorithms for the judgment whether or not the plastic strain lies on the boundary of the complete stagnation region of isotropic hardening and for pulling back the plastic strain to the boundary are necessary in the other models for the cyclic stagnation of isotropic hardening [3,4,44,51].
It is assumed that the isotropic hardening variable H evolves under the following rules.
1. Isotropic hardening is induced if and only if the sub-isotropic hardening surface expands, i.e. the plastic strain rate is directed outward of the sub-isotropic hardening surface. Consequently, it holds that 2. The isotropic hardening rate increases as the plastic strain approaches the normal-isotropic hardening surface, i.e. as the normal-isotropic hardening ratio increases. Then, • H is the monotonically increasing function of R. 3. In order for the isotropic hardening to develop continuously, its rate must be zero, i.e. Then, let the following evolution rule of isotropic hardening be assumed by extending Eq. (53.2): where υ(≥1) is the material constant. Employing the extended isotropic hardening rule in Eq. (76), the rate of similarity-center in Eq. (56) and the plastic modulus in Eq. (57) are modified as follows: On the other hand, the isotropic hardening is restored suddenly when the plastic strain reaches the boundary of the non-hardening region so that the nonsmooth stress-strain curve is predicted in the past formulations for the cyclic stagnation of isotropic hardening [3,44,51,52]. In order to avoid this deficiency, it was proposed by Ohno [44] that the blunting of isotropic hardening is supplemented by the acceleration of kinematic hardening in the state that the plastic strain lies inside the isotropic hardening surface. However, without resorting to such a method, this defect is avoided, always predicting the smooth stress-strain curve in the present formulation based on the notion of the sub-isotropic hardening surface.

Summary of constitutive equations
Summarizing the equations formulated in the foregoing, the set of constitutive equations required for the calculation of the stress-strain relation are listed below.
Subloading surface: Relations of stress, back stress, and similarity-center: Stress loading function: Isotropic hardening function: Stress rate versus strain rate relation: Normalized outward-normal of subloading surface: Loading criterion: Equation for the calculation of normal-yield ratio: Evolution rule of kinematic hardening: Evolution rule of normal-yield ratio Material function in evolution rule of normal-yield ratio: Normal-yield similarity-center ratio: Measure of direction of deviatoric stress observed from similarity-center surface: Normalized outward-normal of similarity-center surface: Sub-isotropic hardening surface: Stagnation function of plastic strain: Evolution rule of size of normal-isotropic hardening surface: Evolution rule of center of normal-isotropic hardening surface: Normalized outward-normal of normal-isotropic hardening surface: Evolution of isotropic hardening: Evolution of similarity-center: Plastic modulus:

Comparisons with test results
The capability of the present model for describing the deformation behavior of metals is verified through comparisons with several basic test data in this Section. The capability of the unconventional plasticity model aimed at describing the plastic strain rate induced by a rate of stress inside the yield surface must be evaluated by a degree to which cyclic loading behavior can be described appropriately. Then, various cyclic loading test data in uniaxial loading are first simulated, and thereafter, a circular strain path test datum is simulated to verify the capability for describing non-proportional loading behavior. The material parameters involved in the present model are shown collectively below, where 16 material constants and five initial values are included to describe the cyclic loading behavior precisely.
The determination of these material parameters is explained below in brief.
1. Young's modulus E and Poisson's ratio ν are determined from the slope and the ratio of lateral to axial strains in the initial part of the stress-strain curve. 2. h 1 , h 2 and F 0 for the isotropic hardening and a α , r α and α 0 for the kinematic hardening are determined from stress-strain curves in the initial and the inverse loadings. 3. R e ,ū and u s are determined from the stress-strain curve in the subyield state, i.e. the elastic-plastic transitional state. 4. c, χ and s 0 are determined from the stress-strain curves in cyclic loading. 5. ξ and τ are determined by the difference of the strain in the non-proportional loading from that in the proportional loading. 6. C, ζ, υ, K 0 and α 0 are determined from the stress-strain curves in cyclic loading with a constant strain amplitude.
All of these material parameters except for ξ and τ for the tangential-inelastic strain rate can be determined only by the stress-strain curves in the uniaxial loading for initial isotropic materials. One can put α 0 = s 0 = α 0 = 0 for the initial isotropy, which is assumed in all the following simulations. Tangential inelasticity is not considered so that ξ and τ for the tangential inelasticity are ignored for the test data in the uniaxial loading. The cyclic loading behavior under the stress amplitude to both positive and negative sides can be predicted to some extent by any models, including even the conventional plasticity model. On the other hand, the prediction of the cyclic loading behavior under the stress amplitude in positive or negative one side, i.e., the pulsating loading inducing the so-called mechanical ratcheting effect requires a high ability for the description of plastic strain rate induced by the rate of stress inside the yield surface. Furthermore, it is noteworthy that we often encounter the pulsating loading phenomena in the boundary value problems in engineering practice, e.g., railways and gears. The comparison with the test data for the 1070 steel under the cyclic loading of axial stress between 0 and +830 MPa after [37] is depicted in Fig. 5, where the material parameters are selected as follows: Material constants: The relation of the axial stress and the axial components of back stress and similarity-center versus the axial strain and the relation of the axial strain versus the number of cycles are depicted in Fig. 5a, where the axial components are designated by () a . The accumulation of axial strain is simulated closely by the present model. The calculation is controlled automatically such that the stress and the plastic strain are attracted to the normal-yield and the normal-isotropic hardening surfaces, respectively, as known from the variations of the normal-yield ratio R and the normal-isotropic hardening ratio R depicted in Fig. 5b. Accumulation of axial strain is overestimated as depicted in Fig. 5c if the reloading behavior is not improved setting u s = 0. Despite the improvement for reloading behavior, however, hysteresis loops are simulated as narrower than those in the test result in order to fit the strain accumulation in the test result. A further improvement is desirable for this insufficiency.
Next, examine the uniaxial cyclic loading behavior under the constant stress amplitude to both positive and negative sides with different magnitudes. Comparison with the test data for the 304L steel under the cyclic loading of axial stress between +250 and −150 MPa after [33] is depicted in Fig. 6 where the material parameters are selected as shown below. Test result ( Xia and Ellyin, 1994 ) Model simulation

Fig. 7
Uniaxial cyclic loading behavior under the constant stress amplitude between −182 and +182 MPa of 304L steel (test data after [50]): a test result and simulation, b variations of normal-yield ratio and normal-isotropic hardening ratio and c simulation by Chaboche model and the normal-isotropic hardening ratio R depicted in Fig. 6b. The relation of the axial stress versus the axial strain and the relation of the axial strain versus the number of cycles simulated using the modified Chaboche model [4] are also depicted in Fig. 6c in which the strain is simulated as larger than the test result, and the hysteresis loops are simulated as narrower than the test data. The prediction of this steel deformation behavior will be improved by incorporating the rate-dependence. Further, we examine the uniaxial cyclic loading behavior for constant symmetric stress amplitude to both positive and negative sides. Comparison with test data of the 304 steel under the cyclic loading of axial stress between +182 and −182 MPa under the constant hoop stress 80 MPa after [50] is depicted in Fig. 7 where the material parameters are selected as follows: Material constants: The relation of the axial stress and the axial components of back stress and similarity-center versus the axial strain and the circumferential strain ε l with the number of cycles are shown in Fig. 7a, while the back stress is induced quite slightly so that it is invisible in this figure. The simulations for the accumulation of axial strain and the hysteresis loops agree well with the test result, except for the prediction of hysteresis loops as narrower than the test result in the initial stage. Here, the axial strain and the lateral strain are accumulated to the compression side and the extension side, respectively, by the application of the hoop stress 80 MPa. The calculation is automatically controlled such that the stress and the plastic strain are attracted to the normal-yield and the normal-isotropic hardening surfaces, respectively, as known from the variations of the normal-yield ratio R and the normal-isotropic hardening ratio R depicted in Fig. 7b. The relations of the axial stress versus the axial and lateral strains and the relation with the number of cycles simulated by Xia and Ellyin [50] (cf. also [11]) are also depicted in Fig. 7c where both the axial and the circumferential strains are overestimated.
Furthermore, examine the uniaxial cyclic loading behavior under the constant symmetric strain amplitudes to both positive and negative sides. Comparison with the test data of the 316 steel under the cyclic loading with the increasing axial strain amplitudes ±1.0, ± 1.5, ± 2.0, ± 2.5, ± 3.0% after [3] is depicted in Fig. 8 where the material parameters are selected as follows: Material The relation of the axial stress and the axial components of back stress and similarity-center versus the axial strain is shown in Fig. 8a. The hysteresis loops and the stagnation of isotropic hardening are simulated closely by the present model. On the other hand, the calculated result without the cyclic stagnation of isotropic hardening overestimates the hardening behavior as shown in Fig. 8b. The calculation is controlled automatically such that the stress and the plastic strain are attracted to the normal-yield and the normal-isotropic hardening surfaces, respectively, as known from the variations of the normal-yield ratio R and the normal-isotropic hardening ratio R depicted in Fig. 8c. The relations of the axial stress and the axial strain simulated by Chaboche [4] and Ellyin and Xia [10] are depicted in Fig. 8d and e, respectively. The strain in the initial stage is simulated as larger than the test result by the former, and the curves predicted by the latter are not smooth but piecewise linear. Finally, we examine the non-proportional loading behavior. Comparison with the test data of the austenitic 17-12 Mo SPH carbon stainless steel for the circular strain path ε a = 0.004 cos α and γ aθ = 0.0036 sin α under σ θ = 50 MPa during 40 cycles after the uniaxial loading to ε a = 0.004 after [8] is depicted in Fig. 9, where σ θ is the circumferential normal strain, and γ aθ is the axial-circumferential engineering shear strain, and α is the angle measured from the axis of ε a in the strain plane (ε a , γ aθ ). The material parameters are selected as follows: Material constants:   The strain path (ε a , ε θ ) (ε θ : circumferential normal strain) and the stress path (σ a , √ 3σ aθ ) (σ aθ : axial-circumferential shear stress) are shown for the test result and the model simulation in Fig. 9a and b, respectively. The simulation of the stress path and the accumulation of lateral strain are in good agreement with the test result. However, the initial stress loop is simulated as a spiral. This point will have to be improved hereinafter. The stress and the plastic strain are attracted to the normal-yield and the normal-isotropic hardening surfaces, respectively, as known from the variations of the normal-yield ratio R and the normal-isotropic hardening ratio R depicted in Fig. 9b. On the other hand, the accuracy of simulation decreases if the tangential-inelastic strain rate is not incorporated by setting ξ = 0 as depicted in Fig. 10. The axial component (D t ) a and the axial-circumferential component (D t ) aθ of the tangential-inelastic strain rate are induced considerably, but the circumferential component (D t ) θ is not induced to a marked degree because of ( • σ t ) θ ∼ = 0 in this test. Therefore, D a /D θ is calculated as smaller and thus ε θ /ε a is calculated as larger if we put ξ = 0.

Concluding remarks
The constitutive equation of metals is formulated refining the past formulation of the subloading surface model by avoiding the indetermination of the subloading surface, by inducing the reloading curve to recover promptly to the preceding loading curve and by introducing the cyclic stagnation of isotropic hardening [3,44] with the important modification by the concept of subloading surface. The applicability of the present model to the description of real deformation behavior of metals is verified through the comparisons with various test data for cyclic loadings. It is attained by virtue of the following distinguished capabilities of the present model.
i. The smoothness condition [15,16,18,21] is fulfilled. Therefore, a continuous variation of tangent modulus for the continuous change of stress state leading to smooth stress-strain curve is predicted always by the present model. ii. The plastic strain rate develops continuously as the stress approaches the normal-yield surface in the present model. iii. The stress always lies on the subloading surface which plays the role of the loading surface. Therefore, only the judgment for the sign of the proportionality factor is necessary for the loading criterion in the present model. iv. The calculation is controlled to fulfill the consistency condition. Therefore, the stress is attracted automatically to the normal-yield surface in the plastic loading process even if it goes out from that surface by virtue of the inclusion of the inequality • R < 0 for R > 1 in Eq. (24) with Eq. (25). Then, it is not necessary to incorporate a computer algorithm for pulling back the stress to the normal-yield surface. Therefore, a stable and robust calculation can be executed even by rather large loading steps in the present model. v. The formulation of the subloading surface model is refined such that the plastic modulus increases in the reloading process, but it decreases in the inverse loading process. Therefore, the cyclic loading behavior is described pertinently suppressing the strain accumulation. vi. The tangential-inelastic strain rate induced by the stress rate tangential to the subloading surface [19,27] is incorporated fulfilling the continuity and the smoothness conditions, which cannot be ignored in the non-proportional loading process. vii. The pertinent formulation of the cyclic stagnation of isotropic hardening is obtained by modifying the formulations by Ohno [44] and Chaboche et al. [3]. First, the rigorous calculation method of the basic variable, called the accumulation of plastic strain and abbreviated as plastic strain, for the cyclic isotropic hardening is given in Eq. (61). Therefore, the formulation for the cyclic stagnation of isotropic hardening becomes applicable to the deformation process subjected to the finite rigid-body rotation. Furthermore, the cyclic stagnation of isotropic hardening is formulated pertinently introducing the concept of a subloading surface. viii. The plastic strain always lies on the sub-isotropic hardening surface. Therefore, only a judgment based on the sign of the quantity tr( ND p ) is necessary in relation to the generation of isotropic hardening in the present model. Then, the judgment of whether or not the plastic strain reaches the normal-isotropic hardening surface is possible. Furthermore, the calculation is controlled automatically to fulfill the consistency condition. Therefore, the plastic strain is attracted to that surface even if it goes out from that surface by virtue of the formulation R • < 0 for R > 1 in Eq. (75). Therefore, a stable and robust calculation can be executed even by rather large loading steps in the present model.
Eventually, it may be concluded that the elastoplastic constitutive equation of metals formulated in this article provides a pertinence in both aspects of the pertinence in the description of plastic deformation and high efficiency and robustness in numerical calculations.
Each particular experimental situation is represented by a different special set of material constants in this paper. Simulations of test data under various loading conditions by means of one single set of material constants would be more difficult. They will be shown in a future paper by finding a pertinent set of test data.