This post is a condensed tour of my doctoral dissertation, Application of High Fidelity Methods to Nuclear Engineering Thermal-Hydraulics Problems: Transient Systems, Method Scalability, and Closure, completed at The George Washington University in 2026 under the direction of Professor Elias Balaras. I've tried to distill the better part of a decade of research into something a technically curious reader can absorb in a single sitting — with enough equations to keep it honest and enough narrative to keep it human.
The throughline is this: direct numerical simulation (DNS) — the practice of solving the Navier-Stokes equations without any turbulence modeling — can serve as a credible numerical laboratory for generating the physical knowledge that nuclear reactor safety analysis desperately needs. But getting there required tearing down legacy software, confronting a surprising numerical failure, and building new physics from first principles.
A Love Letter to Fortran
Let me start with what might be an unpopular opinion in 2026: Fortran is the right language for this work, and I don't think it's close.
The in-house DNS tools I inherited were written in Fortran, relying on software design practices from the 1990s — or worse, the 1970s. Common blocks everywhere: globally visible, mutable state accessible from anywhere in the program, with the footprint of potentially dead code that could interfere with active memory spanning practically the entire codebase. Two separate legacy tools — one a bespoke solver built on the FLASH4 framework, the other a monolithic in-house code I'll call Eddy v1.0 — had fractured into large numbers of incompatible, orphaned branches in version control. The codes worked, in the sense that they produced physically meaningful results, but neither provided a foundation that could be restored to a singular, verified and validated state within any practical scope of work.
The obvious question was: why not rewrite in C++ or Python or Rust? The answer comes down to what we actually need these tools to do. The core computation is the time integration of the incompressible Navier-Stokes equations on a staggered grid:
$$\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} = -\nabla P + \frac{1}{\sqrt{Ra/Pr}} \nabla^2 \vec{u} - \hat{g} T$$
$$\frac{\partial T}{\partial t} + \vec{u} \cdot \nabla T = \frac{1}{\sqrt{Ra \cdot Pr}} \left( \nabla^2 T + S \right)$$
These equations are discretized with second-order central finite differences on structured grids and time-advanced with a semi-implicit projection method — third-order Runge-Kutta for the explicit part, second-order Crank-Nicolson for the implicit part. The staggered grid arrangement with central differences ensures the discrete momentum equation is kinetic-energy-preserving in the inviscid limit, which means minimal numerical diffusion even across the $10^6$ to $10^8$ integration steps observed in this work.
Fortran expresses these mathematical operations on multidimensional arrays transparently — possibly uniquely among compiled languages, it states the necessary constructs in their simplest form and lets the compiler generate the optimized machine code. C++ can certainly match it for raw speed, but subtle complexities around pointer aliasing (the language lacks C's restrict keyword) and the absence of native multidimensional array facilities until std::mdspan in C++23 mean that achieving equivalent optimizations requires libraries like Kokkos and non-trivial template machinery that most researchers (read: graduate students) will not maintain.
The Clean-Sheet Refactor
Rather than attempting to reconcile the two divergent legacy codebases, I performed a clean-sheet refactoring of all validated methods and algorithms into a new tool, Eddy v1.5 — Fortran 2018 standard compliant source, leaning on submodules from the 2008 standard and type-bound procedures from 2003. The key design decisions were:
Modules and submodules for scoped memory. Legacy tools used common blocks — globally visible, mutable state accessible from anywhere. The refactored tool uses Fortran module and submodule features to colocate data with the algorithms that operate on it. Member data are not modifiable by external code but have well-defined, wide exposure. The one exception is the physical vector and scalar fields themselves (velocity, pressure, temperature), where the added complexity of full encapsulation wasn't worth the trade-off for a tool used primarily by researchers.
Upfront memory allocation. All dynamic allocation occurs during problem setup, before the time-integration "hot loop" begins. The legacy tools allocated memory on every call to the Poisson solver — invoked once per Runge-Kutta sub-step, three times per time step — and to the boundary condition routines, invoked dozens of times per time step. Refactoring these allocations out of the critical path and into persistent, setup-time initialization produced measurable performance gains.
A plugin system without polymorphism. New physics capabilities are injected through a linker-based plugin architecture. The main codebase and a default plugin library are compiled as separate object libraries; user-authored Fortran plugins can replace default implementations at link time without modifying the core solver. This avoids the build-time polymorphism of one legacy tool (which used Python to assemble compilation units from a hierarchy of source versions) and the runtime polymorphism that would require object-oriented abstractions most researchers won't maintain.
On code reduction: The refactored Eddy v1.5 codebase achieved near feature parity with the legacy tools — deliberately excepting multiphase flow, which is being rebuilt on the new architecture — at an $\sim 88\%$ reduction in source lines. The legacy Eddy v1.0 was approximately $1.4 \times 10^5$ lines; FLASH4, the other legacy tool, measured $6.1 \times 10^5$ lines at rest and $9.7 \times 10^4$ lines as assembled for compilation. Eddy v1.5 accomplishes this in $1.65 \times 10^4$ lines. Fewer lines means fewer places for bugs to hide, faster onboarding for new students, and more efficient maintenance.
The build system uses GNU Make with a structured hierarchy of Makefile fragments. A selection of build (debug, optimized, release), stack (GNU+OpenBLAS, Intel+MKL), and patch (user plugins) concretizes a single compilation instance. The refactored build system compiles in roughly half the wall time of the legacy tool, and by the time 8 parallel build processes are in play with full optimizations, the legacy tool is a factor of 16 slower due to its inability to parallelize the build effectively.
Runtime performance improved comparably. Across test problems — Rayleigh-Bénard convection, an immersed cylinder in channel flow, turbulent channel flow with an immersed cavity — the refactored tool demonstrates a parallel fraction of ~80–85% (in the sense of Amdahl's law), where the legacy tool varies and is generally limited to less than ~75%. Measured wall-time ratios range from 1× to 15× faster depending on the problem and build configuration, with a single contrived exception: a debug build running the legacy FFT solver, an improbable research configuration.
The Engineering Problem
Now let me explain why any of this matters.
Nuclear reactor safety analysis depends on predicting what happens when things go wrong. The most severe class of accident — a core melt — results in the nuclear fuel rods overheating, melting, and relocating downward through the reactor vessel. At temperatures around 2100 K, the zirconium cladding melts and core material begins relocating. Near 3000 K, the oxidic fuel mixture (roughly 65% $UO_2$, 35% $ZrO_2$) itself liquefies. This degraded core material — a complex mixture of predominantly $U$-$Zr$-$O$-$Fe$-$C$-$B$ — eventually collects in the reactor vessel lower plenum and forms a contiguous, internally heated, naturally circulating molten pool.
The question that reactor safety analysts must answer is: will the reactor vessel survive?
In-Vessel Retention
The strategy of in-vessel retention relies on flooding the cavity outside the reactor vessel with water, providing external cooling to the vessel wall while the molten core attacks it from the inside. The feasibility of this strategy depends entirely on whether the heat flux from the molten pool to the vessel wall remains below the critical heat flux that would cause the external cooling to fail.
The heat transfer from the pool to the vessel is governed by natural convection — buoyancy-driven flow resulting from internal heating by fission product decay. The governing correlations take the form:
$$Nu = C \cdot Ra'^{n} \cdot Pr^{m}$$
where $Nu = hL/k$ is the Nusselt number (dimensionless heat transfer coefficient), $Ra' = \beta g \dot{q} L^5 / \nu \alpha^2$ is the modified Rayleigh number for internally heated flow, and $Pr = \nu / \alpha$ is the Prandtl number. These correlations come from decades of experimental programs — COPO in Finland and France, UCLA and ACOPO in the United States, BALI in France, RASPLAV in Russia, SIMECO in Sweden, LIVE in Germany, SIGMA in South Korea — using simulant materials in scaled geometries.
The experiments are good. The question is whether they're applicable.
The Transient Problem
The critical gap is this: experimental correlations are derived from steady, fully developed flow in idealized geometries. But a real reactor accident is anything but steady. The complete relocation of the reactor core involves between $10^2$ and $10^3$ discrete relocation events over a period of hours. Each event perturbs the molten pool — injecting momentum, changing the pool height, altering the Rayleigh number.
The time for the molten pool to reach steady convection can be estimated from the large-scale convective time constant:
$$t^* = \frac{L^2}{\alpha \sqrt{Ra' \cdot Pr}}$$
Using prototypic parameters — pool height $L \approx 0.8$ m, modified Rayleigh number $Ra' \sim 10^{13}$–$10^{16}$, and a Prandtl number $Pr \sim 0.1$ for a melt between metallic and oxidic composition — the time constant falls in the range of 0.005 to 25 seconds. Transient DNS show that full development requires 3–7 large-scale convective rolls, each taking 1–3 dimensionless time units, which puts the dimensional time to reach steady, fully developed flow between about 0.25 seconds and 3 minutes.
This means the molten pool likely spends a significant fraction of its life recovering from perturbations rather than existing in the steady state that aligns with the experimental correlations. If heat transfer during these transient phases exceeds the steady-state predictions, current safety margins could be non-conservative.
Why this matters for regulation: Several severe accident management strategies for both existing and proposed reactor designs — including the GE Hitachi ESBWR and NuScale VOYGR — credit in-vessel retention as a means of limiting accident consequences. The confidence in these strategies depends directly on the applicability of steady-state heat transfer correlations to inherently transient conditions.
DNS as a Numerical Laboratory
Direct numerical simulation resolves all turbulent scales from the largest energy-carrying eddies down to the Kolmogorov microscale — no turbulence modeling, no empirical closures for the fluid dynamics. For the class of problems I've described, this is not a luxury; it's a necessity.
The alternative high-fidelity methods — RANS (Reynolds-Averaged Navier-Stokes) and LES (Large Eddy Simulation) — both require closure models for unresolved scales. The applicability of these closures in complex multiphase flows is often uncertain, and their validation requires the same kind of prototypic experimental data that is difficult to obtain. DNS sidesteps this entirely: if the numerical procedure is properly verified against benchmark problems, the results are as trustworthy as the governing equations themselves.
Why Immersed Boundary Methods?
The choice of numerical method matters enormously for computational efficiency. I selected structured-grid, finite-difference methods with immersed boundary treatment of complex geometries, based on a clear quantitative advantage: for the same accuracy in integral quantities, immersed boundary methods require conservatively one-third the computational cost of higher-order and general-purpose finite-volume methods in purpose-built benchmarks — even for a stationary body using legacy academic tooling.
The advantages are structural:
- Fast Poisson solvers. The structured grid exploits the block-structured property of the Laplace operator, enabling direct solvers based on Fast Fourier Transforms and cyclic reduction — significantly faster than iterative Krylov-subspace or multigrid methods at scale.
- Kinetic-energy preservation. The staggered grid with central differences ensures minimal numerical diffusion, critical for maintaining accuracy over the millions of time steps required.
- Geometric flexibility. The immersed boundary decouples the fluid grid from the geometry of the body, avoiding the expensive re-meshing required by body-conformal approaches at every sub-iteration of the time advancement — particularly advantageous for moving or deforming boundaries.
The resolution requirements for DNS are stringent:
$$h \lesssim \pi \eta_k = \pi L \left( \frac{Pr^2}{Ra \cdot Nu} \right)^{1/4}, \quad Pr \leq 1$$
where $\eta_k$ is the Kolmogorov length scale. Extrapolated to conditions prototypic of contiguous molten corium pools, this translates to approximately $10^{13}$ computational grid points, with $\sim 10^2$ points across the boundary layer alone. We are not there yet — but the path is clear, and the scaling is favorable.
Transient Convection and Safety Margins
This is the central scientific result of the dissertation, and it carries a surprising twist.
I conducted extensive series of DNS — both two-dimensional and three-dimensional — of Rayleigh-Bénard convection and its internally heated analog, studying transitions between steady convective states across four orders of magnitude in the controlling Rayleigh number ($Ra = 10^6$ through $10^{10}$, with two-dimensional statistical sets extending to $10^{11}$).
The Overshoot in Rayleigh-Bénard Convection
For canonical Rayleigh-Bénard convection (heated from below, cooled from above), transitions from a lower to a higher $Ra$ produced a statistically significant overshoot in the boundary heat transfer. Consider an instantaneous transition from steady three-dimensional convection at $Ra = 10^8$ to $Ra = 10^9$. The naive expectation is that $Nu$ would monotonically approach the new steady value. Instead, the transient $Nu$ peaked at approximately 170% of the final steady-state value — far outside the roughly ±20% band of the steady state's 95th percentile. Across the full set of transition experiments, peak transient $Nu$ reached 1.5–3.0 times the final steady value depending on the correlating factors.
This overshoot is not a numerical artifact. I verified this through multiple lines of evidence:
-
Statistical significance. The overshoot persists across roughly eight statistically independent initial flow field realizations at each transition. In three dimensions, the 5th-percentile transient peak exceeds the 95th-percentile steady state by a factor of ~1.5 (in two dimensions, ~1.5–1.8), and the enhancement appeared in every two-dimensional transition set from $Ra = 10^6$ up through $10^{11}$.
-
Impulse sensitivity. Varying the transition window width over five orders of magnitude reveals behavior characteristic of a dynamical system — underdamped, critically damped, and overdamped responses, with critical damping at a window of roughly six convective time constants. The transient $Nu$ becomes bounded by the adjoining steady states only when the perturbation is imposed over a timescale about an order of magnitude longer than the large-scale convective period.
-
Resonance. Cyclic transitions between two $Ra$ states at a period near ~4× the convective time constant produce a resonance-like amplification, with peak $Nu$ stabilizing above 200% of the final steady value. No other cyclic period — and no single realization that happened to initiate at the right phase — reproduced this excess, implying a genuine resonance rather than a statistical aberration.
The Hysteresis
A clear directional asymmetry exists. Forward transitions (low $Ra$ to high $Ra$) produce significant overshoots; backward transitions produce only an insignificant undershoot. The physical explanation traces to the transient volume of thermal plumes: upon a rapid increase in $Ra$, the boundary layer plumes from the previous state are "oversized" relative to the new equilibrium — steady convection at the higher $Ra$ would have ejected them before they grew so large. They carry excess buoyant force into the bulk flow before they can shrink, producing a transient enhancement in convective heat transfer. Plumes in the reverse direction begin too small, and their growth takes time during which there is no analogous coherent acceleration.
The Good News: Internally Heated Convection is Stable
Here is the result that matters most for nuclear safety: internally heated convection does not exhibit these overshoots.
Despite a wide permutation assessment — varying the method of transition, the phase of initiation within the oscillatory steady $Nu$, cavity aspect ratios, and boundary conditions — the transient $Nu$ for internally heated flows remained bounded by the adjoining steady states. No enhancement was identified.
Implications for reactor safety: The internally heated analog is the physically representative model for a molten fuel pool with fission product decay heat. The absence of transient overshoots suggests that current steady-state correlations used in severe accident codes like MELCOR may be more resilient than previously feared — at least for internally heated pools. The Rayleigh-Bénard overshoot remains relevant for pools where internal heat generation is minor compared to boundary heating, such as the light metallic layer.
Beyond the physics, the study distills its results into an engineering closure: an empirical model for the transient $Nu$ — built from the adjoining steady states, a transition-window factor, and a cyclic-period factor — offered as concrete modeling guidance for the molten pool treatments in system codes.
This finding does not close the book. The internally heated simulations were limited to two dimensions by the available compute allocation, and three-dimensional effects cannot be ruled out. Geometric perturbations — the impact of solid debris structures within the pool — remain unexplored. But the absence of overshoots in internally heated convection is a meaningful reduction in uncertainty for in-vessel retention safety analysis.
It is worth pausing on cost, because it makes the "numerical laboratory" argument concrete: the entire simulation campaign consumed roughly 8.5 million core-hours and under 30 TB of storage — an equivalent capital cost on the order of USD 100K at typical allocation rates. That is economical compared to the facilities cost of an experimental program asking the same questions. Extending the full transition study to $Ra = 10^{11}$, however, was estimated at another five million CPU-hours and fifteen months of wall time — which is what motivated the next chapter of the work.
The Poisson Solver Challenge
Every incompressible flow solver must solve a pressure Poisson equation at every sub-step to project the velocity field into a divergence-free space:
$$\nabla^2 P = \frac{1}{\Delta t} \nabla \cdot \vec{u}^*$$
In our fully resolved Rayleigh-Bénard simulations, this solve accounts for 30–45% of total wall time. On uniform grids, Fast Fourier Transform–based solvers are unbeatable: the transforms cost $\mathcal{O}(N^3 \log N)$, where evaluating the discrete transforms without a fast algorithm would require $\mathcal{O}(N^4)$ dense matrix-vector products. But practical DNS wants stretched grids: fine resolution near walls where boundary layers are thin, coarser resolution in the bulk. Uniform grids waste enormous numbers of grid points.
The Promise of NuFFTs
Non-uniform Fast Fourier Transforms (NuFFTs) have matured significantly in recent years, driven by applications in image processing and signal analysis. The idea is straightforward: if we can perform fast Fourier transforms on non-equidistantly sampled data, we can solve the Poisson equation on stretched grids while retaining the performance characteristics of FFT-based methods. The appeal extends to hardware: FFT libraries come vendor-tuned for GPUs, which is exactly the portability story a solver needs on accelerator-dominated machines.
The approach requires scaling the source distribution by a weight function — the quadrature weights from the local Jacobian of the grid transformation:
$$\phi(x,y,z) = \mathcal{F}\left\{ \frac{\mathcal{H}\{\xi(x,y,z) \rho(x,y,z)\}}{i^2(\omega_p^2 + \omega_q^2 + \omega_r^2)} \right\}$$
where $\mathcal{F}$ is the forward NuFFT, $\mathcal{H}$ is the adjoint NuFFT (playing the role of the inverse), and $\xi$ is the weight function. For a regular grid with orthogonal stretching, $\xi_{l,m,n} \propto \Delta x_l \Delta y_m \Delta z_n$.
Performance results were tantalizing: the NuFFT-based solver computed solutions approximately 20× faster than an equivalent uniform-grid FFT solver while keeping the relative error under 0.25% (the speedup falls to about 3.25× if the error is held under 0.025%). The gain comes from two competing factors: the NuFFT itself is 1.5–4× slower than a standard FFT on the same grid, but stretching allows the same accuracy with vastly fewer grid points — in the benchmark case, $4^3 = 64$× fewer.
The Fundamental Failure
Here is the cautionary tale. Despite the accuracy and performance of the NuFFT-based solver for the solution to the Poisson equation, it fails as a projection operator in the time-integration scheme.
The issue is the residual. For the projected velocity correction to produce a divergence-free field, the solver must be consistent with the finite-difference operators in both the spatial and frequency domains. This consistency is enforced through modified wave numbers that match the finite-difference stencil. On a uniform grid with consistent modified wave numbers, the residual sits at machine precision (~$10^{-12}$%).
On a non-uniform grid, this consistency breaks down — and the residual does not converge as the stretching is reduced. Even at a maximum cell-over-cell stretching ratio of $10^{-12}$%, vanishingly close to uniform, the residual remains above $10^{-10}$% — more than 100× larger than the equidistant grid. In practice, instability of the flow solver is observable within a single time step at even the most modest stretching.
The root cause is a loss of separability. On a non-uniform grid, the finite-difference Laplacian cannot be exactly diagonalized by the NuFFT basis functions: the product of the Fourier matrix and the difference operator becomes poorly conditioned, and its inverse cannot be computed accurately in double precision. The unresolved frequencies do not alias to integer multiples of the base frequency (as they do on uniform grids) but instead spread across many frequency components — producing a characteristic "fringe-like" oscillatory residual pattern.
The lesson: accuracy and residual are different properties, and CFD needs both. Discretization error measures whether the approximate solution is close to the true one; the residual measures whether it satisfies the discrete system. A Poisson solver can be accurate and fast, but the divergence-free projection is a mathematical constraint, not an approximation — when the residual is numerically non-trivial, the corrected velocity retains divergent characteristics, and the time evolution blows up. Numerical efficiency cannot be pursued at the expense of fundamental mathematical properties.
The dissertation's own summary of this chapter is bleakly honest: success was limited to "an irrelevant, accurate, and performant solver; a relevant, accurate, and nonperforming analog; and innumerable inaccurate methods regardless of performance or relevance." Within finite-precision arithmetic, the impediment appears fundamental — unlikely to yield to algebraic trickery. The practical state of the art remains the traditional compromise: Fourier transforms along uniform coordinate directions, direct block-tridiagonal solves along the stretched ones.
Sub-Grid Closure for Nucleate Boiling
The final contribution shifts from the macro-scale convection problem to the micro-scale: developing a mechanistic model for the vapor bubble nucleation cycle in nucleate boiling.
Interface-resolving multiphase DNS can now simulate individual bubble growth, departure, and interaction with remarkable fidelity. But there is a fundamental scale gap: the micro-cavities on a heater surface from which vapor bubbles nucleate are roughly $10^9$ times smaller than the macroscopic fluid domain of interest. No practical DNS can resolve both scales simultaneously.
The standard approach in the literature is to apply an empirically determined, globally constant nucleation frequency to all activated nucleation sites. But no systematic investigation exists into whether a global average frequency produces the same macroscopic heat transfer as the emergent average of locally varying frequencies. This uncertainty is not just academic — nucleation frequency directly affects the predicted surface heat flux, which in turn determines critical heat flux margins.
Bubble Waiting Time
I developed a physically motivated analytical model for the vapor bubble waiting time, building on the classical work of Han and Griffith (1965) with a critical modification: incorporating the static contact angle at the fluid-surface triple contact line.
The classical model assumes a hemispherical bubble nucleus, yielding a simple mechanical equilibrium $\Delta P = 2\sigma / R_c$. This neglects the effect of surface wettability. The modified model uses the Young-Laplace equation with geometric relationships for the contact angle on an idealized elliptical conical cavity:
$$\Delta P = \sigma \left( \frac{\cos \varphi}{R_c} + \frac{\cos \varphi}{R_c / \Gamma} \right)$$
where $\varphi$ is the spherical-cap contact angle (set by the fluid-surface wetting angle and the cavity geometry), $R_c$ is the cavity radius, and $\Gamma$ is the cavity aspect ratio. Combined with the Clausius-Clapeyron relation and a thermal boundary layer growth model, the waiting time criterion for a constant wall temperature (Dirichlet condition) becomes:
$$t_w \geq \frac{9}{4\pi\alpha_l} \left( \frac{(T_w - T_\infty) R_c}{T_w - T_{sat}\left(1 + \frac{\sigma}{\rho_v h_{lv}} \frac{\cos\varphi}{R_c}(1+\Gamma)\right)} \right)^2$$
I also extended this to fixed heat flux (Neumann) boundary conditions, replacing the exponential integral form of the wall conduction solution with its Hermite polynomial analog:
$$t_w \geq \frac{\pi}{4\alpha_l} \left( \frac{3}{2}R_c + \frac{k_l}{q_w''} \left( \frac{\sigma T_{sat}}{\rho_v h_{lv}} \frac{\cos\varphi}{R_c}(1+\Gamma) - (T_\infty - T_{sat}) \right) \right)^2$$
Contact Angle as the Dominant Parameter
Validation against single-bubble boiling experiments in water and FC-72 showed the predicted waiting time is a very strong function of the surface contact angle and only a weak function of the cavity and surface aspect ratios. Fitting the reported cavity sizes and waiting times recovered contact angles of 82–85° for water — squarely within the 80–90° reported for water on the experiments' heater materials — and a plausible ~10° for FC-72, for which no measured value was reported. Identifying the contact angle as the dominant parameter is a useful result in itself: it is a measurable property of a fluid-surface pair, whereas cavity geometry is often poorly characterized in experiments.
On practical applicability: The model is intentionally analytical rather than computational. It must be evaluated every time step across dozens of active nucleation sites for practical simulations. An empirical correlation tuned to a specific fluid-surface pair would be more accurate but less general; a resolved micro-scale simulation of each cavity would be more fundamental but computationally prohibitive. The analytical model occupies the useful middle ground — general, efficient, and stable.
Bubble Growth
For the initial inertially-controlled growth phase, a simplified control volume energy balance relates the growth rate to the heat transfer from the superheated boundary layer and the wall:
$$\frac{\partial V}{\partial t} = \frac{q_{tbl}'' A_{tbl} + q_w'' A_w}{\rho_v h_{lv}}$$
Under the approximation that the bubble base radius is much smaller than the bubble radius ($a \ll r$), this yields a separable differential equation with a transcendental solution, whose roots are found efficiently and stably with simple iterative algorithms — preserving the runtime efficiency needed across many nucleation sites. Satisfyingly, the Neumann formulation reproduced measured growth histories for water without being told the heater surface temperature, where the Dirichlet formulation under-predicted the growth rate even with it.
Future Directions
The work presented here establishes a foundation, but several important extensions remain.
Multiphase integration. Eddy v1.5 deliberately shipped without the legacy multiphase capability, and the framework is architecturally ready for the full three-component model (liquid, vapor, non-condensable gas) needed to attack most of the thermal-hydraulic challenge problems identified in the dissertation. Interface-resolving methods based on level-set approaches are well-validated and ready to implement; the engineering work is in robustness and communication cost — the third-order WENO scheme that keeps the advected level-set sharp demands three guard cells per decomposed domain boundary, and that halo exchange already accounts for up to ~30% of wall time at scale.
Hardware accelerators. Leadership-class machines like Aurora and Frontier derive the majority of their compute from GPUs, and the structured nature of this computation ports straightforwardly — with two exceptions: the fast Poisson solver (an algorithmic challenge) and the parallel decomposition of the Lagrangian body (a memory and communication challenge). The NuFFT approach was the GPU-portable answer to the first, and its failure as a projection operator sends us back to adapting the traditional methods.
Prototypic geometry. The molten pool studies in this work used the canonical parallel-plate geometry. Extension to the curved lower-head geometries of real reactor vessels — and to pools laced with solid debris structures — is exactly what the immersed boundary capabilities of the refactored solver are for. Favorably, the ratio of boundary-layer to bulk resolution requirements approaches only ~5× at prototypic $Ra'$ for internally heated flows, versus ~200× for Rayleigh-Bénard convection at comparably extreme conditions.
Boiling closure refinement. The current sub-grid model over-predicts late-phase growth rates for some fluids and needs a more complete representation of the superheated microlayer beneath the growing bubble — its finite altitude of adherence, its finite volume, and a closure relating the apparent to the microscopic contact angle along the triple line. The existing model is also numerically stiff for fluids with contact angles near 90°, such as water.
The regulatory horizon. New reactor designs increasingly credit passive physical phenomena as the first line of defense, and their designers have less tolerance for the conservatism that semi-empirical closures in system codes like TRACE, RELAP, and MELCOR impose on margins. A DNS campaign that answers a targeted closure question for the cost of a mid-six-figure allocation is a compelling alternative to a prototypic experimental program when a regulatory threshold of confidence must be met.
Closing Thoughts
The transition of nuclear engineering from a discipline of prototypic mockups to one of numerical experimentation is well underway. The computational costs of DNS remain high — the molten pool campaign alone consumed some 8.5 million core-hours on the Expanse and Stampede2 systems, through an XSEDE (now ACCESS) allocation — but the physical insights gained from resolving the smallest scales of turbulence and heat transfer are invaluable.
By building maintainable software, understanding the limits of our numerical tools, and bridging the gap between microscale surface physics and macroscale fluid dynamics, we move closer to a future where nuclear safety is underpinned by the most fundamental laws of physics rather than empirical approximation alone.
If you're interested in the full details, the dissertation is available through The George Washington University library. For questions about the work, the refactored solver, or potential collaborations, don't hesitate to reach out.
Aaron Lentner holds a B.S. in Nuclear Engineering from Purdue University, an M.E. in Nuclear Engineering from Penn State, and a Ph.D. from The George Washington University. His research focuses on high-fidelity computational methods for nuclear reactor thermal-hydraulics.