Full-Core Coupled Neutronic, Thermal-Hydraulic, and Thermo-Mechanical Analysis of Low-Enriched Uranium Nuclear Thermal Propulsion Reactors

Abstract

Nuclear thermal propulsion is an enabling technology for future space missions, such as crew-operated Mars missions. Nuclear thermal propulsion technology provides a performance benefit over chemical propulsion systems by operating with light propellants (e.g., hydrogen) at elevated engine chamber conditions. Therefore, nuclear thermal propulsion reactor cores exhibit high propellant velocities and elevated propellant and fuel temperatures, subsequently leading to relatively high thermal stresses and geometrical deformation. This paper details the numerical approach to solve the thermo-elastic equations, which was implemented into the recently developed ntpThermo code. In addition, this paper demonstrates the extension of the Basilisk multiphysics framework to perform full-core coupled neutronic, thermal-hydraulic, and thermo-mechanical analysis of nuclear thermal propulsion reactors. The analyses demonstrate and quantify thermo-mechanical feedback, which for the investigated cases, acted to reduce maximum fuel temperatures and pressure drop across the fuel element channels. Thermo-mechanical feedback had a significant impact on the mass flow distribution within the reactor core and, thus, a substantial impact on solid-material temperatures and stresses, but only a minor impact on reactivity and local power distributions. Sensitivity studies revealed that the friction factor correlation applied to perform the analysis has a significant impact on the pressure drop across the fuel element channels. The most important observation of this research is the importance of incorporating the thermo-mechanical feedback within an integrated multiphysics solution sequence to enable the consistent design of future nuclear thermal propulsion systems.

1. Introduction

Nuclear thermal propulsion (NTP) systems are a candidate technology for crew-operated Mars missions by the National Aerospace and Space Administration (NASA) and Cis-Lunar operations by the Defense Advanced Research Projects Agency (DARPA). NTP engines use the energy generated by fission to heat the working fluid, typically hydrogen, to achieve extremely high temperatures (~2700 K) and considerable chamber pressures (~7 MPa). This allows NTP engines to produce twice the specific impulse (I_sp) of an ideal chemical rocket engine while also providing thrust levels on the order of 10,000 pounds (lb_f). NASA is currently advancing state of the art low-enriched uranium (LEU) NTP designs with the objective of using NTP engines for human Mars missions in 2039 [1,2,3]. Recently, DARPA announced that the Demonstration Rocket for Agile Cislunar Operations (DRACO) program is targeting a 2026 integrated engine in-space flight demonstration mission to mature NTP technology for near-term applications [4]. Modern LEU designs for Mars transit target thrust values in the range of 12,500–25,000 lb_f with an I_sp of 900 s [1]. In order to achieve an I_sp of 900 s, the average exit bulk coolant temperature must be at least 2700 degrees kelvin, which implies even greater fuel temperatures (~2900 K). An additional important performance requirement is a high thrust-to-weight ratio (≥1.8). Satisfying both the I_sp and thrust-to-weight ratio requirements leads to an extremely high power-density (≥1 MW/L) reactor design, thus, presenting a significant design challenge. Multiphysics modeling and simulation (M&S) tools capable of capturing the relevant feedback mechanisms are required to drive the development of a reactor design that satisfies engine performance requirements while maintaining adequate thermal and mechanical safety margins.

The need to develop advanced M&S capabilities for the design of NTP reactors has been highlighted by a recent National Academy of Science, Engineering, and Medicine (NASEM) study [5], which recommended that “NASA should rely on extensive investments in modeling and simulation”, in addition to ground testing and cargo missions to qualify NTP systems for crewed missions. Recent research [3,6] has demonstrated that a multiphysics solution approach is required for the design and analysis of NTP systems; however, the significance of various feedback mechanisms within the multiphysics solution approach is not well quantified. Typically, neutronics is coupled with thermal-hydraulics (T/H) to account for the impact of variations in solid-material temperatures and coolant density on reactivity and power distributions. Currently, sparse literature exists on the importance of including other feedback mechanisms, such as thermo-mechanics (T/M), within a multiphysics solution scheme to analyze LEU NTP reactors. Previous research has leveraged higher-order codes to perform T/H and T/M analysis of a single supercell [7,8]. However, the previously mentioned analyses have not considered the impact of T/M feedback on the T/H solution. Several papers have recommended the inclusion of T/M feedback within a multiphysics solution approach [3,9]. This recommendation was substantiated by recent research, which demonstrated that T/M feedback impacts the mass-flow distribution and, thus, solid-fuel temperatures of multichannel simulations of an NTP fuel element (FE) [10].

This paper details the T/M solution approach implemented into the ntpThermo code [11] and the implementation of a coupled neutronic, T/H, T/M (N+T/H-T/M) solution sequence into the Basilisk multiphysics framework [6]. The primary objective of this paper is to demonstrate the importance of T/M feedback in coupled N+T/H-T/M simulations. The results presented here demonstrate that the coupled T/M feedback reduces maximum fuel temperatures, pressure drops across the FE channels, and increases maximum von Mises (VM) stresses within the FE channels. The thermal expansion of coolant channels significantly impacts the pressure losses and heat transfer characteristics within the FE channels and, thus, the mass flow and fuel temperature distributions within the core. Neglecting the T/M feedback can lead to inconsistencies in the fuel temperature and VM stress spatial distributions, which are important for systematic design studies. The sensitivity studies performed in this paper reveal that the pressure drop across the FE channels is dependent on the friction factor correlation applied to solve the problem. To conclude, this paper demonstrates the importance of T/M feedback within a multiphysics solution approach in predicting consistent safety margins.

2. Codes and Methods

Section 2.1 contains a brief overview of the ntpThermo code. Section 2.2 contains a description of the thermoelastic equations solution approach. Finally, Section 2.3 describes the iterative T/H-T/M solution scheme and the N+T/H-T/M solution sequence.

2.1. Overview of the ntpThermo Code

The ntpThermo code was developed to fill the gap between expensive high-order computational fluid dynamics (CFD) codes and unvalidated reduced-order codes for the analysis of NTP reactors. Higher-order codes are non-ideal for full-core multiphysics analysis due to their computational cost. While ntpThermo is not intended to replace higher-order codes, it is an ideal code for a coupled full-core N+T/H analysis because of reduced computational cost and adequate accuracy [11]. The code relies on an equivalent model with a 1.5 D conduction–convection solution. A numerical solution scheme was implemented to solve each axial layer with the convective solution relying on Newton’s law of cooling. An iterative solution scheme was developed to solve the counterflow within the moderator element (ME) and capture the heat conducted from the fuel to the moderator for various supercell configurations. This numerical solution scheme was verified against a series of higher-order OpenFOAM CFD benchmarks [11]. The code contains a mass flow convergence scheme that balances the mass flow rate to achieve an equal pressure drop in each channel, where the user provides the inlet pressure of each channel [6]. In this work, each channel is assumed to have the same inlet pressure obtained from the engine cycle analysis summarized in Section 3.1. This solution scheme was applied to converge the mass flow distribution within the FEs and MEs for the analysis detailed in Section 4. Additionally, a comprehensive material properties interface was developed to accommodate temperature and pressure-dependent coolant thermophysical properties as well as user-defined two variable-dependent solid-material properties.

2.2. Thermoelastic Equations Solution Methodology

The T/M solver implemented in ntpThermo leverages the approach developed by Peng and Li [12], where the thermoelastic problem is cast as a second-order Fredholm integral equation. The primary advantage of this approach is the ability to capture temperature-dependent properties in the solution of the thermoelastic equations. This capability is necessary to model an NTP system as the power densities inherent in NTP systems induce significant temperature gradients within the fuel. Additionally, the thermophysical properties of the fuel and cladding materials are a strong function of temperature.

The solution method assumes an axisymmetric plane stress problem where the axial stress (${\sigma }_z$) and strain (${\epsilon }_z$) are assumed to be zero, which is consistent with previous industry studies [1,13]. On the basis of this assumption, the equations for the radial (${\sigma }_r$) and tangential (${\sigma }_{\theta }$) stresses can be detailed by Equations (1) and (2), respectively. The radial (${\epsilon }_r$) and hoop strain (${\epsilon }_{\theta }$) are detailed by Equations (3) and (4), respectively.

$${\sigma }_r\left(r\right)=\frac{{\rm E}\left(r\right)}{1-v{\left(r\right)}^2}\left({\epsilon }_r+v\left(r\right){\epsilon }_{\theta }-\alpha \left(r\right)\left[1+v\left(r\right)\right]T\left(r\right)\right)$$

(1)

$${\sigma }_{\theta }\left(r\right)=\frac{{\rm E}\left(r\right)}{1-v{\left(r\right)}^2}\left({\epsilon }_{\theta }+v\left(r\right){\epsilon }_r-\alpha \left(r\right)\left[1+v\left(r\right)\right]T\left(r\right)\right)$$

(2)

$${\epsilon }_r=\frac{d{u}_r}{dr}$$

(3)

$${\epsilon }_{\theta }=\frac{u_r}{r}$$

(4)

where E is Young’s modulus, $v$ is Poisson’s ratio, and $\alpha$ is the coefficient of thermal expansion, which are all spatially dependent. It is important to note that $T\left(r\right)$ is the difference between the absolute temperature ($T_{abs}$) obtained by the thermal solver and the reference temperature zero stress temperature ($T_{ref}^{\alpha }$); this relationship is detailed in Equation (5).

$$T\left(r\right)=T_{abs}\left(r\right)-T_{ref}^{\alpha }$$

(5)

On the basis of the equations detailed above, the radial displacement ($u_r$) and tangential stress (${\sigma }_{\theta }$) can be represented by Equations (6) and (7), respectively. The equivalent von Mises stress (${\sigma }_{VM}$) is calculated via Equation (8).

$$u_r\left(r\right)=\frac{A}{\chi \left(r\right)}+\frac{1}{\chi \left(r\right)}{{\displaystyle \int }}_a^r\chi \left(r\right)\left\{\frac{1-v{\left(s\right)}^2}{{\rm E}\left(s\right)}{\sigma }_r\left(s\right)+\alpha \left(s\right)\left[1+v\left(s\right)\right]T\left(s\right)\right\}ds$$

(6)

$${\sigma }_{\theta }\left(r\right)=v\left(r\right){\sigma }_r\left(r\right)+E\left(r\right)\frac{u_r\left(r\right)}{r}-E\left(r\right)\alpha \left(r\right)T\left(r\right)$$

(7)

$${\sigma }_{VM}=\sqrt{{\sigma }_r^2-{\sigma }_r{\sigma }_{\theta }+{\sigma }_{\theta }^2}$$

(8)

where $A$ in Equation (6) is a constant to be determined through appropriate boundary conditions, and $\chi \left(r\right)$ can be calculated via Equation (9).

$$\chi \left(r\right)=\exp \left[{{\displaystyle \int }}_a^r\frac{v\left(s\right)}{s}ds\right]$$

(9)

When Equations (3) and (4) are substituted into the equilibrium equation, Equation (10) is produced.

$$\frac{d{\sigma }_r}{dr}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=0$$

(10)

Substituting Equation (6) into Equation (7) and then into Equation (10) yields Equation (11)

$$\frac{d{\sigma }_r}{dr}+\frac{1-v\left(r\right)}{r}{\sigma }_r-\frac{E\left(r\right)}{r^2\chi \left(r\right)}{{\displaystyle \int }}_a^r\frac{\left(1-v^2\left(\rho \right)\right)\chi \left(\rho \right)}{{\rm E}\left(\rho \right)}{\sigma }_r\left(\rho \right)d\rho =F\left(r\right)+A\frac{E\left(r\right)}{r^2\chi \left(r\right)}$$

(11)

where $a$ is the inner radius of the equivalent annulus, and $F\left(r\right)$ is:

$$F\left(r\right)=\frac{E\left(r\right)}{r^2\chi \left(r\right)}{{\displaystyle \int }}_a^r\chi \left(\rho \right)\alpha \left(\rho \right)\left[1+v\left(\rho \right)\right]T\left(\rho \right)d\rho -\frac{E\left(r\right)\alpha \left(r\right)T\left(r\right)}{r}$$

(12)

Integrating Equation (11) with respect to r yields Equation (13)

$$\begin{gathered} {\sigma }_r\left(r\right)+{{\displaystyle \int }}_a^{r}K\left(r,\rho \right){\sigma }_r\left(\rho \right)d\rho =F_1\left(r\right)+Af\left(r\right)+B \\ K\left(r,\rho \right)=\frac{1-v\left(\rho \right)}{v\left(\rho \right)}-\frac{\left[1-v^2\left(\rho \right)\right]\chi \left(\rho \right)}{E\left(\rho \right)} {{\displaystyle \int }}_{\rho }^r\frac{E\left(s\right)}{s^2\chi \left(s\right)}ds \\ {F}_1\left(r\right)={{\displaystyle \int }}_a^{r}F\left(\rho \right)d\rho \\ f\left(r\right)={{\displaystyle \int }}_a^r\frac{E\left(s\right)}{s^2\chi \left(s\right)}ds \end{gathered}$$

(13)

The constants $A$ and $B$ can be found via Equations (14) and (15), respectively, where ${\sigma }_r\left(a\right)$ is the radial stress at the inner radius of the equivalent annulus and ${\sigma }_r\left(b\right)$ is the radial stress at the outer radius of the equivalent annulus.

$$A=\frac{{\sigma }_r\left(b\right)-{\sigma }_r\left(a\right)}{f\left(b\right)}-\frac{F_1\left(b\right)}{f\left(b\right)}+\frac{1}{f\left(b\right)}{{\displaystyle \int }}_a^{b}K\left(b,\rho \right){\sigma }_r\left(\rho \right)d\rho$$

(14)

$$B=-{\sigma }_r\left(a\right)$$

(15)

Now we insert the constants $A$ and $B$ back into the integrodifferential equation for the radial stress to obtain Equation (16) which has the form of a second-order Fredholm integral equation.

$${\sigma }_r={{\displaystyle \int }}_a^{b}L\left(r,\rho \right){\sigma }_r\left(\rho \right)d\rho =h\left(r\right)$$

(16)

The kernel function $L\left(r,\rho \right)$ is detailed in Equation (17) and $h\left(r\right)$ is detailed in Equation (18).

$$L\left(r,\rho \right)=\left\{\begin{array}{cc}K\left(r,\rho \right)-\frac{f\left(r\right)}{f\left(b\right)}K\left(b,\rho \right)& \rho <r\\ -\frac{f\left(r\right)}{f\left(b\right)}K\left(b,\rho \right)& \rho >r\end{array}\right.$$

(17)

$$h\left(r\right)=F_1\left(r\right)+\frac{f\left(r\right)}{f\left(b\right)}\left[-{\sigma }_r\left(a\right)+{\sigma }_r\left(b\right)-F_1\left(b\right)\right]+{\sigma }_r\left(a\right)$$

(18)

The T/M solution is obtained using the finite element data provided by the T/H heat transfer solution. More specifically, the finite element solid-material temperatures are used to obtain the temperature-dependent mechanical properties (e.g., $\alpha$). These finite element temperatures, properties, and radial dimensions are then used to solve the different integrals appearing in Equations (6)–(18) via numerical integration. The latter enables the calculation of stresses and thermal expansion in each finite element radial region. This can be used to update the finite element mesh within an iterative T/H-T/M solution sequence. The described solution relies only on a 1 dim radial approach, also referred to as the equivalent annulus approximation. A detailed description of radial and axial discretization of the mesh can be found in previous research [11,14]. The numerical solution technique implemented in ntpThermo leverages the previously developed and verified approach by Atkinson and Shampine [15]. The algorithm provides a numerical solution for integral equations in the form:

$$\lambda x\left(s\right)-\begin{array}{cc}{{\displaystyle \int }}_a^{b}K\left(s,t\right)x\left(t\right)dt=f\left(s\right)& a\le s\le b\end{array}$$

(19)

which can be written symbolically as:

$$\left(\lambda -K\right)x=f$$

(20)

The routine implemented solves integral equations with kernel functions $K\left(s,t\right)$ that are relatively smooth on $R=\left[a,b\right]\times \left[a,b\right]$ for the interval $\left[a,b\right]$. The procedure relies on Simpson’s rule to discretize the equation to get an approximate solution on all of $\left[a,b\right]$. The solution begins by approximating the integral operator using a quadrature scheme:

$${{\displaystyle \int }}_a^{b}g\left(t\right)dt\cong {{\displaystyle \sum }}_{j=0}^n{w}_{j,n}g\left(t_{j,n}\right)$$

(21)

For smooth kernels, the composite Simpson’s rule for irregularly uneven spaced data with an even number of subintervals, each with a subinterval thickness of $h,$ is detailed in Equation (22).

$${{\displaystyle \int }}_a^{b}f\left(t\right)dt=\frac{h}{3}\left[f\left(x_0\right)+2{{\displaystyle \sum }}_{j=1}^{N/2-1}f\left(x_{2j}\right)+4{{\displaystyle \sum }}_{j=1}^{N/2}f\left(x_{2j-1}\right)+f\left(x_n\right)\right]$$

(22)

It must be noted that the current implementation is limited to only cases where the radial mesh is divided into an even number of regions. However, the individual thickness of each subregion may be different. Simpson’s rule is then used to generate the weight coefficients ($w_{j,n}$) such that $\left(\lambda -K\right)x=f$ can be approximated by Equation (23). Combing the relations detailed above leads to a linear system of equations that can be solved via Equation (24), which is ultimately used to obtain ${\sigma }_r\left(r\right)$. The numerical implementation within ntpThermo has been previously verified against analytical benchmarks [16].

$$\lambda x_n\left(s\right)-{{\displaystyle \sum }}_{j=0}^n{w}_{j,n}K\left(s,t_{j,n}\right)x\left(t_{j,n}\right)=f\left(s\right)$$

(23)

$$x_n\left(s\right)=\frac{1}{\lambda }f\left(s\right)+\left[{{\displaystyle \sum }}_{j=0}^n{w}_{j,n}K\left(s,t_{j,n}\right)x\left(t_{j,n}\right)\right]$$

(24)

2.3. Coupled Neutronic, Thermal-Hydraulic, Thermo-Mechanical Solution Sequence

The T/M solution utilizes equivalent annulus approximation and, therefore, the same radial and axial mesh as the T/H solution. The temperature distributions obtained from the T/H solution can be directly used by the T/M solver in conjunction with temperature-dependent mechanical properties. The coupled T/H-T/M solution sequence, detailed in Figure 1, begins by iteratively solving the multichannel T/H problem to obtain a converged mass flow rate distribution $\left(\overrightarrow{\dot{m}}\right)$. The multichannel T/H solution temperature distribution, dimensions, and temperature-dependent mechanical properties are then input to the T/M solver. The T/M solver returns the radial stress and thermal expansion distributions for each axial layer of each coolant channel specified in the problem. The T/M solution is used to update the radial mesh of each axial layer of each channel, including the node volumes, power densities, and fluid flow areas. The updated geometry is then used by the next T/H solution iteration, where the iteration number is represented by n, to reconverge the mass flow distribution. The outlet pressure is chosen as the convergence criteria for the coupled T/M-T/H sequence, as the T/M feedback primarily affects the pressure drop of each channel.

The coupling framework within the Basilisk code [6] was extended to account for T/M feedback applied to full-core N+T/H-T/M simulations; the solution sequence is detailed in Figure 2. The Serpent code [17] is used as the neutronic solver within the Basilisk code. The methodology for coupled N+T/H simulations has been detailed in previous research [6]. The T/M feedback is accounted for in the coupled N+T/H-T/M simulation via the impact of thermal expansion on the T/H solution; however, the dimensions in the Serpent model are not updated. The coupled sequence utilizes the Serpent multiphysics interface [18] to update the solid-material temperatures and coolant gas densities for each FE and ME in the reactor core throughout each iteration. The coupled analysis sequence converges on k_eff and the fuel power distribution [6]. The excess reactivity $\left({\epsilon }_{keff}\right)$ and fuel power distribution convergence criteria were set at 10 pcm and 1E-5, respectively, for the analysis in this paper.

The energy deposition in the active core is tallied as a function of axial height on an element-by-element basis where each FE and ME is individually tracked. Additionally, each material region within each element is uniquely tallied. However, the energy deposition within that region is assumed to be uniform (i.e., intra-element power peaking is not accounted for); however, future research may expand the T/H-T/M solution methodology to account for the intra-element power distribution within the FEs. The T/H fields are mapped to the neutronic model via the approach presented in Figure 3, where each color represents a unique material region within the ME and FE. The ntpThermo code relies on the equivalent annulus approximation and numerical conduction–convection solver [11]. As such, the volume average of the radial temperature profile obtained from the ntpThermo model of material region $i$ for axial layer $j$ of element $k$ (${\overline{T}}_i^{k,j}$) is applied to the respective region in the neutronic model. Additionally, the bulk coolant densities in each coolant channel for every axial layer are updated.

As previously stated, the Serpent model mesh was not modified during the N+T/H-T/M solution sequence as its impact is relatively small. A two-dimensional x-y Serpent model was used to investigate the impact of explicitly modeling the thermal expansion of FE coolant channels. The T/H fields at the core mid-plane (z = 60 cm) obtained from the full-core N+T/H-T/M analysis in Section 4.1 were used in the analysis. This specific axial height was chosen because of the elevated fuel temperatures. Two models were constructed; the first explicitly modeled the thermal expansion of each FE within the core. The total fuel mass was preserved in the expanded model by modifying each fuel material’s density. The second model did not account for any thermal expansion (i.e., all the coolant channels are identical). Only a 23 pcm difference in reactivity was observed between the two models, where the statistical uncertainty was 13 pcm. The observed RMS error between the two radial power distributions was 0.04%, where the average statistical uncertainty in each power detector was 0.11%. This analysis demonstrates that the explicit modeling of coolant channel thermal expansion has a negligible impact on reactivity and local power distributions.

3. Reactor Design Description

Section 3.1 contains a description of the reactor design considered. Section 3.2 contains a brief description of the temperature-dependent mechanical material properties used in the T/M solution of the FEs.

3.1. Reactor Design Description

The reactor design considered in this paper is a Nuclear Engine for Rocket Vehicle Applications (NERVA) derived design, where the active core region is composed of hexagonal internally cooled FEs and dual-pass tie-tube MEs. The geometry of the fuel and MEs is identical to those previously studied, and several key parameters are detailed in

Table 1. Key design parameters.

Parameter	Value
Number of Fuel Elements	127
Number of Moderator Elements	390
Number of Control Drums	18
Poison Vane Thickness	0.5 cm
Poison Vane B4C Enrichment	95 wt%
Active Core Diameter	31.0 cm
Active Core Height	120.0 cm
Radial Reflector Thickness	11.0 cm
Axial Reflector Thickness	10.0 cm
Active Core H:235U Ratio	334:1
Total HALEU Mass	40.3 kg
235U Enrichment	19.75 wt%
Serpent Model Dry Mass	2508 kg

[6]. The radial and axial core cross-sections are detailed in Figure 4. A detailed design description of the reactor configuration and engine expander cycle can be found in previous research [14,19]. The design uses a Ceramic–Ceramic (CERCER) fuel form, where Zirconium Carbide (ZrC) coated High Assay Low-Enriched Uranium (HALEU) Uranium Nitride (UN) fuel kernels are embedded in a ZrC matrix [1]. This fuel form is currently favored over previously investigated Ceramic–Metallic (CERMET) fuel forms because of its reduced weight and parasitic neutron absorption. The moderator material within the ME is Zirconium Hydride (ZrH₁.₈₉), which is consistent with previous industry studies [1]. An inverse pewee lattice configuration was chosen to ensure that the active core has a similar Hydrogen-to-Uranium-235 atom ratio (H:²³⁵U) as the baseline industry design [1]. The previous calculated steady-state boundary conditions [6] obtained from the NTP-Transient code [19] are applied in this paper.

The previously considered active core lattice pattern [6] was slightly modified to eliminate 26 kg of mass from the reactor subsystem and reduce the required HALEU mass by 11.1 kg while preserving a hot full-power critical drum position of 60 degrees. The Basilisk code contains a fuel loading search algorithm that adjusts the fuel volume loading of each FE to flatten the radial power distribution for a given reactor design [20]. The final fuel loading pattern detailed in Figure 5 contains only 16 unique fuel loadings, and the maximum fuel particle volume loading is 38.9 vol% particles, which is well below the recommended maximum limit of 65 vol% particles [1]. Flattening the radial power profile is necessary to minimize peak fuel temperatures and maximize the overall performance of the engine [20].

3.2. Temperature-Dependent Mechanical Material Properties

A literature review was conducted to obtain the temperature-dependent mechanical properties required to solve the thermoelastic equations. Temperature-dependent correlations for each material are detailed in

Table 2. Temperature-dependent fuel material mechanical properties.

Material	Reference	Property
ZrC	[21]	$E\left(T\right) \left[\mathrm{Pa}\right]=\begin{array}{cc}4.16\times {10}^{11}-1.3\times {10}^{7 }T-12958 T& 295K\le T\le 2253K\end{array}$
	[22]	$\alpha \left(T\right) \left[{\mathrm{K}}^{-1}\right]=\begin{array}{cc}1.7792\times {10}^{-9} T+3.5258\times {10}^{-6}& 1200K\le T\le 2800K\end{array}$
	[23]	$v=0.213$
UN	[24]	$E\left(T\right) \left[\mathrm{Pa}\right]=\begin{array}{cc}25.8\times {10}^6\left(1-2.375\times {10}^{-5} T\right)& 298K\le T\le 1473K\end{array}$
	[25]	$\alpha \left(T\right) \left[{\mathrm{K}}^{-1}\right]=\begin{array}{cc}7.096\times {10}^{-6}+1.409\times {10}^{-9 }T& 298K\le T\le 2523K\end{array}$
		$v=0.281$

. The thermophysical properties detailed in previous research were used in the T/H solution [6]. The zero-stress reference temperature was assumed to be 1500 K to remain consistent with previous studies [7]. The effective CERCER mechanical properties were obtained using the method of mixtures approach [7].

4. Coupled Multiphysics Results

This section presents the coupled N+T/H–T/M analysis. Section 4.1 demonstrates the importance of T/M feedback on full-core N+T/H–T/M analysis. Section 4.2 contains a sensitivity study that explores the impact of various friction factor correlations on the coupled N+T/H–T/M solution.

4.1. Importance of Thermo-Mechanical Feedback in Modeling NTP Systems

The objective of this section is to demonstrate the importance of T/M feedback by comparing N+T/H analysis to N+T/H-T/M analysis. In this work, the T/M feedback solution is applied only to the FEs as they experience the most considerable deformation and stress. The T/M solution of the MEs will be the focus of future studies, including the flow distribution balance analysis when the supply and return channels are expanded or contracted. The analysis presented in this section was obtained with the Churchill friction factor correlation applied to solve the FEs on the basis of the recommendations of previous research [26]. No orificing was applied to the FEs to facilitate comparisons between N+T/H and N+T/H-T/M analysis. The primary feedback that T/M has on the T/H solution is decreasing the pressure loss of higher-power channels, thus decreasing the redistribution of mass flow within the reactor core. Higher-power channels experience greater thermal expansion (increase in flow area) and thus reduced pressure losses.

To quantify the convergence of the multiphysics solution, the L-2 norm error defined by Equation (25) was calculated for multiple spatially distributed T/H fields, such as fuel temperature ($T_{fuel}$), ZrH temperature ($T_{ZrH}$), fuel gas density (${\rho }_{fuel}$), moderator supply gas density (${\rho }_{supply}$), moderator return gas density (${\rho }_{return}$), and fuel power ($Q_{fuel}$). These T/H fields are presented in Figure 6 because they have a significant impact on the spatial power distribution within the reactor [6].

$$l_2^n=\frac{\Vert {\overrightarrow{X}}_n-{\overrightarrow{X}}_{n-1}\Vert {}_2}{\Vert {\overrightarrow{X}}_{n-1}\Vert {}_2}$$

(25)

In Equation (25), $\overrightarrow{X}$ is the spatial distribution of interest for consecutive iterations $n-1$ and $n$ [27]. The trends in Figure 6 demonstrate that temperature, density, and power spatial distributions converged after the sixth iteration. Figure 7 details the excess reactivity for each iteration; the difference in reactivity between the fifth and sixth iteration is within the statistical uncertainty (±7 pcm) of the Monte Carlo solution.

The radial power peaking of each FE in the core with and without T/M feedback is detailed in Figure 8. When comparing Figure 8a,b, it becomes apparent that T/M feedback slightly increases (+0.9%) the maximum radial power peaking. The percent difference between the radial power distributions is detailed in Figure 8c, where the percent difference ($\%\Delta$) is calculated via $\%\Delta =100\times \left(X_{w/TM}-X_{w/oTM}\right)/X_{w/oTM}$, where $X$ is the quantity of interest. A positive %∆ represents where N+T/H analysis underpredicts a given value when compared with N+T/H-T/M analysis. The maximum fuel temperature of each FE is detailed in Figure 9. T/M feedback decreases the maximum fuel temperature by 129 K. Despite the minor increase in radial power peaking, the hot channel maximum fuel temperature decreases as a result of the significant increase in mass flow to the hot channel when T/M feedback is considered. The reduction in solid-fuel temperatures reduces the amount of heat conducted from the central higher-power FE channels to their surrounding ME channels. This, in turn, reduces the local ZrH temperature of the surrounding MEs. Previous research [6] has demonstrated that lower ZrH temperatures increase the local fission rate. Since the local ZrH in the center of the core is reduced when T/M feedback is considered, the radial power peaking increases. The local ZrH temperature has a significant effect on neutron moderation as the thermal peak is shifted with respect to neutron energy [6].

The mass flow distribution of the FEs with and without T/M feedback is compared in Figure 10. It becomes clear when comparing Figure 10a,b that T/M feedback significantly impacts the mass flow distribution within the FE channels as the hot channel experiences a 9.6% increase in mass flow when T/M feedback is considered.

The difference in mass flow distribution is caused by the thermal expansion of coolant channels. The thermal expansion of coolant channels increases the flow area within the FE, resulting in reduced pressure drops, thus, requiring additional mass flow to equalize the pressure drop across the FE channels. The percent difference in the flow area of each FE channel at the core mid-plane (z = 60 cm) with and without thermal expansion is detailed in Figure 11. The mass flow and flow area percent difference distributions illustrate this effect as they have nearly identical shapes. As expected, the hottest channel experiences the greatest increase (+22%) in flow area because of its elevated solid-fuel temperatures.

The outlet bulk coolant temperature distributions with and without T/M feedback are detailed in Figure 12.

Figure 13 compares each component of the friction pressure losses ($\Delta P_{fric}$) for the hot channel with and without T/M feedback. The current implementation does not consider the impact of sudden contraction/expansion of coolant channels on pressure losses within the FE channels. However, future research will investigate the importance of such effects. The with/without ratio of each parameter is presented as a function of axial height in Figure 13. The $\Delta P_{fric}$ and acceleration pressure losses ($\Delta P_{accel}$) equations are detailed in Equations (26) and (27), respectively, where $f$ is the friction factor, $G$ is the mass velocity, $L$ is the length of the axial layer, $D_h$ is the hydraulic diameter of the channel, and $\rho$ is the density of the fluid.

$$\Delta P_{fric}=2f\frac{L}{D_h}\frac{G^2}{\overline{\rho }}$$

(26)

$$\Delta P_{accel}=G^2\left(\frac{1}{{\rho }_{in}}-\frac{1}{{\rho }_{out}}\right)$$

(27)

The mass velocity significantly decreases as the flow area along the channel increases, detailed by the red line in Figure 13. The hot channel pressure losses and bulk coolant pressure obtained via N+T/H and N+T/H–T/M analysis are compared in Figure 14. The decreased mass velocity affects both the $\Delta P_{fric}$ and $\Delta P_{accel}$ as they are both a function of $G^2$. The total pressure drop across the FEs is overpredicted by 150 kPa when T/M feedback is neglected. The results in Figure 13 and Figure 14 illustrate the mechanism causing the mass flow distribution differences in Figure 10c.

T/M feedback has an overall beneficial impact on the T/H performance for the design considered in this paper, as the thermal expansion of coolant channels reduces the pressure drop across the FE channels and redistributes mass flow to higher-power channels. When comparing Figure 9c and Figure 10c, it becomes clear that the heat transfer within the hot channel is reduced as the mass flow increases by 9.6%, but the maximum fuel temperature is only reduced by 4.0%. The hot channel heat transfer coefficient, as a function of axial height, is detailed in Figure 15a and demonstrates that T/M feedback has a negative impact on the heat transfer coefficient. Several competing effects are responsible for this behavior. Similar to Figure 13, the with/without ratios for several parameters which influence the heat transfer coefficient and are significantly impacted by T/M feedback are detailed in Figure 15b. Thermal expansion reduces the coolant velocity ($u$) within the channel, which is partially offset by the increase in the hydraulic diameter ($D_h$) and density ($\rho$). These competing effects result in a minor reduction in Reynolds number near the exit (z > 80 cm) of the hot channel. Hydrogen’s viscosity ($\mu$) and thermal conductivity ($t_c$) increase with temperature [28]. The bulk coolant temperature within the hot channel decreases because of the mass flow rate distribution, thus resulting in a decrease in both the $\mu$ and $t_c$ along the channel. All these factors combined result in the reduction in the heat transfer coefficient within the hot channel.

The maximum VM stress in each FE channel with and without T/M feedback is detailed in Figure 16. When T/M feedback is neglected, the maximum VM stress occurs in the hot channel (Location A). However, when T/M feedback is considered, the maximum VM stress occurs in the FE channel with the highest fuel particle packing fraction (Location B).

The VM stress is commonly used in yielding/fracturing criteria [29] and thermal creep rate correlations [30]. Therefore, the difference in the maximum VM stress and local VM stress distributions should be considered. The FE radial temperature and VM stress profiles for the FE channel at the maximum stress axial layer (z = 82 cm) for location B are detailed in Figure 16. Dashed lines in Figure 17 are used to depict the outer radii of the cladding and fuel meat with T/M feedback (indicated by red lines) and without T/M feedback (indicated by black lines). When T/M feedback is accounted for, the radial mesh is deformed, resulting in the dashed lines moving to the right. Location B’s max VM stress experiences an 11% increase when T/M feedback is considered, which is primarily caused by the local fuel temperature increasing by ~200 K, as detailed in Figure 17a. Figure 17b shows that the greatest VM stress occurs on the exterior of the coolant channel cladding, which is consistent with previous studies [7]. The coefficient of thermal expansion (CTE) mismatch between the fuel matrix and cladding materials has been previously identified as a primary source of thermal stress [7] and has been linked to increased fuel mass loss rates [31]. The FE in location B has the largest CTE mismatch because of the channel-elevated fuel loading (see Figure 4). UN’s thermal expansion coefficient is nearly double that of ZrC and is the primary cause of location B’s elevated VM stresses.

A comparison of several key global parameters is detailed in

Table 3. Comparison of global parameters.

Parameter	w/o T/M	w/ T/M	Difference
Excess reactivity, pcm	23 ± 7	38 ± 7	+15
Max. fuel temperature, K	3204	3075	−129
Max. ZrH temperature, K	491	479	−12
Max. fuel von Mises stress, MPa	273	296	+12
Max. radial power peaking	1.10	1.11	+0.01
Fuel pressure-drop	506	356	−150
Fuel max. outlet bulk coolant temperature, K	3118	2964	−154
Fuel avg. outlet bulk coolant temperature, K	2746	2749	+3

. The analysis in this section demonstrates that T/M feedback has a beneficial impact on the T/H performance of the reactor resulting in reduced fuel temperatures and pressure losses; however, it has a slightly negative impact on the T/M performance.

4.2. Friction Factor Sensitivity Study

The analysis detailed in Section 4.1 used the Churchill friction factor correlation to obtain the pressure losses within the FE channels. Previous research has compared various Nusselt number and friction factor correlations and concluded that the variation in the predicted pressure drop between legacy friction factor correlations was acceptable [26]. However, this study did not consider how variations in the pressure drop between correlations would propagate in multichannel problems or be exacerbated by T/M feedback. Additionally, the study did not consider several legacy friction factor correlations which rely on the Reynolds number evaluated at the coolant channel wall conditions ($R{e}_w$), such as the Taylor-II and modified Koo correlation. The Taylor-II correlation was used to analyze the NRX reactor series during the NERVA program [32]. This correlation used the corrected $R{e}_w$ to ensure that predicted friction factors were within ±10 percent of the experimental data [33].

Coupled N+T/H–T/M analyses were performed with several friction factor correlations applied to the FE T/H solution. The considered friction factor correlations were detailed in previous research [6]. No orificing was applied to the FEs in order to allow the mass flow to redistribute naturally within the core. A summary of the analysis is detailed in

Table 4. Friction factor sensitivity study summary.

	With Wall Correction		Without Wall Correction
Parameter	Taylor-II	Modified Koo	Taylor-I	Churchill	McAdams
Max. fuel temperature, K	3049	3052	3091	3075	3075
Max. ZrH temperature, K	481	483	479	479	478
Max. fuel von Mises stress, MPa	300	300	293	296	295
Max. radial power peaking	1.114	1.116	1.114	1.114	1.115
Fuel pressure-drop, kPa	871	868	150	356	353
Excess reactivity, pcm	47 ± 7	69 ± 7	23 ± 7	38 ± 7	39 ± 7

. The only global parameter that varies considerably based on the applied friction factor correlation is the pressure drop across the FE channels. Therefore, the authors recommend further research to assess the importance of rebalancing the engine cycle within a multiphysics framework to quantify the impact of the pressure drop on the FE and ME inlet boundary conditions.

5. Summary and Conclusions

Nuclear thermal propulsion systems are being actively investigated by NASA and DARPA as an enabling technology for Mars transit and cislunar operations. The design of such systems will require advanced M&S tools capable of capturing the complex multiphysics feedback present within the reactor core. Previous research has recommended [3,9] and demonstrated [10] that T/M feedback should be accounted for in the multiphysics analysis of NTP reactors. This paper details the thermo-elastic equations solution approach implemented into the ntpThermo code. The approach casts the thermo-elastic equations in the form of the second-order Fredholm integral equation. The primary advantage of this approach is that temperature-dependent properties can be applied within the numerical solution scheme to capture the spatial variation of mechanical properties within each radial node and axial layer. Additionally, the implementation of a coupled N+T/H–T/M sequence into the Basilisk code is described. The coupled N-T/H–T/M sequence was applied to analyze a realistic reactor design with a nonuniform fuel loading pattern and temperature-dependent material properties.

Full-core N+T/H analysis was compared to full-core N+T/H–T/M analysis to demonstrate the importance of T/M feedback. When T/M feedback is neglected, the maximum fuel temperature was overestimated by 129 K, pressure losses across the FE channels were overestimated by 150 kPa, and the maximum VM stress was underestimated by 12 MPa while having only a minor impact on reactivity and local power distributions. The thermal expansion of coolant channels increases the flow area within the FE channels. This significantly reduces the mass velocity and, thus the pressure drop across the FE channels. This reduces the pressure losses within the FE channels. Sensitivity studies revealed that the pressure drop across the FE channels varied considerably on the basis of the friction factor correlation applied to solve the FE channels. Correlations that rely on a corrected Reynolds number, which were used to design the NERVA NRX reactors, predict significantly higher pressure drops. The importance of T/M feedback within a multiphysics solution scheme is case-dependent and can be more or less sensitive in some cases (e.g., higher-power peaking). This paper demonstrates that the T/M feedback should be incorporated into the coupled sequence to accurately capture detailed temperature and stress distributions within the analyzed reactor core.

Future research should include the thermal stress/expansion of MEs and its impact on mass flow distribution within the ME channels.