11.3. Metal Fuel Element Models

Although the FPIN2 code was originally developed for the purpose of analyzing oxide fuels, it has been modified for the analysis of metallic fuels. The metal fuel version of the code includes numerous improvements to the basic mechanical analysis models that reflect the experience gained from the in-pile TREAT tests on EBR-II irradiated prototypic fuel of the IFR concept. In addition, models for molten cavity formation, the large gas plenum, molten fuel extrusion, and fuel-cladding eutectic formation are provided to complement the fuel element mechanics calculation.

There are many differences between metal and oxide fuels that affect the transient response of the fuel elements. The existence of a large fission gas plenum in metal-alloy fueled elements plays a major role in determining the internal pin pressure during a transient. The sodium bond and comparatively large radial fuel-cladding gap requires consideration of expelling sodium into the plenum as the fuel expands and the gap closes. Metal fuel thermal conductivity is an order of magnitude larger than oxide fuel, and that leads to a flatter radial temperature profile and a more rapid spreading of the region of molten fuel. Since the metal fuel solidus temperature is well below the cladding melting temperature, in many cases fuel can be expected to melt entirely before the cladding fails. Due to the high fuel thermal conductivity, melting usually begins at or near the top of the fuel column with the axial profile of the fuel centerline temperature more closely following the coolant temperature profile than that for oxide fuel. This prevents large cavity pressurization effects in metallic fuels. In addition, a low melting temperature eutectic alloy is expected to form in metallic fuels between the fuel and cladding at temperatures below the anticipated cladding failure. This eutectic formation allows the fuel and cladding to slip freely in the axial direction and can lead to an accelerated cladding failure at elevated temperatures. Cracking is not a significant phenomenon in metallic fuel pins, but metal fuel swelling can be important.

The current version of the FPIN2 coupled with SASSYS/SAS4A includes various model modifications and additions that address these differences between oxide and metal fuels. The following sections contain discussions of these additions and modifications made to FPIN2 for the analysis of metal-fueled elements. The most significant extension to the code is the inclusion of a model for the plenum region above the fuel column as discussed in Section 11.3.1. Extrusion of molten fuel into the plenum as an additional axial expansion mechanism is discussed in Section 11.3.2. Basic mechanical properties of metal-alloy fueled elements are outlined in Section 11.3.3. The constitutive equations for fuel creep and a non-equilibrium approximation to fuel swelling are presented in Section 11.3.4 and Section 11.3.5, respectively. The plastic flow behavior of three types of cladding materials commonly used in metallic fuel elements of the IFR concept is discussed in Section 11.3.6. Fuel-cladding eutectic formation and its impact on the mechanical analysis of fuel elements are summarized in Section 11.3.7. And finally, the fuel element failure formulation is outlined in Section 11.3.8.

11.3.1. Central Cavity and Plenum Models

Internal pin pressure is calculated in FPIN2 using straightforward volume accounting models of the central cavity and the gas plenum. The radial boundary of the cavity for each axial segment is assumed to be at the point where the fuel has reached the solidus temperature. Elements inside the cavity boundary are assumed to be in a hydrostatic state of stress equal to the cavity pressure and are dropped out of stress-strain calculation. These elements occupy a volume in proportion to their density and contribute to the volume and mass balance iteration that determines the cavity pressure.

Metallic fuel elements lead to relatively low transient induced internal pin pressure since (1) fuel porosity easily accommodates the fuel volume expansion upon melting, (2) fuel melting begins at or near the top of fuel column, and (3) metallic fuels have less restrictive flow path for fission gas passage to the plenum. Thus, the large gas plenum is very effective in mitigating any large pressure buildup inside the fuel elements. The pin plenum model of FPIN2 is based on the following assumptions:

  1. The internal pin pressure, \(p_{\text{g}}\), is considered uniform throughout the interior of a metallic fuel element.

  2. The plenum region is assumed to be at the uniform temperature. 1

  3. The thermal expansion of the fill sodium and its expulsion from the fuel-cladding gap as the gap closes is accounted for. However, sodium compressibility is assumed to be negligible at the expected pressure range.

The first level of iteration, or outer loop, in the mechanical analysis section of the FPIN2 is the search for the internal pin pressure \(p_{\text{g}}\). According to the calculation sequence, the temperatures of the cavity and plenum materials are known at the beginning of a new time step before the pin pressure is found. Using these temperatures, the following can be calculated:

  1. Cavity boundary changes due to fuel melting,

  1. Mass of the gas in the cavity, including the gas released upon melting,

  2. Mass-average temperature of the gas in the cavity,

  3. Volume available from porosity released after fuel melting,

  4. The thermal expansion of the sodium in the plenum.

The algorithm for calculating the plenum and cavity pressures are outlined below for the case where the two volumes are connected and share a common pressure, \(p_{\text{g}}\). The pin pressure, \(p_{\text{g}}\), is calculated by considering the two relationships between the gas volume and pressure that must be satisfied: the ideal gas law, and the interaction between pressure and cavity boundary displacements. Pin molten-cavity and plenum gases are considered as separate entities with individual compositions and temperatures. Since the masses of the gases in the cavity, \(m_{\text{g}}^{C}\), and the plenum, \(m_{\text{g}}^{P}\), as well as the temperatures of the cavity gas, \(T_{\text{g}}^{C}\), and the plenum gas, \(T_{\text{g}}^{P}\), are known, the ideal gas law reduces to an inverse relationship between the total gas volume, \(V_{\text{g}}^{\text{IGL}}\), and the pressure as follows

(11.3‑1)

\[V_{\text{g}}^{\text{IGL}} = \frac{m_{\text{g}}^{C} R^{C} T_{\text{g}}^{C} + m_{\text{g}}^{P} R^{P} T_{\text{g}}^{P}}{p_{\text{g}}}\]

where \(R^{C}\) and \(R^{P}\) are the gas constants for cavity and plenum gases, respectively. The second relationship involves the geometrical calculation of the volume available to the gas that depends on displacement of the solid fuel, cladding, and plenum tube. The liquid fuel and sodium are assumed to be incompressible. Since the gas in the cavity and plenum are considered separately, the total volume has two components

(11.3‑2)

\[V_{\text{g}}^{\text{MECH}}\left( p_{\text{g}} \right) = V_{\text{g}}^{C}\left( p_{\text{g}} \right) + V_{\text{g}}^{P}\left( p_{\text{g}} \right)\]

where \(V_{\text{g}}^{C}\) is the available cavity volume that is calculated using the finite element analysis of solid fuel and cladding, and \(V_{\text{g}}^{P}\) is the volume available for the gas in the plenum. A negative value of \(V_{\text{g}}^{C}\) means that the molten fuel material is extruded upward into the plenum.

The two values for the total volume, \(V_{\text{g}}^{\text{IGL}}\) and \(V_{\text{g}}^{\text{MECH}}\), are equal when the correct value of \(p_{\text{g}}\) is found. In FPIN2, Newton’s method is used to fine the value of \(p_{\text{g}}\) such that

(11.3‑3)

\[f\left( p_{\text{g}} \right) = V_{\text{g}}^{\text{IGL}} - V_{\text{g}}^{\text{MECH}} = 0\]

The Newton iteration equation

(11.3‑4)

\[p_{\text{g}}^{i + 1} = p_{\text{g}}^{i} - \frac{f\left( p_{\text{g}}^{i} \right)}{f'\left( p_{\text{g}}^{i} \right)}\]

requires the derivative of \(V_{\text{g}}^{\text{IGL}}\) and \(V_{\text{g}}^{\text{MECH}}\) with respect to \(p_{\text{g}}\) (\(i\) is the iteration counter). The derivative of \(V_{\text{g}}^{\text{IGL}}\) is easily calculated from the ideal gas law and the derivative of \(V_{\text{g}}^{\text{MECH}}\) is found by forming a finite difference quotient from two calculations of \(V_{\text{g}}^{\text{MECH}}\) at two neighboring values of \(p_{\text{g}}\).

The total mass of sodium, \(m_{\text{Na}}^{T}\), and mass of the gas in the plenum, \(m_{\text{g}}^{P}\), are obtained from the initial conditions. The inventory of sodium within the fuel element consists of sodium in the fuel-cladding gap, \(m_{\text{Na}}^{G}\), and in the plenum, \(m_{\text{Na}}^{P}\) :

(11.3‑5)

\[m_{\text{Na}}^{T} = m_{\text{Na}}^{P} + m_{\text{Na}}^{G}\]

where

(11.3‑6)

\[m_{\text{Na}}^{P} = \rho_{\text{Na}}\left( T^{P} \right) \pi r_{\text{P}}^{2} h_{\text{Na}}^{P}\]

(11.3‑7)

\[m_{\text{Na}}^{G} = \sum_{\text{j}} \rho_{\text{Na}}\left( T_{\text{j}} \right) \pi \left( r_{\text{ci}}^{2} - r_{\text{fo}}^{2} \right)_{\text{j}} \Delta z_{\text{j}}\]

The symbols in Eqs. 11.3-6 and 11.3-7 are defined as follows:

\(\rho_{\text{Na}}\) = sodium density,

\(T^{P}\) = plenum temperature,

\(r_{\text{p}}\) = cladding inner radius in plenum region,

\(h_{\text{NA}}^{P}\) = height of sodium column in plenum,

\(r_{\text{ci}}\) = cladding inner radius for axial segment j,

\(r_{\text{fo}}\) = fuel outer radius for axial segment j,

\(\Delta z_{\text{j}}\) = height of fuel axial segment j.

The mass of the plenum gas is determined by using the ideal gas law:

(11.3‑8)

\[m_{\text{g}}^{P} = \frac{P^{P} \pi r_{\text{P}}^{2}\left( h^{P} - h_{\text{Na}}^{P} \right)}{R^{P}T^{P}}\]

where \(P^{P}\) is the plenum pressure and \(h^{P}\) is initial plenum height. The plenum pressure is obtained from a user specified reference pressure at a reference temperature as follows

(11.3‑9)

\[P^{P} = P_{\text{ref}}\frac{T^{P}}{T_{\text{ref}}}\]

During the transient analysis, the cavity pressure algorithm requires the volume in the plenum available to the plenum fission gas and to the material extruded from the molten cavity. This value is calculated by subtracting the volume of the sodium in the plenum from the total volume of the deformed plenum tube. Thus,

(11.3‑10)

\[V_{\text{g}}^{P} = V^{P} - V_{\text{Na}}^{P}\]

The total plenum volume, \(V^{P}\), is found by assuming that the plenum tube can be treated as a thermoelastic thin shell. The volume of sodium in the plenum is found using the equation

(11.3‑11)

\[V_{\text{Na}}^{P} = \frac{m_{\text{Na}}^{T} - m_{\text{Na}}^{G}}{\rho_{\text{Na}}\left( T^{P} \right)}\]

11.3.2. Extrusion of Molten Fuel into the Plenum

The procedure in FPIN2 for calculating the behavior of the molten cavity materials uses a simple hydrostatic model based on volume accounting. As discussed in previous section, a volume and mass balance iteration determines the cavity pressure. For transients in which initial fuel melting is observed below the fuel-plenum interface, plenum and cavity pressure equations are solved separately, leading generally to a cavity pressure larger than the plenum pressure. Once melting reaches the top axial segment, the two pressures eventually equilibrate to give a common pin pressure. The temporary imbalance between the cavity and plenum pressures plays an important role in the dynamics of the molten fuel extruded into the plenum.

The extrusion of the molten cavity material into the plenum is calculated at the completion of the pressure iteration and consists of calculating the excess volume of the cavity materials over the volume available in the molten cavity. Both the radial and axial dimensions of the cavity section for each axial segment are based on the deformed geometry. The materials in the cavity include liquid fuel at a temperature above the liquidus temperature, solid fuel at a temperature between the solidus and liquidus, and the fission gas that is assumed to obey the ideal gas law and is released to the cavity in proportion to the melt fraction. Open porosity is available to the fission gas volume when the fuel reaches the solidus temperature, but closed grain boundary porosity is release d in proportion to the metal fraction.

For low-burnup fuel pins, the bond sodium in the fuel-cladding radial gap affects the extrusion of molten fuel into the plenum. The metallic fuel can melt completely across the fuel radial cross section in an axial segment. Although the large gap in fresh metallic fuels is not expected to close by the time of 100% fuel melting, the bond sodium in the gap is expected to escape into the plenum as the molten fuel slumps to fill up the gap due to the large density difference between fuel and sodium. This local slumping lessens the fuel extrusion; however, metal fuels at as low as 0.35 at.% burnup contain sufficient fission gas to result in net upward axial fuel motion [11-12]. Two options are provided in FPIN2 to model this behavior. The first option is to assume that the fuel does not slump and expel the bond sodium. In this case, the fuel outer boundary is assumed to remain at the position it was in at the time of 100% melting. The second (default) option allows the molten fuel to move out to the cladding, with the bond sodium being expelled into the plenum and the gap volume being added to the molten cavity.

To handle the second option, a gap closure model has been added to move the molten fuel out to the cladding over several time steps as the bond sodium is expelled into the plenum and the gap volume becomes available to the molten cavity materials. In order to maintain the stability of the numerical calculation, a varying closure rate is coded in FPIN2 in which the fraction of the gap closed per time step is specified. Somewhat artificial time step dependence of this transition is necessary because the temperature drop across the gap may be as large as 30K. Since the radial temperature profile across the fuel element is relatively flat, an excessively rapid reduction of the gap \(\Delta T\) may cause computational difficulties such as fuel refreezing during a power rise.

11.3.3. Basic Metal Fuel Properties

Thermal Properties: In the interfaced mode, the SAS-FPIN2 integrated model uses SASSYS/SAS4A routines to calculate fuel and cladding thermal material properties as discussed in Chapter 10. In the stand-alone mode, FPIN2 uses its own built-in correlations for metallic IFR fuel that have been developed by the “Metallic Fuel Properties Working Group” [11-12]. Melting related properties for the metallic fuels such as solidus and liquidus temperatures, heat of fusion, and volume change on melting are calculated using identical algorithms in FPIN 2 and SASSYS/SAS4A.

Elastic Properties: Tensile test data at room temperature for uranium and its alloys indicate yielding at very low stresses. The equation used in FPIN2 code for the Young’s modulus is as follows

(11.3‑12)

\[E = 0.12 \cdot 10^{6} \left( 1 - 1.2 p \right) \left( 1 - 0.754 \cdot 10^{-3} \left( T - 588 \right) \right)\]

where \(E\) is Young’s modulus (in Bars), \(p\) is fractional porosity, and \(T\) is temperature \(\left( K \right)\). This equation gives a value of one-tenth the handbook [11-7] value and should more properly be called a tangent modulus. However, for the expected monotonic loading of most FPIN2 calculations, the fuel is treated as a pseudo-elastic-plastic material with plastic straining handled as secondary creep and yielding is included as part of the “elastic” behavior. In the code, Poisson’s ratio is calculated from

(11.3‑13)

\[\nu = 0.27 \left( 1 - 0.8 p \right) \left( 1 + 0.854 \cdot 10^{- 3} \left( T - 588 \right) \right)\]

where \(\nu\) is Poisson’s ratio, \(p\) is fractional porosity, and \(T\) is temperature \(\left( K \right)\).

Linear Thermal Expansion: In the metal fuel version of FPIN2, the thermal expansion is expressed in terms of handbook linear thermal expansion data, \(\Delta L/L_{0}\) [11-7]. Using the data rather than the coefficient of thermal expansion, \(\alpha\), results in a more accurate accounting of the total expansion from room temperature to melting since it automatically includes the expansion at solid-to-solid phase transitions. The thermal expansion data and phase transition expansions have been approximated in the code as three straight line segments with breaks at the two solid-to-solid phase change temperature associated with the beginning and ending of the transformation to pure \(\gamma\) phase solid solution. The equations used for the binary fuel are

(11.3‑14)

\[\begin{split} \frac{\Delta L}{L_{0}} = \begin{cases} 1.695 \cdot 10^{- 5}\left( T - 293 \right) & T < 900 \\ 0.0103 + 7 \cdot 10^{- 5}\left( T - 900 \right) & 900 < T < 1000 \\ 0.0173 + 2.12 \cdot 10^{- 5}\left( T - 1000 \right) & T > 1000 \end{cases}\end{split}\]

And, the equations for the ternary fuel are

(11.3‑15)

\[\begin{split} \frac{\Delta L}{L_{0}} = \begin{cases} 1.67 \cdot 10^{-5} \left( T - 293 \right) & T < 864 \\ 0.0095 + 6.7 \cdot 10^{- 5}\left( T - 864 \right) & 864 < T < 950 \\ 0.0153 + 2.12 \cdot 10^{- 5}\left( T - 950 \right) & T > 950 \\ \end{cases}\end{split}\]

where and \(T\) is the temperature \(\left( K \right)\).

11.3.4. Secondary Creep of Metallic Fuels

Uranium alloys deform plastically under constant load at elevated temperatures. The secondary, or minimum, creep rate is defined as the steady-state rate that is attained under these conditions. The equations added to FPIN2 for metal fuels are obtained from a review of the data and theory for modeling secondary creep of U-Pu-Zr alloys that are of interest to the IFR concept [11-13]. Application of creep data to fuel pin analysis often requires correlating the data using mathematical functions of the governing variables. Such correlations range from purely empirical equations to theoretical equations involving fundamental physical properties. In FPIN2, and intermediate approach is preferred where theoretical models are used to obtain the form of the equation, but the parameters are determined form the creep data rather than the fundamental properties. Such an approach allows an interpolation of the data and gives reasonable confidence in often-required extrapolations beyond the database.

The particular form of the secondary creep equation used in FPIN2 is the form used to represent creep of UO2 [11-12]. The total plastic strain rate, \(\dot{\varepsilon}\), is given by

(11.3‑16)

\[\dot{\varepsilon} = C_{1}\left( \dot{F},d \right)\sigma \exp{\left( - \frac{Q}{RT} \right)} + C_{2}\left( \dot{F} \right) \sigma^{n} \exp{\left( - \frac{Q}{RT} \right)} + C_{3} \left( T \right) \sigma \dot{F}\]

where

\(\dot{\varepsilon}\) = secondary creep rate, s-1

\(\sigma\) = equivalent stress, MPa

\(R\) = universal gas constant, (1.987 cal/g-mole-K)

\(T\) = temperature, K

\(Q\) = activation energy, (cal/g-mole)

\(d\) = grain size, µm

\(\dot{F}\) = fission rate, (fissions/cm3-s)

\(C\) = material functions.

In Eq. 11.3-16, the first term represents diffusional creep, the second term dislocation creep, and the third fission-induced creep. Here, the low temperature deformation mechanisms are neglected such as dislocation glide and twinning that may occur in \(\alpha\) uranium at temperatures less than 675K. On the other hand, the expected dependence of in-reactor creep on the fission rate is considered.

Some of the parameters in Eq. 11.3-16 are estimated from knowledge of the structure of material and the values that these parameters take for similar materials. At low temperatures, it is reasonable to postulate that the creep deformation of uranium alloys is dominated by the deformation of the \(\alpha\) uranium matrix. This is particularly true if the other phases present are coarsely dispersed. When the \(\alpha\) uranium matrix deformation is dominant, the activation energy is expected to be the self-diffusion energy in \(\alpha\) uranium. The creep activation energy \(Q=52000\) cal/g-mole has been determined experimentally for \(\alpha\) uranium. It is shown that the temperature dependence of creep can be deduced using such a single activation energy and that the value of n=4.5 is appropriate in the high stress range [11-13]. Using these values for \(Q\) and \(n\), \(C_{1}\) and \(C_{2}\) in Eq. 11.3-16 are chosen to minimize the least squares error between the calculated and measures values of \(\log\left( z \right)\) where \(z\) is the Zener-Holloman parameter

(11.3‑17)

\[z = \dot{\varepsilon} e^{\frac{Q}{\text{RT}}}\]

The result is

(11.3‑18)

\[\begin{split}&C_{1} = 0.5 \cdot 10^{4} &\ & \text{MPa}^{- 1}\text{s}^{- 1} \\ &C_{2} = 6.0 &\ & \text{MPa}^{- 1}\text{s}^{- 1}\end{split}\]

When the calculated time to 2% strain using these values is compared to data from creep tests of U-Pu-Zr alloys, results are found to be in acceptable agreement. Examination of experimental data suggests that the grain-size dependence of \(C_{1}\) can be included as

(11.3‑19)

\[ \begin{align} C_{1} = \left( \frac{160}{d} \right)^{2} 0.5 \cdot 10^{4} && \text{MPa}^{- 1}\text{s}^{- 1} \end{align}\]

for grain sizes near 160 μm.

At temperatures above about 900K, the structure of the U-Pu-Zr alloys no longer contains an \(\alpha\) uranium matrix. The transformations between 900K and 975K are complex and depend strongly on composition. Beyond 975K, however, the solid phase consists of a solid solution of uranium, plutonium, and zirconium over a wide range of compositions. In FPIN2, it is assumed that the creep rate in the intermediate temperature range falls between the creep rate of \(\alpha\) uranium and solid solution \(\gamma\) phase. It is also assumed that creep of uranium alloys is governed by dislocation glide in the \(\gamma\) regime. The phenomenological creep equation is of the form

(11.3‑20)

\[\dot{\varepsilon} = C_{4} \sigma^{3} e^{\left( - \frac{Q_{\gamma}}{RT} \right)}\]

where \(Q_{\gamma}\) is the creep activation energy in the \(\gamma\) phase, and the parameter \(C_{4}\) is, in general, a function of composition and fission rate. Based on tracer diffusion data and/or additional creep data, \(Q_{\gamma}\) is approximated as the activation energy for creep of the pure \(\gamma\) uranium solvent

(11.3‑21)

\[ \begin{align} Q_{\gamma} = 28500 && \text{cal/g-mole} \end{align}\]

The constant \(C_{4}\) is chosen to fit data at 973K as

(11.3‑22)

\[\begin{align} C_{4} = 8.0 \cdot 10^{- 2} && \text{MPa}^{- 3} \text{s}^{- 1} \end{align}\]

Since the data is somewhat limited and the Eqs. 11.3-16 and 11.3-20 are consistent with the creep of alloys with a high percentage of uranium, FPIN2 uses these expressions for all metal fuels.

Experience with earlier versions of the FPIN2 reveals that the evaluation of Eqs. 11.3-16 and 11.3-20 uses a substantial fraction of the computing time. Since the terms involving the exponential function depend only on the temperature, in the current version such terms are evaluated only at the beginning of the mechanics calculation. Therefore, Eq. 11.3-16 is coded in the following form

(11.3‑23)

\[\dot{\varepsilon} = \sigma \left( A + B \sigma^{n - 1} \right)\]

where \(A\) and \(B\) are previously evaluated constants.

The creep equations given above are the same as those given in the Metallic Fuels Handbook [11-7]. However, fission-induced creep term is not included in the FPIN2 because it does not contribute significantly to the fuel stains on the time scale of the transients that the code has been designed to analyze.

11.3.5. Fission Induced Gas Swelling of Metal Fuels

Irradiated metal fuels experience rapid, large-scale swelling when subjected to overheating. From the safety perspective, this transient swelling could provide an inherent, self-limiting, negative reactivity feedback mechanism. In FPIN2, a transient swelling model has been incorporated to get an estimate for the magnitude of this fission-gas-induced swelling. This non-equilibrium model is based on diffusive growth of grain boundary bubbles. 2

The distribution of fission gas retained in the fuel matrix is specified as input in FPIN2, and the remainder of the gas is assumed to be in solution or in small bubbles within the fuel grains. Fission gas that is retained in the fuel during steady-state irradiation provides a source for expansion of both solid and liquid fuel during overheating. The amount of gas in the pin plenum is also important since the plenum pressure is a major contributor to cladding loading and, therefore, cladding failure.

The following assumptions are the basis of the transient swelling model of FPIN2 for the analysis of the metallic fuels:

  1. All “retained gas” (as measured from experiments on small samples of irradiated fuel) is in solution within the fuel or in the form of small bubbles within the fuel grains or on the grain boundaries.

  1. A certain fraction of this gas (input parameter) is contained in grain boundary bubbles, fixed in number and all of the same initial radius (0.1 microns).

  2. The open porosity contains fission gas in equilibrium with the plenum pressure.

  3. Only the grain boundary bubbles contribute to the swelling of solid fuel.

In order to use this model, a relationship between the swelling strain, \(\varepsilon^{s}\), the mean stress, \(\sigma_{m}\), and the fuel temperature, \(T\), is required as explained in Section 11.2.1.2.3. Studies on gas-bubbles growth mechanisms in metallic fuels indicate [11-14] that the transient fuel swelling is dominated by diffusive growth of grain boundary bubbles that is given by

(11.3‑24)

\[\frac{\text{dr}}{\text{dt}} = \frac{\Omega D_{\text{gb}}w \sin^{3}{\theta}}{2k T r^{2} L \eta} \phi\]

where

\(r\) = bubble radius in the grain boundary,

\(D_{\text{gb}}\) = grain boundary diffusion coefficient,

\(w\) = the boundary thickness,

\(\Omega\) = the atomic (or molecular) volume, and

\(kT\) = thermal energy.

The overpressure, \(\phi\), in Eq. 11.3-24 is defined as

(11.3‑25)

\[\phi = p_{\text{g}} + \sigma_{\text{m}} - \frac{2\gamma}{\rho}\]

and it represents the excess of gas pressure, \(p_{\text{g}}\), over the sum of the mean stress, \(\sigma_{\text{m}}\), and the surface tension restraint. This restraint is defined in terms of the specific surface free energy, \(\gamma\), and the radius of curvature, \(\rho\). For the grain-boundary bubbles,

(11.3‑26)

\[\rho = \frac{r}{\sin{\theta}}\]

where \(\theta\) is the intersection angle between the bubble surface and the grain boundary. In FPIN2, the fuel swelling rate is determined from Eq. 11.3-24 for the bubble growth rate and from the assumed bubble density.

The remaining constants in the above equations are geometric constants related to the bubble volume and the bubble spacing [11-14]. The gas in the bubbles is assumed to satisfy the ideal gas law so that

(11.3‑27)

\[p_{\text{g}} = \frac{mRT}{V_{\text{b}}}\]

where \(m\) is the fixed mass of the gas per bubble and \(V_{\text{b}}\) is the volume of the bubbles. The mass of the gas in a bubble is determined by using the equilibrium conditions at the pretransient time when the bubble radius is assumed to be \(r_{\text{o}}\). Thus

(11.3‑28)

\[m = \frac{V_{\text{bo}}}{RT_{\text{o}}}\left( - \sigma_{\text{mo}} + \frac{2\gamma}{r_{\text{o}}} \right)\]

Since all of the grain boundary gas is assumed to reside in the closed porosity bubbles, the number density of bubbles, \(N\), is given by

(11.3‑29)

\[N = \frac{F_{\text{g}}}{m}\]

where \(F_{\text{g}}\) is the grain boundary fission gas density with units of g per cm3 of fuel. For this simple model, only the swelling of the grain boundary bubbles is considered so that

(11.3‑30)

\[\varepsilon^{s} = \frac{1}{3}\left( \frac{\Delta V}{V_{\text{o}}} \right) = \frac{1}{3} N\left( V_{\text{b}} - V_{\text{bo}} \right)\]

These equations, when combined, give the complete specification of the transient swelling model in the form

(11.3‑31)

\[g\left( \varepsilon^{s}, \sigma_{\text{m}}, T \right) = 0\]

The treatment of this constitutive equation for swelling in finite element mechanics is discussed in Section 11.2.2.

11.3.6. Generalized Plastic Flow Behavior of Cladding Materials

The tensile properties of the cladding materials (HT9, D9, and Type 316 stainless steel) are not frequently used in analyses of fuel element performance during normal operation, as designs are usually such that stresses imposed do not challenge the points where yielding should occur. However, under transient conditions at elevated temperatures the applied stresses can result in general yielding. This behavior is modeled in FPIN2 by using a unified plastic deformation model that includes both rate-independent deformation (classical plasticity) and rate-dependent deformation (creep). The equations describing the generalized plastic flow behavior of Type 316 stainless steel, D9, and HT9 cladding are presented in the following sections. The yield strength and ultimate strength of the materials can be calculated from the plastic flow equations by applying them to the geometry of a tensile test.

11.3.6.1. Type 316 SS and D9 Cladding

The deformation behavior of Type 316 stainless steel and D9 cladding can be described over a broad range of temperature and strain-rate with the set of equations shown below. These equations are used in calculations of the strength and deformation at temperatures ranging from room temperature to 1675 K under a variety of loading conditions. In this formalism, the true equivalent flow stress, \(\sigma_{\text{e}}\), is described by the equation

(11.3‑32)

\[\sigma_{\text{e}} = \sigma_{\text{s}} - \left( \sigma_{\text{s}} - \sigma_{\text{l}} \right) \exp\left( - \frac{\hat{\varepsilon}}{\varepsilon^{*}} \right)\]

where \(\sigma_{\text{l}}\) is yield stress of fully annealed unirradiated material, and \(\sigma_{\text{s}}\) is the saturation value of the flow stress that is asymptotically approached at large values of plastic strain. Ongoing true plastic strain is incorporated in the “hardness” parameter, \(\hat{\varepsilon}\) (input in the model), as this strain is accumulated. The hardness parameter also contains contributions from prior cold work, irradiation hardening, and softening caused by annealing, all scaled as true plastic strain. Consequently, \(\hat{\varepsilon} = 0\) for the fully annealed unirradiated material, and \(\hat{\varepsilon} = 0.223\) for 20% cold-worked material (default value in the code). Finally, the value of the parameter \(\varepsilon^{*}\) is determined by the initial hardening rate \(\theta_{\text{l}}\) of fully annealed material \(\left( \hat{\varepsilon} = 0 \right)\) which is given by

(11.3‑33)

\[\theta_{\text{l}} = \frac{\sigma_{\text{s}} - \sigma_{\text{l}}}{\varepsilon^{*}}\]

The initial hardening rate, \(\theta_{\text{l}}\), is found to be only temperature dependent through the shear modulus; so, the rate and temperature dependencies of \(\varepsilon^{*}\) is obtained from the above equation.

The yield stress and saturation flow stress are rate and temperature dependent in accordance with the equations

(11.3‑34)

\[\frac{\sigma_{\text{s}}}{G} = \frac{\sigma_{\text{so}}}{G}\left\lbrack 1 - \exp\left( - \left( \frac{{\dot{\varepsilon}}_{\text{p}}}{{\dot{\varepsilon}}_{\text{os}}} \right)^{mK} \right) \right\rbrack^{1/K}\]

(11.3‑35)

\[\frac{\sigma_{\text{l}}}{G} = \frac{\sigma_{\text{lo}}}{G}\left\lbrack 1 - \exp\left( - \left( \frac{{\dot{\varepsilon}}_{\text{p}}}{{\dot{\varepsilon}}_{\text{ol}}} \right)^{mK} \right) \right\rbrack^{1/K}\]

where \(G\) is the temperature dependent shear modulus, \({\dot{\varepsilon}}_{\text{p}}\) is the equivalent plastic strain rate (the \(\frac{\text{d}{\overline{\varepsilon}}^{p}}{\text{dt}}\) term in Eq. 11.3-26), \(m\) and \(K\) are constants, and \(\sigma_{\text{so}}\), \(\sigma_{\text{lo}}\), \({\dot{\varepsilon}}_{\text{os}}\) and \({\dot{\varepsilon}}_{\text{ol}}\) are temperature dependent functions. The constant \(m\) is the rate sensitivity (the reciprocal of the stress-exponent n in a power-law creep equation), and \(K\) is a non-physical fitting parameter that governs the sharpness of the transition between rate-dependent flow and rate-independent flow. At a given temperature, in the high strain-rate limit, the equation for the saturation stress reduces to

(11.3‑36)

\[\frac{\sigma_{\text{s}}}{G} = \frac{\sigma_{\text{so}}}{G}\]

whereas in the low strain-rate limit, it becomes

(11.3‑37)

\[\frac{\sigma_{\text{s}}}{G} = \frac{\sigma_{\text{so}}}{G}\left( \frac{{\dot{\varepsilon}}_{\text{p}}}{{\dot{\varepsilon}}_{\text{os}}} \right)^{m}\]

which reflects the typical power-law creep behavior observed for these materials at high temperatures and low strain rates. The equation for the yield stress behaves similarly. The temperature dependencies of \(\sigma_{\text{so}}\) and \(\sigma_{\text{lo}}\) are in large part eliminated by dividing by \(G\). However, \({\dot{\varepsilon}}_{\text{os}}\) and \({\dot{\varepsilon}}_{\text{ol}}\) reflect a temperature dependence of high-temperature creep of the form:

(11.3‑38)

\[{\dot{\varepsilon}}_{\text{os}} = {\dot{\varepsilon}}_{\text{oos}} \exp\left( - Q/RT \right)\]

(11.3‑39)

\[{\dot{\varepsilon}}_{\text{ol}} = {\dot{\varepsilon}}_{\text{ool}} \exp\left( - \frac{Q}{RT} \right)\]

where \({\dot{\varepsilon}}_{\text{oos}}\) and \({\dot{\varepsilon}}_{\text{ool}}\) are constants; Q, R, and T are the creep activation energy, gas constant, and absolute temperature, respectively. The values for the parameters in all of the above equations are:

\(G = 92.0 - 4.02 10^{-2} T\) (GPa)

\(\frac{\theta_{\text{l}}}{G} = 3.66 \times 10^{-2}\)

\(\frac{\sigma_{\text{so}}}{G} = 2.00 \times 10^{-2} - 9.12 \times 10^{-6} T\)

\(\frac{\sigma_{\text{lo}}}{G} = 2.06 times 10^{-3} + 7.12 \times 10^{-1} T\)

\(m = \frac{1}{5.35} = 0.187\)

\(\frac{Q}{R} = 38533\) (K)

\({\dot{\varepsilon}}_{\text{oos}} = 1.062 \times 10^{14}\) (s-1)

\({\dot{\varepsilon}}_{\text{ool}} = 3.794 \times 10^{12}\) (s-1)

\(K = 2.0\)

The combination of the Eq. 11.3-32 and the above strain-rate and temperature laws allow one to generate the complete rate and temperature dependent true-stress/true-strain curves. The yield stress, \(\sigma_{\text{l}}\), in these equations can be compared to the 0.2% offset yield stress. The ultimate strength, while not specifically denoted in the equations, can be determined by differentiating the Eq. 11.3-32 in accordance with the construction of true stress, \(\sigma_{\text{u}}\), that corresponds to the engineering ultimate strength:

(11.3‑40)

\[\frac{\text{d}\sigma}{\text{d}\varepsilon} \bigg\rvert_{\varepsilon_{\text{pu}}} = \theta_{\text{l}} \exp\left( - \frac{\varepsilon_{\text{pu}}}{\varepsilon^{*}} \right) = \sigma_{\text{u}}\]

where \(\varepsilon_{\text{pu}}\) is the true plastic strain at the ultimate strength which corresponds to the uniform elongation. Solving the above equation and Eq. 11.3-32 simultaneously gives \(\sigma_{\text{u}}\) and \(\varepsilon_{\text{pu}}\).

11.3.6.2. HT9 Cladding

The deformation behavior of the martensitic-ferritic stainless steel HT9 cladding can be described over a broad range of temperature and strain-rate with the set of equations shown below. These equations are incorporated in FPIN2 and used in calculations of strength and ductility of this alloy at temperatures ranging from room temperature to 1110 K under a variety of loading conditions. The high temperature creep behavior is adequately described by an equation of the Dorn power-law form:

(11.3‑41)

\[\frac{{\dot{\varepsilon}}_{\text{p}}}{{\dot{\varepsilon}}_{\text{oos}}} = \left( \frac{E}{\sigma_{\text{so}}} \right)^{n} \left( \frac{\sigma_{\text{s}}}{E} \right)^{n} \exp\left( - \frac{Q_{\text{c}}}{\text{kT}} \right)\]

In this equation, \({\dot{\varepsilon}}_{\text{p}}\) is the steady-state equivalent creep rate normalized by the constant \({\dot{\varepsilon}}_{\text{oos}}\), \(\sigma_{\text{s}}\) is the equivalent applied stress in the creep test normalized by the constant \(\sigma_{\text{so}}\), \(E\) is the temperature dependent Young’s modulus, \(Q_{\text{c}}\) is the creep activation energy, \(k\) is Boltzmann’s constant, and \(T\) is absolute temperature. The temperature dependence of \(E\) is expressed as

(11.3‑42)

\[\begin{align} E\left( T \right) = 2.12 \cdot 10^{11}\left\lbrack 1.144 - 4.856 \cdot 10^{- 4} T \right\rbrack && \left( \text{Pa} \right) \end{align}\]

and the remaining parameters in the creep equation are given by

\[\begin{align} {\dot{\varepsilon}}_{\text{oos}} = 5.1966 \cdot 10^{10} && \left( \text{s}^{- 1} \right) \end{align}\]
\[\frac{\sigma_{\text{so}}}{E} = 3.956 \cdot 10^{- 3}\]
\[n = 2.263\]
\[\begin{align} \frac{Q_{\text{c}}}{k} = 36739 && \left( \text{K} \right) \end{align}\]

The flow-stress/strain-hardening behavior can be described by an equation of the type

(11.3‑43)

\[\frac{\sigma_{\text{e}}}{E} = \frac{\sigma_{\text{s}}}{E} - \frac{\sigma_{\text{s}} - \sigma_{\text{l}}}{E} \exp\left( - \frac{{\overline{\varepsilon}}_{\text{p}}}{\varepsilon^{*}} \right)\]

where \(\sigma_{\text{e}}\) is the equivalent flow stress at some value of true equivalent plastic strain \({\overline{\varepsilon}}^{p}\) (defined in Eq. 11.3-26), measured from the reference state of as-heat-treated material \(\left( {\overline{\varepsilon}}^{p} = 0 \right)\). As described in Section 11.3.6.1, \(\theta_{\text{l}}\) is the yield stress of as-heat-treated material, \(\varepsilon^{*}\) is a temperature dependent parameter extracted from the hardening rate at the ultimate strength, and \(\sigma_{\text{s}}\) is the saturation stress, or steady-state flow stress, approached at a large plastic strain and assumed compatible with the stress in the above steady-state creep equation. It has been found that assuming that the yield stress \(\sigma_{\text{l}} = 0.8 \sigma_{\text{s}}\) throughout allows the generation of flow stress-strain curves that agree well with data. The true stress at the ultimate strength (maximum load), after differentiating Eq. 11.3-43, is given by

(11.3‑44)

\[\frac{\sigma_{\text{u}}}{E} = \frac{\sigma_{\text{s}} - \sigma_{\text{l}}}{E} \varepsilon^{*} \exp\left( \frac{{- \varepsilon}_{\text{pu}}}{\varepsilon^{*}} \right)\]

where \(\varepsilon_{\text{pu}}\) is the true strain at the ultimate strength corresponding to the uniform elongation. Simultaneously solving this equation in conjunction with the flow stress equation evaluated at \(\varepsilon_{\text{pu}}\), the following expression is obtained for the quantity \(\varepsilon^{*}\):

(11.3‑45)

\[\varepsilon^{*}\left( T \right) = 0.12733 - 3.5027 \cdot 10^{- 4} T + 2.9934 \cdot 10^{- 7}T^{2}\]

An equation to model the transition between high-rate, low-temperature, rate independent flow behavior and creep behavior has the form

(11.3‑46)

\[\frac{\sigma_{\text{s}}}{E} = \frac{\sigma_{\text{so}}}{E}\left\{ 1 - exp\left\lbrack \left( \frac{{\overline{\varepsilon}}_{\text{p}} \exp\left( \frac{Q_{\text{t}}}{kT} \right)}{\varepsilon_{\text{oos}}} \right)^{K/n} \right\rbrack \right\}^{1/K}\]

Because of microstructural changes that occur during long-time creep testing of HT9, and the effects they have on flow stress, there is not a smooth transition between short term tensile behavior and creep behavior. However, setting the non-physical fitting parameter as \(K = 0.2\), and with all the other parameters set as indicated above, Eq. 11.3-46 follows the tensile data quite well, and reduces to the power-law creep equation at very low strain rates.

For temperatures above the completion of the \(\gamma\) transformation temperature (1233 K), the flow stress model in FPIN2 assumes that the deformation rate for HT9 cladding is the same as that given by the Type 316 SS and D9 cladding equations in Section 11.3.6.1. A simple mixture rule is used to calculate the deformation rates in the \(\alpha\)-\(\gamma\) transition region (1110K-1233K).

11.3.7. Fuel-Cladding Eutectic Formation

An additional complication for metallic fuels is the formation of a low melting point eutectic alloy between the fuel and the cladding that can contribute to fuel element failure during transient overheating events. The eutectic alloy forms due to interdiffusion of fuel and cladding constituents at the fuel-cladding interface. It melts at a temperature lower than the fuel and cladding solidus temperatures.

11.3.7.1. Eutectic Penetration of the Cladding

The primary effect of the liquid eutectic alloy formation on fuel element failure is thinning of the cladding wall. The liquid does not further damage the remaining cladding tendon by mechanisms such as liquid metal embrittlement. Experimental data documenting this effect comes from constant temperature time-to-failure experiments on irradiated EBR-II fuel elements and form out-of-pile dipping tests of penetration of iron by molten uranium-iron eutectic alloy. These data have recently been reviewed and the following correlation has been recommended [11-15]

(11.3‑47)

\[\begin{split} \dot{M} = \begin{cases} 0 && T < 1353 \text{ K} \\ 922 + 2.93 \left( T - 1388 \right) - 0.215 \left( T - 1388 \right)^{2} + 0.001134 \left( T - 1388 \right)^{3} && 1355 \text{ K} \leq T \leq 1506 \text{ K} \\ \exp\left( 22.85 - \frac{27624}{T} \right) && T > 1506 \text{ K} \end{cases}\end{split}\]

where \(\dot{M}\) is the melt rate in μm/s, and \(T\) is the absolute temperature in Kelvins. The major feature of the data and these equations is that the melt rate is very rapid for temperatures above 1353 K. For typical cladding dimensions, this means that cladding will completely melt-through is somewhat less than one second once molten fuel reaches 1353K. On the other hand, for the transient rate of concern for a typical FPIN2 application, negligible cladding melt-through occurs until this “rapid eutectic attack” temperature is reached. 3 The Eq. 11.3-47 is used in cladding failure criterion subroutine of FPIN2 to calculate wall thinning.

11.3.7.2. Eutectic Release of Cladding Axial Restraint

When fuel and cladding are in contact, FPIN2 has an option that allows that fuel and cladding to either slip freely axially or to be locked together. In the latter option, it is assumed that the fuel and cladding remain locked, either by metallurgical bonding or by friction, at temperatures below the assumed eutectic alloy melting threshold, \(T_{\text{th}}\) (an input in the code). At temperatures above \(T_{\text{th}}\), the fuel and cladding are assumed to slip freely because of the liquid phase that forms at the fuel cladding interface. This discontinuity in the cladding axial constraint model of FPIN2 at the assumed eutectic alloy melting temperature causes a sudden drop in stress in both fuel and cladding elements as the fuel expands slightly in the axial direction and contracts in the radial direction due to elastic recovery.

When the fuel-cladding gap is open, the fuel and cladding axial displacements are calculated independently. When the fuel and cladding are locked together, the increments in axial strain in the fuel are assumed to be equal to the increments in axial strain in the cladding. This is handled in mechanics calculation by modifying the terms in the stiffness matrix. Specifically, the finite-element algorithm is solved for the resultant fuel and cladding axial forces by combining the fuel and cladding stiffness matrices in a non-iterative procedure. In this procedure, the two equations for the axial displacements are replaced by a single equation for their joint displacement. The redundant equation is then replaced with an identity so that the size of the stiffness matrix and all the associated coding remains the same. When the sudden release of cladding restraint occurs at eutectic alloy melting point, fairly large changes take place in the axial displacements. However, experience has shown that, due to the robust procedure used, this substantial sudden change in the mechanics causes no computational difficulties.

11.3.8. Fuel Element Failure

The failure of the cladding due to temperature and pressure transients is essentially independent of fuel type except that failure in metal fueled elements is augmented by eutectic formation which penetrates and, in effect, thins the cladding as discussed in the previous section. Also, the loss in cladding stiffness is not considered in calculating the continued deformation of the cladding in the oxide fuel version of FPIN2. For conditions where simultaneous eutectic erosion and plastic deformation are important, oxide fuel version under-predicts the permanent cladding strains at the time of failure. Therefore, the FPIN2 analysis of cladding deformation has been modified to account for this eutectic thinning and to determine the permanent cladding strain remaining after a transient more accurately.

The technique that has been implemented in the FPIN2 cladding deformation model to account for the eutectic attack uses an effective ligament thickness ratio to decease the contribution that the material stresses in a given element make to the overall balance between internal and external forces. Elements that are fully liquid are assumed to have a ligament thickness ratio of zero and those that are solid are assumed to have ligament ratio of one. An element that has partially liquefied then has a ratio equal to the fraction of the element that is still solid. In other words, the effect of eutectic formation is included in FPIN2 by considering only the thickness of unaffected cladding that is available to carry the load. This model has the additional benefit that it can easily be modified to account for other damage mechanisms such as cavitation and void growth during tertiary creep. In the current version of the FPIN2, however, these additional mechanisms are not considered.

For transients with relatively short time scale, the eutectic penetration is a factor only if the cladding temperature exceeds the “rapid eutectic attack” temperature (1353K) where the penetration rate increases by three orders of magnitude. The change in cladding wall thickness above 1353K is calculated using Eq. 11.3-47 for Type 316 SS, D9, and HT9 cladding. At temperatures below 1353K, the rate of eutectic penetration is generally insignificant and assumed to be zero.

Cladding rupture is predicted in the code by using the life fraction criteria. The TCD-2 life fraction criterion [11-16] is used for the fuel elements with D9 and Type 316 stainless steel cladding; and the HEDL transient Dorn parameter correlation [11-17] is used for the fuel elements with HT9-cladding. The fuel adjacency effect in the TDC-2 correlation has been neglected. The life fraction over a time step, \(\text{dt}\), is calculated from the rupture time, \(t_{\text{r}}\), for the instantaneous average cladding temperature, \(T_{\text{c}}\), and hoop stress, \(\sigma_{\text{c}}\). Life fractions are summed in the usual way so that cladding failure is predicted to occur at time \(t_{\text{f}}\) when

(11.3‑48)

\[\int_{0}^{t_{\text{f}}}{\frac{1}{t_{\text{r}}\left( \sigma_{\text{c}}, T_{\text{c}} \right)}} \text{dt} = 1\]

The location of the failure is the location of first axial cladding segment to reach a life fraction of 1.0. The hoop stress in the cladding tendon is calculated by the thin-shell equations consistent with the methodology used in developing the life fraction correlations. The thickness of this tendon is reduced by the amount of eutectic attack.

Footnotes

1

In stand-alone FPIN2 calculation, this plenum temperature is assumed to be equal to the coolant outlet temperature.

2

FPIN2 also provides an equilibrium swelling model as an option.

3

The rate at 1000K is about 10-2 µm/s, which would require 10 hours to completely penetrate the full cladding wall thickness.