Fuel Behaviour — 1-D Radial Thermo-Mechanics

Key Facts 

Read this before modifying the fuel behaviour solver.

1D radial steady-state: fuel pellet + gap + cladding
Temperature: \(-\frac{1}{r}\frac{d}{dr}(r k \frac{dT}{dr}) = q'''\) with temperature-dependent conductivity
Gap conductance: radiation + gas conduction (MATPRO correlations)
Cladding stress: thick-wall cylinder (Lamé equations) with internal pressure + thermal stress
Material properties: data/materials/matpro.py (UO2, Zircaloy, gap gases)
IAPWS viscosity fallback: NaN at very high temperatures fixed via fallback (ERR-010)

Overview 

Module 06 simulates the long-term thermo-mechanical behaviour of a single PWR fuel rod under steady-state irradiation (6 years at 200 W/cm). The model couples:

Radial heat conduction — cylindrical UO₂ pellet (30 nodes) and Zircaloy cladding (5 nodes) with temperature-dependent properties.
Thermo-elastic-plastic deformation — fuel and cladding strains including thermal expansion, elastic, creep, and (fuel only) swelling and (clad only) plasticity.
Fuel-cladding gap — helium-filled annular gap with gas-conduction conductance, evolving geometry, and gap closure detection.
Internal gas pressure — ideal gas law with fission gas release (He, Kr, Xe).
Cladding stress — radial equilibrium, strain compatibility, and axial force balance solved as an algebraic linear system.

The simulation discovers the gap closure event (~2.85 years for the default parameters) and continues with closed-gap boundary conditions to end-of-cycle (6 years).

The solver is implemented in solve_fuel_behaviour() in 06.Fuel.Behaviour/fuel_behaviour.py.

Key Design Decision: DAE-to-ODE Restructuring 

The MATLAB original uses ode15s with a diagonal mass matrix to solve a mixed ODE/AE (differential-algebraic) system. Stress components are algebraic variables (zeros on the mass matrix diagonal). scipy.integrate.solve_ivp does not support mass matrices. The BDF time integration approach is shared with the thermal hydraulics module; see ODE Integration.

Approach adopted: Remove stresses from the state vector entirely. At every RHS evaluation, _solve_stress() constructs and solves the \(3(n_f + n_c)\) linear system for stress unknowns. The ODE state vector contains only time-evolving quantities:

\[\mathbf{y} = \bigl[\, T_{\text{fuel}}^{(n_f)},\; T_{\text{clad}}^{(n_c - 1)},\; F,\; \dot{V}_S^{(n_f)},\; \varepsilon_C^{(3 \times n_f)},\; \varepsilon_C^{(3 \times n_c)},\; \bar\varepsilon_P^{(n_c)},\; \varepsilon_P^{(3 \times n_c)} \,\bigr]\]

Default mesh: \(n_f = 30\) fuel radial nodes, \(n_c = 5\) clad radial nodes. Total ODE variables: \(30 + 4 + 1 + 30 + 90 + 15 + 5 + 15 = 190\). (The MATLAB DAE has ~280 variables including stresses.)

Why this works: The stress equilibrium, strain compatibility, and boundary conditions form a linear system in the stresses when the non-elastic strains (thermal, creep, swelling, plastic) are given. The elastic strains depend linearly on stress through Hooke’s law:

\[\varepsilon_i^E = \frac{\sigma_i - \nu(\sigma_j + \sigma_k)}{E} = C_s \, \sigma_i + C_c (\sigma_j + \sigma_k)\]

where \(C_s = 1/E\) (self compliance) and \(C_c = -\nu/E\) (cross compliance). The total strain \(\varepsilon_i = \varepsilon_i^T + \varepsilon_i^E + \varepsilon_i^C + \varepsilon_i^S/3\) (fuel) or \(\varepsilon_i = \varepsilon_i^T + \varepsilon_i^E + \varepsilon_i^P + \varepsilon_i^C\) (clad) is linear in the stresses. Therefore, the equilibrium and compatibility equations are linear in stress, and the system can be solved with np.linalg.solve — no iteration needed.

Tradeoff: The algebraic approach uses the undeformed (fabrication) geometry for the node positions in the stress system. Since strains are O(1%), the geometry error is O(0.01%), which is negligible. The MATLAB DAE solver has the luxury of iterating geometry and stress simultaneously. A small transient settling difference occurs immediately after gap closure (see FB-20260401-004 in the improvement tracker).

Physics Equations 

Fuel Temperature 

Identical to the thermal hydraulics formulation (§ Thermal Hydraulics — Single-Channel LOCA, fuel temperature section), but with the full active fuel height (\(h_0 = 3\) m) rather than an axial mesh. The fuel operates at a constant linear heat generation rate of 200 W/cm = 20 000 W/m.

(1)\[m_j \, c_p(T_j) \frac{dT_j}{dt} = -(Q_{j+1} - Q_j) + \dot{q}'''_j V_{j,0}\]

The heat flows \(Q_j\) at the \(n_f + 1\) radial interfaces are:

\[\begin{split}Q_0 &= 0 \qquad \text{(adiabatic centre)} \\ Q_j &= -k_{\text{UO}_2}(T_{\text{mid}}, Bu, p) \; \frac{T_j - T_{j-1}}{\Delta r} \; A_{\text{mid},j} \qquad j = 1, \ldots, n_f - 1 \\ Q_{n_f} &= h_{\text{gap}} \, (T_{f,n_f} - T_{c,1}) \, A_{\text{gap}}\end{split}\]

The outer clad surface temperature is fixed at \(T_{\text{out}} = 600\) K (no coolant model in this module — the TH coupling is added in Modules 07–08).

Strain Decomposition 

For each strain component \(i \in \{r, \theta, z\}\):

Fuel:

(2)\[\varepsilon_i^{\text{fuel}} = \varepsilon^T(T) + \varepsilon_i^E(\boldsymbol{\sigma}) + \varepsilon_i^C + \frac{\dot{V}_S}{3}\]

where \(\varepsilon^T = \text{thExp}(T)\) is the isotropic thermal expansion, \(\varepsilon_i^E\) is the elastic strain (Hooke’s law), \(\varepsilon_i^C\) is the directional creep strain, and \(\dot{V}_S/3\) is the isotropic volumetric swelling divided equally among the three directions.

Cladding:

(3)\[\varepsilon_i^{\text{clad}} = \varepsilon^T(T) + \varepsilon_i^E(\boldsymbol{\sigma}) + \varepsilon_i^P + \varepsilon_i^C\]

Fuel Swelling 

The volumetric fuel swelling rate has two contributions (MATPRO [MATPRO2003]):

(4)\[\dot{V}_S = \underbrace{2.5 \times 10^{-29} \dot{F}}_{\text{solid FP}} \;+\; \underbrace{8.8 \times 10^{-56} \dot{F} \, (2800 - T)^{11.73} e^{-0.0162(2800 - T)} e^{-8 \times 10^{-27} F}}_{\text{gaseous (T < 2800\,K)}}\]

where \(\dot{F}\) is the fission rate density (fissions/m³/s) and \(F\) is the cumulative fission density (fissions/m³). The gaseous contribution is strongly temperature-dependent and vanishes above 2800 K (restructuring releases fission gases).

Fuel and Clad Creep 

The effective creep rate follows simplified MATPRO correlations:

\[\begin{split}\dot{\varepsilon}_{\text{eff}}^{\text{fuel}} &= 5 \times 10^5 \, \sigma_{\text{VM}} \, \exp\!\left(-\frac{4 \times 10^5}{8.314 \, T}\right) \\ \dot{\varepsilon}_{\text{eff}}^{\text{clad}} &= 10^5 \, \sigma_{\text{VM}} \, \exp\!\left(-\frac{2 \times 10^5}{8.314 \, T}\right)\end{split}\]

The deviatoric components follow the Prandtl–Reuss flow rule:

\[\dot{\varepsilon}_i^C = \dot{\varepsilon}_{\text{eff}} \frac{\sigma_i - \sigma_{\text{mean}}}{\sigma_{\text{VM}}}\]

where \(\sigma_{\text{mean}} = (\sigma_r + \sigma_\theta + \sigma_z)/3\) and \(\sigma_{\text{VM}} = \sqrt{\frac{1}{2}[(\sigma_r - \sigma_\theta)^2 + (\sigma_r - \sigma_z)^2 + (\sigma_\theta - \sigma_z)^2]}\).

Clad Plastic Strain 

Same Norton power-law as the thermal hydraulics module (§ Thermal Hydraulics — Single-Channel LOCA, equation (8)):

\[\dot{\bar\varepsilon}_P = 10^{-3} \left( \frac{\sigma_{\text{VM}} \times 10^6} {K(T) \, |\bar\varepsilon_P + 10^{-6}|^{n(T)}} \right)^{1/m(T)}\]

with Prandtl–Reuss deviatoric components. The \(K(T), m(T), n(T)\) are temperature-dependent Zircaloy strength parameters from MATPRO.

Algebraic Stress Solver 

For each RHS evaluation, the function _solve_stress() constructs and solves the linear system \(\mathbf{A} \boldsymbol{\sigma} = \mathbf{b}\):

System Structure 

\(3(n_f + n_c)\) unknowns ordered as:

\[\boldsymbol{\sigma} = [\, \sigma_r^{f,1} \ldots \sigma_r^{f,n_f},\; \sigma_\theta^{f,1} \ldots \sigma_\theta^{f,n_f},\; \sigma_z^{f,1} \ldots \sigma_z^{f,n_f},\; \sigma_r^{c,1} \ldots \sigma_r^{c,n_c},\; \sigma_\theta^{c,1} \ldots \sigma_\theta^{c,n_c},\; \sigma_z^{c,1} \ldots \sigma_z^{c,n_c} \,]\]

Equations per region (fuel / clad):

Radial stress equilibrium (\(n - 1\) equations per region):

\[\frac{d\sigma_r}{dr} = \frac{\sigma_\theta - \sigma_r}{r}\]

Discretised at node boundary midpoints:

\[\frac{\sigma_{r,j+1} - \sigma_{r,j}}{\Delta r_j} = \frac{(\sigma_{\theta,j} + \sigma_{\theta,j+1})/2 - (\sigma_{r,j} + \sigma_{r,j+1})/2}{r_{\text{mid},j}}\]
Strain compatibility (\(n - 1\) equations per region):

\[\frac{d\varepsilon_\theta}{dr} = \frac{\varepsilon_r - \varepsilon_\theta}{r}\]

Since \(\varepsilon_\theta = C_s \sigma_\theta + C_c(\sigma_r + \sigma_z) + \varepsilon_\theta^{\text{NE}}\) (where NE = non-elastic), this becomes a linear equation in the stresses with the non-elastic strains on the RHS.
Axial strain uniformity (\(n - 1\) equations per region):

\[\varepsilon_z(r_j) = \varepsilon_z(r_{j+1})\]
Boundary conditions (6 equations):
- BC1 — Inner fuel surface: \(\sigma_r(0) = -p_{\text{gas}}\) (central hole) or \(\sigma_r(0) = \sigma_\theta(0)\) (solid cylinder, no central hole)
- BC2 — Outer clad surface: \(\sigma_r(r_{\text{out}}) = -p_{\text{cool}}\)
- BC3 — Fuel/clad gap (radial stress): open gap: \(\sigma_r^f(r_f) = -p_{\text{gas}}\); closed gap: radial stress continuity \(\sigma_r^c(r_{\text{in}}) = \sigma_r^f(r_f)\)
- BC4 — Fuel/clad gap (hoop strain): open gap: \(\sigma_r^c(r_{\text{in}}) = -p_{\text{gas}}\); closed gap: displacement constraint (see below)
- BC5 — Axial force balance (fuel): open gap: \(\int \sigma_z^f \, r \, dr = 0\); closed gap: axial strain continuity across gap
- BC6 — Axial force balance (clad or fuel+clad): \(\int \sigma_z \, r \, dr = p_{\text{in}} r_{\text{in}}^2 - p_{\text{out}} r_{\text{out}}^2\)

Closed-Gap BC4 — Displacement-Based Constraint 

FB-20260401-002

The original MATLAB formulation for BC4 (closed gap) uses a differential form:

\[\frac{\varepsilon_\theta^c(r_{\text{in}}) - \varepsilon_\theta^f(r_f) - \Delta\varepsilon_h} {\delta_{\text{gap}}} = \frac{(\varepsilon_r^f + \varepsilon_r^c)/2 - (\varepsilon_\theta^f + \varepsilon_\theta^c)/2} {r_{\text{gap}}}\]

This requires dividing by the gap width \(\delta_{\text{gap}}\), which is problematic for the algebraic solver:

The MATLAB DAE solver iterates and converges to the deformed \(\delta\)
The algebraic solver uses fabrication geometry, where \(\delta = 100\,\mu\text{m}\) instead of the physical ~6 μm (roughness)
This causes a 17× error in the stress gradient (ERR-013)

Fix adopted: Replace with a displacement-based gap constraint:

(5)\[r_{c,\text{in}}^0 \, (1 + \varepsilon_\theta^c(r_{\text{in}})) - r_{f,\text{out}}^0 \, (1 + \varepsilon_\theta^f(r_f)) = \varepsilon_{\text{rough}}\]

where \(\varepsilon_{\text{rough}} = 6\,\mu\text{m}\). This formulation:

Is physically transparent: the deformed gap width equals the roughness
Is linear in the stresses (through the Hooke’s law dependence of \(\varepsilon_\theta\) on \(\sigma\))
Avoids the \(1/\delta_{\text{gap}}\) amplification that made the differential form ill-conditioned
Eliminates the dependency on \(\Delta\varepsilon_h\) (the hoop strain offset frozen at the moment of closure)

Result: Contact pressure 40.7 MPa vs MATLAB 39.8 MPa (2.2% match).

Gap Closure Event 

The gap between fuel and cladding closes when the deformed gap width \(\delta = r_{c,\text{in}} - r_{f,\text{out}}\) decreases to the surface roughness \(\varepsilon_{\text{rough}} = 6\,\mu\text{m}\):

\[E(t, \mathbf{y}) = \delta(t) - \varepsilon_{\text{rough}}\]

The event function uses a one-shot cache from the last RHS evaluation for efficiency. When \(E\) crosses zero (direction \(-1\), terminal), solve_ivp stops Phase 1.

After closure, the solver:

Records the strain jumps \(\Delta\varepsilon_h\) and \(\Delta\varepsilon_z\) at the fuel/clad contact interface
Switches to closed-gap BCs (BC3, BC4, BC5 change — see above)
Relaxes tolerances from rtol=1e-6 to rtol=1e-4 (matching MATLAB)
Continues integration from the closure state to end-of-cycle

Internal Gas Pressure 

The fuel rod contains an initial charge of helium at 1 MPa (Module 06; 1.2 MPa for Modules 07–08). Fission gas (He, Kr, Xe) is generated proportionally to fission density and released at a constant fraction (FGR = 6%):

\[\mu_{\text{gen},g} = F \cdot V_{\text{fuel}} / N_A \cdot f_g \qquad \mu_{\text{rel},g} = \mu_{\text{gen},g} \cdot \text{FGR}\]

where \(f_g = (0.01, 0.02, 0.23)\) for (He, Kr, Xe) as fractions of fission products.

The gas pressure is computed from the ideal gas law with the gas distributed between the plenum and the gap/void volumes:

\[p_{\text{gas}} = \frac{(\mu_{\text{He},0} + \sum_g \mu_{\text{rel},g}) \, R} {\displaystyle \frac{V_{\text{plenum}}}{T_{\text{plenum}}} + \frac{V_{\text{void}}}{T_{f,1}} + \frac{V_{\text{gap}}}{T_{\text{gap}}}}\]

The gap gas conductivity is computed using the Prandtl mixing rule for He–Kr–Xe mixtures (MATPRO [MATPRO2003]):

\[k_{\text{mix}} = \sum_g \frac{k_g \, x_g} {\sum_h \psi(k_g, k_h, M_g, M_h) \, x_h}\]

where \(\psi(k_1, k_2, M_1, M_2) = (1 + \sqrt{\sqrt{M_1/M_2} \, k_1/k_2})^2 / \sqrt{8(1 + M_1/M_2)}\), \(x_g\) are mole fractions, and \(M_g\) are molecular weights (He=2, Kr=54, Xe=36).

Validation 

Comparison with MATLAB at t = 1 day (open gap, pre-swelling):

Quantity	Python	MATLAB	Diff
Fuel centre T (°C)	1017.80	1017.86	0.06
Fuel outer T (°C)	543.06	543.03	0.03
Clad inner T (°C)	349.66	349.51	0.15
Clad outer T (°C)	327.00	326.85	0.15
Gap width (μm)	71.19	71.29	0.10
Gas pressure (MPa)	2.4855	2.4846	0.04%
Fuel outer r (mm)	4.15197	4.15200	0.03 μm

Comparison at t = 4.69 years (after gap closure):

Quantity	Python	MATLAB	Diff
Fuel centre T (°C)	1142	1124	1.6%
Contact pressure (MPa)	40.7	39.8	2.2%
Gas pressure (MPa)	6.923	6.930	0.1%
Gap width (μm)	6.0	5.8	3%

The open-gap phase achieves excellent parity (< 0.1%). The closed-gap phase shows 1–2% drift due to the algebraic stress solver’s operator-split treatment of the stiff stress-creep feedback (FB-20260401-004).

Known Limitations 

1-D radial only — no axial mesh, no axial conduction or fuel relocation.
Algebraic stress solver uses undeformed geometry — errors are O(0.01%) per strain level, negligible for PWR fuel but may matter for high-burnup fuels with large swelling.
Closed-gap transient settling — the operator-split approach resolves the stress-creep feedback more abruptly than the MATLAB DAE (~0.2 yr vs ~2 yr to reach positive hoop stresses after closure).
Fixed outer clad temperature — no coolant model. The TH coupling is added in Modules 07 and 08.
Simplified creep — single-term Arrhenius correlations. No irradiation creep, no stress relaxation.

References 

[MATPRO2003] (1,2)

D.L. Hagrman et al., MATPRO — A Library of Materials Properties for Light-Water-Reactor Accident Analysis, NUREG/CR-6150, Idaho National Laboratory, 2003.