EXPLICIT AND IMPLICIT INTEGRATION OF CONSTITUTIVE EQUATIONS OF CHABOCHE ISOTROPIC-KINEMATIC HARDENING MATERIAL MODEL

The proper selection of the material model and its parameters in numerical simulations of the material response under loading is a very important task. In the cyclic load one should assume different phenomena, e.g. ratchetting effect, the mean stress relaxation, creep, etc. The Chaboche kinematic hardening model is mainly applied in order to capture the Bauschinger effect due to its good description of the material behaviour under the cyclic loading. The modified Chaboche isotropic-kinematic hardening model is considered in this paper in order to simulate the strain-controlled and stress-controlled uniaxial cyclic tension-compression test. The implicit and explicit integration procedures are tested here, showing advantages and disadvantages of both methods. It is shown that both approaches provide similar results. The implicit integration method of constitutive equations is unconditionally stable and thus less calculation steps are required in comparison to the explicit procedure. The implicit and explicit integration of the Chaboche model including isotropic and kinematic hardening are then implemented for the 3D case in commercial ABAQUS program in the form of the user material procedure. The correctness of results obtained for the user material model is verified by comparison to results for commercially implemented Chaboche model.


INTRODUCTION
An accurate description of the effect of plastic strain on the material is necessary for a reliable prediction of its response beyond the elastic limit.It is especially essential in contact mechanics problems, fatigue analyses, etc. [1,2].The proper selection of the model in numerical calculations which should simulate the material behaviour under loading is a very important issue.It is essential in numerical simulations of cyclic loading problems in which different phenomena, e.g.Bauschinger effect, ratcheting, mean stress relaxation, elastic or plastic shakedown might occur.So far different models have been tested in numerical analyses of cyclic loading problems.Hashemi et al. [3] used the Prager-Ziegler model in numerical simulations of beams subjected to cyclic loading.In [4], the Frederick-Armstrong kinematic hardening model is applied in order to analyze the behaviour of a beam under the cyclic loading test.The Ohno-Wang [5], Johnson-Cook [6], Dafalias-Popov [7], Drucker and Palgen [8], Tseng and Lee [9] and Chaboche [10] models are also tested.The Chaboche model is often used to model the plasticity of metals and their alloys.It is also commonly applied in numerical calculations of cyclic loading problems.The model was originally introduced by Chaboche as the extension of the Frederick and Armstrong (F-A) one by adding additional backstress components of different properties.The original Chaboche model includes only the kinematic hardening of the material which is characterized by the translation of the yield surface in a stress space.The Chaboche model can be also combined with the Voce one assuming the isotropic hardening associated with the expansion of the yield surface.The combination of a Chaboche kinematic and Voce isotropic hardening models results in the modified Chaboche isotropic kinematic hardening (CIKH) model.
Both the original Chaboche kinematic hardening, as well as, the Chaboche isotropic-kinematic hardening models have been applied for solving different engineering problems, including materials forming processes.In [11], the Chaboche model is applied to simulate of the ratcheting phenomena.Zobec and Klemens [12] used the Chaboche model in the numerical analysis of the stress relaxation.The application of the Chaboche model for modelling the shakedown effect is described in [13,14].The proper selection of kinematic and hardening parameters for the CIKH model is a very difficult and time-consuming task.Although the relatively small number of parameters is necessary for the Chaboche model and its modifications, the identification of back stresses often requires solving the optimization problem in order to fit experimental and numerical hysteresis curves.Different methods of the determination of hardening parameters have been applied over the years.Most of them are based on the experimental data obtained under the uniaxial strain-controlled cyclic loading tests [15][16][17].In the author's previous research [18,19], the procedure of the determination of the parameters for the Chaboche-Lemaitre and Voce kinematic-isotropic hardening models based on the half of the first or the last stabilized cycle from the hysteresis curve is described.The explicit, implicit, and semi-implicit integration methods are mainly used to simulate the material response subjected to the external load.In [20], the implicit radial return method with Newton-Raphson iterations is applied for numerical simulations of the ratcheting with the use of the Chaboche kinematic hardening model.The implicit integration of the Chaboche model is also described in [21,22].The application of the explicit integration procedure for the Chaboche model in solving elasticplastic problems, including cyclic plasticity, is contained in [23,25].
Both explicit and implicit integration approaches are described and tested in this paper.The comparison of the explicit and implicit integration procedures for simulating the stress-strain response of material under the total strain-controlled and stresscontrolled cyclic tension-compression tests are presented here.The results obtained using both implicit and explicit integration procedures are compared.It is noted that the implicit integration method gives similar results to the explicit one using less computation increments.

CONSTITUTIVE EQUATIONS OF CHABOCHE ISO-TROPIC-KINEMATIC HARDENING MODEL
The main constitutive equations for the Chaboche isotropic-kinematic hardening (CIKH) model for a three-dimensional problem are described below.The total strain tensor  can be written as a sum of the elastic and plastic parts in line with Eq. 1. ) where:   and   are elastic and plastic parts of a strain tensor, respectively.
In the incremental form, Eq. 1 can be rewritten as (Eq.2): where: ,   and   are the total, elastic and plastic strain increments, respectively.The Hook's law which relates the stress and elastic strain, written in an incremental form, is as follows (Eq.3). ) where:  is the stress increment and  is the elastic constitutive matrix.
The evolution of the yield surface in line with the von Mises yield criterion, assuming both isotropic and kinematic hardening of the material, is written in line with Eq. 4. where: where   is the deviatoric stress,  is a back stress component defining the current center of the yield surface in the stress space,   is the yield stress, () is the isotropic hardening function describing the increase of the yield surface size.The plastic flow occurs when  ≥ 0, otherwise the response of the material is only elastic.

𝑓 = √
For the one-dimensional case, the stress tensor  might be described (Eq.5): ) The consistency condition used in the explicit integration approach is as follows (Eq.6).where:  and  are the backstress and the effective plastic strain increments, respectively.Effective plastic strain increment is (Eq.7).where:   is the plastic strain increment (Eq.8).

𝑑𝑝 = √
Focusing on only a specimen under uniaxial loading, the incremental plastic strain tensor due to incompressibility condition is given in Eq. 9. ) Assuming that  equals the plastic multiplier , the normality condition might be rewritten as (Eq.10): ) where: defines the normal to the yield surface .The normality condition can be rewritten in line with Eq. 11, therefore.
Using plastic multiplier  and the normality condition, the incremental Hook's law can be also written as (Eq.12): The Chaboche model consists of some back stress components   for which the incremental form is written as follows (Eq.13): ) where:   is kinematic hardening parameter representing the translation rate of the yield surface and   is the kinematic parameter determining the relaxation rate of the yield surface translation due to accumulation of the plastic deformation.In the Chaboche model the superposition of the back stresses is assumed (Eq.14): For an uniaxial case, the back stress might be expresses as follows (Eq.15): where:  0 means the initial yield stress and () is isotropic hardening function written as (Eq.17): ) where:  and  are material constants.The  describes the saturated value of () and  determines the rate at which saturation is achieved.

Explicit integration
The explicit integration procedure of constitutive equations of Chaboche isotropic-kinematic hardening (CIKH) model is described in this section.In explicit procedure, the known solution for the time step   is extrapolated to the time  +1 =   + ∆.The first-order forward Euler explicit integration algorithm is used here.It is worth noting that this procedure is conditionally stable and requires very small time increment.
In each computation step, the strain increment and hardening parameters are known and the plastic strain and stress increments are founds well as the back stress components   , are determined.The main steps of the explicit integration procedure of CIKH model is as follows: For a given strain increment the yield condition (see Eq. 4) is checked.
In line with Eq. 4, the behaviour of a material is checked.If  < 0, the material response is only elastic.In this situation,  = 0 and  = 0.The calculations are finished here and the stress increment is calculated in line with Eq. 18.
=   (18.) For  ≥ 0, the material yields and the plastic multiplier  is computed on the basis of the consistency condition, including both isotropic and kinematic hardening (Eq.19): ) For computed plastic multiplier, the plastic strain increment at time  is found (Eq.20).
=  (20.) For time , the backstress and isotropic hardening increments are calculated on the basis of Eq. 13 and Eq. 21.The main advantage of the first-order forward Euler explicit integration method presented here is its simplicity.However, the procedure is conditionally stable.It means that the load increment should be very small which results in a large number of integration steps, therefore.It should be also mentioned that the plastic multiplier  calculated on the basis of the consistency condition causes the yield condition is satisfied only at time .But it might not be also satisfied at time  + ∆.The solution  +∆ after many time steps might drift away from the yield surface, therefore (Fig. 1).Fig. 1 The interpretation of the solution for the explicit integration method

Implicit integration
The explicit integration is not computationally effective for higher number of loops; thus the implicit predictor-corrector procedure is proposed, therefore.The radial-return mapping algorithm (or predictor-corrector) for the CIKH model is described and tested in this paper.Firstly, the 1-D case is presented for simplicity.The steps for the predictor-corrector method of the integration of CIKH constitutive equations are as follows: For the given strain increment  the trial stress   is calculated using the Hooke's law(Eq.24).
=  +  (24.) It should be mentioned that the indices for time  + ∆ are omitted i.e.  is the stress for the time  + ∆.
If the trial stress is beyond the yield surface (Fig. 2), the plastic corrector is needed in order to translate the solution to the yield surface (Eq.25).
The  isotropic component and  backstress one is calculated on the basis of Eq. 28 and Eq.29.
The rest of the procedure is the same as in the explicit onethe plastic strain, isotropic hardening, back stress increments are calculated and the stress increment is computed.Finally, stress, plastic strain, isotropic hardening variable and the back stress are updated for the  step time.
For the 3-D case, the equations above can be rewritten as follows: The trial stress (Eq.30).
The corrected stress value (Eq.31) where part 2∆  is a plastic corrector.
The yield function might be written for 3-D case as follows (Eq.32).

𝑓 = |𝝈 𝒆
where    is the elastic trial stress (Eq.33): The stress increment is computed in line with Eq. 34: where ∆  is calculated as follows (Eq.35): and ∆  is described in line with Eq. 36: where:   is the effective stress.
The increment of the change of the effective plastic strain in  iteration can be solved with the use of the Newton's iterative solution method, giving (Eq.37-38): ) The main advantage of the implicit procedure is that it is unconditionally stable and less time/load increments are necessary in comparison to the explicit procedure.It leads to much more rapid solutions, therefore.The method is more recommended for cyclic loading curves with higher number of hysteresis loops.

Semi-implicit integration
In the implicit method, both    elastic trial stress and  isotropic hardening need to be updated at every iteration in solving (Eq.37, 38).This might be time-consuming task, especially for complex plasticity hardening laws.In the semi-implicit integration algorithm, the effective plastic strain is solved implicitly but the  normal, as well as, the isotropic and kinematic hardening is calculated using the explicit integration.In the semi-implicit scheme, Newton's method is used to determine the effective plastic strain in line with Eq. 37 and the yield function (see Eq. 32) written at the end of the time increment ensures that drift from the yield surface does not occur.The updates of all other quantities are carried out explicitly as follows (Eq.39): where: ∆ and ∆ are computed using Eq.21 and 13.
The semi-implicit scheme is not unconditionally stable and it is necessary to be concerned about its stability and accuracy when the procedure is used.But this approach can shorten the computational time in comparison to the fully explicit integration procedure.

NUMERICAL SIMULATIONS OF THE UNIAXIAL CY-CLIC LOADING TESTS
Both explicit and implicit integration procedures are implemented and tested for 1-D case in a cyclic loading tests.The numerical stress-strain curves for a strain-controlled symmetrical cyclic tension-compression test of PA6 aluminium for ∆ = 5% obtained both with the use of the explicit and implicit procedures are compared in Fig. 3.The comparison of explicit and implicit integration methods in simulations of a non-symmetrical stresscontrolled cyclic loading test is shown in Fig. 4. The following parameters and material constants were used in the simulation:  1 =8000 MPa,  2 =4000 MPa,  3 =1000 MPa,  1 =66,  2 =34,  3 =12, =154 MPa, =7, =0.62•105MPa, =0.33 and the initial yield stress   =374 MPa.
The results obtained for the explicit integration is presented as a solid line consisting of very dense dots representing integration points.The solution for the implicit integration is shown as discrete red points.It is observed that both procedures provide very similar results.However, the implicit procedure requires less number of integration steps.Constitutive equations for the CIKH model presented in Section 2 are also used in developing user material procedure for the 3-D case which is compiled and linked with the commercial Simulia Abaqus finite element method software.Both explicit and implicit procedures for the Chaboche material model have been written.In the user material procedure for a given strain increments, the plastic strain increments and stress increments should be computed and some user data e.g.backstress components, equivalent plastic strain should be also updated and saved in each computation step.It is worth noting that, in Abaqus software the user material procedure is defined in the corotational framework which rotates with the material.Procedure written for a small strains can be also applied for a large strains analysis, because of some additional computations made by Abaqus for calculated stress increments in order to make them objective [28].A very good agreement in stress and strain distributions is obtained which proves that the user material procedure is written correctly.The advantage of using the user material procedure is the easy access to all material data and user parameters in each ) stress (Eq.40-41) [28].where: ̇ is stress rate tensor in a corotational frame,  is Cauchy stress tensor,  is spin tensor which comprises both the deformation and the rotation and  is the angular velocity tensor resulting from a rigid body rotation.The parts - and  are associated with the rigid body rotation and the ̇ term is associated with the material deformation.
The Jaumann stress rate is used for materials commercially implemented in ABAQUS and for a user material subroutine.However, the application of the Green-Naghdi stress rate gives better the material response for a large shearing associated with a large rotations.It is possible to enforce Green-Naghdi objective stress rate by adding two terms shown in Eq. 42.Another important issue is the need to compute continuum Jacobian  = / (explicit approach) or consistent Jacobian  = ∆/∆ (implicit integration) in order to achieve quadratic convergence.This is the most sophisticated task, thus often for simplicity linearization is applied.This subject is out of scope of this paperauthors have implemented the user material procedure for dynamic explicit approach which not require material Jacobian.

CONCLUSIONS
The main conclusions from the work are as follows: 1.The implicit integration procedure is faster and requires less time steps in calculations.2. The explicit integration algorithm is conditionally stable and large number of time steps is required.It is not recommended for complex tasks, therefore.the results obtained with the use of both integration methods are: The implicit and explicit integration procedure for the Chaboche model were implemented in the ABAQUS software as the user material procedures.The correctness of results obtained can be compared with finding obtained for the material model implemented commercially.

2 3𝑑𝜺
:   (7.) number of   backstresses in the model should be small, but sufficient to map the experiment as accurately as possible.According to Chaboche [26], three back stress components are recommended in numerical simulations of cyclic loading tests.The first  1 component should be higly nonlinear and it starts hardening with a very large modulus indicated the quick stabilization.The second backstress  2 should model the transient nonlinear part of a stable hysteresis curve.It should be also nonlinear.The third  3 backstress component might be linear ( 3 = 0) or almost linear ( 3 → 0) in order to represent the subsequent linear part of a hysteresis curve at large strains.It is assumed that  1 >  2 >  3 .The selection of number linear and non-linear back stress components is arbitrary.The isotropic hardening of the material which results the extension of the yield surface is described using the Voce formula (Eq.16): DOI: 10.36547/ams.29.4.1949   =  0 + () (16.) = ( − ) (21.)The stress increment  at time  is computed by the following formula (Eq.22). = ( −   ) (22.) where:   is normal vector at time .Finally, stress, plastic strain, isotropic hardening variable and the back stress are updated for the  +1 step time (Eq.23).{  +∆ =   +    +∆ =    +    +∆ =   +   +∆ =   +  (23.)

Fig. 3 Fig. 4
Fig.3Comparison of explicit and implicit integration procedures for symmetrical cyclic strain-controlled tension-compression tests (stress, engineering strain)

Fig. 5 6 .Fig. 6
Fig. 5 Comparison of von Mises stressuser material procedure (a), Abaqus (b) The exemplary solutions for the plate with holes subjected to the tension obtained for the CIKH model implemented in the user material procedure and for Chaboche model implemented commercially in Simulia Abaqus software are presented in Figs. 5 and 6.
− ) + ( − ) (42.) step as well as possibility of making own modifications and extensions.It is worth mentioned that the user material subroutine UMAT applied here should consider the rotation of the coordinate system and the constitutive equations are written in the corotational frame, therefore.As a result of this objective stress rates should DOI: 10.36547/ams.29.4.1949 computation