Thermal Hydraulics — Single-Channel LOCA
Key Facts
Read this before modifying the thermal hydraulics solver.
Single-channel LOCA transient: 1D axial, two-phase flow
DAE system solved with BDF (
scipy.integrate.solve_ivp)Conservation: mass, momentum, energy for liquid + vapor
Wall heat transfer: Dittus-Boelter (1-phase), Chen (nucleate boiling), film boiling correlations
Water properties: pyXSteam (IAPWS-IF97) via
data/materials/h2o_properties.pyTwo solver versions exist: ODE (original MATLAB port) and DAE (improved)
Gotcha: MATLAB version has a gap conductance bug that ORPHEUS intentionally does NOT reproduce
Overview
Module 07 simulates a Loss-of-Coolant Accident (LOCA) in a single PWR fuel channel. The model couples:
Radial fuel heat conduction — cylindrical UO2 pellet with volumetric heat generation, burnup-dependent conductivity, and thermal expansion.
Radial clad heat conduction — Zircaloy-4 tube with thermal expansion, creep plasticity, and a burst-stress failure criterion.
Axial coolant energy balance — 1-D upward flow with two-phase correlations (Dittus–Boelter, Thom nucleate boiling, Zuber CHF, film boiling).
Gap mechanics — helium-filled gap with deformable fuel/clad geometry and pressure-dependent contact conductance.
Internal gas pressure — ideal gas law in the plenum + gap volume.
The simulation runs 600 s covering normal operation, depressurisation, and LOCA blowdown. Cladding failure (burst) is detected via an event function, after which the integration continues with failed-clad boundary conditions.
The solver is implemented in solve_thermal_hydraulics() in
07.Thermal.Hydraulics/thermal_hydraulics.py.
Physics Equations
ODE State Vector
The flat 1-D state vector packs:
where:
\(T_{\text{fuel}}\) — fuel nodal temperatures (K), shape
(nz, nf)\(T_{\text{clad}}\) — clad nodal temperatures (K), shape
(nz, nc)\(h_{\text{cool}}\) — coolant specific enthalpy (J/kg), shape
(nz,)\(\bar\varepsilon_p\) — effective (scalar) plastic strain
\(\varepsilon_p^{(r,\theta,z)}\) — plastic strain tensor components
Default mesh: \(n_z = 2\) axial levels, \(n_f = 20\) fuel radial nodes, \(n_c = 5\) clad radial nodes. Total DOFs: \(n_z(n_f + 5 n_c) + n_z = 2(20 + 25) + 2 = 92\).
Fuel Temperature
1-D radial heat conduction in a cylindrical UO2 pellet with volumetric heat generation:
Boundary conditions:
Centre (\(r = 0\)): adiabatic, \(\partial T / \partial r = 0\)
Surface (\(r = r_f\)): gap heat flux, \(q_{\text{gap}} = h_{\text{gap}} (T_{f,\text{surf}} - T_{c,\text{in}})\)
Finite-volume discretisation: The fuel pellet is divided into \(n_f\) concentric annular rings. Heat flows \(Q_j\) (W) are computed at the \(n_f + 1\) radial interfaces:
where \(A_{\text{mid},j} = 2\pi r_{\text{mid},j} \Delta z\) is the interface area at the radial midpoint between nodes \(j-1\) and \(j\), and \(A_{\text{gap}}\) is the deformed gap area. The temperature rate for node \(j\) is:
where \(m_j = \rho_f V_j\) is the ring mass and \(V_j = \pi(r_{j+1}^2 - r_j^2) \Delta z\) is the ring volume.
Material properties (MATPRO correlations [MATPRO2003] in matpro.py):
Clad Temperature
1-D radial heat conduction in cylindrical Zircaloy-4 cladding:
Boundary conditions:
Inner surface: gap heat flux from fuel
Outer surface: wall-to-coolant heat transfer, \(q_w = h_{\text{conv}}(T_{c,\text{out}} - T_{\text{cool}})\)
Finite-volume discretisation: Same structure as fuel, with \(n_c\) concentric rings. Heat flows at the \(n_c + 1\) interfaces:
The clad interface areas \(A_{\text{mid},j}^{\text{def}}\) use the deformed clad geometry (including thermal + elastic + plastic strains):
This is important for LOCA analysis where clad ballooning significantly increases the heat transfer area.
Material properties:
Coolant Enthalpy
1-D axial energy balance with upwind differencing:
where \(\dot{m}\) is the junction mass flow rate (kg/s), \(q_w\) is the regime-dependent wall heat flux (W/m2), and \(A_{\text{HX}} = 2\pi r_{\text{out}} \Delta z\) is the heat exchange area.
Discrete form with upwind advection at \(n_z + 1\) junctions:
The junction enthalpies use upwind values: \(h_{z-1/2} = h_{\text{inlet}}\) at the first junction, \(h_{z-1/2} = h_{z-1}\) at interior junctions. Water/steam properties (\(\rho, T, T_{\text{sat}}, h_L, h_V\), etc.) are evaluated via IAPWS correlations through pyXSteam at each axial node.
Gap Conductance
The gap between fuel and cladding is filled with helium. The conductance model depends on the gap state:
where \(k_{\text{He}}(T) = 2.639{\times}10^{-3} \, T^{0.7085}\) W/m-K, \(\delta_{\text{gap}} = r_{c,\text{in}} - r_{f,\text{out}}\) is the deformed gap width, and \(\varepsilon_{\text{rough}} \approx 6\,\mu\text{m}\) is the effective surface roughness for contact conductance.
The gap geometry is deformable: fuel outer radius changes with thermal expansion \(\varepsilon_T^{\text{UO}_2}(T)\), and clad inner radius changes with thermal + elastic + plastic strains.
Internal Gas Pressure
Ideal gas law for the helium inventory distributed between plenum and gap:
where \(\mu_{\text{He},0}\) is the initial helium mole count (conserved), \(V_{\text{gap},z} = \pi(r_{c,\text{in}}^2 - r_f^2) \Delta z_z\) is the deformed gap volume per axial level, and \(T_{\text{gap},z} = (T_{f,\text{surf}} + T_{c,\text{in}})/2\).
Clad Stress — Algebraic Subsystem
The cladding stress is solved as an algebraic (not differential) system at each RHS evaluation. For each axial node, a linear system of \(3 n_c\) equations is solved for \(\sigma_r, \sigma_\theta, \sigma_z\):
Radial equilibrium (\(n_c - 1\) equations):
\[\frac{\partial \sigma_r}{\partial r} + \frac{\sigma_r - \sigma_\theta}{r} = 0\]Strain compatibility (\(n_c - 1\) equations):
\[\frac{\partial \varepsilon_\theta}{\partial r} = \frac{\varepsilon_\theta - \varepsilon_r}{r}\]Axial strain uniformity (\(n_c - 1\) equations):
\[\varepsilon_z(r_j) = \varepsilon_z(r_{j+1})\]Boundary conditions (3 equations):
Outer surface: \(\sigma_r(r_{\text{out}}) = -p_{\text{cool}}\)
Inner surface (gap open): \(\sigma_r(r_{\text{in}}) = -p_{\text{gas}}\)
Axial force balance: \(\int \sigma_z \, r \, dr = p_{\text{gas}} r_{\text{in}}^2 - p_{\text{cool}} r_{\text{out}}^2\)
Strain decomposition:
Linear system construction: For each axial node \(z\), the \(3 n_c\) unknowns are ordered as \([\sigma_{r,1} \ldots \sigma_{r,n_c},\; \sigma_{\theta,1} \ldots \sigma_{\theta,n_c},\; \sigma_{z,1} \ldots \sigma_{z,n_c}]\). The coefficient matrix \(\mathbf{A}\) and RHS vector \(\mathbf{b}\) encode:
Rows 1 to \(n_c - 1\): equilibrium via central differences on the clad node midpoints, \((\sigma_{r,j+1} - \sigma_{r,j})/\Delta r + (\sigma_{r,\text{avg}} - \sigma_{\theta,\text{avg}})/r_{\text{mid}} = 0\)
Rows \(n_c\) to \(2(n_c - 1)\): strain compatibility relating elastic compliance coefficients \(C_s = 1/E\) and \(C_c = -\nu/E\) with non-elastic (thermal + plastic) strain contributions on the RHS
Rows \(2 n_c - 1\) to \(3(n_c - 1)\): axial strain uniformity
Row \(3(n_c - 1) + 1\): outer BC \(\sigma_r(r_{\text{out}}) = -p_{\text{cool}}\)
Row \(3(n_c - 1) + 2\): inner BC (pressure or strain compatibility)
Row \(3 n_c\): axial force balance
The system is dense \((15 \times 15)\) for \(n_c = 5\) and solved
via np.linalg.solve — no iterative method needed.
Gap-closed inner BC: When the gap is closed (Phase 3 of Module 08), the inner BC changes from \(\sigma_r(r_{\text{in}}) = -p_{\text{gas}}\) to a strain compatibility condition:
where \(\Delta\varepsilon_h\) is the hoop strain jump frozen at the moment of closure. This couples the fuel thermal expansion to the clad inner-surface stress.
Plastic Strain Rate — Norton Creep
The effective plastic strain rate follows a Norton power-law:
where \(\sigma_{\text{VM}}\) is the von Mises stress (MPa), and \(K(T), m(T), n(T)\) are temperature-dependent Zircaloy strength parameters from MATPRO. The deviatoric components follow the Prandtl–Reuss flow rule:
Two-Phase Flow Model
The wall heat transfer coefficient is selected by flow regime:
Regime 1 — Single-phase forced convection (\(T_w < T_{\text{sat}}\)):
Regime 2 — Pre-CHF nucleate boiling (\(T_{\text{sat}} < T_w < T_{\text{CHF}}\)):
Regime 3 — Post-CHF film boiling (\(T_w \ge T_{\text{CHF}}\)):
CHF — Modified Zuber correlation with subcooling and void corrections:
Clad Failure Criterion
Failure is detected when the engineering hoop stress exceeds the burst stress:
The burst stress \(\sigma_B(T)\) follows the Kingery–Hobson correlation from MATPRO. After failure, the internal gas pressure equilibrates with coolant pressure: \(p_{\text{gas}} = p_{\text{cool,outlet}}\).
LOCA Boundary Conditions
The LOCA scenario is defined by four time-dependent boundary condition tables, interpolated linearly at each RHS call:
Inlet temperature \(T_{\text{in}}(t)\) — drops during reflood
Outlet pressure \(p_{\text{out}}(t)\) — depressurises from 15.5 MPa
Inlet velocity \(v_{\text{in}}(t)\) — flow coastdown
Power \(P(t)\) — decay heat curve after scram
The LOCA sequence (default PWR parameters):
Time (s) |
Pressure (MPa) |
Velocity (m/s) |
Power (W) |
Phase |
|---|---|---|---|---|
0 – 200 |
15.5 |
4.8 |
69 141 |
Normal operation |
200 – 201 |
15.5 → 0.5 |
4.8 → 0.1 |
69 141 → 4 500 |
Blowdown initiation |
201 – 600 |
0.5 → 0.1 |
0.1 → 0.01 |
4 500 → 2 000 |
LOCA blowdown |
The exact table values are stored in the THParams dataclass (lines 63–100
of thermal_hydraulics.py).
Numerical Methods
ODE Integration
The system is integrated using scipy.integrate.solve_ivp with
method='BDF' (Backward Differentiation Formula), matching MATLAB’s
ode15s. Tolerances: rtol=1e-6, atol=1e-4, max_step=10.0.
The same BDF approach is used by Fuel Behaviour — 1-D Radial Thermo-Mechanics and
Reactor Kinetics — 0-D Point Kinetics + TH.
Chunked Integration for Event Detection
scipy’s built-in event detection (which uses brentq internally) fails
for this problem because the event function can have the same sign at both
ends of a solver step, causing ValueError.
Solution: The integration is split into chunks of one output step (\(\Delta t = 1\) s). After each chunk:
Evaluate \(E(t) = \sigma_B - \sigma_I\) at the chunk endpoint
If \(E\) changed sign (positive → negative), declare clad failure
The chunk endpoint is used as the failure state (1 s resolution). Dense-output bisection was attempted but rejected — the interpolated states produce unphysical water property values (NaN from pyXSteam) because the interpolant does not enforce thermodynamic consistency.
Post-failure integration also uses chunked integration with
try/except ValueError for graceful degradation if any residual NaN
issues arise.
IAPWS Viscosity Fallback
TH-20260405-001
Problem: During post-failure LOCA blowdown, coolant temperature in the
outlet node exceeds 900 °C. pyXSteam’s viscosity function
(my_AllRegions_ph) has an artificial cutoff:
% XSteam.m line 3410
if T > 900 + 273.15
my_AllRegions_ph = NaN;
return
end
This is a conservative validity bound, not a physical singularity. The
NaN viscosity propagates: nu = mu/rho → NaN → Churchill friction
factor → NaN → _compute_pressure returns NaN pressures → solver crash.
Fix: _iapws_viscosity(T_K, rho) in h2o_properties.py
implements the same IAPWS 2008 correlation without the cutoff.
The formula has two parts:
where \(\bar{T} = T / 647.226\) and \(\bar{\rho} = \rho / 317.763\).
Dilute-gas contribution (kinetic theory):
Finite-density correction:
The \(H_{ij}\) coefficient matrix (7 × 6, mostly sparse) is from the IAPWS release [IAPWS2008]. At low pressures (< 1 MPa) and high temperatures, \(\bar{\rho} \ll 1\) so \(\mu_1 \approx 1\) and the viscosity is dominated by the dilute-gas term.
Validation:
T (°C) |
pyXSteam μ (Pa·s) |
IAPWS μ (Pa·s) |
Ratio |
|---|---|---|---|
500 |
2.858e-5 |
2.858e-5 |
1.000000 |
700 |
3.656e-5 |
3.656e-5 |
1.000000 |
899 |
4.405e-5 |
4.405e-5 |
1.000000 |
900 |
4.409e-5 |
4.409e-5 |
1.000000 |
950 |
NaN |
4.590e-5 |
— (fallback) |
1000 |
NaN |
4.767e-5 |
— (fallback) |
1200 |
NaN |
5.450e-5 |
— (fallback) |
Bit-identical below 900 °C; smooth, physically reasonable above. The function could replace pyXSteam’s viscosity entirely (DA-20260405-005).
Validation
The MATLAB reference is in matlab_archive/09.Thermal.Hydraulics/results.m.
MATLAB Gap Geometry Bug
TH-20260405-002
MATLAB’s funRHS.m line 272 contains a bug:
gap.r_ = (clad.r(1) + fuel.dz)/2;
This adds a cladding inner radius (clad.r(1) ≈ 4.22 mm) to the fuel
axial node height (fuel.dz ≈ 1.5 m), producing gap.r_ ≈ 0.752 m
instead of the correct mid-gap radius ≈ 4.17 mm. The intended code was
likely (clad.r(:,1) + fuel.r)/2.
The gap heat transfer area is then:
versus the correct value:
A factor of 180× difference. This massively enhances gap heat transfer, keeping the fuel 332 °C cooler at steady state (808 °C vs 1140 °C).
Consequence for validation: Direct comparison of Python and MATLAB is
only meaningful using thermal_hydraulics_dae.py which replicates the
MATLAB gap geometry. The “correct physics” version (thermal_hydraulics.py)
has no MATLAB reference — it represents the physically correct solution.
Two Versions
Two solver files are maintained for comparison:
Property |
|
|
|---|---|---|
Gap geometry |
Correct: \(r_{\text{mid}} = (r_{c,\text{in}} + r_f)/2\) |
MATLAB match: \(r_{\text{mid}} = (r_{c,\text{in}} + \Delta z_f)/2\) |
Pressure |
Algebraic (computed in RHS) |
State variable with relaxation |
Fuel centre (SS) |
1140 °C |
808 °C (matches MATLAB) |
Clad failure |
287 s |
379 s (MATLAB: 425 s) |
Runs to t = 600 s |
Yes |
Yes |
MATLAB Parity (DAE Version)
Comparison of thermal_hydraulics_dae.py against
matlab_archive/09.Thermal.Hydraulics/results.m at key time points.
All values for outlet node (axial node 2):
Time (s) |
Python h (kJ/kg) |
MATLAB h (kJ/kg) |
Python T_cool (°C) |
MATLAB T_cool (°C) |
Python fuel_T (°C) |
|---|---|---|---|---|---|
1 |
1299 |
1300 |
292.8 |
292.9 |
426.6 |
200 |
1443 |
1443 |
318.4 |
318.4 |
808.2 |
210 |
2357 |
2357 |
283.0 |
283.0 |
368.1 |
250 |
3356 |
3358 |
438.9 |
440.0 |
462.4 |
287 |
3528 |
3527 |
519.5 |
518.9 |
523.0 |
Pre-LOCA steady state and early blowdown match to within 1–2 kJ/kg in coolant enthalpy and 1–2 °C in temperature.
Correct-Physics Steady-State Check
The correct-physics version (thermal_hydraulics.py) has no MATLAB
reference, but the steady-state fuel centre temperature can be checked
against an analytical estimate. For a solid cylindrical pellet with
uniform volumetric heating:
With LHGR = 69 141 / 3.0 = 23 047 W/m and \(k_{\text{UO}_2}\) ≈ 3.0 W/m-K (average over the 300–1100 °C range), this gives:
Adding the gap + clad + coolant resistance (~320 °C above coolant at 300 °C):
The Python solver gives 1140 °C, which is lower than this simple estimate because the actual \(k(T)\) increases at lower temperatures (the conductivity-averaged ΔT is smaller than the constant-k estimate). The agreement within ~8% confirms the correct-physics version is physically reasonable (TH-20260405-006).
The clad failure timing difference (Python 379 s vs MATLAB 425.3 s) is
traced to MATLAB’s additional indexing quirk in the gap temperature
calculation: gap.T = (fuel.T(fuel.nr) + clad.T(1))/2 with
fuel.nr = 20 returns element 20 of the column-major flattened
(2, 20) array — i.e., axial node 2, radial node 10 — not the fuel
surface temperature at each axial level. This produces a different inner
gas pressure (MATLAB ~2.9 MPa vs Python ~4.5 MPa at t = 400 s), which
changes the engineering stress and delays failure.
Investigation History
The NaN at t ≈ 395 s during post-failure LOCA integration was investigated in three phases over sessions 2026-04-01 and 2026-04-05.
Phase 1 — Initial Diagnosis (2026-04-01)
Identified that post-failure Phase 2 integration crashed with
ValueErrorTraced to pyXSteam returning NaN for water/steam properties
Hypothesis: BDF Jacobian numerical differencing perturbs the state vector, pushing enthalpy/pressure outside pyXSteam’s valid range
Mitigation: Chunked integration with
try/exceptfor graceful stopResult: Simulation runs to t = 395 s, stops cleanly
Phase 2 — Root Cause Chain (2026-04-05)
Systematic investigation revealed the full causal chain:
Added diagnostics to the RHS function to print state values at NaN occurrence:
Coolant enthalpy at failure: h = [3172, 4399] kJ/kg — valid for pyXSteam at this pressure
Pressure at failure:
cool_p = [NaN, NaN]— the NaN comes from_compute_pressure, not from pyXSteam directly
Traced NaN through pressure computation: First-pass water properties showed
nu = [1.93e-5, NaN]. Node 2 kinematic viscosity was NaN. Density and temperature were valid.Identified pyXSteam viscosity cutoff:
p_bar, h_kJ = 3.31, 4399.0 T = st.t_ph(p_bar, h_kJ) # 900.7 °C — works fine rho = st.rho_ph(p_bar, h_kJ) # 0.611 kg/m³ — works fine mu = st.my_ph(p_bar, h_kJ) # NaN — 900 °C cutoff!
The cutoff is exactly 900 °C. All other properties (T, ρ, k, c_p) work to well above 1000 °C.
Why Python reaches 900 °C but MATLAB doesn’t: Compared fuel centre temperatures:
Python: 1140 °C at steady state (correct physics)
MATLAB: 808 °C at steady state
The 332 °C difference comes from the gap geometry bug (see MATLAB Gap Geometry Bug). The hotter fuel stores more energy; during LOCA blowdown this energy heats the coolant past 900 °C.
Phase 3 — Approaches Tried and Outcome
# |
Approach |
Result |
Why |
|---|---|---|---|
1 |
Reduce |
Failed |
Coolant genuinely reaches >900 °C; smaller steps don’t prevent it |
2 |
DAE with pressure as state variable + relaxation |
NaN persisted |
Pressure relaxation doesn’t prevent viscosity NaN |
3 |
Replicate MATLAB gap geometry |
NaN eliminated |
Fuel stays at 808 °C; coolant never reaches 900 °C |
4 |
IAPWS viscosity fallback |
NaN eliminated |
Addresses root cause; both versions run to 600 s |
Approach 4 was adopted as the primary fix because it preserves correct
physics without replicating the MATLAB bug. Approach 3 is retained
in thermal_hydraulics_dae.py for parity comparison.
Known Limitations
Clad failure time resolution: 1 s (chunk size). Could be refined via re-integration bisection but not implemented (TH-20260401-006).
No radiation heat transfer across the gap (gas conduction only).
Axial mesh coarse: Only \(n_z = 2\) nodes; no axial conduction.
DAE clad failure timing:
thermal_hydraulics_dae.pyfails at 379 s vs MATLAB’s 425.3 s. The remaining 46 s gap is due to MATLAB indexing quirks in gap temperature and inner gas pressure (not worth replicating further).
References
D.L. Hagrman et al., MATPRO — A Library of Materials Properties for Light-Water-Reactor Accident Analysis, NUREG/CR-6150, Idaho National Laboratory, 2003.
IAPWS, Release on the IAPWS Formulation 2008 for the
Viscosity of Ordinary Water Substance, 2008. The correlation
implemented in _iapws_viscosity and in XSteam/pyXSteam.