A Two-Surface Viscoplastic Model for the Structural Steel 1

As extension of the previous two-surface model in plasticity, a two-surface model for viscoplasticity is presented herein. In order to validate and investigate the performance of the proposed model, several numerical simulations are undertaken especially for structural steel under monotonic and cyclic loading cases, where experimental results and numerical results from the rate dependent kinematic hardening model are also provided for the reference. For all the cases studied, the proposed model can appropriately account for the rate-effects in both maximum stress and hysteretic shapes.


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 Latin American Journal of Solids and Structures 14 (2017) 1000-1016 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 (1977Chaboche ( , 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 [Chaboche and Rousselier (1983b)].Later, it is also used for the simulation of ratcheting [Chaboche (1991)]; 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 [Chaboche (1986)] with dynamic recovery term.Tanaka andYamada (1993), Abdel-Karim andOhno (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.

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 Eq. (1), Eq. ( 2) can be equivalently written as As in the classical rate-independent plasticity, the evolution of plastic strain ( p ij e ) and that of yield stress function (or yield surface) are decided by flow rules and hardening laws in the twosurface 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 Latin American Journal of Solids and Structures 14 (2017) 1000-1016 resides inside the outer yielding surface (this is usually called transition region or the meta-elastic region), this rate-constitutive relation becomes (5) 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), l is lame's first parameter; m is shear modulus; L y s is the inner yield stress; representing the deviatoric stress minus the back stress; Here, we only summarize main formulation for the rate-independent two-surface plasticity model.Interested readers can refer to Banerjee et al. (1987), Chopra andDargush (1994), andSant (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 x .The location of A B is determined by drawing a vector O'AA B parallel to OA, B ij s .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 a ) and inner yield stress ( L y s ).
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 ( B y s ).The translation of inner surface corresponds to kinematic hardening, while the expansion of outer surface produces isotropic hardening.

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 Latin American Journal of Solids and Structures 14 (2017) 1000-1016 with in ij e  and vp ij e  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 [Bradley and Yuen (1983 ); Tirpitz and Schwesig (1992); Lubliner (2008)].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 s 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: L g < 0 (elasticity); L g > 0 and B g < 0 (metaelasticity); B g > 0 (both kinematic and isotropic hardening rules are effective).Once the stressstate 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 l  can be decided by the consistency condition 0 . Otherwise, the magnitude of the viscoplastic strain increment p  is determined from a potential j ( , , ) that is the function of the stress ij s , the back stress ij a , and the magnitude of bounding surface In Eq. ( 13), C and B represent material parameters associated with viscosity.Also, 2 n is a con- trolling parameter for strain-rate sensitivity [Bodner and Partom (1972), Chaboche (1989)] and ij S  represents the deviatoric stress that can be either ij S or ij S .
Thus, the viscoplastic strain rate is identified as where or 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.

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 Latin American Journal of Solids and Structures 14 (2017) 1000-1016 where eij is deviatoric strain tensor, and n and n+1 represent time-steps.Subtracting 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 ( t U , nodal displacements) and other internal variables such as strain ( t e ) and stress ( t s ) 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 ) are computed by using a gaussian quadrature for each element's integration and collecting them in the following manner.
Latin American Journal of Solids and Structures 14 (2017) 1000-1016 At the global system level, the solution is obtained in terms of the incremental nodal displacements ( i U D ) by solving the following set of equations with t t R +D representing the externally applied nodal force vector to the time-step, t+ t  .In solving the global solution of incremental displacement 1 i i U -D , interative scheme is also utilized to update -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 , while the Newton-Raphson iterative method is used to solve the global system equation with updating

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.

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%.

Latin
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.

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.

Nonlinear Kinematic Hardening Model
The main difference between the developed two-surface model and the nonlinear kinematic hardening model [Chaboche, 2008] 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 ( 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.As shown in Table 1, the yielding stress is specified as 30 ksi (206.84MPa) 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.

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.  8.47; n1 = -10.4Material parameters of the hyperbolic sine function: C = 1.E+2; B= 1.E-8; n2 = 5.0 Table 2: Additional material properties and model parameters adopted in the two-surface model.
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 Latin American Journal of Solids and Structures 14 (2017) 1000-1016 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.

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.

Figure 1 :
Figure 1: A two-surface model in the rate-independent plasticity.

Figure 9 :
Figure 9: Numerical results under cyclic loadings of type 1 and type 2 (nonlinear kinematic hardening model).

Figure 13 :
Figure 13: Numerical results under cyclic loadings of type 1 and type 2 (two-surface model).