A Unified Constitutive Model for Creep and Cyclic Viscoplasticity Behavior Simulation of Steels Based on the Absolute Reaction Rate Theory

In this work, the viscoplasticity and creep behavior for modified 9Cr-1Mo and 316 stainless steels were investigated. Based on the absolute reaction rate theory, a unified constitutive model incorporating internal state variables was proposed to characterize the evolution of the back stress. Also, the model was implemented by the ABAQUS system with the semi-implicit stress integration. Compared to the experimental data, the results demonstrated that the proposed approach could effectively simulate the cyclic softening and hardening behavior for such structural steels.


Introduction
Nowadays, many steels were considered for utilization in high temperature structures such as ultrasonic aircrafts and turbine engines.These steels are usually subjected to cyclic loadings during long-term services.Consequently, the operational cycles cause serious thermomechanical fatigues and creep-fatigue failures.An important source of inelasticity in the material structure is creep and viscoplasticity [1,2], which might be produced by microstructural changes [3,4].Especially within such practical engineering materials, more than one damage mechanism might exist, such as the creepviscoplasticity interaction.Therefore, the effects are proven essential in material mechanical performance predictions.On the other hand, it is commonly known that damage accumulation and failure are combined and complicated processes involving the interaction among multiple damage events.In order for structural reliabilities to be ensured, a unified theory that could quantify the coupling effects between creep and viscoplasticity was important to be developed, furthermore identifying the underlying failure mechanisms [5,6].Indeed, Professor Leen and his coworkers were successful in developing the unified constitutive models for utilized advanced steels in high temperature applications [7][8][9].The results demonstrated that the proposed models could reproduce the material behavior accurately in cyclic plasticity, relaxation, and creep.Also, a dislocation-based model for high temperature cyclic viscoplasticity of the 9-12Cr steel was proposed [10].Valanis and Lalwani [11] proposed a modified version for "The Thermodynamics of Internal Variables," considered "The Theory of Absolute Reaction Rates (TARR)."One internal variable was applied to characterize the response of a perfect crystal to external force fields.Although it was proposed that this model could be extended to metals and other materials with a rate dependent behavior, no specific evolution formula for the internal variables or the back stress for these materials was given.Furthermore, an extensive study mainly concentrated on cyclic inelastic behavior simulation based on the experimentally determined back stress and drag stress was conducted by Taguchi and Uno [12].
In recent years, a significant amount of work has been devoted for the effects of prior deformation (prestrain) investigation on the subsequent high temperature creep of metallic materials [13][14][15].As an example, Wilshire and Willis [16] discovered that the prestrain could increase the creep resistance for the 316H/L stainless steel at room temperature.

Advances in Materials Science and Engineering
A negative effect for the nickel-based C263 superalloy was reported by Zhang and Knowles [17], and the results indicated that the prestrain would reduce the creep resistance.Besides, Ohashi and his coworkers [18] performed an investigation of the effects of prior creep on the subsequent plastic deformation behavior for the 316 stainless steel.It was discovered that the creep resistance was enhanced by the preplastic strain by tensile loading, whereas the plastic resistance was also strengthened when prior creep existed.In the conventional inelastic constitutive model, the total inelastic strain was decomposed into the time-independent plasticity and timedependent creep.Besides, it was evidenced that this method was incapable of an accurate description of the coupling effects between viscoplasticity and creep [19,20].
Therefore, in this work a unified constitutive model based on the absolute reaction rate theory was attempted to be developed.Two internal variables related to the rate dependency were embedded into the TARR, whereas the evolution equations were provided in detail.The equivalent plastic strain was considered to account for the interactive mechanism of viscoplasticity and creep.In addition, a numerical simulation using the UMAT of ABAQUS was executed for the validity of the proposed model to be demonstrated.The obtained results would prove to be useful for the analysis of the inelastic behavior.

Physical Mechanism.
The physical mechanisms behind plastic deformation under cyclic loading are derived from the dislocation evolution and corresponding interaction [21].The evolution equation for the accumulated effective plastic strain follows [20]: where   is the plasticity strain.Specifically, the macroscopic strain   is expressed by the slip amount   in the  slip systems [1]: where   denotes the orientation of the  slip systems.The critical shear stress  denotes the resistance to dislocation slips in a material, and it is a function of the sum of the slip increments in all slip systems as For many metals, the dislocation slip and thermal activation dislocation climb are the main reasons for plasticity and creep deformation [20,22,23].Under high temperatures, plastic deformation leads to the formation of cell structures, whereas creep deformation causes the formation of subgrain structures.The characteristic sizes of dislocation cell structures are higher than the characteristic sizes of dislocation subgrain structures.The free distance of dislocation movement becomes shorter than the formation of subgrain structures.In consequence, the various deformation mechanisms would lead to a significant change in the macroresistance.When creep and plasticity interact, the dislocation movement is characterized by dislocation slips in cell structures along with the thermal activation climb in subgrain boundaries.Therefore, an inelastic deformation as the dislocation slip and thermal activation could be considered, whereas the differences in deformation conditions could be described by internal variables.
Based on the aforementioned analysis, an equivalent plastic strain   and an internal variable  were introduced to describe the effects of the precreep   on the subsequent plastic deformation: where  in =   +   . refers to additional hardening, caused by creep deformation.It was provided as where  and  are material constants.
According to the well-known Orowan equation [24] where  is a constant. mean is the average mobile dislocation density. is the average velocity of dislocation movement.Under the absolute reaction rate theory [25], Valanis and Lalwani [11] proposed that the internal variable could be expressed as the average displacement of a group of atoms subjected to a certain energy barrier.The intrinsic displacement   was assumed to be the average displacement of the atomic crowd overcoming the energy barrier of the slip system  formed by the random distribution of point defects and dislocations.When a stress field was applied, the potential energy surface in the vicinity of the group of atoms associated with the internal variable   would suffer a local tilt, and this quantity was assumed to be linearly proportional to the free energy gradient /  [11].This assumption is schematically described in Figure 1.
The average velocity of the atomic crowd overcoming the energy barrier with a height of   0 depends on the energy tilt   under an external stress.An approximation could be provided if the average energy barrier height and energy tilt were considered as   0 and   : where  and   0 are material constants. = 1/  , where  is an absolute temperature and   is the Boltzmann constant.
It was assumed that the average velocity of dislocation movement was proportional to   /.Consequently, ε in was provided by ( 6) and (7): where  is a constant.  is determined by the applied stress and the internal stress and was assumed to be a power function [1]: where  1 is a material constant. is the applied stress.The back stress  denotes the internal resistance, derived from the dislocation tangle, cell structures, subgrain, and pileups.Therefore, the overstress was represented by ( − ).
The combination of ( 8) and ( 9) produced where  =  1 /( 1/ ) is a material constant.The hyperbolic sine function was provided to describe the relationship between the strain rate and stress in (10).From a view of mathematics, when the strain rate was high, the function retreated to the exponential form, often adopted to describe the constitutive under a high strain rate [26].In contrast, at a high temperature and low strain rate, the hyperbolic sine function retreated to the power form that describes the relationship between the strain rate and the flow stress.Interestingly, this equation had a similar form to the form in the phenomenological approach proposed by Delobelle et al. [27].
It should be pointed out that although the absolute reaction rate theory was adopted from chemistry, the corresponding application to the evolution of material microstructure was widely generalized as the unity of chemistry and physics [25].

Constitutive Equations.
The constitutive equations were expressed by where   ,    , and  in  are the total strain, elastic strain, and inelastic strain, respectively.  is the fourth-order tensor of the Hookean elasticity, and ṗ is the effective plasticity.  and   are the deviatoric stress and back stress, respectively.( 1 ) is the drag stress function.For cyclic softening materials, ( 1 ) =  0 (constant). is a material constant.

Evolution Law.
The total back stress rate Ẋ might be divided into two parts, such as Ẋ(1) and Ẋ(2) [28]: The parameter  is related to the additional hardening, caused by creep deformation.The substitution of  into Ẋ(1) produced where  1 ,  Based on the overstress concept, an inelastic strain rate was assumed to be a function of both overstress and drag stress.The drag stress could be recognized as an important variable for the cycle dependence description of overstress.Indeed, the cyclic strain hardening was caused by the hardening of both drag and back stresses.The back stress was divided into two terms,  (1) and  (2) , where the cyclic hardening was caused by the hardening of  (2) .The evolution equations for  (1) and  (2) were of the Bailey-Orwan type.The last term of ( 17) was the thermal recovery term.Therefore, these rules were utilized for the cyclic hardening behavior description.
Furthermore, the cyclic nonhardening behavior could be described as where   represents the cyclic nonhardening region. () and   denote the center and size of the region, respectively.Γ  and   (m = 1,2) are material constants.It should be mentioned that the cyclic nonhardening behavior was proposed by Ohno [29].It was reported that the nonhardening region was defined in the strain space.Specifically, when the strain state fell into the nonhardening region, no hardening or softening behavior existed.Otherwise, hardening or softening occurred.Consequently, the functions ( 1 ) and ( 2 ) might be defined, respectively, as where  1 ,  1 ,  2 ,  1 ,  2 , and  1 are certain other material parameters.

Numerical Algorithmic Formulation
3.1.Semi-Implicit Stress Integration.By the complex hardening laws and evolution equations utilization, the effective plastic strain Δ was always obtained by a fully implicit integration [30].In this study, a practical approach called semi-implicit stress integration was adopted for the inelastic behavior simulation.This approach integrated the effective plastic strain implicitly and updated the other internal variables explicitly for each time increment.
Regarding the cyclic softening material, the flow of the semi-implicit stress integration in the user material subroutine of ABAQUS (UMAT) is presented as follows.
Step 2 (inelastic iteration).At the beginning of each increment step, the inelastic strain increment was assumed to be Following, the effective inelastic strain increment, the back stress increment, and the applied stress increment were calculated Consequently, the back, applied, and overstresses were updated: The new inelastic strain increment Δ in (+1)

𝑖𝑗
and the new equivalent effective strain increment ‖Δ in (+1)

𝑖𝑗
‖ were obtained by The numerical error was estimated to ensure Step 3. Otherwise, the inelastic strain increment, elastic strain increment, back stress, and applied stress were recalculated: In Figure 2, the aforementioned semi-implicit stress integration as a flowchart is presented.
Regarding the cyclic hardening material, the integration flowchart was similar to the cyclic softening material.Besides, the additional state variables for the cyclic nonhardening region were required to be recovered at Step 1 and become updated at Step 3, as presented in Figure 3.
It should be mentioned that if the strain range was considered, a mixed cyclic hardening-softening behavior might be predicted for structural steels by the proposed model.

Consistent Tangent Modulus.
The tangent modulus or the material Jacobian matrix was essential to be provided when the finite element analysis was implemented [31].In this work, Recover state variables and stress, input strain and time increment, set up elastic (initial) stiffness matrix an implicit function in relation to stress and strain increments was deduced.Indeed: Equation ( 12) was provided by the form of an increment By (26) substitution into (27), It was assumed to be By the derivative rule for implicit functions, the consistent tangent modulus [J] could be expressed by It was noted that the material Jacobian matrix obtained was not generally symmetric, and therefore a longer solution time would be required.For computing time minimization, (30) was applied for the symmetric matrix [ sym ] acquisition.
Indeed, many numerical examples demonstrated that the accuracy of the consistent tangent modulus merely affected the convergence speed of the equilibrium iteration, whereas the results from the solution were not affected [32,33].

Results and Discussion
In this section, certain numerical simulations were executed for the prediction capability of the present model to be demonstrated.The proposed model was coded as a UMAT.Unless otherwise specified, the numerical results presented in this paper were produced for the following set of material properties:

Advances in Materials Science and Engineering
Recover state variables and stress, input strain and time increment, set up elastic (initial) stiffness matrix Inelastic strain iteration step as dotted box in figure Update state variables and stress ( The uniaxial cyclic loading behavior of the materials was calculated by the UMAT of ABAQUS and a 3D eightnoded isoparametric brick element at a strain rate of 0.1%/s was employed.The surfaces A, B, and C of the testing element were applied by the normal boundary constraints and subjected to uniaxial strain cycling on the surface D, as presented in Figure 4. Due to uniaxial loading, the applied stress, back stress, overstress, and strain were denoted by , ,  0 , and  without the subscript "" during the analysis.
In addition, the experimental data for the modified 9Cr-1Mo and 316 stainless steels were provided [34].
In Figure 5, the stress-strain hysteresis loops for the modified 9Cr-1Mo and the 316 stainless steels under strain symmetric cyclic loadings are presented.The comparison of the experimental results with the results derived from modeling demonstrated the predictability of the present theory.
The stable hysteresis loops for various strain amplitudes are presented in Figure 6.It was suggested that the cyclic behavior was dependent on the strain range: the response of the stress amplitudes became significantly exaggerated as the applied strain range increased for both the cyclic softening and cyclic hardening materials.
In Figure 7, the stress relaxation curves during the holding time = 3600 s at first and the half-life cycles are presented.It can be observed that the proposed constitutive equations could describe the cyclic stress relaxation behavior for the cyclic strain softening or hardening.Also, Figure 8 is a display of the back stress evolution with the increment of the inelastic strain.
The relationship between the responded stress amplitudes and cyclic number under the strain amplitude of 0.5% for the modified 9Cr-1Mo steel is presented in Figure 9.The strain rates were 1 × 10 −5 , 1 × 10 −4 , 1 × 10 −3 , and 1 × 10 −2 s −1 , respectively.From this figure, the responded stress amplitudes were discovered to decrease along with the increasing number of cycles with the same strain rate and strain amplitude.In addition, the cyclic softening velocity reduced gradually following each cycle.Subsequently to several cyclic loadings, the stress amplitude tended to be consistent.On the other hand, it was apparent that the stress amplitude decreased along with the strain rate decrease following the same number of cycles.This phenomenon could be attributed to the decrement of both yield and flow stresses, due to a reduction in the strain rate [35][36][37].
The back stress amplitude as a function of the number of cycles at the strain rate 0.001/s is presented in Figure 10.From this figure, the back stress  1 became constant in a short time as the cyclic number increased.Therefore, the circulation stabilization phenomenon for  1 was significantly apparent.With an increment of cyclic numbers, the back stress  2 exhibited the cyclic softening behavior, whereas the cyclic softening velocity gradually decreased.It should be noted that the cyclic softening phenomenon was caused by the back stress  2 .
In Figure 11, the relationship of the overstress amplitude and the number of cycles at the strain rate 0.001/s is presented.It was observed that when the cyclic number increased, the amplitude of overstress remained unchanged under strain symmetric cyclic loadings.The results from Figures 9-11 indicated that the cyclic strain softening behavior was caused mainly by the back stress softening behavior.
From Figure 12, it is presented that the responded stress amplitudes increased on the strain rate of 0.001/s as the number of cycles increased.It was remarkable that the stress amplitude could reach a saturated state following the beginning of the initial cyclic loadings, and the responded stress remains unchanged [38].The hardening velocity was gradually reduced to zero.When the number of cycles increased, the back stress  1 for the 316 stainless steel at the strain rate 0.001/s became constant under strain symmetric cyclic loadings, just as presented in Figure 13. 1 exhibited the circulation stabilization phenomenon.In addition, when the number of cycles increased, the cyclic hardening behavior in the back stress  2 existed, whereas the cyclic hardening velocity rapidly decreased.As presented in Figure 14, when the number of cycles increased, the amplitude of overstress remained unchanged under strain symmetric cyclic loadings when the strain rate was 0.001/s.From Figures 12-14, it is presented that the cyclic strain hardening behavior derived mainly from the hardening behavior of the back stress.
Peng and Zeng [39] proposed the relevant constitutive equations to describe the precreep-plasticity behavior according to the endochronic theory.In the current study, an additional hardening effect caused by subgrain structures was utilized to describe the nonequivalent hardening.From Advances in Materials Science and Engineering Figure 15, the prior creep contributed to the subsequent plastic resistance in a different way from the same amount of plasticity.Due to prior creep, the flow stress appeared to be higher than the stress by the pure plastic deformation under uniaxial loading.In addition, under the same total strain, a higher prior creep strain might lead to a higher flow stress in the subsequent plastic deformation.This increment of flow stress due to prior creep nonlinearly varied with the amount of prior creep, which might have been expected as a result of the nonlinear trend in the microstructural behavior by creep.
In order for the experimental results and the calculated values to be compared, a relative error was defined as follows:  where  sim is the simulated stress value, and  exp is the experimental value by Ohashi et al. [18].All the results for various strains are briefly summarized in Tables 1-3.From these tables, the difference errors indicate that the proposed model be appropriate for the description of the combined creep-plastic deformations, because this  unified theory could yield a reasonable agreement with the experimental data.

Conclusions
(1) By the internal variable thermodynamics utilization, a viscoplasticity constitutive model that described the relationship between inelastic strain rate and overstress was established under the framework of the absolute reaction rate theory.
(2) Two internal variables related to the rate dependency were embedded in the proposed model, and the evolution equations for these variables were provided in order for the interactive mechanism of plasticity and creep to be taken into account.
(3) A new expression of the consistent tangent modulus was derived by the derivative rule for implicit viscoplasticity constitutive relationships.
(4) A semi-implicit stress integration algorithm corresponding to the viscoplasticity constitutive model was introduced into the UMAT of ABAQUS.
(5) By comparison of the simulation results and the experimental data, it was derived that the constitutive model proposed in this work could simulate the shape of the stress-strain hysteresis loops.

Figure 1 :
Figure 1: Energy barrier model of atomic motion.

Figure 2 :
Figure 2: Analysis flow chart of a cyclic softening material.

Figure 6 :Figure 7 :
Figure 6: Stress-strain hysteresis loops at half-life cycles for several strain ranges for (a) a modified 9Cr-1Mo steel and (b) a 316 stainless steel.

Figure 8 :Figure 9 :
Figure 8: Back stress behavior simulated by the proposed constitutive model for (a) a modified 9Cr-1Mo steel and (b) a 316 stainless steel.

Figure 10 :
Figure 10: Number of cycle effects on back stress amplitude for a modified 9Cr-1Mo steel.

Figure 11 :Figure 12 :
Figure 11: Number of cycle effects on overstress amplitude for a modified 9Cr-1Mo steel.

Figure 13 :Figure 14 :
Figure 13: Number of cycle effects on back stress amplitude for a 316 stainless steel.

Figure 15 :
Figure 15: Effects of the prior creep strain on the subsequent plastic deformation for a 316 stainless steel by (a) simulations and (b) experiments [18].

Table 1 :
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 0.5%.

Table 2 :
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 0.75%.

Table 3 :
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 1.0%.