Implementation of ABAQUS User Subroutines for Viscoplasticity of 316 Stainless Steel and Zircaloy-4

: This paper describes the formulations for the viscoplasticity of metals based on the Chaboche and Delobelle model. The implementations of the viscoplastic models were detailed herein and then implemented via user subroutines for material models (UMAT) in ABAQUS. Two typical metals, i.e., 316 Stainless Steel and Zircaloy-4, were chosen as examples and their viscoplastic behaviors were captured. Numerical simulations are compared to reported experiments in order to validate the models and the UMAT codes. The typical viscoplastic behaviors of both metals, such as stress relaxation and creep, were captured well through the available experiments. We have publicized all the data and codes.


Introduction
The stress corrosion cracking (SCC) of cladding tube structures caused by a pellet-clad interaction (PCI) seriously threatens the operation safety of nuclear facilities [1][2][3][4].In order to ensure the structural integrity of nuclear fuel cladding during service, it is necessary to accurately calculate the stress and strain levels of cladding materials in service conditions.
316 Stainless Steel and Zircaloy-4 have become widely used cladding materials due to their excellent high temperature mechanical properties, good oxidation resistance and corrosion resistance [5][6][7].Elastic deformation, plastic deformation and high temperature creep are the main deformation behaviors of cladding materials under service conditions.Current fuel rod simulation codes such as FRAPCON, FPIN2 and FALCON, typically simplify complex cladding deformation problems into one-dimensional or two-dimensional problems based on axisymmetric, plane strain and one-dimensional radial representation [8][9][10].However, this method cannot model the multiaxial mechanical behavior of the cladding.Not only that, empirical formulas are generally used in these codes to describe the plasticity and creep behavior of the cladding [11].This approach greatly limits the generalisability and accuracy of the models despite their simplicity in theory and implementation.Therefore, since 2008, Williamson et al. from the Idaho National Laboratory (INL) have been developing the multidimensional finite element fuel code BISON [12][13][14].At the same time, Williamson [15,16] carried out a multidimensional thermomechanical analysis of the multipellet steady and transient light water reactor fuel rod behavior based on the ABAQUS commercial finite element code.A power-law creep formulation for describing the deformation behavior of Zircaloy cladding is implemented via an ABAQUS user creep subroutine.However, 316 stainless steel and Zircaloy-4 have coupled plasticity and creep behavior at high temperatures [17].The classical rate-independent elastoplastic constitutive model, the Norton creep equation [18], or the widely used Johnson-Cook model cannot realize the coupled calculation of plastic deformation and creep, so it is necessary to use a unified viscoplastic model to model the mechanical behavior of cladding materials.
The expression of the unified viscoplastic constitutive equation usually includes two aspects: I. the choice of the viscosity function, that is, the expression of the relationship between viscous stress and equivalent plastic strain rate; and II. the choice of the hardening equations, that is, the expression of the internal variables [19,20].The viscosity equation and the hardening equation are connected through the yield function including overstress , in order to implement the coupled calculation of plasticity and creep.Chaboche and Rousselier [21,22] succeeded in modeling the rate-dependent mechanical behavior of 316 Stainless Steel based on a unified viscoplastic model including a nonlinear isotropic hardening rule, nonlinear kinematic hardening rule, and viscosity equation in the form of the power function.Since then, this model has been widely used to model the nonlinear mechanical behavior of metals and polymers, and it is called the Chaboche model [23][24][25][26][27]. Unlike 316 Stainless Steel, several experimental studies have been carried out to show that Zircaloy-4 has anisotropic viscoplastic properties due to its microstructural texture [28][29][30][31][32], and that Zircaloy-4 has a more complex kinematic hardening behavior.Delobelle [33][34][35] developed a unified viscoplastic constitutive model for Zircaloy-4 modeling by including a viscosity equation in the form of the hyperbolic sine function and three back-stresses playing roles successively instead of simultaneously, with some complicated coupling criteria, and he implemented an approach to modeling the anisotropic yield and plastic flow behavior of Zircaloy-4 by using four fourth-order tensors.
Unfortunately, the above two unified viscoplastic constitutive models are not built into the ABAQUS materials library, and there is neither a detailed description of the numerical solution procedure for the unified viscoplastic constitutive model used for modeling nuclear materials nor the available UMAT codes in the reported literature.In this paper, the Chaboche and Delobelle viscoplastic models for 316 stainless steel and Zircaloy-4, respectively, are implemented via the ABAQUS user material subroutine (UMAT).The specific details of the implementation of these two complex constitutive models, based on the return mapping algorithm and the semi-implicit integration method, are given, and the numerical simulation results are compared with the reported experimental results in order to validate the model and UMAT code.Compared with the existing in-house codes, the two user subroutines developed based on commercial software in this paper allowed open-source access; this can make full use of the advantages of ABAQUS in structural geometric modeling, multi-axial stress analysis and nonlinear finite element solutions, in order to facilitate a more accurate mechanical analysis of cladding structures based on ABAQUS in academia and industry.

The Framework of the Unified Viscoplastic Constitutive Model
Based on the infinitesimal strain theory, the total strain of the constitutive model can be decomposed into the elastic part and viscoplastic part, and the physical mechanism of this method is described in detail in the textbook of Dunne and Petrinic [36].So, the expression for the total strain ε is as follows: where ε e is the elastic strain tensor and ε vp is the viscoplastic strain tensor.The constitutive relation of the elastic deformation part can be expressed by the generalized Hooke's law as follows: where E is the fourth-order elastic tensor, and σ is the Cauchy stress tensor.The plastic strain rate εvp is expressed based on the associated flow rule as follows: where f is the yield function, which is equal to the plastic potential function in the associated flow rule and is defined as follows: where σ eq is the equivalent stress, which is a certain norm of the Cauchy stress tensor, and the specific expression of the equivalent stress is determined by the yield criterion.σ y is the initial yield stress.X i and X i are internal variables in scalar form and second-order tensor form, respectively, which determine the hardening behavior in the viscoplastic constitutive equation.∂ f /∂σ in Equation ( 3) means that the direction of the plastic strain rate is perpendicular to the normal direction of the yield surface, which is known as the normality hypothesis of plasticity.p is the equivalent plastic strain, also called accumulated plastic strain, which reflects the loading history of the materials.ṗ is the equivalent plastic strain rate, indicating the length of the plastic strain increment path in the plastic strain space.The expression of ṗ under the J2 flow rule is as follows: The equivalent plastic strain rate is determined by the viscosity equation φ: where σ v is the viscous stress, also known as overstress [19], which is the time-dependent part of the total stress in the constitutive equation.The yield function including overstress is as follows: Equation (7) shows that the plastic flow behavior described by the unified viscoplastic constitutive model is jointly determined by the hardening rules and the viscosity equation.
The viscosity equation can also be defined as follows: In summary, to carry out the modeling of material mechanical behavior based on the unified viscoplastic constitutive model, the following three parts need to be determined: I. the yield criterion, that is, the expression of the equivalent stress σ eq ; II. the hardening rule, that is, the expression of the internal variables X i and X i ; and III. the expression of the viscosity equation ϕ.

Isotropic Viscoplastic Constitutive Model for 316 Stainless Steel
The viscoplastic constitutive modeling of 316 Stainless Steel is based on the Chaboche model [21,22], and the yield function f expressions based on the von Mises yield criterion, nonlinear isotropic hardening rule and nonlinear kinematic hardening rule are as follows: where s is the deviatoric stress tensor: where tr(σ) is the trace of the stress tensor, that is, tr(σ) = σ 11 + σ 22 + σ 33 .I is the unit second-order tensor, that is, I = δ ij e i e j .The scalar internal variable R in Equation ( 9) is the drag-stress describing the isotropic hardening behavior of the materials, and the secondorder tensor internal variable α in Equation ( 9) is the back-stress describing the kinematic hardening behavior of the materials.The expressions of drag-stress R and back-stress α evolution are as follows: where b and Q are isotropic hardening parameters, and C 1 , C 2 , a 1 and a 2 are kinematic hardening parameters.Equation ( 13) is the nonlinear kinematic hardening rule proposed by Armstrong and Frederick [37].Chaboche [21] first used the back-stress expression of the addition of multiple Armstrong-Frederick rules to obtain better modeling capabilities of kinematic hardening behavior.The viscosity equation φ is defined by Norton's power law: where K and n are viscoplastic parameters.Substituting Equations ( 9) and ( 14) into Equation ( 3), the expression of viscoplastic strain rate based on the Chaboche model can be obtained as follows: Specific details about the implementation of the isotropic viscoplastic constitutive model for 316 Stainless Steel are detailed in Section 3.2.

Anisotropic Viscoplastic Constitutive Model for Zircaloy-4
The viscoplastic constitutive modeling of Zircaloy-4 is carried out based on the Delobelle model [33][34][35].The Hill yield criterion [38] is used to describe the unique anisotropic yield behavior of Zircaloy-4, and the expression of the yield function f is as follows: where M is the fourth-order anisotropy coefficient tensor.Unlike the 316 Stainless Steel's model, there is no initial yield stress term in the yield function of the Zircaloy-4 model, indicating that there is no initial elastic domain in the Zircaloy-4 model, and the initial Hill yield surface is a point.Considering Voigt notation, we have the following: The modeling is based on the assumption of orthotropy and incompressibility because of the essentially radial orientation of the crystallites in the Zircaloy-4 microstructure [34], which leads to only six independent components of the anisotropy coefficient tensor.The fourth-order coefficient tensors N, Q and R, that describe the anisotropic plastic behavior below also have the same form.The hardening internal variable expression is given by the following: α α : R : α α eq αeq ( 24) where q 1 and q 2 are kinematic hardening parameters, and b iso , b θ , Y 0 , Y sat iso and Y sat θ are scalar internal variable hardening parameters.Equations ( 18)-( 20) are kinematic hardening rules in the form of an Armstrong-Frederick rule.The second term on the right side of Equation ( 18) is the static recovery term, which is used to describe the hardening recovery behavior independent of plastic strain, and r m , r m1 , m 0 and m 1 are temperature-dependent material parameters.The internal variable Y * is the asymptotic value of the back-stress in the evolution process, which has a similar expression to the nonlinear isotropic hardening rule.The viscosity equation φ is defined as follows: where ε 0 , N and k are viscoplastic parameters.Creep test results are well described using the viscosity equation in the form of a hyperbolic sine function [33].Substituting Equations ( 16) and ( 26) into Equation ( 3), the expression of viscoplastic strain rate based on the Delobelle model can be obtained as: Specific details about the implementation of the anisotropic viscoplastic constitutive model for Zircaloy-4 are detailed in Section 3.3.

Computational Algorithm
The unified viscoplastic constitutive model proposed in Section 2 was implemented through the ABAQUS user material subroutine.The stress update algorithms of the 316 Stainless Steel model and Zircaloy-4 model are deduced based on semi-implicit integration and return mapping algorithm.The algorithm derived in this section is based on the discretization of the unified viscoplastic constitutive equation, that is, for a time interval [t n , t n+1 ], the variable values at t n being known and the total strain increment ∆ε being given, the purpose of the algorithm detailed hereafter is to calculate the variable values at time t n+1 .

Return Mapping Algorithm
The return mapping algorithm is widely used in computational plasticity [36,39].In general, this approach consists of two parts called "elastic predictor" and "plastic corrector".The generalized Hooke's law given by Equation ( 2) can also be written as follows: where G is the shear modulus and λ is the Lame's constant.Elastic strain can be decomposed as follows: For simplicity, the subscript t in this section indicates the variable at time t n , and the variable without subscript t indicates the variable at time t n+1 .Substitute Equation (29) into Equation ( 28): Considering the incompressible assumption, that is tr(∆ε vp ) = 0, Equation ( 30) can be expressed as follows: Let: where σ tri is the trial stress, also known as "elastic predictor".2G∆ε vp is called "plastic corrector".Substitute Equation (32) into Equation ( 31): After derivation, it can be obtained from the above equation: According to the above derivation, we can get the basic calculation process of the return mapping algorithm as shown in Figure 1: First, assuming that the total strain increments are all elastic strain increments, the yield function f tri of the trial stress state can be obtained (x-y in Figure 1): Secondly, if there is f tri ≤ 0, it means that the material has not entered the yield state, and the stress σ tri eq corresponding to the trial stress is the updated stress.If there is f tri > 0, it means that the material enters the yield state.At this time, the yield function of the real state should be calculated (y-z in Figure 1).Considering Equation (34), the yield function of the true state is as follows: Therefore, obtaining the equivalent plastic strain increment ∆p that satisfies both the hardening rules and the viscous equation is the key problem of the stress update algorithm.Specific details about the derivation of the equivalent plastic strain increment ∆p of the 316 Stainless Steel model and the Zircaloy-4 model are detailed in Sections 3.2 and 3.3.Finally, if ∆p is obtained, the updated stress can be calculated based on the flow rule and generalized Hooke's law (x-z in Figure 1): where n is the direction of plastic flow, that is, n = ∂ f /∂σ.As shown in Equation ( 38), the algorithm implementation of the constitutive model in this paper is carried out based on the constant stiffness iterative method.The scheme of the return mapping algorithm process is illustrated in Figure 2.

Stress Update Algorithm of 316 Stainless Steel Model
The isotropic viscoplastic constitutive model for 316 Stainless Steel in Section 2.2 is discretized according to the semi-implicit integration scheme as follows: According to Equation ( 14): Since the expressions of α and R both contain ∆p n+1 , Equation ( 47) is a nonlinear equation about ∆p n+1 , which can be solved using the Newton-Raphson iterative method.Therefore, it is necessary to derive the Newton-Raphson iteration scheme, rearrange Equation (47), and let: where: ( ∂φ ∂σ eq Calculate the updated internal variables and equivalent plastic strain increment: The above is the Newton-Raphson iterative scheme for solving the nonlinear equation about ∆p (n+1) .When RES = |∆p − φ( f )∆t| → 0, the iteration ends, that is, the ∆p that satisfies both the hardening rules and the viscosity function has been found within the time incremental step [t n , t n+1 ]: Calculate the updated plastic strain increment: where: Calculate the updated stress:

Stress Update Algorithm of Zircaloy-4 Model
The anisotropic viscoplastic constitutive model for Zircaloy-4 in Section 2.3 is discretized according to the semi-implicit integration scheme as follows: (69) ∆α According to Equation ( 26): Substituting Equations ( 18)-( 20) into Equation (82), eliminate dα and rearrange to give (k is the number of iterations): where: Calculate the updated internal variables and equivalent plastic strain increment: The iteration ends when RES = |∆p − φ( f )∆t| → 0: The plastic strain increment is then determined from: where: Calculate the updated stress using Equations ( 65) and (66).

Verification and Validation of FE Modeling and UMAT Code
In order to verify the accuracy of the unified viscoplastic constitutive model and its algorithm implementation, numerical simulations of the viscoplastic mechanical behavior of materials were carried out using the ABAQUS user material subroutine of the 316 Stainless Steel model and the Zircaloy-4 model, and these were compared with the experimental results.The reported experimental studies are described in the cylindrical coordinate system due to the use of hollow circular cross-section rod specimens [22,34,40].In order to correspond to the Cartesian coordinates system adopted for the theory and modeling in this paper, we have: 11 = x = r, 22 = y = θ, 33 = z = z, 12 = xy = rθ, 13 = xz = rz and 23 = yz = θz.
The user must correctly define the number of "*Depvar" and the parameters of "*User Material" in order for the UMAT subroutine to work properly.The former is the number of solution-dependent state variables in the UMAT subroutine, corresponding to "STATEV(K)" in the UMAT code, and the latter is the material parameter, corresponding to "PROP(K)" in the UMAT code.The above operations can be done in the ".inp" file or in the ABAQUS GUI environment.The finite element analysis carried out in this paper is based on the ABAQUS 2021 version.
Before carrying out the verification of the algorithm and UMAT code, the convergence of the 316 Stainless Steel and Zircaloy-4 subroutines should be verified based on a simple uniaxial tensile simulation.The global convergence of the UMAT code in the FE calculation is guaranteed by the ABAQUS implicit analysis solver, that is, the iteration of the incremental step will only end when the residual force reaches the required precision [41].In order to verify the local convergence of the UMAT subroutine, the RES-iteration curves of the UMAT at the last global iteration step in the different incremental steps of the 316 Stainless Steel and the Zircaloy-4 subroutines in the FE calculation are, respectively, output, as shown in Figure 3: RES is the residual equivalent plastic strain rate: RES = |φ − ∆p/∆t| (see the descriptions in Sections 3.2 and 3.3 for detail).It can be seen from Figure 3 that the two UMAT subroutines can achieve sufficient local convergence in each global iteration step.

316 Stainless Steel
Chaboche et al. [22] used the 316 Stainless Steel hollow circular cross-section rod specimen shown in Figure 4a to carry out uniaxial tensile tests at room temperature for different strain rates, and stress relaxation test at room temperature were also carried out: The gauge length of the specimen is a unidirectional stress state along the axis zdirection, so the finite element model is considered as a simple 1 × 1 × 1 cube, as shown in Figure 4b.The displacement boundary condition is defined in the z-direction of the model to simulate the uniaxial tensile behavior of the materials, and the symmetric boundary condition in the z-direction is applied on the other side.The model has meshed into 1000 C3D8 elements for simulation, as shown in Figure 4c.The material constants of 316 Stainless Steel at room temperature are shown in Table 1.[22,40]. Figure 5a is the stress-strain curve of 316 Stainless Steel for different strain rates at room temperature: The numerical simulation results are in good agreement with the experimental results, and the 316 Stainless Steel modeling and UMAT code are verified.When the stress is less than the yield stress σ y , the stress-strain curves for different strain rates overlap, and the stress-strain behavior of the materials show rate-dependence only after yielding, which is the basic characteristic of the elasto-viscoplastic constitutive model.From Equations ( 7) and ( 14), we have the following: The above equation shows that the stress of the materials is positively correlated with the strain rate, that is, the greater the strain rate when the material is loaded, the greater the stress when the material is loaded to the same strain, which is consistent with the phenomenon described by the stress-strain curve in Figure 5a. Figure 5b shows the stress relaxation curve of 316 Stainless Steel at room temperature, and the numerical simulation results are in good agreement with the experiment.When the materials maintain a certain strain, as time increases, the stress gradually decreases until the stress caused by the viscosity of the material disappears completely.It is worth noting that the viscous stress is caused by the plastic strain, as shown in Equation (97).Therefore, if the materials be held at a strain level corresponding to the stress that does not reach the yield stress, materials will not exhibit stress relaxation behavior.Startup-run-shutdown conditions or load fluctuation conditions are common operating states of nuclear reactor structures, such as cladding and pressure vessels, during service, and these conditions may lead to the fatigue failure of the materials [42].The simulation capabilities of the UMAT subroutine developed in this paper cover the low cycle fatigue behavior of materials.Hyde et al. [40] used the 316 Stainless Steel hollow circular cross-section rod specimen shown in Figure 4a to carry out the low cycle fatigue test (LCF test) at 550 • C. The test is controlled by strain, the loading waveform is a triangular wave, the strain amplitude is 0.006, the strain ratio is −1, and the period is 80 s.Numerical simulations also use the finite element model shown in Figure 4b,c, setting the same displacement boundary conditions in the z-direction as the LCF test.As shown in Figure 6a, the stress-strain loops predicted by the model for the first and fiftieth cycles are in good agreement with the test: Figure 6b shows the relationship between the half stress range and the number of cycles, and the numerical simulation results are in good agreement with the experimental results.The model successfully predicted the cyclic hardening behavior of 316 Stainless Steel at high temperature.The combined isotropic and kinematic hardening rules governing the cyclic hardening behavior in the constitutive equations are shown in Equations ( 11)-( 13).

Zircaloy-4
Delobelle et al. [34] used the Zircaloy-4 circular tube specimen shown in Figure 7a to carry out the uniaxial tensile test for different strain rates and creep test in the axial and hoop directions, the tests were carried out at 350 • C.
In the uniaxial tensile test, the circular tube specimen is subjected to a displacementcontrolled loading in the z-direction, and the gauge length of the specimen is a unidirectional stress state along the z-direction, so the finite element model is considered as a simple 1 × 1 × 1 cube, as shown in Figure 7b.The displacement boundary condition is defined in the z-direction of the model to simulate the uniaxial tensile behavior of the specimen, and the symmetric boundary condition in the z-direction is applied on the other side.The model has meshed into 1000 C3D8 elements for simulation, as shown in Figure 7c.Two different metallurgical states of Zircaloy-4 were used for the test, namely recrystallized (R state) Zircaloy-4 and cold worked stress relieved (CWSR state) Zircaloy-4.The Zircaloy-4 material constants of the two metallurgical states are shown in Table 2.
Compared with Zircaloy-4 in the CWSR state, the material constants of Zircaloy-4 in the R state lack the parameters of the internal variable Y iso shown in Equation ( 22), because the Zircaloy-4 in the R state has cyclic loading stability, and Y iso mainly controls the cyclic hardening behavior of the materials.Therefore, the Y iso term is not considered in the modeling of the viscoplastic mechanical behavior of Zircaloy-4 in R state, and the relevant parameters can be set to zero when using UMAT to carry out the simulation of Zircaloy-4 in R state.Figure 8 shows the stress-strain curves of R state and CWSR state Zircaloy-4 for different strain rates at 350 • C.
Based on the test device in Figure 7a, Delobelle et al. [34] carried out uniaxial creep tests of Zircaloy-4 in different material orientations.The loading method of the uniaxial creep test in the z-direction is the same as that of the uniaxial tensile test in the z-direction, and the force sensor located at the clamping end of the ensures that the Zircaloy-4 round tube is in a constant stress loading state.The hoop creep loading of the Zircaloy-4 circular tube is controlled by the pressurizing system, and a constant hoop stress is applied to the specimen to carry out the creep test.This loading scheme is closer to the service conditions of the materials used in nuclear fuel cladding structures.Since the gauge length of the Zircaloy-4 tube is a unidirectional stress state when the creep tests are carried out in two directions, the finite element model is still considered as a simple 1 × 1 × 1 cube.As shown in Figure 7d, the stress boundary condition is defined in the z-direction of the model to simulate the uniaxial creep behavior of the specimen in the axial direction, and the symmetrical boundary condition in the z-direction is applied on the other side.As shown in Figure 7e, the stress boundary condition is defined in the y-direction of the model to simulate the uniaxial creep behavior of the specimen in the hoop direction, and the symmetric boundary condition in the y-direction is applied on the other side.The model has meshed into 1000 C3D8 elements for simulation, as shown in Figure 7c.The creep stresses of Zircaloy-4 in the R state are 125 MPa, 135 MPa, 150 MPa, and 170 MPa, respectively, and the creep stresses of Zircaloy-4 in the CWSR state are 170 MPa and 250 MPa, and each creep stress level is maintained for 200 hours and then increased to the next stress level.Figure 9 shows the incremental uniaxial creep curves of Zircaloy-4 in R state and CWSR state at 350 • C: The numerical simulation results are basically consistent with the experimental results.Neglecting the uncertainties in the loading control and deformation acquisition of the high-temperature creep test, the strain amplitudes, as well as the flow directions, are, on the whole, correctly modeled, so the simulation results are generally acceptable.The model and UMAT code successfully simulated the anisotropic creep behavior of Zircaloy-4 in the R state and the CWSR state, and the anisotropic creep behavior of the materials was more obvious at higher creep stress.

Multiaxial Stress State Analysis
In order to verify the ability of the developed two subroutines to simulate the multiaxial stress state, the finite element analysis of the thick-walled pipe structure under internal pressure based on the 316 Stainless Steel and Zircaloy-4 subroutines was carried out.The geometric models and boundary conditions were recommended by the Axisymmetric thick cylinder benchmark in the National Agency for Finite Element Methods and Standards (NAFEMS) [41], as shown in Figure 10: The geometric model is to be a 1/4 model of the thick-walled pipe structure on the XY plane.As shown in Figure 10a, symmetrical boundary conditions are applied on the two symmetrical surfaces, and the z-direction displacement is constrained on the upper and lower surfaces to ensure that the analysis results of the model can be arrayed in the z-direction.Uniform internal pressure is applied on the inside surface of the pipe, in which the 316 Stainless Steel model applies the internal pressure (IP) from 0 to 200 MPa with a ramp loading waveform within 20 s, and the Zircaloy-4 model applies the internal pressure from 0 to 300 MPa with a ramp loading waveform within 30 s.The model has meshed into 2300 C3D8 elements for simulation, as shown in Figure 10b.The results of the multiaxial stress state analysis are shown in Figure 11.
Figure 11a,b are the von Mises stress distribution diagrams of the 316 Stainless Steel and Zircaloy-4 thick-walled pipe models respectively.It can be seen that the von Mises stress gradually decreases from the inside surface of the pipe to the outside surface.In the cylindrical coordinate system, the radial stress σ r and hoop stress σ θ are output along path A and path B, and they are compared with the results of NAFEMS benchmark, as shown in Figure 11c,d.Since the NAFEMS benchmark is based on the elastic-perfectly plastic constitutive model, only qualitative comparative analysis is carried out in this paper.It can be seen that the radial stress σ r is the compressive stress near the inside pipe wall, and close to zero near the outside pipe wall.As the IP increases, the compressive stress on the inside pipe wall increases gradually.The hoop stress σ θ is tensile stress, and the distribution of hoop stress σ θ is related to IP.When the IP is small at the initial stage of loading, the σ θ near the inside pipe wall is greater than the σ θ near the outside pipe wall.With the increase of IP, the area of maximum hoop stress gradually moves to the outside pipe wall.When the IP increases to a higher level, the σ θ near the outside wall is larger than the σ θ near the inside wall.The reason for this phenomenon is that with the increase of IP, plastic flow occurs in the area close to the inside pipe wall, and the main load-bearing area of the structure gradually moves to the outside pipe wall.Compared with the NAFEMS benchmark, both the 316 Stainless Steel and Zircaloy-4 UMAT subroutines can complete the modeling of the multiaxial stress state of the thick-walled pipe structure under internal pressure.The simulation results successfully captured the rate-dependent stress-strain behavior, stress relaxation behavior and cyclic hardening behavior of 316 Stainless Steel, and the rate-dependent stress-strain behavior and anisotropic creep behavior of Zircaloy-4 in two metallurgical states.This demonstrates that the simulation capability of the user subroutines we developed covers the various loading conditions of the cladding materials in service.We hope our efforts provide useful material models and available source codes to help for users from academia or industry to carry out more accurate mechanical analyses of cladding structures.

Figure 1 .
Figure 1.The stress-strain update process based on return mapping algorithm, and the corresponding yield surface.

Figure 2 .
Figure 2. Flow diagram of return mapping algorithm.

Figure 5 .
Figure 5.Comparison of 316 Stainless Steel experimental and numerical results at room temperature: (a) uniaxial tensile stress-strain curves for different strain rates, (b) stress relaxation curve.

Figure 6 .
Figure 6.Comparison of model predictions of cyclic hardening behavior of 316 Stainless Steel with experimental data at 550 • C: (a) stress-strain loops of the 1st cycle and the 50th cycle, (b) half stress range (∆σ/2)-cycle number (N) curve.

Figure 7 .
Figure 7. Zircaloy-4 specimen and finite element model: (a) Zircaloy-4 specimen and loading system [30,34], (b) boundary condition for uniaxial loading in the z-direction, (c) meshing model, (d) boundary condition for uniaxial creep in the z-direction, (e) boundary condition for uniaxial creep in the y-direction.

Figure 8 .
Figure 8. Uniaxial tensile stress-strain curves of CWSR state and R state Zircaloy-4 for different strain rates.

Figure 11 .
Figure 11.The results of the thick-walled pipe structure FE analysis: (a) von Mises stress distribution of 316 Stainless Steel model, (b) von Mises stress distribution of Zircaloy-4 model, (c) Stress components along path A compared with NAFEMS benchmark, (d) Stress components along path B compared with NAFEMS benchmark.
Within this study, we present the implementation of the ABAQUS user material subroutines of the Chaboche unified viscoplastic constitutive model and the Delobelle unified viscoplastic constitutive model for 316 Stainless Steel and Zircaloy-4.The Chaboche model includes nonlinear isotropic hardening, nonlinear kinematic hardening, and a viscosity equation in the form of a power function.The Delobelle model includes a viscosity equation in the form of the hyperbolic sine function and three back-stresses playing roles successively instead of simultaneously, and the anisotropic viscoplastic behavior of the materials is considered.The derivation of the stress update algorithm that implements the Chaboche model and the Delobelle model is based on the return mapping algorithm and the semi-implicit integration scheme, and the UMAT user subroutines for 316 Stainless Steel and Zircaloy-4 are developed.The simulation results with the use of subroutines were verified against existing experimental data from the literature.The models and UMAT subroutines successfully simulated the viscoplastic mechanical behavior of 316 Stainless Steel and Zircaloy-4.The numerical implementation scheme of the constitutive model in this paper is not limited to the implementation of the Chaboche model and the Delobelle model but, as a set of general methods can also be used for the numerical implementation of other material rate-dependent plastic behavior modeling.

Table 1 .
Material constants for 316 Stainless Steel at room temperature and 550 • C