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

Authors: Huayan Chen; Xiangguo Zeng; Yang Guo; Fang Wang
Journal: Advances in Materials Science and Engineering
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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 creep-viscoplasticity 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–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–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. 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 time-dependent 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.


## 2. Theoretical Framework

### 2.1. 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]:(1)Q˙=ε˙p,where εp is the plasticity strain.Specifically, the macroscopic strainεp is expressed by the slip amount γn in the n slip systems [1]:(2)dεp=∑nαndγn,where αn denotes the orientation of the n 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(3)τ=f∑n∫dγn.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ε-p and an internal variable Q were introduced to describe the effects of the precreep εc on the subsequent plastic deformation:(4)dε-p=μdεindQ=dε-p=μdεin,where dεin=dεc+dεp. μ refers to additional hardening, caused by creep deformation. It was provided as(5)μ=Dεca+1,where D and a are material constants.According to the well-known Orowan equation [24](6)ε˙in=bPmeanV,where b is a constant. Pmean is the average mobile dislocation density. V 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 qr was assumed to be the average displacement of the atomic crowd overcoming the energy barrier of the slip system r 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 qr would suffer a local tilt, and this quantity was assumed to be linearly proportional to the free energy gradient ∂ψ/∂qr [11]. This assumption is schematically described in Figure 1.Figure 1
Energy barrier model of atomic motion.The average velocity of the atomic crowd overcoming the energy barrier with a height ofε0r depends on the energy tilt ωr under an external stress. An approximation could be provided if the average energy barrier height and energy tilt were considered as ε-0r and ω-r:(7)dqrdt=Ke-βε-0rsinh⁡βω-r,where K and ε-0r are material constants. β=1/KBT, where T is an absolute temperature and KB is the Boltzmann constant.It was assumed that the average velocity of dislocation movement was proportional todqr/dt. Consequently, ε˙in was provided by (6) and (7):(8)ε˙in=BPmeanV=Hsinh⁡βω-r,where H is a constant. ω-r is determined by the applied stress and the internal stress and was assumed to be a power function [1]:(9)ω-r=σ-XN1n,where N1 is a material constant. σ is the applied stress. The back stress X denotes the internal resistance, derived from the dislocation tangle, cell structures, subgrain, and pileups. Therefore, the overstress was represented by (σ-X).The combination of (8) and (9) produced(10)ε˙in=Hsinh⁡σ-XNnsgnσ-X,where N=N1/β1/n 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].


### 2.2. Constitutive Equations

The constitutive equations were expressed by(11)εij=εije+εijin,(12)σij=Eijklεkl-εklin,(13)ε˙ijin=32p˙Sij-XijSij-Xij,(14)p˙=Hsinh⁡Sij-XijNη1n,where εij, εije, and εijin are the total strain, elastic strain, and inelastic strain, respectively. Eijkl is the fourth-order tensor of the Hookean elasticity, and p˙ is the effective plasticity. Sij and Xij are the deviatoric stress and back stress, respectively. Nη1 is the drag stress function. For cyclic softening materials, Nη1=N0constant. H is a material constant.


### 2.3. Evolution Law

The total back stress rateX˙ij might be divided into two parts, such as X˙(1)ij and X˙(2)ij [28]:(15)X˙ij=X˙1ij+X˙2ij.The parameterμ is related to the additional hardening, caused by creep deformation. The substitution ofμ into X˙(1)ij produced(16)X˙1ij=C123μB1ε˙ijin-X1ijε˙ijin,(17)X˙2ij=C223Bη2ε˙ijin-X2ijε˙ijin-C3X2ijH2-1X2ij,where C1, C2, C3, and H2 are material constants. Bη2 refers to the cyclic hardening function. For cyclic softening materials, Bη2=B2L+1-Lexp-Wη2, η˙2=ε˙ijin, and B1, B2,L, and W are material parameters.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,X(1) and X(2), where the cyclic hardening was caused by the hardening of X(2). The evolution equations for X(1) and X(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(18)η˙m=Γmε˙ijin,m=1,2,gm=εijin-αmij-ρm≤0,ρ˙m=ΓmRmε˙ijin,α˙mij=Γm1-Rmε˙ijin,Γm=1;gm=0,∂gm∂εijinε˙ijin>00;gm<0,∂gm∂εijinε˙ijin≤0,where gm represents the cyclic nonhardening region. α(m)ij and ρm denote the center and size of the region, respectively. Γm and Rm (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 Nη1 and Bη2 might be defined, respectively, as(19)Nη1=N11+L1η1M1N1M1,Bη2=B11+L2η2M2B1M2-1,where N1, L1, L2, M1, M2, and B1 are certain other material parameters.


## 2.1. 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]:(1)Q˙=ε˙p,where εp is the plasticity strain.Specifically, the macroscopic strainεp is expressed by the slip amount γn in the n slip systems [1]:(2)dεp=∑nαndγn,where αn denotes the orientation of the n 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(3)τ=f∑n∫dγn.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ε-p and an internal variable Q were introduced to describe the effects of the precreep εc on the subsequent plastic deformation:(4)dε-p=μdεindQ=dε-p=μdεin,where dεin=dεc+dεp. μ refers to additional hardening, caused by creep deformation. It was provided as(5)μ=Dεca+1,where D and a are material constants.According to the well-known Orowan equation [24](6)ε˙in=bPmeanV,where b is a constant. Pmean is the average mobile dislocation density. V 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 qr was assumed to be the average displacement of the atomic crowd overcoming the energy barrier of the slip system r 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 qr would suffer a local tilt, and this quantity was assumed to be linearly proportional to the free energy gradient ∂ψ/∂qr [11]. This assumption is schematically described in Figure 1.Figure 1
Energy barrier model of atomic motion.The average velocity of the atomic crowd overcoming the energy barrier with a height ofε0r depends on the energy tilt ωr under an external stress. An approximation could be provided if the average energy barrier height and energy tilt were considered as ε-0r and ω-r:(7)dqrdt=Ke-βε-0rsinh⁡βω-r,where K and ε-0r are material constants. β=1/KBT, where T is an absolute temperature and KB is the Boltzmann constant.It was assumed that the average velocity of dislocation movement was proportional todqr/dt. Consequently, ε˙in was provided by (6) and (7):(8)ε˙in=BPmeanV=Hsinh⁡βω-r,where H is a constant. ω-r is determined by the applied stress and the internal stress and was assumed to be a power function [1]:(9)ω-r=σ-XN1n,where N1 is a material constant. σ is the applied stress. The back stress X denotes the internal resistance, derived from the dislocation tangle, cell structures, subgrain, and pileups. Therefore, the overstress was represented by (σ-X).The combination of (8) and (9) produced(10)ε˙in=Hsinh⁡σ-XNnsgnσ-X,where N=N1/β1/n 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].


## 2.2. Constitutive Equations

The constitutive equations were expressed by(11)εij=εije+εijin,(12)σij=Eijklεkl-εklin,(13)ε˙ijin=32p˙Sij-XijSij-Xij,(14)p˙=Hsinh⁡Sij-XijNη1n,where εij, εije, and εijin are the total strain, elastic strain, and inelastic strain, respectively. Eijkl is the fourth-order tensor of the Hookean elasticity, and p˙ is the effective plasticity. Sij and Xij are the deviatoric stress and back stress, respectively. Nη1 is the drag stress function. For cyclic softening materials, Nη1=N0constant. H is a material constant.


## 2.3. Evolution Law

The total back stress rateX˙ij might be divided into two parts, such as X˙(1)ij and X˙(2)ij [28]:(15)X˙ij=X˙1ij+X˙2ij.The parameterμ is related to the additional hardening, caused by creep deformation. The substitution ofμ into X˙(1)ij produced(16)X˙1ij=C123μB1ε˙ijin-X1ijε˙ijin,(17)X˙2ij=C223Bη2ε˙ijin-X2ijε˙ijin-C3X2ijH2-1X2ij,where C1, C2, C3, and H2 are material constants. Bη2 refers to the cyclic hardening function. For cyclic softening materials, Bη2=B2L+1-Lexp-Wη2, η˙2=ε˙ijin, and B1, B2,L, and W are material parameters.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,X(1) and X(2), where the cyclic hardening was caused by the hardening of X(2). The evolution equations for X(1) and X(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(18)η˙m=Γmε˙ijin,m=1,2,gm=εijin-αmij-ρm≤0,ρ˙m=ΓmRmε˙ijin,α˙mij=Γm1-Rmε˙ijin,Γm=1;gm=0,∂gm∂εijinε˙ijin>00;gm<0,∂gm∂εijinε˙ijin≤0,where gm represents the cyclic nonhardening region. α(m)ij and ρm denote the center and size of the region, respectively. Γm and Rm (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 Nη1 and Bη2 might be defined, respectively, as(19)Nη1=N11+L1η1M1N1M1,Bη2=B11+L2η2M2B1M2-1,where N1, L1, L2, M1, M2, and B1 are certain other material parameters.


## 3. Numerical Algorithmic Formulation

### 3.1. Semi-Implicit Stress Integration

By the complex hardening laws and evolution equations utilization, the effective plastic strainΔp 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 1.
Input state variablesX1ijt,X2ijt, initial stress σijt, strain increment εijt+Δt, time increment Δt, and elastic (initial) stiffness Eijkl.Step 2 (inelastic iteration).
At the beginning of each increment step, the inelastic strain increment was assumed to be(20)Δεijink=δ0k=0.
Following, the effective inelastic strain increment, the back stress increment, and the applied stress increment were calculated(21)Δεijink=23ΔεijinkΔεijink,ΔX1ijk=C123μB1Δεijink-X1ijtΔεijink,ΔX2ijk=C223Bη2Δεijink-X2ijtΔεijink-C3X2ijtH2-1X2ijtΔt,Δσijk=EijklΔεij-Δεklink.
Consequently, the back, applied, and overstresses were updated:(22)X1ijk=X1ijt+ΔX1ijk,X2ijk=X2ijt+ΔX2ijk,σijk=σijt+Δσijk.
The new inelastic strain incrementΔεijin(k+1) and the new equivalent effective strain increment Δεijin(k+1) were obtained by(23)Δεijink+1=32Hsinh⁡Sijk-XijkNη1nSijk-XijkSijk-XijkΔt,Δεijink+1=23Δεijink+1Δεijink+1.
The numerical error was estimated to ensure(24)absΔεijink+1-ΔεijinkΔεijink<err.Step 3.
Otherwise, the inelastic strain increment, elastic strain increment, back stress, and applied stress were recalculated:(25)X1ijt+Δt=X1ijk,X2ijt+Δt=X2ijk,σijt+Δt=σijk,εijint+Δt=εijint+Δεijink,εijint+Δt=εijint+Δεijink+1.In Figure2, the aforementioned semi-implicit stress integration as a flowchart is presented.Figure 2
Analysis flow chart of a cyclic softening material.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 Step1 and become updated at Step 3, as presented in Figure 3.Figure 3
Analysis flow chart of a cyclic hardening material.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.


### 3.2. 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, an implicit function in relation to stress and strain increments was deduced. Indeed:(26)Δp=23ΔεijinΔεijin=23Δεij-Eijkl-1ΔσklΔεij-Eijkl-1Δσkl.Equation (12) was provided by the form of an increment(27)Δσij=EijklΔεkl-Δεklin.By (26) substitution into (27),(28)Δσij=EijklΔεkl-Δεklin=EijklΔεkl-3223Δεij-ΔεijeΔεij-ΔεijeSkl-XklSkl-Xkl.It was assumed to be(29)FΔσij,Δεij=Δσij-EijklΔεkl-3223Δεij-ΔεijeΔεij-ΔεijeSkl-XklSkl-Xkl=0.By the derivative rule for implicit functions, the consistent tangent modulus [J] could be expressed by(30)J=dΔσijdΔεkl=-∂F∂Δσij-1∂F∂Δεkl.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 [Jsym] acquisition.(31)Jsym=12J+JT.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].


## 3.1. Semi-Implicit Stress Integration

By the complex hardening laws and evolution equations utilization, the effective plastic strainΔp 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 1.
Input state variablesX1ijt,X2ijt, initial stress σijt, strain increment εijt+Δt, time increment Δt, and elastic (initial) stiffness Eijkl.Step 2 (inelastic iteration).
At the beginning of each increment step, the inelastic strain increment was assumed to be(20)Δεijink=δ0k=0.
Following, the effective inelastic strain increment, the back stress increment, and the applied stress increment were calculated(21)Δεijink=23ΔεijinkΔεijink,ΔX1ijk=C123μB1Δεijink-X1ijtΔεijink,ΔX2ijk=C223Bη2Δεijink-X2ijtΔεijink-C3X2ijtH2-1X2ijtΔt,Δσijk=EijklΔεij-Δεklink.
Consequently, the back, applied, and overstresses were updated:(22)X1ijk=X1ijt+ΔX1ijk,X2ijk=X2ijt+ΔX2ijk,σijk=σijt+Δσijk.
The new inelastic strain incrementΔεijin(k+1) and the new equivalent effective strain increment Δεijin(k+1) were obtained by(23)Δεijink+1=32Hsinh⁡Sijk-XijkNη1nSijk-XijkSijk-XijkΔt,Δεijink+1=23Δεijink+1Δεijink+1.
The numerical error was estimated to ensure(24)absΔεijink+1-ΔεijinkΔεijink<err.Step 3.
Otherwise, the inelastic strain increment, elastic strain increment, back stress, and applied stress were recalculated:(25)X1ijt+Δt=X1ijk,X2ijt+Δt=X2ijk,σijt+Δt=σijk,εijint+Δt=εijint+Δεijink,εijint+Δt=εijint+Δεijink+1.In Figure2, the aforementioned semi-implicit stress integration as a flowchart is presented.Figure 2
Analysis flow chart of a cyclic softening material.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 Step1 and become updated at Step 3, as presented in Figure 3.Figure 3
Analysis flow chart of a cyclic hardening material.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.


## 3.2. 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, an implicit function in relation to stress and strain increments was deduced. Indeed:(26)Δp=23ΔεijinΔεijin=23Δεij-Eijkl-1ΔσklΔεij-Eijkl-1Δσkl.Equation (12) was provided by the form of an increment(27)Δσij=EijklΔεkl-Δεklin.By (26) substitution into (27),(28)Δσij=EijklΔεkl-Δεklin=EijklΔεkl-3223Δεij-ΔεijeΔεij-ΔεijeSkl-XklSkl-Xkl.It was assumed to be(29)FΔσij,Δεij=Δσij-EijklΔεkl-3223Δεij-ΔεijeΔεij-ΔεijeSkl-XklSkl-Xkl=0.By the derivative rule for implicit functions, the consistent tangent modulus [J] could be expressed by(30)J=dΔσijdΔεkl=-∂F∂Δσij-1∂F∂Δεkl.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 [Jsym] acquisition.(31)Jsym=12J+JT.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].


## 4. 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:(i)
For the modified 9Cr-1Mo steel at 550°C(32)H=8.777×10-27,H2=11.8,n=0.245,L=2.45×10-1,W=9.13×10-1,N0=0.425,B1=88.7,B2=1.21×102,C1=3.23×102,C2=6.72×103,C3=3.86×10-27.(ii)
For the 316 stainless steel at 600°C(33)H=1.22×10-10,H2=5.45,n=0.426,D=27.933,a=0.528,L1=1.2×103,L2=3.43×103,C1=5.0×103,C2=6.73×102,C3=7.2×10-15,M1=0.138,M2=0.358,R1=0.011,R2=0.023,B1=63.0,N1=4.686×10-2.The uniaxial cyclic loading behavior of the materials was calculated by the UMAT of ABAQUS and a 3D eight-noded isoparametric brick element at a strain rate of 0.1%/s was employed. The surfaces ①, ②, and ③ of the testing element were applied by the normal boundary constraints and subjected to uniaxial strain cycling on the surface ④, as presented in Figure4. Due to uniaxial loading, the applied stress, back stress, overstress, and strain were denoted by σ, X, σ0, and ε without the subscript “ij” during the analysis. In addition, the experimental data for the modified 9Cr-1Mo and 316 stainless steels were provided [34].Figure 4
Testing element.In Figure5, 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.Figure 5
Stress-strain hysteresis loops at first and half-life cycles for 0.01 strain range for (a) a modified 9Cr-1Mo steel and (b) a 316 stainless steel.
(a) (b)The stable hysteresis loops for various strain amplitudes are presented in Figure6. 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.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.
(a) (b)In Figure7, 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.Figure 7
Stress relaxation curves at first and half-life cycles for (a) a modified 9Cr-1Mo steel and (b) a 316 stainless steel.
(a) (b)Figure 8
Back stress behavior simulated by the proposed constitutive model for (a) a modified 9Cr-1Mo steel and (b) a 316 stainless steel.
(a) (b)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 Figure9. 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–37].Figure 9
Number of cycle effects at various strain rates on stress amplitude for a modified 9Cr-1Mo steel.The back stress amplitude as a function of the number of cycles at the strain rate 0.001/s is presented in Figure10. From this figure, the back stress X1 became constant in a short time as the cyclic number increased. Therefore, the circulation stabilization phenomenon for X1 was significantly apparent. With an increment of cyclic numbers, the back stress X2 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 X2.Figure 10
Number of cycle effects on back stress amplitude for a modified 9Cr-1Mo steel.In Figure11, 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.Figure 11
Number of cycle effects on overstress amplitude for a modified 9Cr-1Mo steel.From Figure12, 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.Figure 12
Number of cycle effects on total stress amplitude for a 316 stainless steel.When the number of cycles increased, the back stressX1 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. X1 exhibited the circulation stabilization phenomenon. In addition, when the number of cycles increased, the cyclic hardening behavior in the back stress X2 existed, whereas the cyclic hardening velocity rapidly decreased.Figure 13
Number of cycle effects on back stress amplitude for a 316 stainless steel.As presented in Figure14, 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.Figure 14
Number of cycle effects on overstress amplitude for a 316 stainless steel.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 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.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].
(a) (b)In order for the experimental results and the calculated values to be compared, a relative error was defined as follows:(34)Er=σsim-σexpσexp×100%,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.Table 1
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 0.5%.
Effective strain (%) σ e x p (MPa) σ s i m (MPa) E r (%) 3 271 272 0.37 4 290 283 2.41 5 309 291 5.83Table 2
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 0.75%.
Effective strain (%) σ e x p (MPa) σ s i m (MPa) E r (%) 3 277 283 2.17 4 299 297 0.37 5 314 309 1.59Table 3
Comparison between experimental and simulated stresses at various effective strains with a prior creep strain of 1.0%.
Effective strain (%) σ e x p (MPa) σ s i m (MPa) E r (%) 3 288 314 9.03 4 311 328 5.47 5 328 337 2.74From these tables, the difference errors indicate that the proposed model appeared to be appropriate for the description of the combined creep-plastic deformations, because this unified theory could yield a reasonable agreement with the experimental data.


## 5. 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.



---

*Source: 2129478-2017-02-08.xml*