Finite Element Implementation of a Temperature-Dependent Cyclic Plastic Model for SA508-3 Steel

: A new temperature-dependent cyclic plastic model, combining the nonlinear cyclic softening and kinematic hardening rules is established for a nuclear material of SA508-3 steel. A modiﬁed isotropic hardening rule is proposed to capture the temperature-dependent cyclic softening, and a modiﬁed kinematic hardening rule is established to improve the prediction of the ratcheting behavior by introducing an exponential function related to the accumulated plastic strain. The stress is updated by the radial return mapping algorithm based on the backward Euler integration. A new expression of consistent tangent modulus for the equilibrium iteration is derived, and then the proposed model is implemented into the ﬁnite element software ABAQUS by using the user material subroutine (UMAT) to simulate the temperature-dependent ratcheting behaviors of SA508-3 steel. Finally, the ratcheting evolutions of notched bars at elevated temperature are obtained by uniaxial stress-controlled cyclic tests, and the nonuniform strain ﬁelds on the surface of plates with a center hole is measured by using the digital image correlation (DIC) technology. Comparisons between experimental and simulated results of a material point and structural examples show that the implemented model can provide reasonable predictions for ratcheting behaviors and nonuniform strain ﬁelds of structures at different temperatures for SA508-3 steel.


Introduction
Ratcheting is a progressively accumulated inelastic strain under cyclic stressing with a non-zero mean stress. In the process of life prediction and safety assessment, both ASME and RCC-M codes specify that the ratcheting deformation occurring in structural components of nuclear reactor must be addressed. Numerous experiments were carried out to reveal the evolutions of ratcheting behaviors and many nonlinear kinematic hardening models were established to predict ratcheting deformation at different loading conditions [1][2][3][4]. Commercial finite element software, such as ANSYS and ABAQUS, are widely applied in engineering design and analysis. The nonlinear kinematic hardening rule, such as the Chaboche's rule [5], is available in these commercial software tools, which can capture the cyclic plastic deformation of some metal materials. For example, a Chaboche-type material model was applied in the Low-Cycle Fatigue finite element assessment on 316L pipes for the "Wendelstein 7-X" stellarator-type nuclear fusion experiment [6], but it commonly gives a conservative prediction [7].
To reasonably describe ratcheting behaviors of different materials, suitable modifications in the kinematic hardening rule are necessary. For example, the material parameter m i in the Ohno-Wang II kinematic hardening rule [8,9] was modified as a function of the maximum equivalent stress ε can be written as the sum of three parts, i.e., the elastic strain rate . ε e , plastic strain rate . ε p and thermal strain rate According to the Hooke's law, the elastic constitutive relation can be expressed as following: . ε e = D(T) −1 : . σ (2) where D(T) denotes the temperature-dependent fourth-order elastic tensor, .
σ is the stress rate tensor. The plastic strain rate is formulated based on the flow rule as: where .
p is the accumulated plastic strain rate, s and α represent the deviatoric stress and backstress tensors, respectively. denotes the norm. The thermal strain rate is defined as: where α is the thermal expansion coefficient and .
T is the rate of temperature. 1 is the second-order unit tensor.
The von Mises yield surface can be written in the general form: where Q represents the isotropic deformation resistance and reflects the change in the radius of yield surface. The nonlinear kinematic hardening rule revealing the evolution of backstress makes it possible to describe the progressive ratcheting deformation. Based on the Armstrong-Frederick's rule [21], Chaboche [5] firstly decomposed the backstress of into several components, which can be expressed as: .
where C i and r i are material constants and M is the number of decomposed backstress components. If the number of components increases, a higher prediction accuracy in the accumulation of plastic strain will be obtained, that is, the curve of stress versus plastic strain will be correspondingly divided into more parts (as shown in Figure 1), but it would cost greater computational consumption. It is noted that the plastic hardening moduli C i can be calculated by differentiating the curve of stress versus accumulated plastic strain with fixed r i . Figure 2 shows an exponential function of the plastic hardening modulus C with respect to the accumulated plastic strain at 250 • C and 300 • C, respectively.
where i C and i r are material constants and M is the number of decomposed backstress components. If the number of components increases, a higher prediction accuracy in the accumulation of plastic strain will be obtained, that is, the curve of stress versus plastic strain will be correspondingly divided into more parts (as shown in Figure 1), but it would cost greater computational consumption.  It is noted that the plastic hardening moduli i C can be calculated by differentiating the curve of stress versus accumulated plastic strain with fixed i r . Figure 2 shows an exponential function of the plastic hardening modulus C with respect to the accumulated plastic strain at 250 °C and 300 °C, respectively. Therefore, an exponential relationship can be used to describe the dependence of the plastic hardening moduli i C on the accumulated plastic strain and temperature as: where a , b and c are material constants and can be determined by fitting the corresponding experimental results. Therefore, the evolution of backstress in the kinematic hardening rule can be modified as: where ( ) , C p T represents the function of the accumulated plastic strain p and the temperature T , and ( ) r T is only relevant to T .
Considering the temperature dependence of the isotropic hardening, the temperaturedependent isotropic deformation resistance Q is expressed as:

Finite Element Implementation
The cyclic plastic model proposed in Section 2 is implemented into ABAQUS by using the UMAT. The implicit stress integration algorithm and a new expression of the consistent tangent modulus are derived based on the radial return mapping and backward Euler integration methods, respectively. Therefore, an exponential relationship can be used to describe the dependence of the plastic hardening moduli C i on the accumulated plastic strain and temperature as: where a, b and c are material constants and can be determined by fitting the corresponding experimental results. Therefore, the evolution of backstress in the kinematic hardening rule can be modified as: .
where C(p, T) represents the function of the accumulated plastic strain p and the temperature T, and r(T) is only relevant to T.
Considering the temperature dependence of the isotropic hardening, the temperature-dependent isotropic deformation resistance Q is expressed as: where γ(T) is a material parameter controlling the rate of resistance . Q(p, T); Q sa (T) represents the saturated isotropic deformation resistance only related to temperature. It should be noted that, when the accumulated plastic strain equals to zero, the initial deformation resistance Q 0 can be determined from Equation (10). For cyclic softening or cyclic hardening materials, it is observed that the initial deformation resistance Q 0 in tensile tests is apparently different from that in the subsequent cycles, so the initial value of Q 0 should be changed at the end of tensile curve in the first cycle.

Finite Element Implementation
The cyclic plastic model proposed in Section 2 is implemented into ABAQUS by using the UMAT. The implicit stress integration algorithm and a new expression of the consistent tangent modulus are derived based on the radial return mapping and backward Euler integration methods, respectively.

Discretization of Constitutive Equations
Firstly, it is necessary to discretize constitutive equations by the backward Euler integration method. It is noted that ratcheting behaviors of SA508-3 steel are investigated at different temperatures in this work; however, each temperature is fixed in the process of cyclic loading, and thus the temperature gradient .
T and thermal strain rate . ε T can be neglected in following derivations.
Considering the interval from a state n to n + 1, if the terms with n subindexes are known, then the terms with subindexes n + 1 can be calculated with the following expressions: The discrete evolution rule of kinematic hardening is given by The isotropic hardening rule can be discretized as:

Implicit Stress Integration Algorithm
It is assumed that the variables at time t n are known, then once the values of ∆ε n+1 and ∆t n+1 are given, σ n+1 and the variables at time t n+1 can be calculated subsequently.
Assuming the strain increment as elasticity, the trail stress can be obtained as For the trial stress state, the yield function is given by where s * n+1 = σ * n+1 − 1 3 tr σ * n+1 I. If F * y(n+1) ≤ 0, then the trial stress state is accepted as the true stress state. If F * y(n+1) > 0, then plastic yielding occurs, and then the stress state should be corrected as: where D(T) : ∆ε p n+1 is the plastic corrector. It is seen from Equation (21) that σ n+1 can be obtained once ∆ε p n+1 is solved. A non-linear scalar equation can be derived to solve this increment. The deviatoric part of Equation (21) can be given as: where G(T) denotes the shear modulus. From Equation (17), α n+1 can be rewritten as where θ n+1 is defined by the following form: Eliminating ∆ε p n+1 in Equations (22) and (23) with (14) and (15) yields Substituting Equation (25) into the yield function Equation (20) gives It is noted that Q n+1 (p, T) and C n+1 (p, T) are functions related to ∆p n+1 and T, the above non-linear scalar equation can be solved by the successive substitution method. The variables at time t n+1 will be updated after ∆p n+1 is obtained. The convergence of the successive iteration process had been proved by Kobayashi and Ohno [22]. The scheme of the implicit stress integration process is illustrated in Figure 3.
Substituting Equation (25) into the yield function Equation (20) gives It is noted that 1 ( , ) is obtained. The convergence of the successive iteration process had been proved by Kobayashi and Ohno [22]. The scheme of the implicit stress integration process is illustrated in Figure 3.

Consistent Tangent Modulus
In the process of finite element implementation, the consistent tangent modulus is necessary to ensure the global equilibrium iteration. In this subsection, a new expression of consistent tangent Differentiating Equations (13)- (15) and (18) gives

Consistent Tangent Modulus
In the process of finite element implementation, the consistent tangent modulus is necessary to ensure the global equilibrium iteration. In this subsection, a new expression of consistent tangent modulus d∆σ n+1 d∆ε n+1 is derived. Differentiating Equations (13)-(15) and (18) gives It is noted that D(T) : d∆ε p = 2G(T)d∆ε p , and taking the deviatoric part of Equation (27) gives where is the deviatoric operation of a tensor, ⊗ stands for operation of tensor product and I is the fourth order unit tensors.
Combination of Equations (28) and (30) gives Substituting Equation (29) into Equation (28) gives According to literature [21], it is assumed that the increment of deviatoric back stress has the following differential form: where H n+1 is fourth-order hardening modulus of back stress. Substituting Equations (35) and (32) into Equation (34) to eliminate d∆α n+1 and d∆s n+1 yields where L n+1 is defined as Then elimination of d∆ε Therefore, the expression of consistent tangent modulus d∆σ n+1 d∆ε n+1 is derived as The main differences between the above derived consistent tangent modulus and that done by Kobayashi and Ohno [22] are reflected in two aspects: the temperature effect considering in the variables and the tensor H n+1 (shown in Equations (35) and (37)) resulting in different formation of L n+1 . Though applying different material parameters in the kinematic hardening rule, the general form of Equation (39) seems to be similar to the classical approach of Chaboche model, because the backstress is hided in the tensor L n+1 . In order to employ the successive substitution algorithm, it is necessary to obtain the tensor H n+1 . Detailed derivation of H n+1 are given below.

Ratcheting Behaviors in A Material Point
Stress-controlled cyclic tests are performed to reveal ratcheting behaviors of SA508-3 steel at different temperatures. The experimental temperature is set as 25 • C, 150 • C and 325 • C, the loading condition is set as 50 ± 450 MPa (mean stress = 50 MPa; stress amplitude = 550 MPa), and the stress rate is prescribed as 100 MPa/s. Corresponding material parameters are listed in Table 1.  The experimental stress-strain curves are shown in Figure 4a-c. Only the hysteresis loops in the 1st-5th, 25th and 50th cycles are extracted for more clear observation. It is shown that in the 50th cycle, the strain at 150 • C is greater than that at 25 • C when the temperature increases; but it shows a contrary tendency when the temperature is up to 325 • C. This is due to the dynamic strain aging effect, which can improve the material's ability to resist deformation at elevated temperature. Liu et al. [23] also observed similar dynamic strain aging effect of SA508-3 steel when the temperature exceeded 260 • C at a certain loading rate.

Ratcheting Behaviors of Notched Bars
Uniaxial stress cyclic experiments were performed to investigate ratcheting behaviors of notched bars. The SA508-3 steel was machined into a round notched bar specimen with a section diameter of 10 mm and a gauge length of 30 mm. The finite element model with two-dimensional (2D) axisymmetrical four-node element is established in ABAQUS, as shown in Figure 6. Simulations in a material point by using the finite element software ABAQUS are performed to verify the implemented model. The ratcheting strain is defined as the averaged value of maximum strain and minimum strain in each cycle. Take the case of 25 • C as an example, the simulations are made by Chaboche model and the proposed model respectively. The corresponding material parameters for Chaboche model at 25 • C are: C 1 = 7e3 MPa, C 2 = 2.8e3 MPa, C 3 = 800 MPa, r 1 = 100, r 2 = 25, r 3 = 30. The simulated results are compared in Figure 5a. It can be seen that the Chaboche model gives a very conservative estimation on the ratcheting strain, while the proposed model gives a relatively good prediction. There is no need to repeat predictions by Chaboche model at 150 • C and 325 • C, which are also notably over-estimated. Comparisons between experimental and simulated ratcheting strain evolution by the proposed model are demonstrated in Figure 5b. It can be seen that the simulated results match the experimental ones well. Therefore, the implemented model can provide good simulations to ratcheting behaviors of SA508-3 steel. (c)

Ratcheting Behaviors of Notched Bars
Uniaxial stress cyclic experiments were performed to investigate ratcheting behaviors of notched bars. The SA508-3 steel was machined into a round notched bar specimen with a section diameter of 10 mm and a gauge length of 30 mm. The finite element model with two-dimensional (2D) axisymmetrical four-node element is established in ABAQUS, as shown in Figure 6.

Ratcheting Behaviors of Notched Bars
Uniaxial stress cyclic experiments were performed to investigate ratcheting behaviors of notched bars. The SA508-3 steel was machined into a round notched bar specimen with a section diameter of 10 mm and a gauge length of 30 mm. The finite element model with two-dimensional (2D) axi-symmetrical four-node element is established in ABAQUS, as shown in Figure 6. The stress-controlled loading condition is prescribed as 50 ± 550 MPa, and the stress rate is prescribed as 100 MPa/s. Ratcheting strains near the notch center can be obtained from experiments, as shown in Figures 7a-9a. The notch center enters plastic zone because of the stress concentration. With the increasing of number of cycles, the plastic strain of the notch center accumulates gradually. For the notched bar, the structural strain field is not uniform because of stress concentration effect so that the simulated precision is not as good as that in a material point. The model is modified based on the Chaboche model, though still some efforts should be made to improve it, the conservative prediction has been improved remarkably. Comparisons between experimental and simulated results at 25 °C, 150 °C and 325 °C are shown in Figures 7-9 respectively. It can be seen that the simulated results can give similar evolutions of ratcheting strain to experimental ones. Though the predicted ratcheting strains in the last cycles are higher than those in the initial ones, the proposed model still can give acceptable predictions of ratcheting behaviors of notched bars. The stress-controlled loading condition is prescribed as 50 ± 550 MPa, and the stress rate is prescribed as 100 MPa/s. Ratcheting strains near the notch center can be obtained from experiments, as shown in Figures 7-9. The notch center enters plastic zone because of the stress concentration. With the increasing of number of cycles, the plastic strain of the notch center accumulates gradually. For the notched bar, the structural strain field is not uniform because of stress concentration effect so that the simulated precision is not as good as that in a material point. The model is modified based on the Chaboche model, though still some efforts should be made to improve it, the conservative prediction has been improved remarkably. Comparisons between experimental and simulated results at 25 • C, 150 • C and 325 • C are shown in Figures 7-9 respectively. It can be seen that the simulated results can give similar evolutions of ratcheting strain to experimental ones. Though the predicted ratcheting strains in the last cycles are higher than those in the initial ones, the proposed model still can give acceptable predictions of ratcheting behaviors of notched bars. on the Chaboche model, though still some efforts should be made to improve it, the conservative prediction has been improved remarkably. Comparisons between experimental and simulated results at 25 °C, 150 °C and 325 °C are shown in Figures 7-9 respectively. It can be seen that the simulated results can give similar evolutions of ratcheting strain to experimental ones. Though the predicted ratcheting strains in the last cycles are higher than those in the initial ones, the proposed model still can give acceptable predictions of ratcheting behaviors of notched bars.

Ratcheting Behaviors of Plates with A Center Hole
To accurately simulate the plastic strain accumulation, it is necessary to verify the prediction of strain field distribution. Therefore, the displacement-controlled cyclic tests are performed to obtain the strain field distribution on the surface of plates with a center hole. The thickness of the plate is 2 mm, the gauge length is 10 mm, and the diameter of center hole is 2 mm. A non-contact digital image correlation (DIC) system ARAMIS5M (from GOM mhH Ltd., Brunswick, Germany) is used to measure the full-field strain within the gauged section of specimens. It is difficult to accurately obtain the strain field distribution of plates with a center hole by our DIC device at an elevated temperature since the specimen is heated by a surrounding electric oven and the temperature is also too high to use a transparent temperature controller, which would cut off the field of vision of cameras. Thus, experiments are only carried out at room temperature. The cyclic loading is controlled by displacement with a nominal strain rate of 0.2%/s and the nominal strain amplitude is 0.6%. Material parameters used in simulations are shown in Table 1.
It is noted that the tensile strength and fatigue performance of specimens was reported to be affected by scale and manufacturing mechanisms [24]. All tested specimens in this work were machined by the WEDM (Wire Electrical-Discharge Machining, Xinhe machinery processing Ltd., Chengdu, China) technology, which is considered as a valid alternative to milling in tensile specimen manufacturing [25]. Besides, the ratcheting deformation is a type of secondary deformation superimposing on the primary deformation, it is very difficult to predict accurately. In fact, although the proposed model provided more accurate predictions of ratcheting behaviors than the Chaboche model, some differences between predictions and simulations can be observed in all simulations. That is say, comparing to the differences in predictions of ratcheting behaviors, the effects of scale and manufacturing on ratcheting behaviors are relatively small, and thus they are neglected in finite element simulations. A 2D finite element model copying the specimen size is established in ABAQUS, as shown in Figure 10.
Comparisons between experimental and simulated results at the 1st, 10th and 20th cycle are shown in Figures 11 and 12. It can be seen that the simulated results of the gauged zone have similar strain field distribution with experimental ones measured by the DIC technology. Besides, the strain value of Y direction (along the loading direction) on a specific path is extracted to compare the difference between simulated results and measured ones by the DIC technology. The path location is shown in Figure 13a. It can be seen from Figure 13b that the simulated strain along the specific path is acceptable comparing with the experimental results, which means the strain field of the plate with a center hole can be reasonably simulated by the implemented cyclic plastic model.

Ratcheting Behaviors of Plates with A Center Hole
To accurately simulate the plastic strain accumulation, it is necessary to verify the prediction of strain field distribution. Therefore, the displacement-controlled cyclic tests are performed to obtain the strain field distribution on the surface of plates with a center hole. The thickness of the plate is 2 mm, the gauge length is 10 mm, and the diameter of center hole is 2 mm. A non-contact digital image correlation (DIC) system ARAMIS5M (from GOM mhH Ltd., Brunswick, Germany) is used to measure the full-field strain within the gauged section of specimens. It is difficult to accurately obtain the strain field distribution of plates with a center hole by our DIC device at an elevated temperature since the specimen is heated by a surrounding electric oven and the temperature is also too high to use a transparent temperature controller, which would cut off the field of vision of cameras. Thus, experiments are only carried out at room temperature. The cyclic loading is controlled by displacement with a nominal strain rate of 0.2%/s and the nominal strain amplitude is 0.6%. Material parameters used in simulations are shown in Table 1.
It is noted that the tensile strength and fatigue performance of specimens was reported to be affected by scale and manufacturing mechanisms [24]. All tested specimens in this work were machined by the WEDM (Wire Electrical-Discharge Machining, Xinhe machinery processing Ltd., Chengdu, China) technology, which is considered as a valid alternative to milling in tensile specimen manufacturing [25]. Besides, the ratcheting deformation is a type of secondary deformation superimposing on the primary deformation, it is very difficult to predict accurately. In fact, although the proposed model provided more accurate predictions of ratcheting behaviors than the Chaboche model, some differences between predictions and simulations can be observed in all simulations. That is say, comparing to the differences in predictions of ratcheting behaviors, the effects of scale and manufacturing on ratcheting behaviors are relatively small, and thus they are neglected in finite element simulations. A 2D finite element model copying the specimen size is established in ABAQUS, as shown in Figure 10.  Comparisons between experimental and simulated results at the 1st, 10th and 20th cycle are shown in Figures 11 and 12. It can be seen that the simulated results of the gauged zone have similar strain field distribution with experimental ones measured by the DIC technology. Besides, the strain value of Y direction (along the loading direction) on a specific path is extracted to compare the difference between simulated results and measured ones by the DIC technology. The path location is shown in Figure 13a. It can be seen from Figure 13b that the simulated strain along the specific path is acceptable comparing with the experimental results, which means the strain field of the plate with a center hole can be reasonably simulated by the implemented cyclic plastic model.

Conclusions
A temperature-dependent cyclic plastic model for SA508-3 steel, considering the nonlinear cyclic softening and kinematic hardening rules, is implemented into a commercial finite element code, i.e., ABAQUS.


The discrete material parameters in the developed kinematic hardening rule are considered as an exponential function related to the accumulated plastic strain.  Details of implementation scheme based on the radial return mapping and backward Euler integration methods are presented.  A new expression of consistent tangent modulus for overall equilibrium iteration is derived.
Then the model is implemented into ABAQUS by using the UMAT. Finite element simulations are carried out to predict ratcheting behaviors of both material point and structures, including notched bars and plates with a center hole. Comparisons between experimental and simulated results show that:


The implemented temperature-dependent cyclic plastic model can provide reasonable predictions for ratcheting behaviors of structure components of SA508-3 steel at different temperatures.
This work is helpful to evaluate ratcheting behaviors of nuclear reactor structures at different temperatures in further work. It should be mentioned that the working temperature range of the researched material is from room temperature to 350 °C, which has been covered by the present experiments. The adaptability of the proposed model beyond the working temperature range should be validated by performing more experimental researches in further work. Distance along X direction (mm) Strain in Y direction (%) DIC Simulation 1st: 10th: 20th: Figure 13. Comparison between the simulated results and measured ones by the DIC technology: (a) evaluated path location; (b) curves of strain in Y direction versus distance along X direction at different temperatures.

Conclusions
A temperature-dependent cyclic plastic model for SA508-3 steel, considering the nonlinear cyclic softening and kinematic hardening rules, is implemented into a commercial finite element code, i.e., ABAQUS.

•
The discrete material parameters in the developed kinematic hardening rule are considered as an exponential function related to the accumulated plastic strain. • Details of implementation scheme based on the radial return mapping and backward Euler integration methods are presented. • A new expression of consistent tangent modulus for overall equilibrium iteration is derived.
Then the model is implemented into ABAQUS by using the UMAT. Finite element simulations are carried out to predict ratcheting behaviors of both material point and structures, including notched bars and plates with a center hole. Comparisons between experimental and simulated results show that:

•
The implemented temperature-dependent cyclic plastic model can provide reasonable predictions for ratcheting behaviors of structure components of SA508-3 steel at different temperatures.
This work is helpful to evaluate ratcheting behaviors of nuclear reactor structures at different temperatures in further work. It should be mentioned that the working temperature range of the researched material is from room temperature to 350 • C, which has been covered by the present experiments. The adaptability of the proposed model beyond the working temperature range should be validated by performing more experimental researches in further work.