1 INTRODUCTION

The general description of material nonlinearity primarily concerns a constitutive equation. Essential ingredients to establish this constitutive equation can be characterized by yield surface (or yield stress function), plastic flow rule, normality rule and hardening law- the yield surface divides the elastic and plastic region; the flow rule relates the plastic strain to stress state; the normality rule leads the incremental plastic strain to the normal level of the yield surface; the hardening law describes the plastic evolution. While the classical rate independent plasticity (so called plasticity) is often illustrated by a single yield surface with either isotropic hardening or kinematic hardening, these models can not completely express the hardening behavior of material. In particular, in the single yield surface model, the elastic domain is considered to be too large compared to experiment results and it is difficult to describe sudden changes from elastic to plastic or plastic to elastic region ((^{Chaboche, 2008})). Also, it cannot solely account for cycle-by-cycle accumulation of permanent deformation (ratcheting), which often happen in engineering practice ((Bari and Hassan, 2008)). In order to account for more realistic response of plasticity, previously, two-surface and multiple yield surface models have been developed by (^{Mroz, 1967}), (^{Dafalias and Popov, 1975}), (^{Krieg, 1975}), and (^{Banerjee et al., 1987}), where key ideas reside in combining the kinematic and the isotropic hardening rules. Among these models, (^{Dargush and Soong, 1995}) tested the two-surface model by (^{Banerjee et al., 1987}) with application to metallic plate dampers, which shows excellent agreement to experimental force-displacement data provided by (^{Tsai et el, 1993}).

For much elaborate description of material nonlinearity with consideration of strain/stress rate effects under various loading conditions like monotonic, cyclic, and more complex transient dynamic loadings, there have been numerous researches especially for the development of rate-dependent plastic material models and their practical applications in engineering problems. For examples, (^{Bodner and Partom, 1972}) established a theory of rate-dependent plasticity (or viscoplasticity); (^{Chaboche, 1977}, ^{1989}) developed a viscoplastic constitutive model with nonlinear kinematic hardening- this model is applied to the 316 stainless steel under cyclic loading and creep relaxation ENT#091;(^{Chaboche and Rousselier, 1983b})ENT#093;. Later, it is also used for the simulation of ratcheting ENT#091;(^{Chaboche, 1991})ENT#093;; (^{McDowell, 1992}) developed viscoplastic nonlinear kinematic hardening model under thermomechanical cyclic conditions, while (^{Tanaka, 1994}) developed a viscoplastic constitutive model under non-proportional loading; (^{Ohno and Wang, 1993}) modified Armstrong-Frederick model ENT#091;(^{Chaboche, 1986})ENT#093; with dynamic recovery term. (^{Tanaka and Yamada, 1993}), (^{Abdel-Karim and Ohno, 2000}) continued with this study especially on the nonlinear kinematic hardening model with steady-state ratchetting, while Chaboche-type nonlinear kinematic hardening models for ratchetting are investigated by (^{McDowell, 1995}), (^{Kang et al., 2001}), (^{Kang et al., 2002}), (^{Kang el al., 2004}), (^{Yaguchi and Takahashi, 2000}), (^{Yaguchi and Takahashi, 2005}), and (^{Kang et al., 2006}). Main differences in such various viscoplastic models and the classical rate-independent plasticity can be described by overstress and time-dependent behavior. Thus, while elastic strain and strain hardening rules in the viscoplasticity are the same as those in the classical rate-independent plasticity, the stress state goes beyond the elasticity domain whereas this overstress is not allowed in the classical rate-independent plasticity. Also, the viscoplasticity model can describe creep phenomenon (time dependent irreversible deformation for long term response), and it can account for the rate of loading in strain-stress responses, which may become a major issue for earthquake excitation and high velocity impact.

In this paper, we extend the rate independent two-surface model by (^{Banerjee et al., 1987}) to the rate dependent two-surface model with consideration of rate-effects. The proposed model is then implemented in commercial finite element software (ABAQUS, 2008) by using a user subroutine (UMAT), and some representative examples are considered to elucidate the features of the proposed model.

2 A RATE INDEPENDENT TWO-SURFACE PLASTICITY MODEL

In the plasticity, a total strain at a given stress can be decomposed into two parts, which correspond to an elastic strain and plastic strain. For the multiaxial case, this can be generalized as a rate-form as

with the superposed dot indicating a derivative with respect to the time. Thus, and representing the elastic strain rate and the plastic strain rate.

In Eq. (1), the elastic strain rate is related to a stress rate with a fourth-order elastic constitutive tensor _{
Ceijkl
} as

With Eq. (1), Eq. (2) can be equivalently written as

As in the classical rate-independent plasticity, the evolution of plastic strain (_{
εpij
} ) and that of yield stress function (or yield surface) are decided by flow rules and hardening laws in the two-surface model. However, this time, there are two yield surfaces, where the kinematic yield surface (or loading surface) in principal stress axes resides inside the isotropic yield surface (or bounding surface) as shown in Figure 1.

For the isotropic material, the constitutive relation in this two-surface model of plasticity is identifed as

in a rate-form of elastic response, where a stress point moves until it reaches the inner yield surface and strains are fully recoverable. Also, when the stress state stays on the inner yielding surface but resides inside the outer yielding surface (this is usually called transition region or the meta-elastic region), this rate-constitutive relation becomes

Further, as the stress state becomes larger, the inner yield surface continuously approaches to the bounding surface and it finally touches the bounding surface and make this bounding surface expand. For this case, the rate-constitutive relation becomes

In Eqs. (4)-(6), λ is lame’s first parameter; *μ* is shear modulus; _{
σLy
} is the inner yield stress; _{
σBy
} is the outer yield stress; _{
Sij
} is the deviatoric stress; is _{
Sij
} - _{
αij
} representing the deviatoric stress minus the back stress; _{
Hp
} = _{
hB
} is an isotropic hardening modulus dependent on parameters h^{B} (h^{B}
_{0}+ h^{B}
_{1}
_{
σBy
} ), _{
hB
}
_{0}, _{
hB
}
_{1}, *n*
_{1}, and ξ = Here, we only summarize main formulation for the rate-independent two-surface plasticity model. Interested readers can refer to (^{Banerjee et al., 1987}), (^{Chopra and Dargush, 1994}), and (^{Sant, 2002}).

As described in Figure 1, geometrically, the kinematic hardening of the inner yield surface depends on a vector that joins the stress state to the bounding surface, ξ_{
ij
} . The location of A^{B} is determined by drawing a vector O'AA^{B} parallel to OA, _{
σBij
} . The direction of α_{
ij
} is then determined to be paralled to AA^{B}. The inner surface, which separates the elastic range and inelastic range, is composed of its center and radius expressed by the back stress (α_{
ij
} ) and inner yield stress (_{
σLy
} ). Meanwhile, the outer surface, which always contains the inner yield surface, is located on the center of stress space with a radius represented by the outer yield stress (_{
σBy
} ). The translation of inner surface corresponds to kinematic hardening, while the expansion of outer surface produces isotropic hardening.

3 A RATE DEPENDENT TWO-SURFACE MODEL

3.1 Formulation

Unlike earlier single yield surface models in viscoplasticity that assume all inelastic strain to be rate dependent, in a present two-surface model, the inelastic strain is divided into rate-independent and rate-dependent parts as

with and representing the inelastic strain-rate and the viscoplastic strain-rate, respectively.

In the past, the above approach had been proven to be effective when accounting for the hysteresis loop with the rounded corners in stress-transition from elastic to inelastic or vice versa ENT#091;(^{Bradley and Yuen, 1983}); (^{Tirpitz and Schwesig, 1992}); (^{Lubliner, 2008})ENT#093;. In addition, as in the existing two-surface plasticity model, a present two-surface viscoplasticity model utilizes the following yield functions such as

and

depending on kinematic hardening and isotropic hardening rules. Here, _{
σe
} is the von Mises effective stress and the symbol “:” designates the product contracted twice. Thus, three separate regions in Figure 1 are identified with these yield functions: _{
gL
} < 0 (elasticity); _{
gL
} > 0 and _{
gB
} < 0 (meta-elasticity); _{
gB
} > 0 (both kinematic and isotropic hardening rules are effective). Once the stress-state is identified with these criteria, each portion of the inelastic strain is given by flow rule and nomality hypothesis:

In Eq. (10), the magnitude of the plastic strain increment can be decided by the consistency condition _{
gI
} = _{
ġI = 0
} . Otherwise, the magnitude of the viscoplastic strain increment is determined from a potential *φ*

that is the function of the stress _{
σij
} , the back stress α_{
ij
} , and the magnitude of bounding surface _{
σBy
} .

With adoption of the hyperbolic sine function for *φ* (^{Dunne and Petrinic, 2005}), can be explicitly written as

In Eq. (13), *C* and *B* represent material parameters associated with viscosity. Also, *n*
_{2} is a controlling parameter for strain-rate sensitivity ENT#091;(^{Bodner and Partom, 1972}), (^{Chaboche, 1989})ENT#093; and represents the deviatoric stress that can be either or _{
Sij
} .

Thus, the viscoplastic strain rate is identified as

and consequently, the inelastic strain is given by

Overall, in the present two-surface viscoplasticity model, the consitutive relation of isotropic material results in

for the elastic response, kinematic hardening response, and both kinematic and isotropic hardening response, respectively.

3.2 Finite Element Implementation

The present two-surface model for viscoplasticity was implemented by a user subroutine (UMAT) in the (ABAQUS, 2008). Once the small increment of strain is given, new updated state variables such as stress, back stress, plastic strain, viscoplastic strain are obtained by integrating the constitutive equations. For this numerical simulation, a higher-order adaptive step size Runge-Kutta method is also used to integrate the constitutive equations until a high level of accuracy is obtained. All the equations are expressed by incremental form, which are easily implemented in the Fortran source code. Resulting equations become

where _{
eij
} is deviatoric strain tensor, and *n* and *n+1* represent time-steps.

Subtractig Eq.(21) from Eq.(20) with the strain decomposition, and the relationship between strain and deviatoric strain in Eq.(22), one finds

Using a plastic flow rule and viscoplastic flow rule, Eq.(23) yields

Finally, one can obtain the following incremental formulation.

Eq. (25) is numerically into the commercial finite element code, ABAQUS, where a certain iteration scheme is introduced to solve nonlinear equations at the element and global system level ((^{Bathe and Cimento, 1980})). The following is a brief explanation about how the nonlinear solver works (here, subscript indices are omitted to avoid complexity). First, all the variables are initialized with corresponding time-step. Then, the solution (_{
tU
} , nodal displacements) and other internal variables such as strain (_{
tε
} ) and stress (_{
tσ
} ) are stored at time (*t*), where the iteration counter (*i*) is fixed as 1. Variables at the time-step of (*t*+∆*t*) are calculated and updated through iterations with convergence criteria and/or a maximum number of iterations.

For example, the stress is updated from a known converged solution as

At the element level, the tangent constitutive matrix (^{
t
}
^{+∆}
_{
tCi-
}
^{1}) is updated with a certain iterative scheme, and consequently, the global stiffness matrix (^{
t
}
^{+-∆}
_{
tCi-
}
^{1}) and nodal force vector (^{
t
}
^{+∆}
_{
tFi-
}
^{1}) corresponding to the interanl element stresses (^{
t
}
^{+∆}
_{
tσi-
}
^{1}) are computed by using a gaussian quadrature for each element’s integration and collecting them in the following manner.

At the global system level, the solution is obtained in terms of the incremental nodal displacements (∆_{
Ui
} ) by solving the following set of equations

with ^{
t
}
^{+∆}
_{
tR
} representing the externally applied nodal force vector to the time-step, *t*+∆*t*. In solving the global solution of incremental displacement ^{
t-
}
^{1}∆_{
Ui
} , interative scheme is also utilized to update ^{
t
}
^{+∆}
_{
tKi-
}
^{1} along with convergence criteria and/or a maximum number of iterations. Finally, the numerical solution of displacements at the time-step (*t*+∆*t*) is updated by

and strains are calculated from these updated displacements.

In the present work, we adopt the Cash-Karp method (one of embedded Runge-Kutta methods) as the iterative method for updating ^{
t
}
^{+∆}
_{
tCi-
}
^{1}, while the Newton-Raphson iterative method is used to solve the global system equation with updating ^{
t
}
^{+∆}
_{
tKi-
}
^{1}.

4 NUMERICAL EXAMPLES

In this section, we verify the present two-surface model for viscoplasticity with representative examples including both monotonic and cyclic loading cases, where both experimental results by (^{Chang, 1985}) and numerical simulation results by the kinematic hardening model are provided for the reference.

4.1 Experimental Results

For the monotonic loading case, Chang experimentally tests A36 structural steel with three different strain rates of 10^{-3}/sec, 10^{-4}/sec, and 10^{-5}/sec, in each of which loading was increased until 2% axial strain was achieved. Also, in cyclic tests, there are four different types of loading as shown in Figure 2, where all the specimen were loaded in the axial direction and the loading continued until it was stabilized. In particular, for the case 1 and the case 2, the loading was increased until 0.8% axial strain was attained, while the loading was increased until 0.6%, 1.2%, and 1.5% axial strain were obtained for the case 3 and the case 4.

Figure 3-Figure 6 show experimental results of axial stress and axial strain for each loading case, where results from the loading type 3 and the loading type 4 are separately provided for the sake of clarification.

As shown in Figure 3, three different values of maximum were observed in the monotonic test: 32.5 ksi after initial peak (strain rate 10^{-3}/sec); 30ksi (strain rate: 10^{-4}/sec); 28.5 ksi (strain rate: 10^{-5}/sec). This indicates that the faster strain rate gives the higher stress after yielding. Also, it is observed that the higher strain rate gives the longer plastic plateau and strain hardening effects are negligible in these monotonic tests. With reference to the strain rate of 10^{-5}/sec, the yield stress at a strain rate of 10^{-4}/sec was increased by 5%, and the yield stress at strain rate of 10^{-3}/sec was increased by 14%.

For cyclic tests as shown in Figure 4~Figure 6, the plastic plateau is no longer observed, however, again, the faster strain rate yields the higher stress. In particular, maximum stresses in the loading type 1 (10^{-4}/sec) and the loading type 2 (10^{-2}/sec) are given by 38.5 ksi and 40.8 ksi, respectively. Also, for the loading type 3 (10^{-4}/sec) and the loading type 4 (strain rate: 5 x 10^{-3}/sec), maximum stresses are identified as: 34 ksi and 37.5 ksi at the 0.6% axial strain; 39 kis and 42.1 ksi at the 1.2% axial strain; 39 ksi and 44.7 ksi at the axial 1.5% strain.

4.2 Numerical Simulation Results

With a user subroutine (UMAT) in the ABAQUS, both rate-dependent nonlinear kinematic hardening model two-surface models are implemented by the authors. For fair comparison between these two models, we employ one four-node bilinear axisymmetric element (CAX4) for the cylinder-type specimen in every numerical simulation examples, as shown in Figure 7. Also, in numerical simulation, Newton-Raphson method is adopted with fixing the maximum number of iteration as 1000 and incremental displacement as 10-9 under such conditions, we checked that all the numerical solutions satisfy convergence criteria of residual force (less than 10-5) and displacement (less than 10-8), respectively.

4.2.1 Nonlinear Kinematic Hardening Model

The main difference between the developed two-surface model and the nonlinear kinematic hardening model ENT#091;(^{Chaboche, 2008})ENT#093; for the rate-dependent plasticity is that von Mises yield surface cannot expand in the nonlinear kinematic hardening model, while the developed two-surface model allows both expansion and translation of von Mises yield surface. In addition, in the developed two-surface model, we differentiate inelastic strain as the rate-dependent () and the rate-independent part () as in Eq. (7). Thus, we have the following equations for the nonlinear kinematic hardening model

whereas the developed two-surface model has Eqs. (15)-(18).

In numerical simulation of the nonlinear kinematic hardening model, we take material properties and parameters as shown in Table 1.

Young's modulus(E): 28,500 ksi(196,500 MPa), Poisson's ratio (υ): 0.35 |

Yield stress(σ_{γ}): 30 ksi (206.84 MPa) |

Material parameters of the hyperbolic sine function: C = 1.E+2; B = 1.E-8; n
_{2} = 5.0 |

As shown in Table 1, the yielding stress is specified as 30 ksi (206.84 MPa) for A36 steel, and this value is adopted in numerical simulation with calibration of experimental results from (^{Chang, 1985}).

Figure 8-Figure 10 show numerical simulation results obtained from the kinematic hardening model. In Figure 8, maximum stress is computed as 32.3 ksi (strain rate 10^{-3}/sec), 30ksi (strain rate: 10^{-4}/sec), and 28.6 ksi (strain rate: 10^{-5}/sec) for the monotonic loading. Also, in cyclic tests, we have maximum stress of 45.6 ksi (loading type 1) and 46.4 ksi (loading type 2). For the other loading types, maximum stresses at each axial strain are obtained as: 33.52 ksi (loading type 4) and 33.51 ksi (loading type 3) at the 0.6% axial strain; 41.8 ksi (loading type 4) and 41.7 ksi (loading type 3) at the 1.2% strain; 46.9 ksi (loading type 4) and 42.2 ksi (loading type 3) at the 1.5% strain.Thus, compared to experimental results, one can check that the nonlinear kinematic hardening model yields reasonably good results in terms of maximum stress for every case. However, this model gives almost linear shape of hysteresis in every stress-strain response that is not well matched to experimental results.

4.2.2 Two-Surface Model

For numerical simulation of the two-surface model, basically, the same parameters in Table 1 are utilized. Additional parameters adopted in numerical simulation to account for combined effects of both kinematic nd isotropic hardening are summarized in Table 2. In particular, both inner yielding stress and outer yielding stress are introduced here instead of a single yielding stress. Again, every parameter in Table 2 is calibrated from experimental results by (^{Chang, 1985}), here.

Young's modulus(E): 28,500 ksi(196,500 MPa), Poisson's ratio (υ): 0.35 |

Inner yielding stress (σ_{
Ly
} ) = 30 ksi (206.84 MPa) |

Outer yielding stress (σ_{
By
} ) = 60.91 ksi (420MPa) |

Hardening parameters: _{
HB
}
_{0} = 1650 MPa; _{
hB
}
_{1} = -8.47; n
_{1} = -10.4 |

Material parameters of the hyperbolic sine function: C = 1.E+2; B = 1.E-8; n
_{2} = 5.0 |

Figure 12-Figure 15 show numerical simulation results obtained from the two-surface model. For the monotonic loading test as shown in Figure 12, numerical results from the two-surface model are exactly the same as those from the rate dependent nonlinear kinematic hardening model. This may indicate that magnitude of the loading and rate of the loading are not large enough to make every response go beyond the outer yielding surface. On the other hand, for cyclic tests, where every response occurs beyond the inner and outer yielding surfaces, the two-surface model is more appropriate than the kinematic hardening model in accounting for both shape of hysteresis and maximum stress values. Thus, for the loading type 1 and the loading type 2, we have maximum stress of MPa (42.5ksi) and Mpa (44.9ksi). Also, for the loading type 3 and the loading type 4, maximum stresses at each axial strain are obtained as: 33.5 ksi (loading type 4) and 37.5 ksi (loading type 3) at the 0.6% axial strain; 41.8 ksi (loading type 4) and 42.1 ksi (loading type 3) at the 1.2% strain; 46.9 ksi (loading type 4) and 44.7 ksi (loading type 3) at the 1.5% strain. All these maximum stresses are much closer to experimental results, compared to those from the nonlinear kinematic hardening model. Also, for all the cases, the two-surface model gives reasonably good rounded shape of hysteresis as detected in experiments.

5 CONCLUSIONS

In this paper, we develop a two-surface model for rate-dependent plasticity as an extension of the two-surface model in the rate-independent plasticity by (^{Banerjee et al., 1987}). The model combines both isotropic and kinematic hardening rules, where the constitutive relation is modified to account for rate-effects. In particular, this model has three separate regions with two bounding surfaces. Thus, when the stress remains in the lower bound, the response is elastic. Also, when stress exceeds the upper bound, the response shows both kinematic and isotropic hardening. In the middle range where the stress resides beteween the lower bound and the upper bound, the response only follows the kinematic hardening rule. With use of a subroutine UMAT, the present model is implemted in the commercial finite element software, ABAQUS. Then, this model is validated through both monotic and cycling loading cases with comparison to experiments and the nonlinear kinematic hardening model. The present model shows excellent agreement with experiments in both maximum stress and shape of hysteresis, while the nonlinear kinematic hardening model is not suitable to account for the shape of hysteresis.