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

1. Introduction

The stress corrosion cracking (SCC) of cladding tube structures caused by pellet clad interaction (PCI) seriously threatens the operation safety of nuclear facilities[1,2]. 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[3,4]. 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[5,6,7]. 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[8]. This approach greatly limits the generalizability and accuracy of the models despite their simplicity in theory and implementation. Therefore, since 2008, Williamson et al from Idaho National Laboratory (INL) have been developing the multidimensional finite element fuel code BISON[9,10,11]. At the same time, Williamson[12] carried out the multidimensional thermomechanical analysis of multipellet steady and transient light water reactor fuel rod behavior based on ABAQUS commercial finite element code. A power-law creep formulation for describing the deformation behavior of Zircaloy cladding is implemented via ABAQUS user creep subroutine. However, 316 stainless steel and Zircaloy-4 have coupled plasticity and creep behavior at high temperatures[13], 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; II. the choice of the hardening equations, that is, the expression of the internal variables[14]. The viscosity equation and the hardening equation are connected through the yield function including overstress to implement the coupled calculation of plasticity and creep. Chaboche and Rousselier[15,16] succeeded in modeling the rate-dependent mechanical behavior of 316 Stainless Steel based on a unified viscoplastic model including 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[17,18,19,20,21,22]. 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[23,24,25,26]. Delobelle[27,28,29] 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 role 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. 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 are given, and numerical simulation results are compared with reported experimental results to validate the model and UMAT code. Two user subroutines developed in this work allowed open-source access to facilitate a more accurate mechanical analysis of cladding structures based on ABAQUS in academia and industry.

2. Viscoplastic Constitutive Equation

2.1. The Framework of Unified Viscoplastic Constitutive Model

Based on the Infinitesimal strain theory, the total strain of the constitutive model can be decomposed into elastic part and viscoplastic part, and the physical mechanism of this method is described in detail in the textbook of Dunne and Petrinic[30]. So the expression for the total strain $\mathit{\epsilon }$ is:

$$\mathit{\epsilon }={\mathit{\epsilon }}^e+{\mathit{\epsilon }}^{vp}$$

(1)

where ${\mathit{\epsilon }}^e$ is the elastic strain tensor and ${\mathit{\epsilon }}^{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:

$${\mathit{\epsilon }}^e={\mathit{E}}^{-1}:\mathit{\sigma }$$

(2)

where $\mathit{E}$ is the fourth-order elastic tensor, and $\mathit{\sigma }$ is the Cauchy stress tensor. The plastic strain rate ${\dot{\mathit{\epsilon }}}^{vp}$ is expressed based on the associated flow rule as follows:

$${\dot{\mathit{\epsilon }}}^{vp}=\dot{p}\frac{\partial f}{\partial \mathit{\sigma }}$$

(3)

where f is the yield function, which is equal to the plastic potential function in the associated flow rule and is defined as:

$$f=f({\sigma }_{eq},{\sigma }_y,p,X_i,{\mathit{X}}_i),i=1,2,\cdots$$

(4)

where ${\sigma }_{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. ${\sigma }_y$ is the initial yield stress. $X_i$ and ${\mathit{X}}_i$ are internal variables in scalar form and second-order tensor form, respectively, which determine the hardening behavior in the viscoplastic constitutive equation. $\partial f/\partial \mathit{\sigma }$ 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. $\dot{p}$ is the equivalent plastic strain rate, indicating the length of the plastic strain increment path in the plastic strain space. The expression of $\dot{p}$ under the J2 flow rule is:

$$\dot{p}={\left(\frac{2}{3}{\dot{\mathit{\epsilon }}}^{vp}:{\dot{\mathit{\epsilon }}}^{vp}\right)}^{1/2}$$

(5)

The equivalent plastic strain rate is determined by the viscosity equation $\varphi$:

$$\dot{p}=\varphi \left({\sigma }_v\right)$$

(6)

where ${\sigma }_v$ is the viscous stress, also known as overstress[14], which is the time-dependent part of the total stress in the constitutive equation. The yield function including overstress is as follows:

$$f={\sigma }_v>0$$

(7)

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:

$$\left\{\begin{array}{cc}\dot{p}=0\hfill & ,f\le 0\hfill \\ \dot{p}=\varphi \left(f\right)>0\hfill & ,f>0\hfill \end{array}\right.$$

(8)

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. Yield criterion, that is, the expression of the equivalent stress ${\sigma }_{eq}$; II. Hardening rule, that is, the expression of the internal variables $X_i$ and ${\mathit{X}}_i$; III. The expression of the viscosity equation $\phi$.

2.2. Isotropic Viscoplastic Constitutive Model for 316 Stainless Steel

The viscoplastic constitutive modeling of 316 Stainless Steel is based on the Chaboche model[15,16], and the yield function f expressions based on the von Mises yield criterion, nonlinear isotropic hardening rule and nonlinear kinematic hardening rule are:

$$\begin{array}{cc}\hfill f& ={\sigma }_{eq}-R-{\sigma }_y\hfill \\ & ={\left[\frac{3}{2}(\mathit{s}-\mathit{\alpha }):(\mathit{s}-\mathit{\alpha })\right]}^{1/2}-R-{\sigma }_y\hfill \end{array}$$

(9)

where $\mathit{s}$ is the deviatoric stress tensor:

$$\mathit{s}=\mathit{\sigma }-\frac{1}{3}tr\left(\mathit{\sigma }\right)\mathit{I}$$

(10)

where $tr\left(\mathit{\sigma }\right)$ is the trace of the stress tensor, that is, $tr\left(\mathit{\sigma }\right)={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}$. $\mathit{I}$ is the unit second-order tensor, that is, $\mathit{I}={\delta }_{ij}{\overrightarrow{e}}_i{\overrightarrow{e}}_j$. The scalar internal variable R in Equation (9) is the drag-stress describing the isotropic hardening behavior of the materials, and the second-order tensor internal variable $\mathit{\alpha }$ in Equation (9) is the back-stress describing the kinematic hardening behavior of the materials. The expressions of drag-stress R and back-stress $\mathit{\alpha }$ evolution are:

$$\dot{R}=b(Q-R)\dot{p}$$

(11)

$$\sum _{i=1}^2{\mathit{\alpha }}_i$$

(12)

$${\dot{\mathit{\alpha }}}_i=C_i\left(\frac{2}{3}{a}_i{\dot{\mathit{\epsilon }}}^{vp}-{\mathit{\alpha }}_i\dot{p}\right)$$

(13)

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[31]. Chaboche[15] first used the back-stress expression of the addition of multiple Armstrong-Frederick’s rules to obtain better modeling capabilities of kinematic hardening behavior. The viscosity equation $\varphi$ is defined by Norton’s power law:

$$\dot{p}=\varphi \left(f\right)={\left(\frac{f}{K}\right)}^n$$

(14)

where K and n are viscoplastic parameters. Substituting Equation (9) and (14) into Equation (3), the expression of viscoplastic strain rate based on the Chaboche model can be obtained as:

$${\dot{\mathit{\epsilon }}}^{vp}=\frac{3}{2}{\left(\frac{f}{K}\right)}^n\frac{\mathit{s}-\mathit{\alpha }}{{\sigma }_{eq}},f>0$$

(15)

Specific details about the implementation of the isotropic viscoplastic constitutive model for 316 Stainless Steel are detailed in Section 3.2.

2.3. Anisotropic Viscoplastic Constitutive Model for Zircaloy-4

The viscoplastic constitutive modeling of Zircaloy-4 is carried out based on the Delobelle model[27,28,29]. The Hill yield criterion[32] is used to describe the unique anisotropic yield behavior of Zircaloy-4, and the expression of the yield function f is:

$$f={\sigma }_{eq}={\left[\frac{3}{2}(\mathit{s}-\mathit{\alpha }):\mathit{M}:(\mathit{s}-\mathit{\alpha })\right]}^{1/2}$$

(16)

where $\mathit{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 Zircaloy-4 model, and the initial Hill yield surface is a point. Considering Voigt notation, we have:

$$\mathit{M}=\left(\begin{array}{cccccc}{M}_{11}& M_{12}& M_{13}& 0& 0& 0\\ M_{12}& M_{12}& M_{23}& 0& 0& 0\\ M_{13}& M_{23}& M_{33}& 0& 0& 0\\ 0& 0& 0& M_{44}& 0& 0\\ 0& 0& 0& 0& M_{55}& 0\\ 0& 0& 0& 0& 0& M_{66}\end{array}\right)with\begin{array}{c}\hfill M_{11}+M_{12}+M_{13}=0\\ \hfill M_{12}+M_{22}+M_{23}=0\\ \hfill M_{13}+M_{23}+M_{33}=0\end{array}$$

(17)

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[28], which leads to only six independent components of the anisotropy coefficient tensor. The fourth-order coefficient tensors $\mathit{N}$, $\mathit{Q}$ and $\mathit{R}$ that describe the anisotropic plastic behavior below also have the same form. The hardening internal variable expression is given by:

$$\dot{\mathit{\alpha }}=q\left[\frac{2}{3}{Y}^{*}\mathit{N}:{\dot{\mathit{\epsilon }}}^{vp}-\mathit{Q}:\left(\mathit{\alpha }-{\mathit{\alpha }}^{\left(1\right)}\right)\dot{p}\right]-\left[r_m{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_0}+r_{m1}{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_1}\right]:\mathit{N}:\mathit{R}:\frac{\mathit{\alpha }}{{\alpha }_{eq}}$$

(18)

$${\dot{\mathit{\alpha }}}^{\left(1\right)}=q_1\left[\frac{2}{3}{Y}^{*}\mathit{N}:{\dot{\mathit{\epsilon }}}^{vp}-\mathit{Q}:\left({\mathit{\alpha }}^{\left(1\right)}-{\mathit{\alpha }}^{\left(2\right)}\right)\dot{p}\right]$$

(19)

$${\dot{\mathit{\alpha }}}^{\left(2\right)}=q_2\left[\frac{2}{3}{Y}^{*}\mathit{N}:{\dot{\mathit{\epsilon }}}^{vp}-\mathit{Q}:{\mathit{\alpha }}^{\left(2\right)}\dot{p}\right]$$

(20)

$$Y^{*}=Y_0+Y_{iso}+Y_{\theta }$$

(21)

$${\dot{Y}}_{iso}=b_{iso}\left(Y_{iso}^{sat}-Y_{iso}\right)\dot{p}$$

(22)

$${\dot{Y}}_{\theta }=b_{\theta }\left(Y_{\theta }^{sat}f\left(\theta \right)-Y_{\theta }\right)\dot{p}$$

(23)

$$f\left(\theta \right)=1-|cos\theta |,\theta =arccos\left(\frac{3}{2}\frac{\mathit{\alpha }:\mathit{R}:\dot{\mathit{\alpha }}}{{\alpha }_{eq}{\dot{\alpha }}_{eq}}\right)$$

(24)

$${\alpha }_{eq}={\left(\frac{3}{2}\mathit{\alpha }:\mathit{R}:\mathit{\alpha }\right)}^{1/2},{\dot{\alpha }}_{eq}={\left(\frac{3}{2}\dot{\mathit{\alpha }}:\mathit{R}:\dot{\mathit{\alpha }}\right)}^{1/2}$$

(25)

where $q_1$ and $q_2$ are kinematic hardening parameters, and $b_{iso}$, $b_{\theta }$, $Y_0$, $Y_{iso}^{sat}$ and $Y_{\theta }^{sat}$ are scalar internal variable hardening parameters. Equations (18)(19)(20) are kinematic hardening rules in the form of Armstrong-Frederick’s 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 $\varphi$ is defined as:

$$\dot{p}=\varphi \left(f\right)={\epsilon }_0{\left[sinh\left(\frac{f}{N}\right)\right]}^k$$

(26)

where ${\epsilon }_0$, N and k are viscoplastic parameters. Creep test results are well described using the viscosity equation in the form of hyperbolic sine function[27]. Substituting Equation (16) and (26) into Equation (3), the expression of viscoplastic strain rate based on the Delobelle model can be obtained as:

$${\dot{\mathit{\epsilon }}}^{vp}=\frac{3}{2}{\epsilon }_0{\left[sinh\left(\frac{f}{N}\right)\right]}^k\frac{\mathit{M}:\left(\mathit{s}-\mathit{\alpha }\right)}{{\sigma }_{eq}},f>0$$

(27)

Specific details about the implementation of the anisotropic viscoplastic constitutive model for Zircaloy-4 are detailed in Section 3.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 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 $\Delta \epsilon$ being given, and the purpose of the algorithm detailed hereafter is to calculate the variable values at time $t_{n+1}$.

3.1. Return Mapping Algorithm

The return mapping algorithm is widely used in computational plasticity[30,33]. 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:

$$\mathit{\sigma }=2G{\mathit{\epsilon }}^e+\lambda tr\left({\mathit{\epsilon }}^e\right)\mathit{I}$$

(28)

where G is the shear modulus and $\lambda$ is the Lame’s constant. Elastic strain can be decomposed as follows:

$${\mathit{\epsilon }}^e={\mathit{\epsilon }}_t^e+\Delta {\mathit{\epsilon }}^e={\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }-\Delta {\mathit{\epsilon }}^{vp}$$

(29)

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):

$$\mathit{\sigma }=2G\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }-\Delta {\mathit{\epsilon }}^{vp}\right)+\lambda tr\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }-\Delta {\mathit{\epsilon }}^{vp}\right)\mathit{I}$$

(30)

Considering the incompressible assumption, that is $tr(\Delta {\mathit{\epsilon }}^{vp})=0$, Equation (30) can be expressed as:

$$\mathit{\sigma }=2G\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }\right)+\lambda tr\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }\right)\mathit{I}-2G\Delta {\mathit{\epsilon }}^{vp}$$

(31)

Let:

$${\mathit{\sigma }}^{tri}=2G\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }\right)+\lambda tr\left({\mathit{\epsilon }}_t^e+\Delta \mathit{\epsilon }\right)\mathit{I}$$

(32)

where ${\mathit{\sigma }}^{tri}$ is the trial stress, also known as “elastic predictor”. $2G\Delta {\mathit{\epsilon }}^{vp}$ is called “plastic corrector”. Substitute Equation (32) into Equation (31):

$$\mathit{\sigma }={\mathit{\sigma }}^{tri}-2G\Delta {\mathit{\epsilon }}^{vp}$$

(33)

After derivation, it can be obtained from the above equation:

$${\sigma }_{eq}={\sigma }_{eq}^{tri}-3G\Delta p$$

(34)

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 (①-② in Figure 1):

$$f^{tri}={\sigma }_{eq}^{tri}-{\sigma }_y$$

(35)

Secondly, if there is $f^{tri}\le 0$, it means that the material has not entered the yield state, and the stress ${\sigma }_{eq}^{tri}$ 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 (②-③ in Figure 1). Considering Equation (34), the yield function of the true state is:

$$f={\sigma }_{eq}-{\sigma }_y=\left({\sigma }_{eq}^{tri}-3G\Delta p\right)-{\sigma }_y={\sigma }_v$$

(36)

Therefore, obtaining the equivalent plastic strain increment $\Delta 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 $\Delta p$ of the 316 Stainless Steel model and the Zircaloy-4 model are detailed in Section 3.2 and Section 3.3. Finally, if $\Delta p$ is obtained, the updated stress can be calculated based on the flow rule and generalized Hooke’s law (①-③ in Figure 1):

$$\Delta {\mathit{\epsilon }}^{vp}=\Delta p\mathit{n}$$

(37)

$$\Delta \mathit{\sigma }={\mathit{E}}^{-1}\Delta {\mathit{\epsilon }}^e={\mathit{E}}^{-1}:\left(\Delta \mathit{\epsilon }-\Delta {\mathit{\epsilon }}^{vp}\right)$$

(38)

$$\mathit{\sigma }={\mathit{\sigma }}_t+\Delta \mathit{\sigma }$$

(39)

where $\mathit{n}$ is the direction of plastic flow, that is, $\mathit{n}=\partial f/\partial \mathit{\sigma }$. 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.

3.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:

$$f_{(n+1)}={\sigma }_{eq}-R_{(n+1)}-{\sigma }_y={\sigma }_{eq}^{tri}-3G\Delta p_{(n+1)}-R_{(n+1)}-{\sigma }_y$$

(40)

$${\mathit{\alpha }}_{(n+1)}={\mathit{\alpha }}_{\left(n\right)}+\Delta {\mathit{\alpha }}_{(n+1)}$$

(41)

$$\Delta {\mathit{\alpha }}_{(n+1)}=\Delta {\mathit{\alpha }}_{1(n+1)}+\Delta {\mathit{\alpha }}_{2(n+1)}$$

(42)

$$\Delta {\mathit{\alpha }}_{1(n+1)}=C_1\left(a_1\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-{\mathit{\alpha }}_{1\left(n\right)}\Delta p_{(n+1)}\right)$$

(43)

$$\Delta {\mathit{\alpha }}_{2(n+1)}=C_2\left(a_2\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-{\mathit{\alpha }}_{2\left(n\right)}\Delta p_{(n+1)}\right)$$

(44)

$$R_{(n+1)}=R_{\left(n\right)}+\Delta R_{(n+1)}$$

(45)

$$\Delta R_{(n+1)}=b\left(Q-R_{\left(n\right)}\right)\Delta p_{(n+1)}$$

(46)

According to the Equation 14:

$$\dot{p}=\frac{dp}{dt}=\frac{\Delta p}{\Delta t}=\varphi \left(f\right)=\varphi (\Delta p,\mathit{\alpha },R)$$

(47)

Since the expressions of $\mathit{\alpha }$ and R both contain $\Delta p_{n+1}$, the Equation 47 is a nonlinear equation about $\Delta p_{n+1}$, which can be solved by Newton-Raphson iterative method. Therefore, it is necessary to derive the Newton-Raphson iteration scheme, rearrange the Equation 47, and let:

$$\psi (\Delta p,\mathit{\alpha },R)=\Delta p-\varphi (\Delta p,\mathit{\alpha },R)\Delta t=0$$

(48)

Take the first-order Taylor expansion of $\psi$:

$$\psi +\frac{\partial \psi }{\partial \Delta p}d\Delta p+\frac{\partial \psi }{\partial \mathit{\alpha }}:d\mathit{\alpha }+\frac{\partial \psi }{\partial R}dR=0$$

(49)

Substituting Equation 48 into Equation 49:

$$\Delta p-\varphi \Delta t+\left(1-\frac{\partial \varphi }{\partial \Delta p}\Delta t\right)d\Delta p-\frac{\partial \varphi }{\partial \mathit{\alpha }}:d\mathit{\alpha }\Delta t-\frac{\partial \varphi }{\partial R}dR\Delta t=0$$

(50)

Substituting Equation 11, 12 and 13 into Equation 50, eliminate $d\mathit{\alpha }$, $dR$ and rearrange gives (k is the number of iteration):

$$d\Delta p^{(k+1)}=\frac{{\varphi }^{\left(k\right)}-\Delta p^{\left(k\right)}/\Delta t}{\frac{1}{\Delta t}-{\left(\frac{\partial \varphi }{\partial \Delta p}\right)}^{\left(k\right)}+\frac{3}{2}{\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}\left(C_1{a}_1+C_2{a}_2\right)-{\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}\left(C_1{\mathit{n}}_{\left(n\right)}:{\mathit{\alpha }}_1^{\left(k\right)}-C_2{\mathit{n}}_{\left(n\right)}:{\mathit{\alpha }}_2^{\left(k\right)}\right)-{\left(\frac{\partial \varphi }{\partial R}\right)}^{\left(k\right)}b\left(Q-R^{\left(k\right)}\right)}$$

(51)

where:

$${\varphi }^{\left(k\right)}={\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}-R^{\left(k\right)}-{\sigma }_y}{K}\right)}^n$$

(52)

$$\frac{\partial \varphi }{\partial \mathit{\alpha }}=-\frac{\partial \varphi }{\partial {\sigma }_{eq}}\mathit{n}$$

(53)

$${\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}=\frac{n}{K}{\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}-R^{\left(k\right)}-{\sigma }_y}{K}\right)}^{n-1}$$

(54)

$${\left(\frac{\partial \varphi }{\partial R}\right)}^{\left(k\right)}=-\frac{n}{K}{\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}-R^{\left(k\right)}-{\sigma }_y}{K}\right)}^{n-1}$$

(55)

$${\left(\frac{\partial \varphi }{\partial \Delta p}\right)}^{\left(k\right)}=-\frac{3nG}{K}{\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}-R^{\left(k\right)}-{\sigma }_y}{K}\right)}^{n-1}$$

(56)

Calculate the updated internal variables and equivalent plastic strain increment:

$$\Delta p^{(k+1)}=\Delta p^{\left(k\right)}+d\Delta p^{(k+1)}$$

(57)

$${\mathit{\alpha }}^{(k+1)}={\mathit{\alpha }}_1^{(k+1)}+{\mathit{\alpha }}_2^{(k+1)}$$

(58)

$${\mathit{\alpha }}_1^{(k+1)}={\mathit{\alpha }}_{1\left(n\right)}+\Delta {\mathit{\alpha }}_1^{(k+1)}={\mathit{\alpha }}_{1\left(n\right)}+C_1\left(a_1{\mathit{n}}_{\left(n\right)}-{\mathit{\alpha }}_{1\left(n\right)}\right)\Delta p^{(k+1)}$$

(59)

$${\mathit{\alpha }}_2^{(k+1)}={\mathit{\alpha }}_{2\left(n\right)}+\Delta {\mathit{\alpha }}_2^{(k+1)}={\mathit{\alpha }}_{2\left(n\right)}+C_2\left(a_2{\mathit{n}}_{\left(n\right)}-{\mathit{\alpha }}_{2\left(n\right)}\right)\Delta p^{(k+1)}$$

(60)

$$R^{(k+1)}=R_{\left(n\right)}+\Delta R^{(k+1)}=R_{\left(n\right)}+b\left(Q-R_{\left(n\right)}\right)\Delta p^{(k+1)}$$

(61)

The above is the Newton-Raphson iterative scheme for solving the nonlinear equation about $\Delta p_{(n+1)}$. When $RES=|\Delta p-\varphi (f)\Delta t|\to 0$, the iteration ends, that is, the $\Delta p$ that satisfies both the hardening rules and the viscosity function has been found within the time incremental step $[t_n,t_{n+1}]$:

$$\Delta p_{(n+1)}=\Delta p^{(k+1)}$$

(62)

Calculate the updated plastic strain increment:

$$\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}=\Delta p_{(n+1)}{\mathit{n}}_{\left(n\right)}$$

(63)

where:

$${\mathit{n}}_{\left(n\right)}={\left(\frac{3}{2}\frac{\mathit{s}-\mathit{\alpha }}{{\sigma }_{eq}}\right)}_{\left(n\right)}$$

(64)

Calculate the updated stress:

$$\Delta {\mathit{\sigma }}_{(n+1)}={\mathit{E}}^{-1}:\left(\Delta {\mathit{\epsilon }}_{(n+1)}-\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}\right)$$

(65)

$${\mathit{\sigma }}_{(n+1)}={\mathit{\sigma }}_{\left(n\right)}+\Delta {\mathit{\sigma }}_{(n+1)}$$

(66)

3.3. 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:

$$f_{(n+1)}={\sigma }_{eq}={\sigma }_{eq}^{tri}-3G\Delta p_{(n+1)}$$

(67)

$${\mathit{\alpha }}_{(n+1)}={\mathit{\alpha }}_{\left(n\right)}+\Delta {\mathit{\alpha }}_{(n+1)}$$

(68)

$${\mathit{\alpha }}_{(n+1)}^{\left(1\right)}={\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}+\Delta {\mathit{\alpha }}_{(n+1)}^{\left(1\right)}$$

(69)

$${\mathit{\alpha }}_{(n+1)}^{\left(2\right)}={\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}+\Delta {\mathit{\alpha }}_{(n+1)}^{\left(2\right)}$$

(70)

$$\begin{array}{cc}\hfill \Delta {\mathit{\alpha }}_{(n+1)}=& q\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:\left[{\mathit{\alpha }}_{\left(n\right)}-{\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}\right]\Delta p_{(n+1)}\right\}\hfill \\ & -\Delta t\left[r_m{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_0}+r_{m1}{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_1}\right]\mathit{N}:\mathit{R}:\frac{{\mathit{\alpha }}^{\left(n\right)}}{{\alpha }_{eq}}\hfill \end{array}$$

(71)

$$\Delta {\mathit{\alpha }}_{(n+1)}^{\left(1\right)}=q_1\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:\left[{\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}-{\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}\right]\Delta p_{(n+1)}\right\}$$

(72)

$$\Delta {\mathit{\alpha }}_{(n+1)}^{\left(2\right)}=q_2\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:{\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}\Delta p_{(n+1)}\right\}$$

(73)

$$Y_{(n+1)}^{*}=Y_0+Y_{iso(n+1)}+Y_{\theta (n+1)}$$

(74)

$$Y_{iso(n+1)}=Y_{iso\left(n\right)}+\Delta Y_{iso(n+1)}$$

(75)

$$Y_{\theta (n+1)}=Y_{\theta \left(n\right)}+\Delta Y_{\theta (n+1)}$$

(76)

$$\Delta Y_{iso(n+1)}=b_{iso}\left(Y_{iso}^{sat}-Y_{iso\left(n\right)}\right)\Delta p_{(n+1)}$$

(77)

$$\Delta Y_{\theta (n+1)}=b_{\theta }\left(Y_{\theta }^{sat}f{\left(\theta \right)}_{\left(n\right)}-Y_{\theta \left(n\right)}\right)\Delta p_{(n+1)}$$

(78)

According to the Equation 26:

$$\dot{p}=\frac{\Delta p}{\Delta t}=\varphi (\Delta p,\mathit{\alpha })$$

(79)

Similar to the 316 Stainless Steel model, the Newton-Raphson iterative scheme for solving the equivalent plastic strain increment $\Delta p_{(n+1)}$ of the Zircaloy-4 model is derived, rearrange the Equation 79, and let:

$$\psi (\Delta p,\mathit{\alpha })=\Delta p-\varphi (\Delta p,\mathit{\alpha })\Delta t=0$$

(80)

Take the first-order Taylor expansion of $\psi$:

$$\psi +\frac{\partial \psi }{\partial \Delta p}d\Delta p+\frac{\partial \psi }{\partial \mathit{\alpha }}:d\mathit{\alpha }=0$$

(81)

Substituting Equation 80 into Equation 81:

$$\Delta p-\varphi \Delta t+\left(1-\frac{\partial \varphi }{\partial \Delta p}\Delta t\right)d\Delta p-\frac{\partial \varphi }{\partial \mathit{\alpha }}:d\mathit{\alpha }\Delta t=0$$

(82)

Substituting Equation 18, 19 and 20 into Equation 82, eliminate $d\mathit{\alpha }$ and rearrange gives (k is the number of iteration):

$$d\Delta p^{(k+1)}=\frac{{\varphi }^{\left(k\right)}-\frac{\Delta p^{\left(k\right)}}{\Delta t}+{\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}\Delta t\left[r_m{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_0}+r_{m1}{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_1}\right]{\mathit{n}}_{\left(n\right)}:\mathit{N}:\mathit{R}:\frac{{\mathit{\alpha }}^{\left(\mathbf{k}\right)}}{{\alpha }_{eq}}}{\frac{1}{\Delta t}-{\left(\frac{\partial \varphi }{\partial \Delta p}\right)}^{\left(k\right)}+q{\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}\left\{\frac{2}{3}{Y}^{*\left(k\right)}{\mathit{n}}_{\left(n\right)}:\mathit{N}:{\mathit{n}}_{\left(n\right)}-{\mathit{n}}_{\left(n\right)}:\mathit{Q}:\left[{\mathit{\alpha }}^{\left(k\right)}-{\mathit{\alpha }}^{\left(1\right)\left(k\right)}\right]\right\}}$$

(83)

where:

$${\varphi }^{\left(k\right)}={\dot{\epsilon }}_0{\left\{sinh\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}}{N}\right)\right\}}^k$$

(84)

$${\left(\frac{\partial \varphi }{\partial {\sigma }_{eq}}\right)}^{\left(k\right)}=\frac{{\dot{\epsilon }}_{0}k}{N}{\left\{sinh\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}}{N}\right)\right\}}^{k-1}cosh\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}}{N}\right)$$

(85)

$${\left(\frac{\partial \varphi }{\partial \Delta p}\right)}^{\left(k\right)}=-\frac{3G{\dot{\epsilon }}_{0}k}{N}{\left\{sinh\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}}{N}\right)\right\}}^{k-1}cosh\left(\frac{{\sigma }_{eq}^{tri}-3G\Delta p^{\left(k\right)}}{N}\right)$$

(86)

Calculate the updated internal variables and equivalent plastic strain increment:

$$\Delta p^{(k+1)}=\Delta p^{\left(k\right)}+d\Delta p^{(k+1)}$$

(87)

$$\begin{array}{cc}\hfill {\mathit{\alpha }}^{(k+1)}=& {\mathit{\alpha }}_{\left(n\right)}+\Delta {\mathit{\alpha }}^{(k+1)}\hfill \\ \hfill =& {\mathit{\alpha }}_{\left(n\right)}+q\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:\left[{\mathit{\alpha }}_{\left(n\right)}-{\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}\right]\Delta p_{(n+1)}\right\}\hfill \\ & -\Delta t\left[r_m{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_0}+r_{m1}{\left(\frac{{\alpha }_{eq}}{{\alpha }_0}\right)}^{m_1}\right]\mathit{N}:\mathit{R}:\frac{{\mathit{\alpha }}^{\left(n\right)}}{{\alpha }_{eq}}\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill {\mathit{\alpha }}^{\left(1\right)(k+1)}=& {\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}+\Delta {\mathit{\alpha }}^{\left(1\right)(k+1)}\hfill \\ \hfill =& {\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}+q_1\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:\left[{\mathit{\alpha }}_{\left(n\right)}^{\left(1\right)}-{\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}\right]\Delta p_{(n+1)}\right\}\hfill \end{array}$$

(89)

$$\begin{array}{cc}\hfill {\mathit{\alpha }}^{\left(2\right)(k+1)}=& {\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}+\Delta {\mathit{\alpha }}^{\left(2\right)(k+1)}\hfill \\ \hfill =& {\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}+q_2\left\{\frac{2}{3}{Y}_{(n+1)}^{*}\mathit{N}:\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}-\mathit{Q}:{\mathit{\alpha }}_{\left(n\right)}^{\left(2\right)}\Delta p_{(n+1)}\right\}\hfill \end{array}$$

(90)

$$Y^{*(k+1)}=Y_0+Y_{iso}^{(k+1)}+Y_{\theta }^{(k+1)}$$

(91)

$$Y_{iso}^{(k+1)}=Y_{iso\left(n\right)}+\Delta Y_{iso}^{(k+1)}=Y_{iso\left(n\right)}+b_{iso}\left(Y_{iso}^{sat}-Y_{iso\left(n\right)}\right)\Delta p^{(k+1)}$$

(92)

$$Y_{\theta }^{(k+1)}=Y_{\theta \left(n\right)}+\Delta Y_{\theta }^{(k+1)}=Y_{\theta \left(n\right)}+b_{\theta }\left(Y_{\theta }^{sat}f{\left(\theta \right)}_{\left(n\right)}-Y_{\theta \left(n\right)}\right)\Delta p^{k+1}$$

(93)

The iteration ends when $RES=|\Delta p-\varphi (f)\Delta t|\to 0$:

$$\Delta p_{(n+1)}=\Delta p^{(k+1)}$$

(94)

The plastic strain increment is then determined from:

$$\Delta {\mathit{\epsilon }}_{(n+1)}^{vp}=\Delta p_{(n+1)}{\mathit{n}}_{\left(n\right)}$$

(95)

where:

$${\mathit{N}}_{\left(n\right)}={\left(\frac{3}{2}\frac{\mathit{M}:\left(\mathit{s}-\mathit{\alpha }\right)}{{\sigma }_{eq}}\right)}_{\left(n\right)}$$

(96)

Calculate the updated stress as Equation 65 and 66.

4. 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 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[16,18,28]. In order to correspond to the Cartesian coordinates system adopted for the theory and modeling in this paper, so we have: $11=x=r$, $22=y=\theta$, $33=z=z$, $12=xy=r\theta$, $13=xz=rz$ and $23=yz=\theta z$.

4.1. 316 Stainless Steel

Chaboche et al[16] used the 316 Stainless Steel hollow circular cross-section rod specimen shown in Figure 3a 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 z-direction, so the finite element model is considered as a simple $1\times 1\times 1$ cube, as shown in Figure 3b. 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 3c. The material constants of 316 Stainless Steel at room temperature are shown in Table 1:

Figure 4a 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 ${\sigma }_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 Equation 7 and 14, we have:

$${\sigma }_v=K{\dot{p}}^{1/n}$$

(97)

that is:

$$\mathit{\sigma }\propto \dot{\mathit{\epsilon }}$$

(98)

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 4a. Figure 4b 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 stress that does not reach the yield stress, materials will not exhibit stress relaxation behavior.

Hyde et al[18] used the 316 Stainless Steel hollow circular cross-section rod specimen shown in Figure 3a to carry out the low cycle fatigue test (LCF test) at 550 °C. The test is controlled by strain, the loading waveform is triangular wave, the strain amplitude is $0.006$, the strain ratio is $-1$, and the period is $80s$. Numerical simulations also use the finite element model shown in Figure 3b and c, setting the same displacement boundary conditions in the z-direction as the LCF test. As shown in Figure 5a, the stress-strain loops predicted by the model for the first and fiftieth cycles are in good agreement with the test:

Figure 5b 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 Equation 11, 12 and 13.

4.2. Zircaloy-4

Delobelle et al[28] used the Zircaloy-4 circular tube specimen shown in Figure 6a 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 displacement-controlled 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\times 1\times 1$ cube, as shown in Figure 6b. 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 6c. 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 the 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 7 shows the stress-strain curves of R state and CWSR state Zircaloy-4 for different strain rates at 350°C:

The numerical simulation results are in good agreement with the experimental results, which verifies the correctness of Zircaloy-4 modeling and UMAT code. According to Equation 7 and 26, we have:

$${\sigma }_v=N{sinh}^{-1}\left[{\left(\frac{\dot{p}}{{\dot{\epsilon }}_0}\right)}^{1/k}\right]$$

(99)

Similar to 316 Stainless Steel, the stress of Zircaloy-4 is positively correlated with the strain rate (Equation 98), 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 7.

Based on the test device in Figure 6a, Delobelle et al[28] 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 specimen 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 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\times 1\times 1$ cube. As shown in Figure 6d, 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 6e, 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 6c. The creep stresses of Zircaloy-4 in the R state are $125MPa$, $135MPa$, $150MPa$, and $170MPa$, respectively, and the creep stresses of Zircaloy-4 in the CWSR state are $170MPa$ and $250MPa$, and each creep stress level is maintained for 200 hours and then increased to the next stress level. Figure 8 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.

5. Conclusions

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. 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 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 various loading conditions of the cladding materials in service. We hope our efforts provide useful material models and available source codes for users from academia or industry to carry out more accurate mechanical analysis of cladding structures.