Inspiration

The fuel in a particle bed nuclear thermal rocket engine consists of a bed of fine particles through which the hydrogen passes directly. This allows for a huge heat transfer surface in a very small core; however, it has the disadvantage of being subject to a certain kind of failure, called thermal-hydraulic instability. When one part of the bed runs a bit too hot, the hydrogen around it becomes less dense and more viscous and thus more resistant to flowing. Less hydrogen flows, making that part even hotter.

This is supposed to be prevented by the use of an open wall structure known as a cold frit. As per the NASA analysis in 1993, a bed without a cold frit was considered unstable in low power density conditions but further analysis required more detailed calculation based on specific designs. Since we have similar interest in nuclear engineering field, we decided to do just that by simulating one published fuel element and checking the stability using only the published data.

How we built it

Physical model. The fuel element consists of 20 separate flows that go through the cold frit, the particle bed, and the hot frit. Each flow passes through a common inlet plenum and a common outlet channel; however, there is variation in pressure in both channels, and therefore each slice has its own pressure drop:

$$ \Delta p_i(\dot m_i, T_i) = p_{\text{plenum}}(z_i) - p_{\text{channel}}(z_i), \qquad \sum_{i=1}^{20} \dot m_i = \dot m_{\text{total}} $$

The above differences increase the flow rate in certain slices relative to others. This phenomenon was called the flow shift and became the most important factor for our solution.

Hydrogen properties. Hydrogen used in this project varies between 200 K and above 3000 K, so there isn’t any single source providing all needed properties. In our implementation we mixed CoolProp and NASA (Svehla) transport properties and NIST-JANAF tables between 200 and 5000 K ensuring continuity of values and slopes between every seam.

Solver. Coupled solver solves nonlinear problems regarding flow split, pressure drop, and heat transfer ensuring mass, momentum, and energy conservation. If the direct method fails to converge, it tries continuation technique described before. Case is considered unsolved when all options are explored and the continuation technique fails.

Two independent stability tests.

Ledinegg (static) test. The slice is considered statically unstable if its pressure drop decreases when its flow increases:

$$ \frac{\partial\, \Delta p_i}{\partial \dot m_i} < 0 \;\Rightarrow\; \text{unstable} $$

Eigenvalue test. The linearization of the fuel-temperature dynamics about the steady-state condition yields a Jacobian whose eigenvalues determine stability:

$$ \frac{d\,\delta \mathbf{T}}{dt} = \mathbf{J}\,\delta \mathbf{T}, \qquad \text{stable} \iff \max_k \operatorname{Re}(\lambda_k) < 0 $$

The Jacobians for both tests are validated using brute-force finite differences.

Sweeps and interface. In order to identify the conditions when stability is lost, sweeps are performed using reactor power (from 5% to 100%), cold-frit resistance, power distribution and single slice temperature disturbance. Using a Streamlit app, you can launch the solver in real time with any input parameters and visualize validation results, stability plots, a runaway animation, and the comparison of the solver with the surrogate model generated by training.

Verification and validation. 101 tests (pytest), grid convergence (20×20 is within 0.02% of 40×40), and comparison with published results. All pass/fail criteria and all sources of input are committed into the repository prior to the first validation. All later decisions are logged as amendments prior to the run in question. All result files carry the Git commit ID associated with them.

Results

Validation (Tier 1): 7 out of 8 pre-registrations checks pass for two published fuel elements. Tsinghua element pressure drop was calculated as 1060.8 and 1225.3 kPa (outlet and inlet based) compared to 1189.2 kPa from the published value within ±25% rule. Ji et al. element pressure drop was calculated as 512.0 kPa compared to 529.8 kPa.

The single fail, V7. With no cold frit, our model will be in the steady state regime only at the lowest power level (1 MW/L), and it will be unstable at it. In all higher power regimes, the flow redistribution will cause the starvation of the bottom fuel slice up to exceeding the maximum hydrogen temperature limit (5000 K). Thus, the instability regime will start already at 1 MW/L that is below the 1.33 – 12 MW/L interval established previously. We included the fail. Both stability checks agreed on it.

Our temperatures are higher than those from the papers. Prior to coding, we verified the power and flows against NIST-JANAF hydrogen data and find a temperature of ~3310 K at the output (not 3000 K), suggesting an error of ~13%-15% in the calculated hydrogen enthalpy of the paper. Our solver agreed with that hand calculation to within 3 K. We use the reference data and compare trends..

How cold frit can be enough? In the simplified model of Kerrebrock (without considering plenum and channel), we found that a cold frit with 20-30% of the element pressure drop eliminates the instability. Including plenum and exit channel, 75.5% of the element pressure drop would be needed. The design in the paper uses 90%.

For stability, power density doesn’t decide which works. Cold-frit-controlled reactors were found to be stable between 1 MW/L and 30 MW/L while hot-frit-controlled reactors were not, even for low power density. This is contrary to some aspects of the published theory explaining when each technology should apply.

An optimized frit sacrifices margin. Ji et al.'s optimized cold frit maintained stability from 5% to 100% power, yet only had one-half to two-thirds the published frit's stability margin. It also operates at a higher temperature compared to the published frit, particularly at low power, exceeding the 3800 K fuel limit at 20-22.5% power.

Flow shift costs 600 K, which can be compensated by a U-shaped flow pattern. Flow shift can be measured as the coefficient of variation of the slice flows:

$$ \text{CoV} = \frac{\sigma(\dot m_i)}{\overline{\dot m}} $$

For the conventional Z configuration (both plenum and channel flows going the same direction), CoV = 0.150; the most fed slice receives 1.6 times as much flow as the least fed one. This increases peak fuel temperature already from 3064 K to 3671 K, very close to 3800 K limit. If we use the opposite end of the plenum to feed the flow (U configuration), then CoV drops to 0.030 and the peak fuel temperature returns to 3131 K.

The Tsinghua large reactor (Tier 2). The four checks here are not evaluable: each needs a case that goes beyond the 5000 K limit for our property data (peaks reached in paper itself go up to 7121 K). Before running the check, we decided that such a check will never be considered as passing. The cases that did solve are listed separately, as observations, not counted. At core heights 80 and 100 cm, our pressure drops agree within 0.3% and 1.3% to the published values, without any tuning of inputs.

The trained surrogate and why it does not satisfy our criteria. Ji et al. mention an efficient surrogate as future work, which we provided in the form of gradient-boosted surrogates based on 2000 solver runs. It is 180,000 times faster than the solver, and its mean peak-temperature error of 41.5 K satisfies our criterion. However, we required that the surrogate have a balanced accuracy of at least 0.95, where:

$$ \text{balanced accuracy} = \frac{1}{2}\left(\frac{TP}{TP + FN} + \frac{TN}{TN + FP}\right) $$

it received a score of 0.915. It fails to detect 8 out of 48 designs which do not have any steady states, erroneously suggesting there is one. Accuracy is 0.978 and seems high only because 88% of designs converge. The strict criterion is what it fails, hence it cannot substitute the solver by our criterion.

Runaway demonstration. Take one slice 20 K warmer and observe. In the presence of the cold frit, the perturbation decays exponentially to zero, with characteristic time of 0.39 s. In its absence, the hot slice stops receiving coolant flow and its outlet gas exceeds the 3800 K fuel limit at 7.6 s.

Issues we encountered

Data integration problems. The coolprop database, Svehla and NIST-JANAF differ somewhat when overlapping. Combining them without making kinks which could trigger the solver required significant effort.
Convergence on the critical point. Solver finds it difficult to converge precisely around the instability line – the point of interest to us – and this is what the continuation procedure is needed for.
Disagreeing published papers relative to reference data. We decided to report the disagreement rather than adjust the model to match.
Failing one validation test. Given that the criteria have been pre-established, we present the failure and provide an explanation.
Results that cannot be evaluated. The large reactor operates at a temperature higher than that which can be obtained from the property data. Pre-defining what "not evaluable" means prevented us from quietly choosing a different problem.
A well-known weakness in the model. The plenum pressure recovery term employs $\rho u^2$, while our assumption calls for $\tfrac{1}{2} \rho u^2$. The fix had been entered, but was subsequently overwritten in a re-write of the solver; hence all results presented herein employ $\rho u^2$. Results using Kerrebrock's configuration do not employ the plenum.

What we learned

Having pass/fail criteria written out beforehand in light of no results altered the process we followed: all became a finding, from the test that failed to the surrogate that did not make the cut. We learned to report inconsistencies between sources rather than force them to agree, and that simple accuracy can hide exactly the cases that matter. We also had an experience related to version control of our own: a change to the physics code that was committed but quietly reverted later on, which we discovered via the history.

What's next

Reintroduce the momentum correction for plenum as a new model version and repeat all results using it. Afterwards, include neglected factors such as hydrogen dissociation and cross flow between slices and re-train the surrogate so that “OK” also means that fuel must be below 3800 K, based on a new pre-registered set.