Nonlinear dynamic analysis of nuclear fuels considering thermo–elastic–plastic–creep phenomena

Abstract

Recently, research on the behavior of nuclear fuel rods in accident situations was performed. The research conducted was very sophisticated and involved complex tasks related to thermal, elastic, plastic, and creep phenomena. Previously, a nonlinear iterative static analysis for the behavior of a nuclear fuel rod was performed taking into consideration the thermal, elastic, plastic, and creep effects. However, the analysis cannot be applied to an unstable situation such as a bulging phenomenon of a nuclear fuel rod when a loss of coolant accident occurs. In this study, a nonlinear iterative dynamic analysis is performed considering large strain and thermo-mechanical effects, and the results of this dynamic analysis were obtained. The analysis is similar to those of static analysis; however, it differs from the static analysis after loss of coolant accident because of the inertia term of the dynamic equation.

Introduction

Energy consumption of human beings is drastically increasing till date and will increase further more in the future. The fossil energy will cease to serve sooner or later, and the alternative energy is still far less than the necessity. As the importance of nuclear energy is expected to be increased more and more, the safe control of a nuclear power system including more precise prediction of the accidental behavior of nuclear fuels is required because of the formidable situation that would be encountered in case of an accident, which might happen without warning. Recently, researches on the behavior of nuclear fuel rods in accidents are being performed; the research conducted are very sophisticated and involved complex task related to thermal, elastic, plastic, and creep phenomena. Moreover, the strain may not be small; thus, we need a large strain plasticity analysis.

Thermo–elastic–plastic–creep analysis was fully developed using an equivalent stress function algorithm for small strain problems. ¹ Large strain plasticity analyses were studied for isotropic and anisotropic problems.^2,3 Thermo–elastic–plastic–creep analysis considering large strain plasticity has been applied to nuclear fuel in loss of coolant accident (LOCA) situation. ⁴ However, the quasi-static analysis could not handle the bulging phenomenon at the final stage of deformation.

The bulging phenomenon is a type of stability problem, and neither the quasi-static approach nor the post-buckling approach can be used to deal with this phenomenon, as these approaches are not suitable for problems where properties depend upon time and temperature. The only appropriate approach is a dynamic analysis considering the thermo–elastic–plastic–creep effects. This dynamic analysis will be developed and applied to the LOCA situation to show the thermo-mechanical behavior of fuels at the final stage of deformation. The results will show the exact time when the nuclear fuel rods will fail in a LOCA situation after the bulging phenomenon, which is a critical factor in the design and recovery of nuclear power plants.

Finite element analysis of nuclear fuels including the bulging phenomena

For the analysis of the bulging phenomenon of a nuclear fuel rod, the heat transfer equation and the dynamic equation of motion will be presented. However, when the inertia term is neglected, the dynamic equation may be reduced to an equilibrium equation.

Nonlinear transient heat transfer equation

For the thermo–elastic–plastic–creep analysis of a nuclear fuel rod, it is necessary to obtain time-dependent temperatures by solving the following well-known transient equation of heat transfer

$$k({\mathrm{T}}_{,\mathrm{xx}}+{\mathrm{T}}_{,\mathrm{yy}})+\mathrm{Q}-c\rho \stackrel{·}{\mathrm{T}}=0$$

$$\mathrm{q}={\mathrm{f}}_B\mathrm{on}{\Gamma }_f$$

$$\mathrm{q}=h(\mathrm{T}-{\mathrm{T}}_{\infty })\mathrm{on}{\Gamma }_h$$

(1)

where T is the temperature, k is the coefficient of the thermal conductivity, ρ is the density, c is the specific heat, q is the heat flow rate, f _B is the heat flux, and h is the convective heat transfer coefficient. The partial differential equation of heat transfer (equation (1)) is expressed in other functional forms to be minimized as in equation (2). ⁵

Consider

$${\displaystyle \prod ={\displaystyle \int \left(\frac{1}{2}{T}_{\partial }^T\kappa T_{\partial }-QT+c\rho T\stackrel{.}{T}\right)}}dV+{\displaystyle \int \left(f_{B}T+h{T}_{\infty }T-\frac{1}{2}h{T}^2\right)dS}$$

(2)

where $\kappa$ is the conductivity matrix and T_∂ is the temperature gradient. After we introduce the finite elements with proper shape functions and minimize the functional parameters, the following finite element equation is derived

$$\begin{array}{c}{{}^{t+\Delta t}\mathbf{C}}^{(i)t+\Delta t}{\stackrel{·}{\mathbf{T}}}^{(i-1)}+({{}^t\mathbf{K}}^K+{{}^t\mathbf{K}}^C)\Delta {\mathbf{T}}^{(i)}\hfill \\ ={}^{t+\Delta t}\mathbf{Q}+{{}^{t+\Delta t}\mathbf{Q}}^{C(i-1)}-{{}^{t+\Delta t}\mathbf{Q}}^{K(i-1)}\hfill \end{array}$$

(3)

$${{}^{t+\Delta t}\mathbf{T}}^{(i)}={{}^{t+\Delta t}\mathbf{T}}^{(i-1)}+\Delta {\mathbf{T}}^{(i)}$$

(4)

Superscript i indicates the quantity of iterative state i at time $t+\Delta t$ ; therefore, ${}^{t+\Delta t}C$ is the heat capacity matrix at time $t+\Delta t$ , ${}^{t+\Delta t}\stackrel{·}{\mathbf{T}}$ is the rate of temperature at time $t+\Delta t$ , ${{}^t\mathbf{K}}^K$ is the conductivity matrix, and ${{}^t\mathbf{K}}^C$ is the convection matrix at time t. $\Delta \mathbf{T}$ is the incremental nodal point temperature, ${}^{t+\Delta t}\mathbf{Q}$ is the nodal point heat input vector, ${{}^{t+\Delta t}\mathbf{Q}}^C$ is the nodal point heat contribution due to convection, and ${{}^{t+\Delta t}\mathbf{Q}}^K$ is the equivalent nodal point heat contribution. Accordingly, to solve the unknown nodal point temperature ${}^{t+\Delta t}\mathbf{T}$ , the equilibrium iteration will be performed repeatedly until the difference of the right-hand terms of equation (3) is smaller than the convergence tolerance.

Iterative nonlinear equation of motion

After obtaining the temperature distributions of the fuel, the incremental equilibrium equation for the motion of the fuel can be solved to obtain the incremental displacements. Generally, the method for analyzing nonlinear problems includes the total Lagrangian and updated Lagrangian formulation. The updated Lagrangian formulation, which calculates the state at time $t+\Delta t$ referring to the configuration at time t, is the method for nonlinear problems using the principle of the virtual work and the Cauchy stress tensor.

The principle of virtual work at time t may be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{{}^{t}R}& ={\int }_{0V}{}_0{}^t{S}_{ij}{\delta }_0^t{\epsilon }_{\mathrm{ij}}{d}^{0}V\hfill \\ \multicolumn{1}{c}{}& ={\int }_{0V}{}_0{}^t{C}_{ijrs}({}_0{}^t{e}_{rs}+{}_0{}^t{\eta }_{rs}){\delta }_0^t{e}_{\mathrm{ij}}{d}^{0}V\hfill \\ \multicolumn{1}{c}{}& +{\int }_{0V}{}_0{}^t{S}_{ij}\delta {\eta }_{\mathrm{ij}}{d}^{0}V\hfill \end{array}$$

(5)

where $C_{\mathrm{ijrs}}$ is the material property tensor of stress–strain.

The relation between the external virtual work ${}^{t}R$ and the external force vector ${}^t\mathbf{P}$ is as follows ⁶

$$\begin{array}{cc}\multicolumn{1}{c}{{}^{t}R}& ={\delta }^t{\mathbf{U}}^T{}^t\mathbf{P}\hfill \\ \multicolumn{1}{c}{}& ={\int }_{0V}{\delta }^t{\mathbf{U}}^T{}_0{}^t{\mathbf{B}}_{L0}{}^t\mathbf{D}{\mathbf{B}}_L{}^t\mathbf{U}{d}^{0}V\hfill \\ \multicolumn{1}{c}{}& +{\int }_{0V}{\delta }^t{\mathbf{U}}^T{}_0{}^t{\mathbf{B}}_{NL}^T{}_0{}^t{\mathbf{S}}_0^t{\mathbf{B}}_{\mathrm{NL}}{}^t\mathbf{U}{d}^{0}V\hfill \end{array}$$

(6)

Applying the quotient law in order to obtain the equivalent nodal forces from equation (6)

$${}^t\mathbf{P}={\displaystyle {\int }_{{}^{0}V}{}_0{}^t{\mathbf{B}}_{L0}^T\mathbf{D}{\mathbf{B}}_L{}^t\mathbf{U}d{}^{0}V}+{\displaystyle {\int }_{{}^{0}V}{}_0{}^t{\mathbf{B}}_{NL}^T{}_0{}^t{\mathbf{S}}_0^t{\mathbf{B}}_{NL}{}^t\mathbf{U}d{}^{0}V}={}_0{}^t\mathbf{F}$$

(7)

where ${}_0{}^t\mathbf{F}$ is the equivalent nodal force vector corresponding to the internal stresses at time t. It is the equilibrium state when ${}^t\mathbf{P}={}_0{}^t\mathbf{F}$ .

Because it establishes the equilibrium state when the external force is equal to the equivalent nodal force at time $t+\Delta t$

$${}^{t+\Delta t}\mathbf{P}={}_0{}^{t+\Delta t}\mathbf{F}$$

(8)

Equation (8) can be written in the increment form as

$${}^{t+\Delta t}\mathbf{P}={}_0{}^{t+\Delta t}\mathbf{F}={}_0{}^t\mathbf{F}+{}_0\mathbf{F}$$

(9)

where ${}_0\mathbf{F}$ is the incremental nodal force vector, which can be obtained approximately from the tangent stiffness matrix as

$${}_0\mathbf{F}\approx {}_0{}^t\mathbf{K}\Delta \mathbf{U}$$

(10)

$\Delta \mathbf{U}$ denotes the increment of the nodal displacement vector ${}^{t+\Delta t}\mathbf{U}$ . Substituting equation (10) to equation (9)

$${}_0{}^t\mathbf{K}\Delta \mathbf{U}={}^{t+\Delta t}\mathbf{P}-{}_0{}^t\mathbf{F}$$

(11)

is obtained, and because the inertia force is added to body force, in dynamic problems, equation (11) can be written as

$${}_0{}^t\mathbf{K}\Delta \mathbf{U}={}^{t+\Delta t}\mathbf{P}-{}_0{}^t\mathbf{F}-{\mathbf{M}}^{t+\Delta t}\stackrel{··}{\mathbf{U}}$$

(12)

where $\mathbf{M}={\int }_{0V}{}_0{}^{\rho }{\mathbf{H}}^T\mathbf{H}{d}^{0}V$ and $\mathbf{H}$ is the shape function matrix. The stiffness matrix in equation (12) can be decomposed into ${}_0{}^t\mathbf{K}={}_0{}^t{\mathbf{K}}_L+{}_0{}^t{\mathbf{K}}_{NL}$ , where the latter term of stiffness is known as geometric nonlinear stiffness with regard to large displacement, buckling, and instability behaviors

$$\mathbf{K}={\mathbf{K}}_L+{\mathbf{K}}_{\mathrm{NL}}={\int }_{0V}{({\mathbf{B}}_{L0}+{\mathbf{B}}_{L1})}^{T}D({\mathbf{B}}_{L0}+{\mathbf{B}}_{L1})d^{0}V$$

(13)

As equations (11) and (12) are linearized equations, they have some errors in comparison to the exact solutions. Therefore, iteration is required to correct the errors. Equation (12) can be expressed with the equilibrium by the Newton–Raphson iteration and Newmark method as follows ⁶

$$\begin{array}{l}{}_0{}^{t+\Delta t}{\mathbf{K}}^{(i-1)}\Delta {\mathbf{U}}^{(i)}={}^{t+\Delta t}\mathbf{P}-{}^{t+\Delta t}{\mathbf{F}}^{(i-1)}\\ -\mathbf{M}\left\{\frac{4}{{(\Delta t)}^2}\left({}^{t+\Delta t}{\mathbf{U}}^{(i-1)}-{}^t\mathbf{U}\right)-\frac{4}{\Delta t}{}^t\stackrel{.}{\mathbf{U}}-{}^t\stackrel{··}{\mathbf{U}}\right\}\end{array}$$

(14)

Superscript i indicates the quantity of iterative state i. In order to simplify the displacement, equation (14) can be written as equation (15)

$${}_0{}^{t+\Delta t}{\mathbf{K}}^{(i-1)}\Delta {\mathbf{U}}^{(i)}={}^{t+\Delta t}\mathbf{P}-{}^{t+\Delta t}{\mathbf{F}}^{(i-1)}-{\mathbf{M}}^{t+\Delta t}{{}_0\stackrel{··}{\mathbf{U}}}^{(i)}$$

(15)

As an iterative equilibrium equation which combines the Newton–Raphson iteration method and the Newmark integration method is used in the numerical calculation, equation (14) analyzes the nonlinear dynamic problems in this article. Ignoring the mass term in equation (15), the nonlinear incremental equilibrium equation for the quasi-static motion of the fuel can be obtained.

Thermo–elastic–plastic–creep analysis

The deviatoric stress is computed from the deviatoric strain while regarding the plastic and creep strain^1,4 adopting the same symbols in Kojic and Bathe ¹ as

$${}^{t+\Delta t}\mathbf{S}=\frac{{}^{t+\Delta t}E}{1+{}^{t+\Delta t}\mathbf{\nu }}({}^{t+\Delta t}\mathbf{e}\prime -{}^{t+\Delta t}{\mathbf{e}}^P-{}^{t+\Delta t}{\mathbf{e}}^C)$$

(16)

The mean stress is computed from the mechanical mean strain as equation (17)

$${}^{t+\Delta t}{\sigma }_m=\frac{{}^{t+\Delta t}E}{1-2^{t+\Delta t}\nu }({}^{t+\Delta t}{e}_m-{}^{t+\Delta t}{e}_{\mathrm{th}})$$

(17)

where ${}^{t+\Delta t}\mathbf{S}$ is the deviatoric stress tensor, ${}^{t+\Delta t}\mathbf{e}\prime$ is the deviatoric strain tensor, ${}^{t+\Delta t}{\mathbf{e}}^p$ is the plastic strain tensor, ${}^{t+\Delta t}{\mathbf{e}}^c$ is the creep strain tensor, ${}^{t+\Delta t}{\sigma }_m$ is the mean stress, ${}^{t+\Delta t}{e}_m$ is the mean strain, ${}^{t+\Delta t}{e}_{\mathrm{th}}$ is the thermal strain, ${}^{t+\Delta t}\mathbf{E}$ is the Young’s modulus and ${}^{t+\Delta t}\nu$ is the Poisson’s ratio at time $t+\Delta t$ .

If we split the plastic and creep stains at $t+\Delta t$ into the terms at time t and incremental terms, equation (17) then becomes

$${}^{t+\Delta t}\mathbf{S}=\frac{{}^{t+\Delta t}E}{1+{}^{t+\Delta t}\nu }({}^{t+\Delta t}\mathbf{e}″-\Delta {\mathbf{e}}^p-\Delta {\mathbf{e}}^c)$$

(18)

where

$${}^{t+\Delta t}\mathbf{e}″={}^{t+\Delta t}\mathbf{e}\prime -{}^t{\mathbf{e}}^p-{}^t{\mathbf{e}}^c$$

(19)

where ${}^t{\mathbf{e}}^p$ and ${}^t{\mathbf{e}}^c$ are the plastic and creep strain at the start of the current time step.

Creep strain may be represented by various creep laws such as the power creep law or the exponential creep law. To confirm plastic stress based on the given displacement, ${{}^{t+\Delta t}u}^{(i)}$ , it must satisfy one of the yield criteria. Adopting the von Mises yield criterion

$$f_y={{}^{t+\Delta t}\overline{\sigma }}^2-{{}^{t+\Delta t}\sigma }_y^2=0$$

(20)

where

$${}^{t+\Delta t}\overline{\sigma }=\sqrt{1.5{}^{t+\Delta t}\mathbf{S}\cdot {}^{t+\Delta t}\mathbf{S}}:\mathrm{von\; Mises\; equivalent\; stress}$$

(21)

$${{}^{t+\Delta t}\sigma }_y=\mathrm{yield\; stress}$$

(22)

We have to solve the equation to obtain the deviatoric stress tensor. By applying following flow rules of plastic and creep deformations

$$\Delta {\mathbf{e}}^p=\Delta \lambda {}^{t+\Delta t}\mathbf{S}$$

(23)

$$\Delta {\mathbf{e}}^c=\Delta t^{\tau }{\gamma }^{\tau }\mathbf{S}$$

(24)

$$\tau =(1-\alpha )t+\alpha (t+\Delta t):\mathrm{relaxed\; time}$$

(25)

where $\alpha$ is the relaxation parameter.

We obtain

$${}^{t+\Delta t}\mathbf{S}=\frac{1}{{}^{t+\Delta t}{a}_E+\alpha \Delta t{}^{\tau }\gamma +\Delta \lambda }\left[{}^{t+\Delta t}({\mathbf{e}}^{″}-(1-\alpha )\Delta t{}^{\tau }\gamma {}^t\mathbf{S})\right]$$

(26)

where

$$a_E=\frac{1+\nu }{E}$$

(27)

Taking the product of both sides of the equation, we obtain the following equation

$$f_g=a^{2t+\Delta t}{\overline{\sigma }}^2+b^{\tau }\gamma -c^{2\tau }{\gamma }^2-d^2=0$$

(28)

where $a={}^{t+\Delta t}{a}_E+\alpha \Delta t^{\tau }\gamma +\Delta \lambda =(1/2G\prime )$ , $b=3(1-\alpha )\Delta t^{t+\Delta t}\mathbf{e}″·{}^t\mathbf{S}$ , $c=(1-\alpha )\Delta t^t\overline{\sigma }$ , and $d^2=1.5^{t+\Delta t}\mathbf{e}″·{}^{t+\Delta t}\mathbf{e}″$ .

The function on the left side of equation (28) is known as the effective-stress function, and it has to satisfy the condition of zero. Because $\overline{\sigma }$ is a function of the plastic strain, and $\Delta \lambda$ is also a function of the plastic strain and yield stress; equation (28) must be solved simultaneously for $\overline{\sigma }$ and $\Delta \lambda$ using a root finding algorithm. Both the Newton method and the bisection method were used to find the proper roots. The latter method is very reliable without being affected by the derivative of the function. For the assumed ${}^{t+\Delta t}\overline{\sigma }(={{}^{t+\Delta t}\sigma }_y)$ , $\Delta \lambda$ was computed first and then checked whether the equation is satisfied. If the equation cannot be satisfied, then the usual root finding algorithm can be used to find ${}^{t+\Delta t}\overline{\sigma }(={{}^{t+\Delta t}\sigma }_y)$ and $\Delta \lambda$ simultaneously such that the chosen type of yield criterion is satisfied. The tangent constitutive matrix ${\mathbf{C}}^{\mathrm{EPC}}$ can be obtained using the perturbation method.^1,7

Here, we propose a first stage constitutive modulus as

$$C″=\frac{1}{{}^{t+\Delta t}{a}_E+\lambda +\Delta t^{\tau }\gamma }=2G″$$

(29)

where $C″$ is approximately closer to the tangent modulus than $C\prime$ . We could determine fast convergence with $C″$ than $C\prime$ in several examples.

Large strain (Hencky logarithmic) measure

For large strain analysis, some basic relations are required which are briefly introduced here.^2,3,8–11 The deformation gradient is represented by

$${}_0{}^t\mathbf{X}=\frac{{\partial }^t\mathbf{x}}{{\partial }^0\mathbf{x}}={}^t{\mathbf{J}}^T({}^0{\mathbf{J}}^{-1}{)}^T$$

(30)

where $\mathbf{J}$ is the transformation between the Cartesian coordinate and natural coordinate systems. It can be decomposed into rotation tensor ${}_0{}^t\mathbf{R}$ and right stretch tensor ${}_0{}^t\mathbf{U}$

$${}_0{}^t\mathbf{X}={}_0{}^t\mathbf{R}{}_0{}^t\mathbf{U}$$

(31)

The right Cauchy–Green gradients are expressed using the deformation gradient

$${}_0{}^t\mathbf{C}={}_0{}^t{\mathbf{X}}^T{}_0{}^t\mathbf{X}$$

(32)

Solving the following equivalent problem

$${}_0{}^t\mathbf{C}{}_0{}^t\mathbf{P}={}_0{}^t{\mathbf{\Lambda }}^2{}_0{}^t\mathbf{P}$$

(33)

A diagonal matrix ${}_0{}^t\mathbf{\Lambda }$ of principal stretches and their direction matrix is then obtained. The right stretch matrix is

$${}_0{}^t\mathbf{U}={}_0{}^t\mathbf{P}{}_0{}^t\mathbf{\Lambda }{}_0{}^t{\mathbf{P}}^T$$

(34)

Therefore, the rotation matrix is given by

$${}_0{}^t\mathbf{R}={}_0{}^t\mathbf{X}{}_0{}^t{\mathbf{U}}^{-1}$$

(35)

Then, the Hencky logarithmic strain tensor is given by

$${}_0{}^t{\mathbf{E}}_H=\frac{1}{2}{}_0{}^t\mathbf{P}{}_0{}^t\mathbf{H}{}_0{}^t{\mathbf{P}}^T$$

(36)

where

$${}_0{}^t\mathbf{H}=\ln {\mathbf{\Lambda }}^2$$

(37)

The Green–Lagrange strain tensor can be represented as

$${}_0{}^t{\mathbf{E}}_G=\frac{1}{2}({}_0{}^t\mathbf{C}-\mathbf{I})$$

(38)

and the Almansi strain tensor as

$${}_0{}^t{\mathbf{E}}_A=\frac{1}{2}(\mathbf{I}-{}_0{}^t\mathbf{C})$$

(39)

The Hencky logarithmic strain tensor can also be written as

$${}_0{}^t{\mathbf{E}}_H=\ln {}_0{}^t\mathbf{U}=\frac{1}{2}\ln {}_0{}^t\mathbf{C}$$

(40)

Using the equations above and comparing them (three kinds of marks) with theoretical equations (equations (41)–(44)) in Figure 1 then

$${}_0{}^t{e}_{11}=\frac{1}{2}({}_0{}^t{\lambda }^2-1):\mathrm{Green}-\mathrm{Lagrange\; strain}$$

(41)

$${}^t{e}_{11}={}_0{}^t\lambda :\mathrm{small\; strain}$$

(42)

$${}_0{}^t{E}_{11}^{(H)}=\ln {}_0{}^t\lambda :\mathrm{logarithmic\; strain}$$

(43)

$${}_t{}^t{e}_{11}^A=\frac{1}{2}(1-{}_0{}^t{\lambda }^{-1}):\mathrm{Almansi\; strain}$$

(44)

Figure 1.

Strain–stretch relations for various strain measures.

Using the procedure of thermo–elastic–plastic–creep with logarithmic strain measure and additive decomposition of plastic strain, the developed Fortran program with additive decomposition method follows the procedures in Bathe and colleagues,^2,8,9 except that plastic deformation gradient is not decomposed multiplicatively as

$${\mathbf{X}}_A^E={}_0{}^{t+\Delta t}\mathbf{X}{\left({}_0{}^{t+\Delta t}{\mathbf{X}}^p\right)}^{-1}$$

(45)

at the trial elastic state first, but considered additively as in equations (16) and (18). Therefore, continuation between later large strain procedure and former small strain procedure including creep deformation is not lost.

Verification of developed program

For the analysis of the bulging phenomenon of a nuclear fuel rod, analysis program is developed. The developed FORTRAN program was tested with several cases.

Thermal analysis

Thermal analysis ⁴ was performed using the present program and commercial package ADINA. The boundary conditions are shown in Figure 2.

Figure 2.

Boundary condition for thermal analysis.

Temperatures are analyzed and the results of three nodes are compared in Table 1 showing a good agreement.

Table 1.

Results of thermal analysis.

	Temperature (°K)
	Node 9	Node 11	Node 13
Present program	366.6	354.3	331.93
ADINA	366.6	354.3	331.93

Plasticity analysis

Plasticity analysis ⁴ was performed using the present program and commercial package ADINA. The boundary conditions are shown in Figure 3.

Figure 3.

Boundary condition for plasticity analysis.

The von Mises, radial, and hoop stresses are analyzed for three cases of perfect, bilinear, and general plasticities. The initial value of the yield stress is 25 MPa, and the mechanical pressure load increases up to 9 MPa with 90 steps. As shown in Figure 4, the present program is proven to have a verified accuracy for the plasticity problems.

Figure 4.

Results of plasticity analysis: (a) perfect plasticity, (b) bilinear plasticity, and (c) general plasticity.

Creep analysis

Creep analysis was performed using the present program and commercial package ADINA. ⁴ The boundary conditions are shown in Figure 5.

Figure 5.

Boundary condition for creep analysis.

The von Mises, radial, and hoop stresses of the present program and ADINA are compared each other for creep problems as shown in Figure 6. The present program is reliable in the analysis of elasto-creep as well as in the thermal and plasticity analyses.

Figure 6.

Results of creep analysis.

Large strain plasticity

For the case without creep deformation, the following example compares the two procedures of large strain problems: additive decomposition (developed present program) and multiplicative decomposition. ² The first one is uniaxial isochoric extension with the following conditions

$$\begin{array}{c}\mathrm{Young}\prime \mathrm{s\; modulus}:E=200\mathrm{MPa}\hfill \\ \mathrm{Shear\; modulus}:\mu =76.92\mathrm{MPa}\hfill \\ \mathrm{Bulk\; modulus}:\kappa =166.7\mathrm{MPa}\hfill \\ \mathrm{Poisson}\prime \mathrm{s\; ratio}:\nu =0.3\hfill \\ \mathrm{Yield\; stress}:{\sigma }_y=0.75\mathrm{MPa}\hfill \\ \mathrm{Hardening\; modulus}:H=2\mathrm{MPa}\hfill \end{array}$$

An in-plane isochoric deformation with no rotation of the principal stretch directions, given by the deformation gradient

$$[X]=\left[\begin{array}{ccc}\hfill \lambda \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 1/\lambda \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill 1\hfill \end{array}\right]$$

(46)

where λ ∈ [1, 1.5].

The results are compared in Figure 7 showing a good agreement.

Figure 7.

Normal stress σ_x along stretch λ.

The second example is shear deformation given by the deformation gradient

$$[X]=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \gamma \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 1\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill 1\hfill \end{array}\right]$$

(47)

where γ ∈ [0, 0.5].

The result compared to the literature results (solid square) obtained by the multiplicative decomposition ² in Figure 8.

Figure 8.

Shear stress τ₁₂ along shear strain γ.

Behavior of nuclear fuel rods after LOCAs considering the bulging effect

After a LOCA during the behavior of nuclear fuels, an unstable situation develops where the cladding stiffness becomes zero because of a local uncontrollable deformation. ⁴ In this case, an analysis using the post-buckling method applied to quasi-static problems is not a proper method for analyzing this transient problem because properties of the materials, which composed the nuclear fuel rods, are not constant. The material properties are considered functions of time and temperature. In order to analyze this type of phenomenon without divergence, the dynamic analysis method was adopted. Consequently, the bulging motion of the inflated nuclear fuel rods during a LOCA situation can be analyzed using the dynamic equation.

In the usual case of LOCA, the control rods are inserted into the fuel rods to shut down the reactor, and the power is reduced successfully. However, during events, such as the Fukushima accident, the control rods could not be inserted into the fuel rods because the insertion system, which is located at the lower part of the reactor, is malfunctioned. Therefore, it was assumed that the heat generation of the fuels was partly reduced.

The material properties of the fuel components are the same as in the literature.^12,13 Mesh generation of the fuel rod was conducted using two-dimensional (2D) axisymmetric and quadrilateral eight-node elements, and the boundary conditions were applied as shown in Figure 9. The boundary conditions include the application of heat generation to the pellet as well as coolant pressure and convection heat transfer from the cladding to coolant. In addition, the gas pressure of the gap between the pellet and the cladding is also applied. The actual nuclear fuel rods are activated with mechanical and thermal properties shown in Table 2. In this article, the gas pressure of 34.5 MPa is applied to the analyses to assume the severe situation of the nuclear fuel rods. Figure 10 shows the temperature changes at the pellet center and cladding outer surface at steady-state conditions prior to 100 s.

Figure 9.

(a) Solid nuclear fuel model and (b) boundary conditions.

Table 2.

Mechanical and thermal properties.

	Pellet	Gap	Cladding	Coolant
T ∞	–	–	600°K	–
Pressure (MPa)	–	3.45	–	15.5
Density (mg/mm3)	1.0412e−08	1.0412e−08	6.5500e−09	–
Heat generation (mJ/mm3 s)	459	–	–	–
Convection (mJ/s mm2 °K)	–	–	100	–
Specific heat (mJ/mg °K)	0.313889e+09	0.313889e+07	0.325896e+09	–
Thermal conductivity (mJ/s mm °K)	0.344808e+01	0.5	0.129064e+02	–

Figure 10.

Temperatures changes at pellet center and cladding outer surface along time. ⁴

The LOCA event was assumed to occur at 70 s resulting in a loss of coolant pressure and change in the ambient heat transfer coefficient from 7500 to 75 mJ/s mm² °K. In addition, the control rod is also malfunctioned. It was assumed that one of four conditions of heat generation (complete shutdown = 0%, 55%, 60% and full heat generation = 100%) was maintained after LOCA. Temperature changes at the cladding outer surface are shown in Figure 11. As expected, there were rapid changes in temperature at both the outer surface of the cladding because of the shutdown. It was detected that the outer surface of the cladding shows a rapid increase in temperature and then the temperature increases slowly afterward at both 60% and 100% while temperature of the reactor continued to decrease at complete shutdown after 70 s. After LOCA, there were sharp increases in temperature especially for full heat generation. The cladding withstood the internal pressure load from the fuel rod. However, Young’s modulus and strength of the clad decreased since the cladding temperature beyond the steady-state temperature condition of 600 °K. As a result, the fuel rod cannot resist the pressure load and deformation begins to occur, which causes larger stresses leading to fuel failure. This phenomenon is similar to buckling instability. The asterisk marks denote the divergence time at each heat generation in the quasi-static analysis of the fuel rods.

Figure 11.

Temperatures at cladding outer surface with the LOCA at 70 s.

Figure 12 shows the radial displacements of the cladding outer surface during operation with considering dynamic effects. The heat generation above 55% reactor power caused cladding instability, which then induced continuous expansion of the fuel rod ultimately leading to burst fuel rods. However, adopting the dynamic equation, the analyses was not been stopped due to the instability conditions. The variations of hoop stresses at the cladding surface are shown in Figure 13.

Figure 12.

Radial displacements of cladding outer surface with the LOCA at 70 s.

Figure 13.

Hoop stresses at cladding outer surface with the LOCA at 70 s.

In the case of 60% and 100% heat generation after a LOCA, the behavior of the fuel rods could be analyzed using the dynamic equation without considering the cladding instability. This approach differed from Kim et al. ⁴ where the analyses do not consider the dynamic effects.

Conclusion

For the safe design and operation of a nuclear power plant, a precise thermo-mechanical analysis regarding the behavior of nuclear fuel rods in accidental situations is required as well as the analysis for controlled operations. Thermo-mechanical analysis, including the thermal, elastic, plastic (large strain), and creep effect, on nuclear fuel comprehensively is rare. Moreover, the bulging phenomenon, which is the last stage of deformation followed by the bursting failure, cannot be analyzed by a conventional quasi-static approach. We presented a nonlinear dynamic analysis of nuclear fuel rods considering the thermal, elastic, plastic (large strain), and creep effects for analyzing the bulging phenomenon in LOCA cases. It can be summarized as follows:

• A comprehensive nonlinear iterative static analysis of nuclear fuel rods was made including the thermal, elastic, large strain plastic, and creep effects.

• A nonlinear iterative dynamic analysis procedure was formulated including the thermo-mechanical effects to handle the bulging phenomenon in LOCA situation.

• The results of the dynamic analysis are very close to those of the static approach before the bulging. During the bulging phenomenon, dynamic results could be obtained while static analysis failed to proceed beyond the critical point.