Subloading-Overstress Model: Unified Constitutive Equation for Elasto-Plastic and Elasto-Viscoplastic Deformations Under Monotonic and Cyclic Loadings

Various elasto-plastic models for the rate-independent deformation, various elasto-viscoplastic models for the rate-dependent deformation and their combinations have been proposed during a long time more than one or more centuries. Firstly, the history of the development of the elastoplasticiy and the elasto-viscoplasticity is reviewed comprehensively. Unfortunately, each of these models possesses their own drawbacks and limitations. The unified constitutive equation of the elasto-plastic and the elasto-viscoplastic deformations is provided by incorporating the subloading surface model into the overstress model in this article, which is capable of describing the monotonic and the cyclic loadings at the general rate ranging from the quasi-static to the impact loading. The validity of the unified model is verified by the comparison with various test data of metals under various loading conditions. Consequently, the elastoplastic constitutive equation can be disused hereinafter, while it is expressed by the cumbersome formulation including the complicated plastic modulus but limited to the description of the purely static deformation which is not induced actually.


Introduction
The researches on the elastoplasticity for the quasi-static deformation behavior and the viscoplasticity for the dynamic deformation behavior have the long history in one or more centuries. In this context, the main goal of these researches would be the formulation to describe them in a unified equation for the monotonic and the cyclic loadings at the general rate of deformation, while the solids and structures in engineering practice are subjected to monotonic and the cyclic loadings at variable rate of deformation from the quasi-static to the impact loading.
Firstly, the history of the development of the elastoplastic and the elasto-viscoplastic constitutive equations will be reviewed comprehensively in the following.
The elasto-plastic and the elasto-viscoplastic deformations have been often described separately by the elastoplastic model for the quasi-static deformation and the creep model for the dynamic deformation, respectively, Norton [97], Odqvist [98], Robotnov [120], Rice [118,119], Bodner and Partom [10], Lemaitre and Chaboche [86], Abdel-Karim and Ohno [1], Saleeb et al. [122], Betten [7], Asaro and Lubarda [5], Chaboche [17], Ohno et al. [101], de Sauza Neto et al. [25], Peirce et al. [111,112], etc. However, it should be noticed that the creep strain rate in the creep model is induced even in any low stress level since it is formulated by the stress or the ratio of the stress to the yield stress and thus the response of the creep model is not reduced to that of the elastoplastic model in a quasi-static rate of deformation. Therefore, the deformation behavior in the intermediate rate cannot be described pertinently by them, because different responses are described by the creep 1 3 and the elastoplastic constitutive equations. Consequently, the combination of the creep model and the elastoplastic model cannot be used for deformation analysis at variable rate. Besides, various empirical models, Johnson and Cook [77], Steinberg et al. [125], Zerilli and Armstrong [144], Follansbee and Kocks [34], etc. have been proposed and used to the engineering in practice but they are ad. hoc models for particular phenomena without generality. Further, the fundamentally-irrational models, Bailey [6], Garofalo [35], etc. involving the time itself in order to describe the creep behavior have been proposed, which result in the loss of the objectivity as the constitutive relation Oldroyd [104].
On the other hand, the response of the overstress model, Bingham [9], Hohenemser and Prager [70], Prager [116], Perzyna [114,115], Duvaut and Lion [28], Krempl et al. [82], Chaboche [12,14], Lubliner [93], Perić [113], Arnold and Saleeb [4], Simo [123], Simo and Hughes [124], Kaneko and Oyamada [78], Ho and Krempl [69], Lubarda [92], Ottosen and Ristinmaa [105], Chaboche [16,17], Saleeb and Arnold [121], Mayama et al. [95], de Sauza Neto et al. [25], Ho [68], Chaboche et al. [19], Guo et al. [36], Chen and Feng [20], Chen, et al. [21,22], etc. is reduced to that of the elastoplastic constitutive equation in the quasi-static rate of deformation. Therefore, the overstress model possesses the capability for the description of the deformation at the variable rate. However, the existing overstress model is incapable of describing the cyclic loading behavior under the general stress amplitude, since it adopts the conventional elastoplastic constitutive equation with the yield surface enclosing the purely-elastic domain, and it is also incapable of describing the impact loading behavior since it expresses the elastic deformation behavior leading to the infinitely high stress in the infinitely high rate of deformation. The subloading surface model, Hashiguchi [37,38,40,48,50] is capable of describing the elastoplastic deformation under the monotonic and the cyclic loading behaviors, while it does not incorporate a purely-elastic domain so that the tangent stiffness modulus changes continuously always satisfying the continuity condition in the large, i.e. the smoothness condition, Hashiguchi [41] and the mechanical ratchetting phenomenon is described for any low stress level. The subloading surface model has been tried to describe the ratedependent behavior by incorporating the overstress model by Hashiguchi [46,48,50], Hashiguchi et al. [53], Moranha et al. [94], Fincato and Tsutsumi [31,32] and Anjiki and Hashiguchi [2]. However, the peculiar stress vs. strain curve exhibiting the softening is described by the formulations by Hashiguchi [46,48,50], Hashiguchi et al. [53] and Moranha et al. [94]. The formulation by Fincato and Tsutsumi [31,32] exhibits the discontinuous tangent stiffness modulus such that the inflexion point is induced in the stress-strain curve. The formulation by Anjiki and Hashiguchi [2] is quite complex such that the similarity-center of the normal-yield and the subloading surfaces moves departing from the elasticcore and thus the quite cumbersome numerical calculation is required.
In this article, the unified constitutive equation for the elasto-plastic and the elasto-viscoplastic deformations is provided, which is capable of describing the deformation behaviors during the monotonic and the cyclic loadings at the general rate ranging from the quasi-static to the impact loading by incorporating the subloading surface model into the overstress model. Further, the isotropic stagnation phenomenon is incorporated into the model. Besides, the physical interpretations of all the formulations are given comprehensively. Then, the explicit functions included in the constitutive equation are given for metals. The validity of the present model is verified by the simulations of test data for metals under the various loadings at several strain rates including variable strain rate and the stress relaxations at fixed strains, Inoue et al. [72], Abel-Karim and Ohno [1], Saleeb and Arnold [121], Lemaitre and Chaboche [86]. Further, the test data reported as the rate-independent elastoplastic deformations at room temperature, Jiang and Zhang [74] and Chaboche et al. [18] are simulated accurately by the present model.

Subloading Surface Model for Elastoplastic Deformation
The complete formulation of the subloading surface model for the elastoplastic deformation modifying the past one, Hashiguchi [37,38,40,48,50] is described concisely in this section, while the subloading surface model for metals with the simplified evolution rule of the elastic-core limited to the nonhardening state and without the isotropic hardening stagnation is installed in the commercial FEM software Marc. By combining the subloading surface model and the overstress model, the subloading-overstress model will be formulated for the description of the monotonic and the cyclic loading behaviors in the general rate from the quasistatic to the impact loading in the next section. The distinguished feature of the subloading surface model, which is not furnished in the other elastoplastic models, e.g. the superposed kinematic hardening model, Chaboche et al. [18], Ohno and Wang [102], Hassan et al. [66], etc., the multi surface model, Mróz [96], Iwan [73] and the two surface model, Dafalias and Popov [23], Krieg [83], Hashiguchi [39], Yoshida and Uemori [140], etc., is the exclusion of the purely elastic domain so that the plastic strain rate due to the variation of stress inside the yield surface can be described. Therefore, the continuity condition in the large, i.e. the smoothness condition, Hashiguchi [41,42,44] leading to the continuous variation of the tangent stiffness modulus tensor is fulfilled always. It is of crucial where c k is the material constant with the stress dimension and b k (< 1) is the non-dimensional material constant. Equation (9) is the slight modification of the nonlinear kinematic hardening rule of Armstrong-Frederick [3] in which b k F is defined to be the material constant with the stress dimension, while it should be noticed that the translational range of the yield surface would not be independent of the size of the yield surface. Now, let the subloading surface be incorporated, which is similar to the yield surface, renamed as the normal-yield surface, with respect to the similarity-center and thus it is described as follows (see Fig. 1).
where R(0 ≤R ≤ 1) is the normal-yield ratio, i.e. the ratio of the size of the subloading surface to that of the normal-yield surface and stands for the conjugate (similar) point in the subloading surface to the point representing the back-stress, i.e. the kinematic hardening variable in the normal-yield surface. Here, note that the purely-elastic behavior is induced when the stress coincides with the similarity-center and thus the normal-yield ratio becomes zero (R = 0) . Then, let the similarity-center be called the elastic-core. Here, the following relation holds by virtue of the similarity (see Fig. 1). which yields Besides, y in Fig. 1 is the conjugate stress on the normalyield surface to the current stress on the normal-yield surface, fulfilling − = R( y − ) . The time-derivative of is described from Eq. (13) as The subloading surface in Eq. (11) is rewritten by Eq. (15) as follows: R is calculated by solving Eq. (18) in the elastic loading process ̇ e ≠ ,̇ p = by substituting the updated values of , , and F . Equation (18) is expressed by the quadratic equation and thus the normal-yield ratio can be expressed analytically for the Mises material as will be shown in Sect. 2.7.
The evolution rule of the normal-yield ratio R is given by where the function U(R) of R satisfies the conditions The material parameter R e is interpreted to be the ratio of the fatigue limit stress to the yield stress y . Note here that the incorporation of the material parameter R e does not mean the incorporation of the surface enclosing a purely-elastic domain as known from the fact: The plastic strain rate is predicted for the cyclic loading with a small stress amplitude under a high average stress by the subloading surface model with the incorporation of R e . Incidentally, the tangent stiffness modulus changes continuously satisfying the continuity condition in the large, i.e. the smoothness condition always in the monotonic loading process. However, the normalyield ratio R must be calculated from Eq. (18) for the elastic loading process (̇ p = 0,̇ e ≠ 0) . It can be expressed explicitly by Eq. (84) for metals as will be shown in Sect. 2.7.
for R e < R < 1 (subyield state) = 0 for R = 1 (normal-yield state) < 0 for R > 1 (over normal-yield state) The explicit equation satisfying Eq. (20) is given by the cotangent function where u is the material parameter and ⟨ ⟩ is the Macaulay's bracket defined by ⟨s⟩ = (s+�s�)∕2 , i.e. s < 0 ∶ ⟨s⟩ = 0 and s ≥ 0 ∶ ⟨s⟩ = s ( s : arbitrary scalar variable). Equation (21) fulfills the continuity condition in the large, i.e. the smoothness condition (Hashiguchi [42]) since the function U(R) decreases continuously from the infinite value U(R e )( → ∞) . If u is fixed to be constant, Eqs. (19) with (21) can be integrated analytically as where p ≡ ∫ ||̇ p ||dt and p 0 is the initial value of p , whilst one must set R 0 = R e for R 0 < R e .

Evolution Rule of Elastic-Core
The elastic-core is not allowed to approach the normal-yield surface closely, because the purely elastic deformation is induced when the stress lies on the elastic-core but the apparent plastic deformation is induced when the stress lies on the normal-yield surface. In addition, it should be noted that a smooth elastic-plastic transition leading to the continuous variation of the tangent stiffness modulus tensor is not described violating the continuity condition in the large, i.e. the smoothness condition Hashiguchi [41], if the elastic-core lies on the normal-yield surface. On the other hand, the other plasticity models, i.e. the cylindrical superposed kinematic hardening model, Chaboche et al. [18], Ohno and Wang [102], Hassan et al. [66], etc., the multi surface model, Mróz [96], Iwan [73] and the two surface model, Dafalias and Popov [23], Krieg [83], Hashiguchi [39], Yoshida and Uemori [140], etc. violate the smoothness condition, since the purely-elastic domain contacts to the yield surface directly.
Then, it is assumed that the elastic-core surface which passes through the elastic-core is similar to the normal-yield surface with respect to the back-stress is smaller than the normal-yield surface. Therefore, the following inequality holds.
where ℜ c designates the ratio of the size of the elastic-core surface to the normal-yield surface (see Fig. 1) and (< 1) is the material constant designating the maximum value of ℜ c . The material-time derivative of Eq. (23) at the limit state that lies on the limit elastic-core surface f (̂ ) = F(H) yields: Here, noting on account of the Euler's homogeneous function f (̂ ) in degree-one for the variable ̂ , the substitution of Eq. (25) to Eq. (24) leads to Now, assume the equation where c e is the material constant and is the conjugate stress on the limit elastic-core surfaces to the current stress on the subloading surface, i.e.
which is based on the relation where y is the conjugate stress on the normal-yield surface. Equation (27) means that the elastic-core translates so as to approach the conjugate stress in the nonhardening state: It follows for Eqs. (26) with (27) that as far as the normal-yield surface satisfies the convexity condition (cf. Hashiguchi [48]), noting that f (̂ )∕ which is the outward-normal of the elastic-core surface at the current elastic-core and − makes an obtuse angle with f (̂ )∕ when lies on the limit elastic-core surface, while lies on the limit elastic-core surface. ̂ c is the normalized outwardnormal of the elastic-core surface (see Fig. 1), i.e. Therefore, Eq. (27) satisfies the enclosing condition of the elastic-core so that the elastic-core can never go out from the limit elastic-core surface.
Then, the evolution rule of the elastic-core is given from Eq. (27) with Eqs. (8) and (9) by The translation of the elastic-core would have to be influenced by the variation of the expansion/contraction (isotropic hardening) and the translation (kinematic hardening) of the normal-yield surface in general, since its movement is limited to the interior of the normal-yield surface which expands and translates with the plastic deformation. This aspect is different from the evolution rule of the kinematic hardening which causes the translation of the normal-yield surface.

Plastic Strain Rate
Adopt the associated flow rule for the subloading surface: where The rates of the isotropic hardening in Eq. (8) and the rate of kinematic hardening in Eq. (9) for the state in which the current stress lies on the normal-yield surface are extended to for the state in which the current stress lies on the subloading surface.
Equation (32) with Eqs. (35) and (36) leads to The material-time derivative of Eq. (11) leads to the consistency condition of the subloading surface as follows: Here, one has The substitutions of Eq. (33) with Eqs. (19), (35), (36) and (37) into Eq. (45) leads to from which the plastic multiplier • and the plastic strain rate ̇ p are given as follows: The stress rate is given from Eq. (5) with Eq. (50) as follows: The loading criterion is given as follows (Hashiguchi [43,45,47]): where the judgment whether or not the stress reaches the yield surface is not required since the plastic strain rate is induced continuously as the stress approaches the normalyield surface. Equation (52) is applicable not only to the hardening state but also to the perfectly-plastic and softening state.

Improvement of Inverse and Reloading Behavior
The is the function of only the normal-yield ratio R . Therefore, p induced during a certain change of R in the monotonic loading process is identical irrespective of the difference of loading processes, e.g. the initial loading, the reloading and the inverse loading and of the proportional and non-proportional loadings. This property causes 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 in the extended subloading surface model.
As the fact observed in test data, the plastic strain rate is suppressed in the loading direction in which the elastic-core exists observing from the center of the normal-yield surface, i.e. the back stress and inversely the large plastic strain rate is induced in its opposite direction leading to the Masing effect. This tendency is more remarkable when the elastic-core is closer to the normal-yield surface leading to the larger value of ℜ c in Eq. (23). The degree of the loading/unloading is expressed by the scalar variable ̂ c ∶ where ̂ c is the normalized outward-normal of the elastic-core surface as defined in Eq. (31) and is the normalized outward-normal of the subloading surface, i.e. the direction of the plastic strain rate in Eq. (34). The magnitude of the plastic strain rate is controlled by the material parameter u in the function U(R) in Eq. (21) for the evolution rule of the normal-yield ratio R in Eq. (19).
Then, let the material parameter u be extended to the following equation for the extended subloading surface model in which the elastic-core, i.e. the similarity-center moves with the plastic deformation, introducing the variables ℜ c and C n .

Cyclic Stagnation of Isotropic Hardening
As a plastic deformation induced by the mutual slips of solid particles proceeds, the isotropic hardening proceeds with the loss of the dislocation mobility caused by the piling-up of obstacles or the formation of various networks. However, when a load reversal occurs, the isotropic hardening ceases temporarily by the recovery of the dislocation mobility. This phenomenon considerably affects the cyclic loading behavior in which the reverse loading is repeated. To describe this phenomenon, the concept of the cyclic stagnation of isotropic hardening, i.e. the existence of nonhardening region was proposed by Chaboche et al. [18] (see also Chaboche [13,14]). The concept insists that the isotropic hardening does not proceed when the plastic strain given by the time-integration of plastic strain rate lies inside a certain region, called the nonhardening region, in the plastic strain space.
Assuming that the isotropic hardening stagnates when the plastic strain p lies inside a certain region, let the following surface, called the normal-isotropic hardening surface, be introduced.
where g(̃ p ) is chosen to be the homogeneous degree-one of the variable ̃ p . The variables K and designate the size and the center, respectively, of the normal-isotropic hardening surface, the evolution rules of which will be formulated later. Furthermore, we introduce the surface, called the sub-isotropic hardening surface, which always passes through the current plastic strain p and which has a similar shape and a same orientation to the normal-isotropic hardening surface (see Fig. 2). Here, the plastic strain p is Fig. 2 Normal-and sub-isotropic hardening surfaces regarded to be the internal variable as • = c k̇ p is used in the Prager's linear kinematic hardening rule. On the other hand, the kinematic hardening variable, i.e. back-stress is used instead of p by Yoshida and Uemori [140]. However, would not be relevant but p would be relevant as the isotropic stagnation variable, since the mechanism of the isotropic hardening would be different from that of the kinematic hardening.
The subloading-isotropic stagnation surface is expressed by the following equation.
where R (0 ≤R ≤ 1) is the ratio of the size of sub-isotropic hardening surface to that of the normal-isotropic hardening surface. It plays the role as the measure for the approaching degree of the 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 = g(̃ p )∕K in terms of the known values of p , and K .
The consistency condition of the subloading isotropic hardening surface is given by Equation (60) is reduced to the following relation which must be fulfilled when the plastic strain just lies on the normalisotropic hardening surface. i.e.
Then, we assume the following equations so as to fulfill Eq. (61).
where 0 ≤ C ≤ 1 and ( ≥ 1) are the material constants and Substituting Eqs. (62) and (63) for the evolution rules of K and into Eq. (60), the rate of the normal-isotropic hardening ratio is given by which is the monotonically-decreasing function of R fulfilling Therefore, the normal-isotropic hardening ratio R increases when the plastic strain moves to the outward of the sub-isotropic hardening surface but it decreases such that the normal-isotropic hardening surface encloses the plastic strain when the plastic strain tends to go out from the normal-isotropic hardening surface by virtue of the inequality • R < 0 for R > 1 as shown in Eq. (66). Furthermore, needless to say, the judgment of whether the plastic strain reaches the normal-isotropic hardening surface is not necessary in the present formulation consistently based on the subloading concept in all aspects. In contrast, the judgment whether the plastic strain or the back-stress reaches the isotropic hardening (stagnation) surface is required and the input loading increments must be infinitesimal such that the isotropic stagnation surface encloses the plastic strain in the other models, Chaboche et al. [18], Chaboche [13][14][15], Ohno [99], Ohno and Kachi, [100], Ohno et. al. [101], Yoshida and Uemori [140,141].
Eventually, let the evolution rule of isotropic hardening be assumed as follows: where ( ≥ 1) is the material constant and Employing the extended isotropic hardening rule in Eq. (67) instead of Eq. (35) into Eq. (48), the plastic modulus is modified as follows: The function g(̃ p ) is given in the simplest form as follows: In the subloading surface model, the stress is pulled-back to the normal-yield surface and the isotropic hardening surface translates so as to involve the plastic strain automatically, while these mechanical functions are not furnished in the other plasticity and the viscoplasticity models. Therefore, the rather large incremental steps can be input even in the Euler forward method in the numerical calculation by the subloading surface model.

Mises Metals
The material functions for metals are given in this section, while the subloading surface model has been applied to various materials other than metals, i.e. soils, Hashiguchi The isotropic hardening function is given by where s r and c H are the material constants. The rates of the kinematic hardening variable and the elastic-core are given by Eqs. (36) and (37) as The plastic modulus is given by substituting Eqs. (73), (74), (78), (79) and (80)  i.e.
The normal-yield ratio R is derived from the quadratic Eq. (83) as follows:

Subloading-Overstress Model: Extension to Description of Elasto-Viscoplastic Deformation in General Rate
The subloading surface model will be extended to the subloading-overstress model which is capable of describing the monotonic and the cyclic deformations in general rate from quasi-static (elastoplastic) to the impact loading. The strain is additively decomposed into the elastic strain e and the viscoplastic strain vp instead of Eq. (2) as follows:

Creep Model and Overstress Model for Viscoplastic Deformation
The viscoplastic models are classified into the creep model (Norton [97] and the overstress model whose rheology models are shown in Fig. 3. The creep model is regarded as the nonlinearization of the Maxwell model for the viscoelastic deformation so that the creep strain rate is induced in any low stress level. On the other hand, the overstress model, Bingham [9] is the incorporation of the slider exhibiting the yield condition into the dashpot in parallel extending the creep model such that there exists the threshold value for the generation of the viscoplastic deformation. Therefore, the viscoplastic deformation is induced by the overstress from the yield stress but it is not induced for the stress lower than the yield stress. Therefore, the response of the creep model is not reduced to that of the • =̇ e +̇ vp elastoplastic constitutive equation but that of the overstress model is reduced to the elastoplastic constitutive equation for the quasi-static deformation process. Therefore, the creep model is not acceptable, although it is used widely (e.g. Lemaitre and Chaboche [86], Abdel-Karim and Ohno [1], Betten [7], Chaboche [17], Ohno et al. [103], de Sauza Neto et al. [25]) without regard to this serious defect but the overstress model is acceptable for the description of the viscoplastic deformation. Consequently, the overstress model should be adopted for the description of the viscoplastic deformation.

Static and Limit Subloading Surfaces
The static subloading and the limit subloading surfaces (see Fig. 4) are incorporated in addition to the subloading, the normal-yield, the elastic-core and the limit elastic-core surfaces in the subloading surface model for the rate-independent elastoplastic deformation. Here, the subloading surface on which the current stress lies may become larger than the normal-yield surface in general. The normal-yield ratio R which may become larger than unity in general is calculated by the identical equation for the rate-independent elastoplastic deformation, i.e. Equation (18) in general and by Eq. (92) for metals.
The viscoplastic strain rate is induced by the overstress from the static subloading surface expressed in the following equation which is given by replacing the current stress , the center and the normal-yield ratio R to their conjugate points s , s and the static normal-yield ratio R s ( ≤ 1) , respectively, in the subloading surface in Eq. (11) Fig. 3 Rheology models for creep model and overstress model R s ( ≤ 1) is called the static normal-yield ratio because the quasi-static deformation proceeds in the state R = R s . It can be analytically calculated by Eq. (22) with the replacement of the plastic strain rate to the viscoplastic strain rate for ℜ c C n = const. as follows: under the initial condition R s = R s0 for vp = vp 0 , where vp = ∫ ||̇ vp ||dt.

Viscoplastic Strain Rate
Now, let the flow rule of the viscoplastic strain rate be given by incorporating the crucially important variable in the subloading surface model, i.e. the normal-yield ratio R into the overstress model as follows: where is the positive viscoplastic multiplier given by v (viscoplastic coefficient) and n( > 1) are the material constants. The smooth elastic-viscoplastic transition is described by adopting the overstress due to R−R s instead of R−1 , because R s increases gradually up to 1 with the (90) viscoplastic deformation, while R s was fixed as R s = 1 in the past overstress model. The function U for the evolution rule of the normal-yield ratio R and the positive proportionality factor (plastic multiplier) • are modified to Eqs. (54) and (56), respectively, by incorporating the variables ℜ c ( ≡ f (̂ )∕F) and C n (≡̂ c ∶ ) in the extended subloading surface model for the rate-independent elastoplastic deformation. Analogously, for the ratedependent deformation, the function U s and the viscoplastic multiplier Γ be modified to Eq. (89) and where u c is the material constant. Thus, the rising and the lowering of the stress due to the decrease and the increase of strain rate are enforced in the state ℜ c C n > 0 and ℜ c C n < 0 , respectively. It leads to the description of the global Masing effect in the viscoplastic deformation process.
Further, let Eq. (93) be extended such that the infinite viscoplastic strain rate is induced for R → c m R s ( c m ( ≫ 1) : material constant) leading to → ∞ in the impact loading process as follows: The variable c m R s (R < c m R s ≤ c m ) is called the limit normal-yield ratio. Here, the surface to which the stress can reach at most is given by setting R = c m R s in the subloading surface in Eq. (11) is called the limit subloading surface. Fig. 4 Limit subloading, subloading, normal-yield, staticsubloading, limit elastic-core, elastic-core and static-subloading surfaces in subloadingoverstress model 1 3 On the other hand, the elastic response is described unrealistically for the impact loading in the past overstress models, Perzyna [114,115], Chaboche [14,17], de Souza Neto et al. [25], Lemaitre and Chaboche [86], etc. The rational description of deformation in the general rate ranging from the quasi-static to the impact loading is attained by the above-mentioned subloading-overstress model in which the crucially-important variables, i.e. the normal-yield ratio R , the static normal-yield ratio R s , ℜ c , C n and the limit normalyield ratio c m R s are incorporated. The rates of the internal variables H , and are given from Eqs. (35), (36) and (37) with the replacement of ̇ p to ̇ vp (= ) as follows: The isotropic hardening stagnation can be incorporated by the identical equations formulated in Sect. 2.6 with the replacement of the plastic strain rate ̇ p to the viscoplastic strain rate ̇ vp .

Strain Rate vs. Stress Rate Relation
The strain rate vs. stress rate relations are given from Eq. (85) with Eqs. (3) and (91) as follows: which is represented in the incremental form as follows: Then, it follows from Eqs. (99) with (92) in the quasistatic deformation that so that the stress changes along the static subloading surface given by Eq. (86). Consequently, the response of the subloading−overstress mode is reduced to that of the original subloading surface model for the rate−independent elastoplastic deformation behavior described in Sect. 2. The subloading−overstress model is no more than the generalization of the subloading surface model to the description of the elastoplastic deformation The stress-strain curve described by the subloadingoverstress model is illustrated in Fig. 5 The smooth elasticviscoplastic transition and the cyclic loading behavior can be described appropriately by the unified equation.

Verification of Subloading-Overstress Model by Simulations of Test Data
The validity of the subloading-overstress model will be verified by the comparisons with the test data of metals under various loading conditions in this section.

Material Parameters
The following 18 material constants and 5 initial values of the internal variables at most are involved in the subloadingoverstress model. The determinations of these material parameters are 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 stress-strain curve. (2) s r , c H and F 0 for the isotropic hardening and c k , b k and 0 for the kinematic hardening are determined from the stress-strain curves in the initial and the inverse loading process.
(3) u , u c and R e ( < 1) for the evolution of the normal-yield ratio are determined from the stress-strain curve in the subyield state, i.e. the elastic-plastic transitional state. (4) c e , and 0 for the elastic-core are determined from the stress-strain curves in cyclic loading. (5) C , , , K 0 and 0 for the isotropic hardening stagnation are determined from the stress-strain curves in loading-unloading process. (6) v , u c , n( > 1) and c m ( ≫ 1) are determined by the viscoplastic deformation characteristics.
All of these material parameters can be determined only by the stress-strain curves in the uniaxial loading for initial isotropic materials. One can put 0 = 0 = 0 = for the initial isotropy, which is assumed in all the subsequent simulations. We calculate K 0 by K 0 = ||̃ p 0 || leading to R 0 = 1 by inputting infinitesimal values of p 0 and 0 in order that the isotropic hardening rule in Eq. (78) without the isotropic hardening stagnation holds in the initial loading process for all the present simulations. Besides, the three material constants C(0 ≤ C ≤ 1), ( ≥ 1) and Elastic moduli: E(Young � s modulus), (Poisson � s ratio) Hardening isotropic: s r (limit isotropic hardening ratio), c H kinematic: c k , b k (limit kinematic hardening ratio) ( < 1) Evolution of normal -yield ratio ∶ u, u c , R e (fatigue limit normal -yield ratio)( < 1) Translationon of elastic -core: c e , (elastic -core limit normal -yield ratio)(< 1) Stagnation of isotropic hardening: Viscoplastic deformation: v (viscoplastic viscous coefficient), u c , n( > 1), c m (limit normal -yield ratio)(≫ 1) Normal -yield surface size: F 0 center: 0 Elastic -core ∶ 0 Normal -isotropic hardening surface size:K 0 center: 0 ( ≥ 1) for the isotropic hardening stagnation can be fixed as C = 0.5 , = 3 and = 3 in all the simulations shown later in this section. Consequently, the number of material constants which must be determined in each test datum can be reduced to be 15 as a total. The automatic determination function of the material parameters for the simplified version of the subloading surface model for the elastoplastic deformation of metals with the evolution rule of the elastic-core limited to the nonhardening state and without the isotropic hardening stagnation is installed in the commercial software Marc which can be used by all the Marc users. The efficient and accurate determinations of the material parameters for the simplified version can be referred to Liu et al. [91] adopting the artificial neural network.

Dynamic Loading Process Inducing Elasto-Viscoplastic Deformation
The uniaxial monotonic loading behavior at various constant strain rates at 600 °C in the test data for 21/4Cr-1Mo steel (SA 387, Gr.22) after Inoue et al. [72] is shown in Fig. 6, where the material parameters are chosen as follows: Material constants: Initial values: The close simulations are shown in Fig. 6, while the simulated stress for the strain rate 0.000001∕s is slightly lower than the stress in the test data.
The monotonic loading behavior at various constant strain rates at 550 • C in the test data for Modified 9Cr-1 Mo steel after Abel-Karim and Ohno [1] is shown in Fig. 7, where the material parameters are chosen as follows: Material constants: Initial values: The quite close simulations are also seen in Fig. 7. The accuracy of the present simulations are same levels of the simulations by the experimenters Abdel-Karim and Ohno [1] based the creep model, who have provided the test data, while the creep model is physically irrational lacking the generality, i.e. incapable of describing the rate-independent elastoplastic deformation behavior as explained in Introduction.
The stress relaxation behavior with same prior strain rate 0.0005∕s but the three levels of fixed strain, i.e. 0.139, 0.400 and 0.602% in the test data of TIMETAL21S at 650 °C after Saleeb and Arnold [121] is shown in Fig. 8, where the material parameters are chosen as follows: Material constants: Initial values: The quite close simulations are shown in Fig. 8 for the three fixed strain levels.
The deformation behavior at the variable strain rate between 0.001/s and 0.000001/s for 304 stainless steel at room temperature 20 °C in the test data after Krempl, cf. Lemaitre and Chaboche [86] is shown in Fig. 9, where the material parameters are chosen as follows: Material constants: Initial values: The quite close simulations are shown in Fig. 9a. The accuracies of the simulations hardly decrease by ignoring the isotropic hardening stagnation as seen in Fig. 9b.

Quasi-static Loading Process Inducing Elastoplastic Deformation Behaviors
It was verified in the previous section that the subloadingoverstress model is capable of describing the elasto-viscoplastic deformations. Further, it will be verified below that the subloading-overstress model is capable of describing the elastoplastic deformations by the simulations of test data for the quasi-static rate of deformation at the room temperature.
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 Drucker [27] with the yield surface enclosing the purely elastic domain, e.g. Chaboche model, Chaboche et al. [18], Chaboche [13,14], Ohno and Wang [102], etc. 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 Isotropic hardening function: F 0 = 160 MPa. 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 +830MPa at room temperature after Jiang and Zhang [74] is depicted in Fig. 10a, where the material parameters and the strain rate are selected as shown in the following, while the axial strain rate is chosen to be quite low, i.e. ̇ a = 10 −7 ∕s , although the strain rate is not written in the paper of Jiang and Zhang [74]. Material constants: The accuracy of simulation hardly decreases by ignoring the isotropic hardening stagnation as seen in Fig. 12b.
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 Chaboche et al. [18] is depicted in Fig. 11 where the material parameters and the strain rate are selected as shown in the following, while the axial strain rate is chosen to be quite low, i.e. ̇ a = 10 −8 ∕s , although the strain rate is not written in the paper of Chaboche et al. [18].
Material constants: Initial values: The quite close simulation is shown in Fig. 11a. However, the accuracy of simulation decreases seriously by ignoring the isotropic hardening stagnation as shown in Fig. 11b, differing from the simulations of the test data shown in Figs. 9 and 10. It can be inferred from these simulation results that the influence of the isotropic hardening stagnation is induced in the loading process in which the full inverse loading is repeated many times.
The calculated results for the axial strain rates (a) 10 −10 ∕s,(b) 10 −9 ∕s, (c) 10 −7 ∕s , (d) 10 −6 ∕s and (e)10 -5 /s are shown in Fig. 12. It is observed from this figures that the calculated results are identical for the strain rates slower than 10 −7 ∕s but the stress increases with the increase of the strain rate faster than 10 −7 ∕s . In other words, the quasi-static deformation leading to the elastoplastic deformation is induced for the strain rate slower than 10 −7 ∕s . Consequently, the deformation behavior of the 316 steel at room temperature under the strain rate faster than the quite slow strain rate 10 −7 ∕s may not be calculated accurately by the elastoplastic constitutive equation but must be calculated by the subloading-overstress model. The accuracies of simulations of the test data in Figs. 10 and 11 by Hashiguchi et al. [64] and Hashiguchi and Ueno [61] using the rate-independent subloading surface model for the elastoplastic deformation are improved to the simulations depicted in these figures by using the subloadingoverstress model.

Concluding Remarks and Future Works
The subloading-overstress model is provided and it is verified that the model is capable of describing not only the dynamic but also the quasi-static deformation behaviors by the comparison with various test data of metals under the monotonic and the cyclic loading processes at the various strain rates.
The elastoplastic deformation behaviors has been simulated by using the elastoplastic constitutive equation. However, it should be noticed that the elastoplastic constitutive equation holds only in the purely quasi-static deformation which does not occur actually. Therefore, an infinitely long time is required in order to perform a purely quasi-static deformation test and thus test data reported as the elastoplastic deformation behaviors should not be regarded as test data for a quasi-static deformation. In facts, the test data reported as the elastoplastic deformation behaviors can be closely simulated at the low strain rate by the subloadingoverstress model as verified in Sect. 4.3. Incidentally, it should be noted that the purely quasi-static deformation is hardly induced actually in test data and engineering practice. Therefore, the subloading-overstress model should be used for the deformation analyses in any rate of deformation including the quasi-static deformation, disusing the elastoplastic constitutive equation involving the plastic modulus M p which is rather complex as shown in Eq. (81). Besides, it is revealed that the isotropic hardening stagnation can be ignored in the deformation behavior as far as the full inverse loadings are not repeated many times.
The distinguished advantages of the subloading-overstress model are summarized as follows: i. The plastic/viscoplastic strain rate can be described in any low stress level, while it cannot be described in all the other plastic/viscoplastic models because they assume the purely-elastic domain inside which the stress lies always. Incidentally, the incorporation of the subloading-overstress model would be inevita-1 3 ble for the analysis of fatigue phenomenon without resorting to use of the unnatural method by the twoscale damage model premised on the plasticity model assuming the purely-elastic domain, Lemaitre [84,85], Lemaitre and Doghri [87], Lemaitre et al. [88], Lemaitre and Desmorat [89], etc. based on the law of micromechanics, Eshelby [29], in which the inclusion of the microscopic voids or cracks is incorporated in a mesoscale representative volume element. ii. The continuity condition in the large, i.e. the smoothness condition, Hashiguchi [41,42,44] is fulfilled always leading to the continuous variations of the tangent stiffness tensor and the isotropic hardening by virtue of the subloading concept, while it is violated leading to the discontinuous variations of the tangent stiffness modu-lus tensor and the isotropic hardening in all the other plasticity and the viscoplasticity models because of the incorporation of the purely-elastic domain. iii. The elastoplastic and the elasto-viscoplastic deformations under the monotonic and the cyclic loading at the general deformation rate ranging from the quasi-static to the impact loading can be described in the unified formulation. Therefore, the elastoplastic constitutive equation can be disused.
It should be emphasized that these distinguished advantages in the subloading-overstress model is provided by the incorporation of the basic variable incorporated uniquely in the subloading surface variable, i.e. the normal-yield ratio R. The following extensions of the subloading-overstress model are desirable to describe a more general deformation behavior, while the formulations and the verifications with test data will have to be performed in future works.
(1) The elastoplastic deformation behavior is influenced by the temperature. However, the dependence of material functions and constants in the elastoplastic constitutive equation on the temperature would be complex so that they are specified for each temperature resulting in the piecewise-linear relation in the deformation analysis as seen in the research papers, Kou [81], Li et al. [90], Ohno et al. [103], Pandya et al. [109], Yang et al. [138], etc. and the commercial software (e.g. Marc and Abaqus). Such method for incorporation of the piecewise-linear relation is easy-going way but requires the large computational load. Therefore, it is desirable to formulate material functions and constants as functions of temperature. (2) It is desirable to perform the simulation of test data in an impact loading behavior. However, it would be difficult to obtain an accurate test data of stress-strain relation since an intense vibrational deformation is induced in test specimen at an impact loading so that it would have to be analyzed as the dynamic boundary-value problem. It will be executed in the future work. (3) The present formulation will have to be extended to the multiplicative hyperelastic-based viscoplasticity for the description of the monotonic and the cyclic loading behaviors by improving the past formulations, Hashiguchi [48][49][50], Hashiguchi and Yamakawa [65], Fincato and Tsutsumi [33] which are incapable of describing the viscoplastic deformation pertinently predicting the peculiar stress-strain curve with an inflection point even in the monotonic loading process as delineated in the introduction. Here, it should be noticed that the unconventional elastoplasticity models, i.e. the multi surface model, Mróz [96], Iwan [73] and the two surface model, Dafalias and Popov [23], Krieg [83] other than the subloading surface model cannot be extended to the hyperelastic-based plasticity model. In fact, only the conventional elastoplasticity model with the yield surface enclosing the purely-elastic domain is adopted to the formulation of the multiplicative hyperelasticbased plastic constitutive equations, Wallin et al. [131,132], Dettmer and Reese [26], Wallin and Ristinmaa [131], Vladimirov et al. [129,130], etc. which are capable of describing the finite elastoplastic deformation only during the monotonic loading process. Besides, it should be noted that the hypoelastic-based formulations, Truesdell [126], Truesdell and Noll [127], Hill [67], etc. are capable of describing the finite plastic deformation but they are limited to the description of the infinitesimal elastic deformation as explained in detail by Hashiguchi [48,50] and shown by the numerical examples by Brepopols et al. [11] and Jiao and Fish [75,76] and required the cumbersome time-integration of stress rate in numerical calculation (cf. Hughes and Winget [71], de Souza Neto et al. [25], etc.). In addition, it should be noticed that the multiplicative hyperelastic-based plastic constitutive equation is irrelevant to the hypoelastic-based plastic and viscoplastic constitutive equations but it can be derived from the infinitesimal hyperelastic-based plastic and viscoplastic constitutive equations formulated in this article (cf. Hashiguchi [48][49][50]) by the additive to the multiplicative transformation of the deformation measures.
Funding No funding was received to assist with the preparation of this manuscript.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.