A non-unified viscoplastic constitutive model based on irreversible thermodynamics and creep-fatigue life prediction for Type 316 stainless

In this work, to describe the cycle behavior considering fatigue-creep interaction, a non-unified viscoplastic constitutive model for 316 stainless steel is derived within the irreversible thermodynamic framework. The internal variables considering kinematic and isotropic hardening properties are selected to construct the evolution equation of visco-plastic and creep terms. The proposed constitutive model was validated by the comparison with the existing literature. It was manifested that this constitutive model could successfully predict the hardening behavior and stress relaxation process under the cyclic loading. During the dwell period, the increment of the inelastic strain is decomposed into the viscoplastic and creep term. The viscoplastic deformation dominates first stage of the stress relaxation, while the stable stage is controlled by the creep term. Finally, the predicted values of mean stress are taken into the Manson-Coffin law, the low cycle fatigue life prediction are carried out based on the numerical model, which showed robust correlation with experimental results.


Introduction
316 stainless steel is widespread applicated in pressure-bearing parts in the fields of chemical industry, nuclear power, and aerospace due to its reliable safety and high temperature stability [1][2][3][4][5]. Generally, during the operation of the equipment, the pressure and temperature fluctuation produces cyclic loading which will cause the fatigue damage to the equipment and the inherent stress relaxation during steady operation induces creep deformation [6,7]. For example, the main container of the ITER tritium storage system bed is made of 316 stainless steel [8,9]. The cycle absorption and desorption operation lead to a potential risk of structural failure because of the interaction of creep and fatigue. The interaction between fatigue and creep has become an important factor limiting the service life of such equipment. Instead of a linear summation, this interaction behavior exhibits that the creep deformation and fatigue damage promote each other, and eventually accelerate the failure of the structure Therefore, the prediction of fatigue-creep life of 316 stainless steel is generally valued by many research in recent years. To consider this minus point, numerous method and standard have been carried out to forecast the lifetime under creep-fatigue conditions like LDS, SPR, and DE [10,11]. Although these methods can predict the life of 316 stainless steel, their prediction results are relatively conservative. With the development of computer technology, artificial intelligence neural network algorithms [12,13] are applied to calculate the fatigue life of 316 stainless steel. In addition, as a traditional fatigue life prediction method, the Coffin-Manson criterion is still used in various studies proposed by Tirbonod [14] and Yuan et al [15]. For the above evaluation method, the creep-fatigue interaction deformation is always the most crucial part for life prediction.
The most common way to get this combined mechanical behavior is low-cycle fatigue (LCF) experiments with incorporation of hold time at constant strain or stress in the laboratory. Shankar et al [16] proposed an experimental study about the time dependent effects on the LCF behavior of 316 stainless steel weld joint. The different parameters influence (including strain rate, temperature, strain range, and hold time) are discussed based on their experiment results and microstructure observation. Another experimental creep-fatigue interaction research for austenitic 316 stainless steel weld joint was given by Kumar et al [17] The interaction of creep deformation with isothermal low cycle fatigue cycling was assessed by incorporating dwell time. They found that the most sensitive factor for this structure response belong to the temperature range and hold period under isothermal and thermomechanical fatigue loading. Carroll et al [18] studied the creep-fatigue deformation of an ultra-fine precipitate strengthened advanced austenitic alloy and the corresponding cyclic deformation response is directly compared to 316 stainless steel. For both of them, the similar failure value could be obtained from the experiments. A number of uniaxial creep-fatigue tests of austenitic stainless steel were conducted by Takahashi et al [19]. Based on their results, they developed a creep-fatigue method which could be successfully used in 316 stainless steel life prediction. In the past decades, many similar experiments have been carried out to analysis the creep-fatigue interaction influence of 316 stainless steel. It was found that various factors such as temperature, strain rate, hold period and strain range had significant influences on LCF resistance. However, these experiments were not suitable for the engineering applications due to the complex experimental conditions and huge expenses.
The other way to obtain the cycle deformation of material is based on the constitutive model. In order to describe the mechanical response of materials better under cyclic loading, a lot of research work are devoted to the development and application of constitutive models [20]. One is to use a yield surface function to define an inelastic region. The stress state outside the yield surface can lead to the evolution of inelastic deformation. This concept was first proposed by Perzyna [21,22], developed and applied by many researchers especially Chaboche et al [23,24]. Due to the accuracy and advancement of the Chaboche unified viscoplasticity model, it has been used to simulate many cycle experiments with different materials. Hyde et al [25] used Chaboche model to reproduce the experimental cycle behavior of 316 stainless steel under high temperature conditions. By introducing an interpolation function, the anisothermal cycle loading behavior could be simulated in the Chaboche model, the comparison results with the experimental data showed reasonably good correlation. Roy et al [26] also used the Chaboche model to simulate the cycle behavior for 316 stainless steel at room temperature. The simulation plastic strain energy dissipation was carried out as the parameters to give a successful life prediction. Chaboche model contains a lot of material constants, which makes it inconvenient to use this model. To solve this problem, Palma [27] proposed an iterative optimization method to determine the parameters in Chaboche model. Besides the austenitic stainless steels, Chaboche model was successfully used to simulate stress state in other materials like 42CrMo4 steel [28], Aluminium Alloy 7075-T6 [29], and 9Cr-1Mo steel [30]. The other method to construct the governing equations of constitutive model is based on the utilize of internal variables without considering the yield surface like Bodner-Partom model [31]. Bodner-Partom model was extended by Sands et al [32] to simulate the response of materials with extreme kinematic hardening, the results showed that the modification of this model can better describe the hardening process of the material. In the past decades, many studies attempted to establish and derive new constitutive models based on the theory of irreversible thermodynamics [33][34][35][36][37]. This method is developed from the laws of thermodynamic, which is assumed that the deformation behavior of the material is required to satisfy the basic conservation relation in thermodynamics. It provides the basis and framework for the development of constitutive relation. Based on this, the various constitutive model for different materials has been conducted in past decades like asphalt mixtures [35], NiTi shape memory alloy [38], and wrought Mg alloys [39].
However, the current constitutive research lacks a unified description of the interaction between fatigue and creep. The existing constitutive models can be divided into unified one and non-unified one according to their different descriptions of inelastic strain. For the unified constitutive model, the increment of the inelastic strain could be described by a flow equation and a set of discrete internal variables. The inelastic strain could be obtained by the unified equations. However, in the non-unified constitutive model, the inelastic strain should be decomposed into time-independent plastic strain and time-dependent creep strain. The creep and plastic behavior are solved by the different evolution equations. Compared with the unified constitutive model, the non-unified model could more accurate stress-strain response description when the high temperature creep deformation is the dominant factor [40]. In this work, to describe the mechanical behavior of fatigue-creep interaction, a non-unified creep-viscoplastic constitutive model has been proposed based on irreversible thermodynamic theory, where inelastic strain is decomposed into creep strain and visco-plastic strain. The creep evolution term has been conducted by the back stress and drag stress consistent with the irreversible thermodynamic framework. The simulation results agree well with other existing reference. Finally, the predicted value of mean stress was taken into the Coffin-Manson law to give a life prediction for 316 stainless steel.

Constitutive model 2.1. Preliminary assumptions
Some assumptions made in the present work for visco-plastic constitutive model are as follows: 1. Plastic flow of material is incompressible; 2. Stress tensor could be identified by strain tensor and the internal variable, the evolution equation of internal variables follows a set of ordinary differential equations: The study of microscopic mechanism shows that visco-plastic flow is the result of crystal slip and diffusion. Therefore, we can define a scalar generalized viscous potential function from the perspective of irreversible thermodynamics and dislocation dynamics:

Irreversible thermodynamic framework
To fulfill the conservation requirement of thermal and deformation energy. The following equation has to be satisfied: where r is specific heat supply rate, q is heat flux gradient, y is the Helmholtz specific free energy. Meanwhile, the Clausius-Duhem inequality also need to be satisfied when the inelastic deformation is considered to ensure that the entropy of the system cannot decrease: Substituting equation (4) into equation (5) to eliminate the specific heat supply rate, the following inequality can be obtained: For isothermal conditions, equation (6) can be written as Clausius-Duhem form: Assuming the Helmholtz free energy function of unit can be written as: Taking the time derivative of y in equation (8): Substituting equations (2) and (9) into equation (6): (10), the following equations can be obtained:

Due to the dissipation inequality shown in equation
Due to the independence of e  , ij e e  ij in and V , i equation (10) can be simplified as following form under isothermal and adiabatic conditions.
Where A i is the generalized force associated with the internal variable V . i Inequality (12) indicates that the free energy is dissipative when the internal variable changes. In order to satisfy the Clausius-Duhem inequality, potential function can be written as: Substituting equations (13b) and (13c) into equation (12): Convex outward and nonnegative restrictive condition need to be followed by the construction of s Q = Q( ) A , ij i and l  respectively.

Construction of helmholtz free energy function
According to equation (2), the unit Helmholtz function can be written as the sum of elastic part y , e inelastic part y , in and drag stress part y : The elastic term y e is given as: The visco-plastic energy function can be constructed by the internal variables a ij related to the back stress W ij vp and satisfy the following forms: The drag stress free energy can be constructed by the internal variable b: The equation (15) can be rewritten as following form by Equations (16) Substituting equations (16) and (15) into equation (11), the evolution equation of internal variable can be obtained as: In order to satisfy Clausius-Duhem inequality, the potential function is constructed as: where S ij is the deviator of stress, a b c , , are material constants.
Substituting equation (19) into equation (13): Three main forms are used to construct l  , including power function (Ax n ), exponential function x sinh m n ). For the most metals, power functions are applicable. The traditional yield surface equation can be expressed as follows: The l l   , 1 2 can be expressed as follows: The total governing equations are summarized as follows:  The derivation process of constitutive model is shown in figure 1. The parameters of constitutive model for the 316 stainless steel are shown in table 1.

Incremental form and process of constitutive model
Step 1. An assumed initial value is given at the first step. Step 3. Calculating the drag stress term.
Step 4. The stress is determined.
Step 5. The new back stress and stress term are obtained.
, , Step 6. The yield surface is defined based on the back stress value.  Step 7. The new equivalent visco-plastic strain increment is as follows:

Uniaxial tension
The qualitative analysis for strain rate hardening and temperature softening effect are performed in this part. Figure 3(a) presents the tensile response curve of the 316 stainless steel under different strain rates with a constant temperature value (873 K). It can be seen from the figure that the material exhibits obvious ratedependent characteristics. The rate-dependent characteristics are mainly manifested in two aspects. Firstly, the yield point of the material increases with the strain rate. When the loading rate increased from 10 −3 to 10 −1 , its yield point increased about 30 MPa. Secondly, the rate of plastic evolution will decelerate remarkably under the high strain rate value which leads to the larger slope of the stress-strain curve in the plastic phase. At 0.6% engineering strain, the stress value is about 330 MPa when the strain rate is 10 −3 , and the stress value reaches about 380 MPa when the strain rate raises to 10 −1 . The temperature softening effect exhibits the opposite trend compared with the strain rate, the increasing of temperature will lead to the decline of elastic modulus and yield point which is shown in figure 3(b). Moreover, the high temperature will promote the evolution of plasticity which will result in a significant reduction of stress amplitude under the same engineering strain. Therefore, the constitutive model simulation data could reflect the rate and temperature sensitivity of the stainless steel. Figure 4 shows uniaxial tension and unloading stress-strain curves under different strain amplitudes with strain rate:10 −3 and temperature:873 K. The predicted curve with strain amplitude of 0.3% agreed well with the experimental results proposed by Hyde et al [25], which illustrates that the model is reasonable to describe the loading and unloading process of 316 stainless steel. Meanwhile, it can be seen from the figure, as the strain amplitude increases, due to the improvement of isotropic and kinematic hardening effects, the curve in the unloading stage shows an obvious nonlinearity tendency. For example, when the strain amplitude is 0.15%, the unloading curve is almost a straight line due to the isotropic and kinematic hardening effect. However, for the case of 0.6%, the loading strain beyond the hardening range shown as the curve slope can't keep constant value during the unloading process. In addition, the final stress value after tension and unloading is also affected by loading strain amplitude. Compare the situation of 0.15% and 0.5%, the final stress value drops from −49 MPa to −198 MPa. The reason for this phenomenon can be explained that the large deformation promotes the accumulation of inelastic strain in the internal of the material, and the greater stress is required for the material to return the initial state.

Cyclic hysteresis curve
The validation for the developed model under the cycle loading is conducted by comparing the prediction value to the experimental data proposed by Hyde et al [25]. The figure 5 gives the comparison results with the experimental temperature at 873 K, (a) Strain amplitude=0.6%, (b) Strain amplitude=1.0%. It can be seen from the above two figures that the comparison with the experimental data reflect that the numerical model can predict the cyclic results of 316 stainless steel well. Although the errors of strain amplitude=0.6% is larger than the situation of strain amplitude=1.0%, the maximum value is controlled within 10%. As shown in figure 5, at the initial stage of cyclic loading, the mechanical response against the cyclic loading is an unstable process which is shown as the peak value of stress increases with the number of cycles and the corresponding hysteresis curve cannot be closed (cycle hardening effect). However,this unstable hardening effect will be steady after more than 10 cycles, and finally a closed hysteretic loop is formed when the cyclic hardening effect reaches saturation. After that, the energy dissipation of material tends to a fixed value in each cycle loading which could be used to predict low cycle fatigue life.
The influence of strain amplitudes for hysteresis loop under 20 cycles with strain ration=0 is shown in figure 6. As shown in this figure, the height of the stress-strain hysteresis loop (i.e., the stress amplitude) exhibits  a decreasing tendency with the cyclic loading proceeds, indicating that the peak stress is reduced throughout the whole process, and the material undergoes a cyclic stress relaxation process. Under the loading condition of strain control, the total strain amplitude remains constant and the cyclic stress relaxation phenomenon implies that the reduction of elastic strain and the accumulation of plastic strain evidently. Under the cyclic loading process, the material will dissipate energy repeatedly due to plastic deformation, and its macroscopic manifestation is the hysteretic effect. In general, it is considered that the dissipation energy is basically equivalent to the plastic strain energy. During the uniaxial tension cyclic loading process, the plastic strain energy can be identified based on hysteresis loop area. The calculation results showed that the dissipation energy will be enhanced by the larger deformation loading which is shown as the larger steady-state hysteresis loop area for the high strain amplitude situation. Furthermore, the above two figures exhibit different trends with the continuous cyclic loading. When the strain ratio R=0 and the strain amplitude is 0.15%, the steady-state hysteresis loop is extremely narrow, and almost all experimental points are on a straight line shown in figure 6(a), indicating that the stress and strain corresponding to the cyclic load are mainly in the elastic range, and the material has a longer fatigue life. However, when the strain amplitude is 0.3%, the steady-state hysteresis loop is relatively extended, and the hysteresis curve area is larger than the former one presented in figure 6(b), indicating that the corresponding deformation state are in the plastic range, and its service life is weakened by the influence of the cyclic loading.
The figure 7 shows the comparison of hysteresis curves under different strain loading history, that is, the material is first loaded to a strain value and then cycled. It can be seen from the figure 7, with the value of the strain loading history gradually increases, the plasticity strain in the material gradually accumulates, so that the  stress-strain gap of the hysteresis curve quickly narrows, and finally reaches a relatively stable state. When the material is loaded to a high strain value, the mechanical properties of the hardening part have been rapidly strengthened, and its plastic dissipation capacity under cyclic loading has been gradually weakened, which is shown as the strong brittleness.

Cyclic loading with holding time
In general, the cyclic stress response of the 316 stainless steel is greatly affected by the dwell time. For the sake of the validation of the developed model, the material parameters are replaced to be the same with the Shimada et al experiments [40] and compared with the cycle test under the condition of holding load. The figure 8 and figure 9 show the comparison of the hysteresis curve with holding time and the experimental comparison of the stress relaxation phenomenon caused by viscoplastic and creep under the holding load condition respectively. The prediction value by the developed model agrees well with the data proposed in Shimada et al [40] experimental results, which validates the accuracy of the developed model. During the holding strain period, the visco-plastic and creep term will continue to evolve, which leads to grown of inelastic strain proportion. This phenomenon is shown as the stress value gradually decreases with the time observed from figure 8. The stress relaxation effect depicted against the dwell time is shown in figure 9. Although the simulation results of the stress value after 300 s is higher than the experimental data, the errors of the stress are controlled within 10 MPa.
The figure 10 reflects the cycle-hold time stress response under different strain amplitudes with the physical parameters shown in table 1. The temperature and strain rate are 873 K and 10 −3 respectively. Figure 10(a) shows the stable hysteresis loop with strain amplitudes of 0.4%, 0.7%, and 1.0% under load retention. It can be found from the figure that the higher strain amplitude results in a more obvious stress relaxation effect. It is observed as the longer descent line in the hysteresis loop. Figure 10(b) present the corresponding relaxation  curve during hold time. The stress relaxation effect could be roughly divided into two stages. In the first stage, inelastic strain increases with a high rate, which leads to the stress value drop rapidly. Then, the rate of stress relaxation will gradually decrease to a constant value in the second stage. The comparison of stress relaxation value under different time stage is shown in figure 11. For the different strain amplitude cases, the obvious contrast is observed in the first 200 s. In the first 10 s, the stress has fallen by roughly 21 MPa when the strain amplitude is 1.0%, while the stress decreases by approximately 13 MPa when the strain amplitude is 0.4%. In this constitutive model, the inelastic strain is composed of visco-plastic and creep term. The figure 12 present the relationship between the inelastic strain rate and the stress. The visco-plastic rate, and creep rate are also plotted for the convenience of observation. It can be observed that visco-plastic line and creep line intersect at the rate about 1.0×10 −6 . A similar point value was obtained in the experimental investigation proposed by Shimada et al [40]. Due to the hardening effect, after 10 cycles, the material will have greater stress value under the same strain. For the first stage of the stress relaxation process, viscoplastic deformation dominates the inelastic strain increment, the reduction of the stress exhibits a higher rate. In the second stage, the stress decreases with a minor constant rate compared with the first stage. The evolution of internal variables tends to be stable at this stage, the magnitude of stress relaxation is determined by the creep deformation. In general, the stress relaxation process is roughly divided into three stages. In the last stage, due to the accumulation of creep damage, the inelastic strain will increase rapidly until reaching the point of rupture. However, it is worth noting that the damage effect is not considered in this constitutive model. Therefore, only the first two stages of stress relaxation process can be described in this work.

Mean stress
The mean stress is expressed as the average of the maximum and minimum stress during a single cycle. It is generally agreed that mean stress is the decisive factor of life prediction. The longer prediction life value corresponds to the larger mean stress due to the less generation of inelastic deformation. Therefore, we will discuss the effect of different parameters on the mean stress including strain amplitude, strain rate and temperature. For the convenience of comparison, the strain ratio is set to 0 in this section. Figure 13(a) reflects the influence of different strain amplitudes on the mean stress (temperature:873 K, strain rate:10 −3 ) with strain amplitudes of 0.15%, 0.2%, and 0.3%. The mean stress in the initial part is found to depend on the strain amplitude. It could be observed from figure 13(a) that the initial value is increased by approximately 20 MPa when strain amplitude increases from 0.15% to 0.3%. The mean stress relaxation phenomenon could also be found in the figure 13 which is shown as the reduction of the mean stress along with the cycle numbers. The stain amplitude promotes this phenomenon, it is noticed that the 0.3% case has the fastest declining rate shown in figure 13(a). However, the curve will reach a stable value after several cycles. This stable value will decrease as the improving of the strain amplitude. Figure 13(b) shows the variation of the mean stress curve of 316 stainless steel with the number of cycles under different strain rates. The strain rate hardening effect has a strong positive influence on the initial value, the initial mean stress will be improved significantly by the increasing of the strain rate. Due to the strain rate strengthening effect, the yield point of the material is increased, which slows down the relaxation effect of the mean stress. A similar tendency could be found between the figure 13(a) (the strain amplitude is 0.15% and the strain rate is 10 −3 ) case and figure 13(b) case (the strain amplitude is 0.2% and the strain rate is 10 −1 ). This could be explained that the increasing strain rate results in the higher yield point of the material, which leads to the cyclic stress response are mainly distributed in the elastic range. Meanwhile, the stable value of mean stress is also improved by the strain rate hardening effect. Figure 13(c) presents the mean stress curves at different temperatures. It can be seen from the figure that, compared with strain rate and strain amplitude, the influence of temperature on the mean stress change is relatively inapparent. When the temperature rises from 773 K to 873 K, the initial and stable values of the mean stress only decrease by about 5 MPa. Furthermore, due to the temperature softening effect will reduce the yield point of the material, the high temperature will promote the inelastic strain accumulation in the material under the same strain amplitude. This process is shown as the mean stress decreasing slightly with the higher temperature. In general, the higher stable value of mean stress means the longer fatigue life. Therefore, the improvement of strain rate has a positive impact on the service life of materials, while the improvement of temperature and strain has a negative effect.

Fatigue-creep life prediction
It is mostly recognized that the service life of a material is positively related to its stable mean stress under cyclic loading. I.e., the longer life prediction value could be found under the higher mean stress case. In this section, the mean stress value obtained from the cyclic behavior modeling based on the developed constitutive equations are The experimental service life value could be found in the reference proposed by Srinivasan et al [12]. Figure 14(a) shows the fatigue life prediction of different strain amplitudes. It can be seen that the prediction results are relatively ideal, and all 4 points fall within the range of 2. The prediction results of different hold time given in figure 14(b) are not as good as those shown in figure 14(a). However, all predicted points are within the range of the double zone, and the predicted results are acceptable.

Conclusion
To describe the 316 stainless steel mechanical behavior of fatigue-creep interaction under cycle loading, a nonunified viscoplastic constitutive model is derived consistent with the irreversible thermodynamic theory. The main results and conclusions are summarized as: 1. The developed non-unified viscoplastic constitutive model could describe cyclic stress response and strain state precisely under the cyclic loading without dwell loading. Meanwhile, due to the conducted evolution equations of viscoplasticity and creep, the good agreement of the simulation results and experimental data in other reference proved that this model could also predict the cyclic hysteresis loops with the dwell stress relaxation effect.  2. Viscoplasticity and creep dominate the two stages of stress relaxation respectively. For the first stage of the stress relaxation process, viscoplastic deformation dominates the inelastic strain, the reduction of the stress exhibits a higher rate. During the dwell period, the demarcation point is about 1.0×10 −6 for visco-plastic rate and the creep rate. After that, the stress relaxation is determined by the creep deformation shown as the stress decreases with a minor constant rate.
3. Based on the calculation results, this work discusses the relationship between the mean stress of the material under cyclic loading under different strain amplitude, strain rate and temperature. The improvement of strain rate has a positive impact on the mean stress of materials while the improvement of temperature and strain has a negative effect.
4. The mean stress calculated by the constitutive model are taken into the Coffin-Manson law, some life prediction results are given in the article. The prediction life value agrees well with the experimental data proposed in other reference.