Integrated Steady-State System Package for Nuclear Thermal Propulsion Analysis Using Multi-Dimensional Thermal Hydraulics and Dimensionless Turbopump Treatment

: Nuclear thermal propulsion is an evolving technology that can be utilized for long-distance space travel. This technology yields the advantage of a high thrust and specific impulse, but requires an examination of the potential design adjustments necessary to enhance its feasibility. The development of nuclear thermal propulsion requires a comprehensive understanding of the system-level behavior during transient and steady-state operation. This paper extends our previous research by including the proper handling of turbomachinery with multi-channel thermal hydraulic simulations only for steady-state solutions. The system-level approach presented here enables the treatment of the turbopump components through non-dimensional analysis that eliminates the assumption of constant efficiencies. All the other components within the system (e.g., reflector and core) can be discretized to multiple channels and layers, in which the full thermal hydraulic solution is established. The approach chosen here enables the realistic modeling of the propellant flow within the expander cycle by capturing the pressure losses, mass flow rate splits, and enthalpy gain for various operational conditions. The verification of the package is completed through point comparisons of previous investigations into similar system designs. Furthermore, sensitivity studies are used to benchmark the capabilities of the package and investigate solution variations due to the perturbation of operational conditions and regimes. The sensitivity studies performed here are important to capture variation in flow characteristics (e.g., temperature, pressure, mass flow rates) for different design objectives such as the thrust and specific impulse. This work demonstrates that system-level simulations lacking multi-channel capability and proper turbomachinery treatment may yield higher uncertainties in understanding the engine’s response and characteristics to changing various requirements. This is extremely important when screening the design space of such propulsion systems and when transient simulations are required.


Introduction
Nuclear thermal propulsion (NTP) has several advantages over traditional propulsion systems for the future of space travel given its high thrust and specific impulse capability [1].The common practice of modern space rocket engines includes variations of chemical propulsion systems yielding specific impulses (I sp ) of ~450 s.For worthwhile implementation, alternatives to this must offer significantly higher I sp given the established nature of this technology [2].Nuclear thermal propulsion devices are targeting specific impulses of about 900 s [3].NTP systems can also provide a reasonably high thrust, which is technology-dependent, but is generally up to 100 klb f , thus providing a boon for their development.The current state of NTP development is focused on the Demonstration Rocket for Agile Cislunar Operations (DRACO) project under an agreement between the Defense Advanced Research Projects Agency (DARPA) and the National Aeronautics and Space Administration (NASA), which is proposed to be an unmanned nuclear thermal rocket that is planned to be demonstrated in 2027 [4].Given the rapid pace at which this technology must be developed, accurate and reliable NTP simulation and modeling is imperative [5].This need was emphasized by the National Academy of Science, Engineering, and Medicine in a study [6] that focused on the requirements to mature NTP technology.
To accurately model an NTP, all components of the system must be coupled along with the neutronics and thermal hydraulics (T/H).In previous related research, system modeling was demonstrated using the Numerical Propulsion System Simulation (NPSS) package [7].In addition, reduced-order system T/H with static turbomachinery (i.e., constant efficiencies) have been used to perform transient simulations of startup and shutdown sequences using the previous iterations of the ntpSystem [8], which is the focus of this paper.These simulation tools have previously been compared to a single case that relied on the Nuclear Engine for Rocket Vehicle Application (NERVA) study [9,10].The NPSS analysis utilized full turbomachinery modeling, including hydrogen mass flow rates, appropriate chamber conditions, and full-system implementation [7].The strength of the NPSS system lies in the treatment of the turbomachinery components.However, the system implementation utilizes 0/1-dimensional analysis via treating each component in the system as single and independent T/H channels.In addition, NPSS lacks the flexibility to define the kinetics of different reactor cores, including temperature reactivity coefficients, delayed neutron fractions, and precursors' decay constants.This led to the development of kinetics capabilities in the ntpSystem code [8], with the emphasis on robust neutronic capabilities including the flexibility to input thermal feedback coefficients, reactivity control, and neutron precursor concentrations throughout various operation scenarios [11].However, that version of the ntpSystem treated each component as 0-dimensional, defined only by the inlet and outlet conditions rather than a proper multi-channel T/H solution.Furthermore, the code lacked the comprehensive nature of turbomachinery implementation and relied on simple models, which assumed constant efficiency values.A deficiency of such a method is that the turbine and pump efficiency values remain constant for various operational conditions.Furthermore, future research focusing on transient simulations [12] requires linking the turbomachinery rotational speed and efficiency with the system solution.
The goal of the current work is to build a comprehensive package, i.e., an upgrade to the ntpSystem, with the addition of an accurate turbomachinery treatment that will remove static efficiency assumptions.In addition, a multi-channel T/H solver is embedded within the framework, which allows the user to define how to solve the various components within the engine, including the nozzle, reflector, moderator, solid nuclear fuel, and piping.The main advantage of the current implementation is the ability to define each component in an arbitrary manner (e.g., having 0-, 1-, or 3-dimensions) and control how these components are coupled and numerically solved.To do so, the T/H solutions are obtained though integrating the ntpThermo package [13] directly into the ntpSystem.The ntpThermo package was developed as an engineering tool to perform multi-channel full core T/H analysis specifically tailored for coupled neutronic and T/H simulations.The package was previously validated [14] and code-to-code verified [13].The integration of turbomachinery together with the flexibility to perform multi-channel T/H simulations of the various components provides a comprehensive full system analysis [15], enabling an accurate steady-state solution that can be extended for transient reactor operations.In addition, having accurate boundary conditions can assist other higher fidelity simulation tools to capture the realistic behavior of the NTP engine.A noteworthy example is the utilization of the Multiphysics Object Oriented Simulation Environment's (MOOSE) tools for high-fidelity multiphysics simulations [16].Specifically, that work used Griffin [17] and RELAP-7 [18] to model the expander cycle.However, the latter research did not focus on the equivalent turbomachinery fidelity.
Energies 2024, 17,3068 3 of 26 The steady-state system implementation and solution is verified to ensure validity through a comparison with results obtained using NPSSs utilizing the large NERVA model as well as thorough sensitivity studies.NTP development has previously utilized assumptions derived from previous NERVA-related efforts [19].For example, current regulation requirements and trends favor the use of High-Assay Low-Enriched Uranium (HALEU [19]) as opposed to the previous NERVA designs that relied on High-Enriched Uranium (HEU) fuel.This transition has a strong impact on the core design [20] and requires the use of highly efficient moderating elements [21] that may alter the characteristic of the system's flow.More specifically, the boundary conditions of the different components within the flow cycle are strongly correlated with the ratio between these moderating elements and the nuclear fuel [22].Among the most visible characteristic changes are the component-wise heat deposition (including its spatial variation) and the pressure drops across the different components, which determines the mass flow rate splits.This paper aims to unfold some of the most important dependencies and clarify the validity of common assumptions through performing multiple sensitivity studies.Therefore, investigations of varying pump pressure ratios, the specific impulse, the thrust, the neutronic power deposition through moderating element/fuel element ratios, and multi-channel T/H comparisons are carried out.These analysis methods are vital not only for the system validity, but also for design optimization given that every component within the system experiences property alterations as a result.Furthermore, clarity on design optimization is provided from these studies through presenting the resulting system behavior.For example, the pump pressure ratio has previously utilized acceptable values to achieve a solution, but sensitivity studies examine the optimization to provide design guidance.The revised ntpSystem package provides valuable high-resolution information that yields realistic boundary conditions which are vital in full core steady-state simulations.In addition, multiphysics investigations can be conducted to determine how the thermomechanical feedback of various materials may hinder system performance.The finalized steady-state solution can then be implemented to initialize transient simulations of reactor startup, operation, and shut-down scenarios in order to provide a comprehensive NTP simulation.

Materials and Methods
This section presents a description of the implemented methodology to obtain the steady-state solution.The breakdown consists of a discussion on the propellant flow pathway (denoted here as flow circuits) in Section 2.1, beginning with the hydrogen tank and concluding with the nozzle chamber.Then, the solution scheme for evaluating the properties throughout the system is presented in Section 2.2 in a step-by-step manner.Finally, an explanation of the treatment of the turbomachinery components is provided in Section 2.3.

Circuit Flow
The circuit described in Figure 1 is an expander cycle in which the cryogenic tank conditions are assumed to be known.Together with the objective chamber conditions (e.g., pressure and temperature or the specific impulse and thrust), the operational conditions (e.g., total system mass flow rate and power) are determined so that the chamber and tank boundary conditions are matched.The expander cycle consists of a defined number of turbines and pumps that are connected through a shaft with a modifiable gear ratio in the form of turbopumps.This enables the pump and turbine operation to feed into each other, and, for the state point calculation, to match the shaft power.The hydrogen also acts as a coolant for the reactor components and is ejected as a propellant, thus maximizing the thrust while minimizing the material requirements.The circuit flow begins with cryogenic hydrogen that is passed from the pressurized tank at State Point (SP) 1 and through the pumps, where the flow is split evenly based on the number of pumps (SP2), and the flow is then split into two separate flow paths.One of the paths leads into the nozzle and reflector circuit (SP3a), while the other leads into the moderator element circuit (SP3b).The nozzle circuit continues from the nozzle (SP3a) to the reflector (SP4a) before reconvening with the other flow split.The other side of the split flow continues to the moderating element, with a separate calculation being made between the supply channel (SP3b) and return channel (SP4b) of the moderator to distinguish the separate state points.Here, the supply is representative of the channel in which the coolant is provided to the moderator in the center of the element, and the return is the counterflow in the outer radius of the moderator.The two circuits then converge into an outlet to develop a combined head (SP6).The flow then splits to match the number of desired turbines, which, in this case, is two, and in each that flow is split, leading to the Turbine Bypass Control Valve (TBCV) (SP7a/b) and the turbines themselves (SP8a/b and SP9a/b).The turbine side drives the shaft that connects to the pump at the beginning of the entire circuit, thus providing additional head to the pump.The flow then reconvenes (SP10) to feed into the reactor fuel element coolant channels (SP11).Finally, the flow is heated by the fuel elements and enters the nozzle chamber (SP12).
Energies 2024, 17, x FOR PEER REVIEW 4 of 27 cryogenic hydrogen that is passed from the pressurized tank at State Point (SP) 1 and through the pumps, where the flow is split evenly based on the number of pumps (SP2), and the flow is then split into two separate flow paths.One of the paths leads into the nozzle and reflector circuit (SP3a), while the other leads into the moderator element circuit (SP3b).The nozzle circuit continues from the nozzle (SP3a) to the reflector (SP4a) before reconvening with the other flow split.The other side of the split flow continues to the moderating element, with a separate calculation being made between the supply channel (SP3b) and return channel (SP4b) of the moderator to distinguish the separate state points.
Here, the supply is representative of the channel in which the coolant is provided to the moderator in the center of the element, and the return is the counterflow in the outer radius of the moderator.The two circuits then converge into an outlet to develop a combined head (SP6).The flow then splits to match the number of desired turbines, which, in this case, is two, and in each that flow is split, leading to the Turbine Bypass Control Valve (TBCV) (SP7a/b) and the turbines themselves (SP8a/b and SP9a/b).The turbine side drives the shaft that connects to the pump at the beginning of the entire circuit, thus providing additional head to the pump.The flow then reconvenes (SP10) to feed into the reactor fuel element coolant channels (SP11).Finally, the flow is heated by the fuel elements and enters the nozzle chamber (SP12).

Steady-State Solution Scheme
The state solution is obtained through numerical iteration to balance the engine's turbomachinery to satisfy the cryogenic tank and chamber boundary conditions.The rocket's geometry, including the number of fuel and moderating elements, along with the relevant dimensions (e.g., engine component lengths and diameters) must be defined.There are various sources that describe how to obtain the complete solution of the expander cycle (or similar cycles), but some steps are often omitted or only partially described [23][24][25].A complete and self-contained description is provided in this chapter, mainly for reproducibility purposes.The flowchart presented in Figure 2 and the following brief procedure are presented to guide through the critical steps that are described in the following subsections.It is important to note that the current implementation performs system-level simulations in a standalone manner through requiring the user to define all components necessary for a solution, such as the fuel channels and moderating elements.The data for each component are also provided by the user and include the geometry description (e.g., coolant channels, pitch), materials used (including properties), as well as the power or heat flux distributions.

Steady-State Solution Scheme
The state solution is obtained through numerical iteration to balance the engine's turbomachinery to satisfy the cryogenic tank and chamber boundary conditions.The rocket's geometry, including the number of fuel and moderating elements, along with the relevant dimensions (e.g., engine component lengths and diameters) must be defined.There are various sources that describe how to obtain the complete solution of the expander cycle (or similar cycles), but some steps are often omitted or only partially described [23][24][25].A complete and self-contained description is provided in this chapter, mainly for reproducibility purposes.The flowchart presented in Figure 2 and the following brief procedure are presented to guide through the critical steps that are described in the following subsections.It is important to note that the current implementation performs system-level simulations in a standalone manner through requiring the user to define all components necessary for a solution, such as the fuel channels and moderating elements.The data for each component are also provided by the user and include the geometry description (e.g., coolant channels, pitch), materials used (including properties), as well as the power or heat flux distributions. 1.The boundary conditions and performance characteristics are defined (Section 2.2).
(a) The engine performance is defined by the specific impulse and thrust or directly using chamber pressure and temperature.(b) The system code calculates the required system mass flow rate and power (Section 2.2.1). 2. The front and back circuits are defined according to Figure 1, and both are solved separately.
(a) Nozzle reflector and moderator circuit (SP1-SP7a/b).This circuit follows the propellant flow through the pump (Section 2.3.1),nozzle, reflector, piping, and moderator elements.The operational conditions at SP7a/b and SP8a/b are known, but the mass flow rate split between the turbine ( ) and the bypass  valve is not known at this stage.(b) Solid fuel elements circuit (SP9a-SP13).The solution starts with the known chamber conditions, and the problem is solved in reverse order for the solid fuel core, as well as the fuel pipes connecting the core and the turbine.After solving this circuit, the operational conditions at SP10 are known.The later state point is a mixing junction between the mass flow rate from the turbine and the bypass (Equation (1)). 1.
The boundary conditions and performance characteristics are defined (Section 2.2).

(a)
The engine performance is defined by the specific impulse and thrust or directly using chamber pressure and temperature.(b) The system code calculates the required system mass flow rate and power (Section 2.2.1).

2.
The front and back circuits are defined according to Figure 1, and both are solved separately.
(a) Nozzle reflector and moderator circuit (SP1-SP7a/b).This circuit follows the propellant flow through the pump (Section 2. Solid fuel elements circuit (SP9a-SP13).The solution starts with the known chamber conditions, and the problem is solved in reverse order for the solid fuel core, as well as the fuel pipes connecting the core and the turbine.After solving this circuit, the operational conditions at SP10 are known.The later state point is a mixing junction between the mass flow rate from the turbine and the bypass (Equation (1)).
where h SP10, mix , h SP9a , and h bypass are the values for the mixing junction (SP10), outlet turbine (SP9a/b), and bypass branch.At SP10, the mass flow rate recombines to the system mass flow rate.
(c) As the turbine and bypass mass flow rates are not known in advance, an iterative solution is applied to vary the turbine outlet enthalpy (h SP9a ).Given the inlet and iteratively changed outlet conditions to the turbine, the solution described in Section 2.3.2 yields .m t and .m bypass .The iterative procedure continues until the right-hand side of Equation ( 20) matches the known mixing enthalpy at SP10 (left-hand side of Equation ( 20)).
The desired chamber temperature (T c ) and pressure (P c ) are used to calculate the overall initial system mass flow rate [22,25].The solution sequence requires the user to input the required tank and chamber conditions to evaluate the system conditions.
m sys is the mass flow rate of the system, Q sys is the system thermal power, h c , is chamber enthalpy, and h tank is the hydrogen tank enthalpy.The user can specify either the required power or the mass flow rate.Alternatively, the specific impulse and thrust can be used [21] to calculate the flow conditions [22], such as the required power and mass flow.The throat pressure is calculated via setting the Mach number to unity, and the chamber conditions are evaluated accordingly.A brief description of this methodology is provided in Section 2.2.1.The use of the thrust and specific impulse requires knowing the nozzle throat area/exit area and the nozzle expansion ratio.These specific elements are vital because the system is evaluated at both ends of the expander cycle to converge on unknown properties in evaluating the turbopump solution.As shown in Figure 1, the mass flow rate is split between the nozzle and moderator circuits.The ntpSystem package enables the user to choose between a pre-determined constant split ratio (denoted by χ) or converge on the mass flow rate distribution by requiring the same pressure losses across the shared plenums.This capability is described in more detail in Section 2.4.For generality, the outlet enthalpy of each component (j) is evaluated via Equation (3).
where h in and h out denote the inlet and outlet enthalpies, respectively, and PF and .
m are the power fraction and mass flow rate of the specific component j.In the case where a constant split ratio is used, the mass flow rate in the nozzle circuit is simply χ × .m sys .The flow through the nozzle continues to the reflector (SP4a5a) as it gains additional enthalpy.Regardless of if a user-defined split ratio is used or calculated, the mixed enthalpy between the moderator and nozzle/reflector circuits can be found using Equation (4).

Nozzle Flow and Throat Conditions
The methodology discussed relies on the assumption that the flow at the nozzle is choked [25].The flow is established at the throat of the nozzle, as characterized in Equation (5). .
where A t is the throat area, P c is the stagnation chamber pressure, and T c is the chamber temperature.γ is the ratio of the specific heats and R is the universal gas constant, both of which are presented in Equation (6).
Energies 2024, 17, 3068 where c p and c v represent the specific heat capacity at a constant pressure and volume, respectively.MW represents the molecular weight of the propellant (H 2 ) and R ′ is the universal gas constant.The velocity at the throat is sonic, meaning that the Mach number (velocity relative to the velocity of the sound) is 1.0.Satisfying the choked flow conditions allows the calculation of the needed mass flow, temperature, and pressure at the chamber.

Turbomachinery
The mathematical relations pertaining to the turbopump methodology are provided in refs.[26,27].Although complete, the methodology is cumbersome to implement as the values for many constants are either not provided or not justified.A more practical notation is given in [23], followed by some specific useful implementation tips which are provided in [24].For completeness, consistency, and reproducibility, the upcoming sections describe the complete set of procedures needed to implement the treatment of the turbomachinery in conjunction with the fuel, reflector, nozzle, moderator, and pipe components.

Steady-State Pump
The pump provides the necessary chamber pressure and thrust while also enabling the operation of the turbines via surmounting the frictional and acceleration pressure losses throughout the engine channels.The steady-state methodology pump and turbine performances are implemented in the current research with a main reliance on refs.[23,24].The solution adopts the use of dimensionless parameters.The main terms needed to describe a turbopump are known as the similarity parameters, which are represented by the specific speed, n s , and specific diameter, d s , and are described in Equation (7) and Equation (8), respectively.Later, the concept of the suction specific speed s s also appears (as shown in Equation ( 12)).
where ω is the rotational speed of the rotor in rad/s, .
V is the volumetric flow rate (= .m p /ρ in ), and the pump head rise is H p .Additionally, the general form of the specific diameter is represented by the impeller diameter, D.
The inlet conditions, among which are the inlet pressure, P in , and propellant density, ρ in , together with the specified pressure ratio or outlet pump pressure are used to calculate the pump head (in meters) as follows: where g is the gravitational constant.The inlet enthalpy and density for the pump are evaluated from the known inlet pressure and temperature values.Thereafter, the netpositive suction head (NPSH), H sr , that is required by the pump to prevent cavitation is calculated using Equation (10).
H sr = P in min − P v gρ in (10) where P v is the coolant vaporization pressure in Pascals and P in min = P in × 1 − P marg is the product of the pump inlet pressure (P in ), with a user-defined margin (P marg ) of typically 10% to account for various uncertainties (e.g., thermophysical properties).The actual pump head and NPSH are used to define the Thoma cavitation parameter (σ cav ) according to Energies 2024, 17, 3068 8 of 26 An issue that can arise within the pump is cavitation, particularly when the pressure at the pump inlet is lower than the vapor pressure of the fluid.An important parameter that can be used to determine if cavitation will occur is the suction specific speed (Equation ( 12)), which, if greater than 5.0, indicates that cavitation will likely occur.
Equation ( 12) is only presented for completeness, but is not actually used in the analysis.In practice, the user is responsible for defining the maximum suction specific speed (s s ).Combined with the Thoma parameter evaluated by Equation ( 11), the maximum specific speed (n sp,max ) is interpolated using pre-generated tabulated data [23], as depicted in Figure 3.The figure shows all the possible solutions tabulated as a function of the dimensionless specific speed (n s ).The Thoma cavitation parameter can be set as the upper constraint, drawn as the dashed horizontal line in Figure 3.Then, a feasible solution that satisfies the maximum suction specific speed is found (illustrated by the vertical dashed line in Figure 3).The solution must satisfy the required suction specific speed, and, in some cases, no solution can be obtained.
the product of the pump inlet pressure ( ), with a user-defined margin ( ) of typically 10% to account for various uncertainties (e.g., thermophysical properties).
The actual pump head and NPSH are used to define the Thoma cavitation parameter ( ) according to An issue that can arise within the pump is cavitation, particularly when the pressure at the pump inlet is lower than the vapor pressure of the fluid.An important parameter that can be used to determine if cavitation will occur is the suction specific speed (Equation ( 12)), which, if greater than 5.0, indicates that cavitation will likely occur.
Equation ( 12) is only presented for completeness, but is not actually used in the analysis.
In practice, the user is responsible for defining the maximum suction specific speed ( ).
Combined with the Thoma parameter evaluated by Equation ( 11), the maximum specific speed ( , ) is interpolated using pre-generated tabulated data [23], as depicted in Figure 3.The figure shows all the possible solutions tabulated as a function of the dimensionless specific speed ( ).The Thoma cavitation parameter can be set as the upper constraint, drawn as the dashed horizontal line in Figure 3.Then, a feasible solution that satisfies the maximum suction specific speed is found (illustrated by the vertical dashed line in Figure 3).The solution must satisfy the required suction specific speed, and, in some cases, no solution can be obtained.The next interpolation step is to obtain the efficiency for a set of specific speeds and specific diameters.In the current research, pre-generated tabulated data taken from Ref. [23] are used.The data are shown as a diagram of contours for different efficiency values (Figure 4).The maximum specific speed ( , ) obtained in the previous step is used to set the rightmost value for the specific speed (denoted by the red dashed line in Figure 4).Then, the implemented procedure iteratively seeks the optimum efficiency to the left of  , . This iterative procedure uses a log-log interpolation scheme between the The next interpolation step is to obtain the efficiency for a set of specific speeds and specific diameters.In the current research, pre-generated tabulated data taken from ref. [23] are used.The data are shown as a diagram of contours for different efficiency values (Figure 4).The maximum specific speed (n sp,max ) obtained in the previous step is used to set the rightmost value for the specific speed (denoted by the red dashed line in Figure 4).Then, the implemented procedure iteratively seeks the optimum efficiency to the left of n sp,max .This iterative procedure uses a log-log interpolation scheme between the contours, and the outcome is a set of values that represent the specific speed and diameter as well as the pump efficiency (η p ).
The specific speed can be used in Equation (7) to extract the rotor speed, ω, in radians per second.Additionally, the impeller diameter, D, in meters is found through knowing the specific diameter, d sp , and inverting Equation ( 8).It must be pointed out that the package includes multiple pump types, such as radial and axial, with the control parameters not described here but provided in ref. [23].The specific speed can be used in Equation ( 7) to extract the rotor speed, ω, in radians per second.Additionally, the impeller diameter, D, in meters is found through knowing the specific diameter, d , and inverting Equation ( 8).It must be pointed out that the package includes multiple pump types, such as radial and axial, with the control parameters not described here but provided in Ref. [23].
Finally, the work performed by the pump (W ) in Watts can be found from Equation ( 13), while the shaft work necessary to drive the pump can be found from Equation (14).
The outlet enthalpy found after passing through the pump is then found using the pump work, which can then be used to determine the properties at the pump outlet.

Steady-State Turbine
Many similarities reside across the methodology between the turbine and the pump.The treatment of the turbine is largely taken from Refs.[23,24].The solution begins by finding the head across the turbine as follows: Please note that both the inlet and outlet enthalpies are assumed to be known.The description about how to iterate on the turbine's outlet enthalpy is provided at the beginning of Section 2.2.
The system has the capability of utilizing (or not) a gear ratio (), which represents the ratio between the pump and turbine.The turbine velocity, in radians per second, is simply  =  × .When the gear ratio is unity ( = 1), both the turbine and pump have identical angular velocities ( =  ).
The next step is to evaluate the maximum allowable turbine diameter (in meters) via using the material-dependent yield stress ( ) as follows: Finally, the work performed by the pump (W p ) in Watts can be found from Equation ( 13), while the shaft work necessary to drive the pump can be found from Equation (14).
The outlet enthalpy found after passing through the pump is then found using the pump work, which can then be used to determine the properties at the pump outlet.

Steady-State Turbine
Many similarities reside across the methodology between the turbine and the pump.The treatment of the turbine is largely taken from refs.[23,24].The solution begins by finding the head across the turbine as follows: Please note that both the inlet and outlet enthalpies are assumed to be known.The description about how to iterate on the turbine's outlet enthalpy is provided at the beginning of Section 2.2.
The system has the capability of utilizing (or not) a gear ratio (Gr), which represents the ratio between the pump and turbine.The turbine velocity, in radians per second, is simply ω t = ω p × Gr.When the gear ratio is unity (Gr = 1), both the turbine and pump have identical angular velocities (ω t = ω p ).
The next step is to evaluate the maximum allowable turbine diameter (in meters) via using the material-dependent yield stress (σ y ) as follows: where ξ is the safety factor (fixed to 1.2 here), S is the rotor shape/geometric factor (fixed to 0.2 here), ω is the rotational speed, and D is the diameter of the turbine.The same analysis can also be applied to the pump.Additionally, the rotor's yield stress, σ y , and material density, ρ r , are presented in Table 1 (the data rely on ref. [24]).The package also enables the user to input their own material with a unique density and yield stress.As the objective is to maximize the turbine's efficiency, the volumetric flow rate of the propellant through the turbine should be higher than that in the pump assembly.High propellant volume flow rates will require high rotor speeds, resulting in high rotor stress levels due to the centrifugal forces.Therefore, it may be necessary to operate the turbine at suboptimal rotor speeds to ensure the rotor disk and blade root stresses are within the acceptable thresholds [23].
The solution approach here aims to match the shaft power obtained in Section 2.3.1.This requires an iterative solution procedure being applied to the volumetric flow rate (and hence the mass flow rate) admitted through the turbine.In practice, the volumetric flow rate .V = .m t /ρ is iterated until it satisfies the same shaft power.The density can be evaluated at the inlet, outlet, or as an average, as the inlet and outlet turbine operational conditions are assumed to be known.The mass flow rate admitted through each turbine, .m t , must be smaller or equal to the system mass flow rate scaled by the number of turbines in the expander cycle.For consistency, we will denote a solution at a specific iteration as .
Similar to the pump treatment, the optimum efficiency, η (i) p , is obtained along with the specific diameter from the pre-generated data taken from ref. [23]; these are presented in Figure 5.The specific speed, n (i) st , from Equation ( 18) is represented by the vertical dashed red line.The upper constraint for the specific diameter, d st,max , is calculated using Equation (19) and is represented by the horizontal dashed line in Figure 5.The maximum turbine efficiency, η (i) t , and specific diameter, d (i) st , are found based on a log-log interpolation.
The iteration-wise specific diameter, d (i) st , and speed, n st , together with the turbine efficiency, η (i) t , are then used to calculate the predicted shaft work as follows: The iterative procedure continues until W (i) sha f t matches the value of W sha f t found in Equation (14).Upon convergence, the volumetric flow rate through the turbine is used to calculate the mass flow rate through the turbine and the bypass valve as follows: m sys are the mass flow rates (kg/s) in a single turbine, bypass valve, and system flow, respectively, while N represents the number of turbopumps.The turbine impeller diameter is calculated from the following: It is important to note that the user can also provide the impeller diameter, in which case, the gear ratio will be used to calculate the turbine rotational speed, which in turn is used to find the specific speed and diameter.The specific speed and diameter are then used directly to extract the efficiency.The procedure described still requires an iterative approach to match the shaft work.The iteration-wise specific diameter,  ( ) , and speed,  ( ) , together with the turbine efficiency,  ( ) , are then used to calculate the predicted shaft work as follows: The iterative procedure continues until  ( ) matches the value of  found in Equation ( 14).Upon convergence, the volumetric flow rate through the turbine is used to calculate the mass flow rate through the turbine and the bypass valve as follows: where  ,  , and  are the mass flow rates (kg/s) in a single turbine, bypass valve, and system flow, respectively, while  represents the number of turbopumps.The turbine impeller diameter is calculated from the following: It is important to note that the user can also provide the impeller diameter, in which case, the gear ratio will be used to calculate the turbine rotational speed, which in turn is used to find the specific speed and diameter.The specific speed and diameter are then used directly to extract the efficiency.The procedure described still requires an iterative approach to match the shaft work.

ntpThermo Background and Integration
The ntpThermo code is a reduced order code specifically designed to perform efficient coupled simulations [28,29].The code was developed using the python programming language due to its open-source licensing/libraries and flexible data structures.This code uses a modular framework for the construction of complex solution schemes from relatively simple classes.The code was tailored to address the complexities of a general NTP system, such as the propellant counter flow or the heat transfer between adjacent solid regions (e.g., moderating and fuel elements).A conceptual illustration of a complete setup model is presented in Figure 6.
An object-oriented environment allows the user to define custom material properties and general correlations.Alternatively, the database includes various pressure-and temperature-dependent properties as well as the friction factor and heat transfer correlations.The problem definition begins with the most basic constitute, denoted as the Element.It contains information about the element geometry and materials, such as the dimensions, number of coolant channels, geometry type (e.g., hexagonal, rectangle, circular), etc.Several element definitions can be used to define a thermal hydraulic channel (denoted as the Channel), each having a unique number of layers and boundary conditions.Then, multiple channels can be used to define the final model (denoted as Core), which is the largest structure that can be constructed.Once the full core is constructed, it can be solved using a variety of convergence schemes (e.g., the mass flow).An object-oriented environment allows the user to define custom material properties and general correlations.Alternatively, the database includes various pressure-and temperature-dependent properties as well as the friction factor and heat transfer correlations.The problem definition begins with the most basic constitute, denoted as the Element.It contains information about the element geometry and materials, such as the dimensions, number of coolant channels, geometry type (e.g., hexagonal, rectangle, circular), etc.Several element definitions can be used to define a thermal hydraulic channel (denoted as the Channel), each having a unique number of layers and boundary conditions.Then, multiple channels can be used to define the final model (denoted as Core), which is the largest structure that can be constructed.Once the full core is constructed, it can be solved using a variety of convergence schemes (e.g., the mass flow).
The ntpThermo code is directly integrated into the system level package.All the system components except the turbine and pump are defined separately using the ntpThermo input deck.These components include pipes, the reflector, nozzle, moderating elements, and fuel elements.This flexible definition allows the user to define the different components and solution schemes in the most generic manner.For example, the user can treat the core as a single T/H channel or as multiple channels.Furthermore, the user can choose from various convergence schemes (e.g., mass flow) to realistically model the problem at hand.For example, if all the channels share the same inlet and outlet plenum, the solution scheme will adjust the mass flow rate distribution until the pressure drop is identical within all channels.An additional available convergence scheme that relies on adjusting the orifice can be used to control the mass flow rates such that the outlet propellant temperatures are uniform.This in turn typically leads to a reduced fuel temperature, but also increased pressure drops.Figure 7 presents the setup for the main components used in the current analysis.The nozzle reflector and moderating channels (Figure 7a) were defined using the same core object, and the mass flow rate within each circuit was obtained via enforcing identical pressure drops.Thereafter, the outlet conditions were obtained by weighting the outlet temperatures from the reflector and moderating the elements with their corresponding mass flow rates ( ,  ).The fuel channels can be defined as a separate core object (Figure 7b), and all channels within the fuel are solved in synergy.The latter approach will yield a more realistic pressure drop compared to using a single T/H channel.It must be noted that the connecting pipes (e.g., pump to the nozzle and moderator circuits, and turbine to fuel) are also set using the ntpThermo input deck.The ntpThermo code is directly integrated into the system level package.All the system components except the turbine and pump are defined separately using the ntpThermo input deck.These components include pipes, the reflector, nozzle, moderating elements, and fuel elements.This flexible definition allows the user to define the different components and solution schemes in the most generic manner.For example, the user can treat the core as a single T/H channel or as multiple channels.Furthermore, the user can choose from various convergence schemes (e.g., mass flow) to realistically model the problem at hand.For example, if all the channels share the same inlet and outlet plenum, the solution scheme will adjust the mass flow rate distribution until the pressure drop is identical within all channels.An additional available convergence scheme that relies on adjusting the orifice can be used to control the mass flow rates such that the outlet propellant temperatures are uniform.This in turn typically leads to a reduced fuel temperature, but also increased pressure drops.Figure 7 presents the setup for the main components used in the current analysis.The nozzle reflector and moderating channels (Figure 7a) were defined using the same core object, and the mass flow rate within each circuit was obtained via enforcing identical pressure drops.Thereafter, the outlet conditions were obtained by weighting the outlet temperatures from the reflector and moderating the elements with their corresponding mass flow rates ( . m ME ).The fuel channels can be defined as a separate core object (Figure 7b), and all channels within the fuel are solved in synergy.The latter approach will yield a more realistic pressure drop compared to using a single T/H channel.It must be noted that the connecting pipes (e.g., pump to the nozzle and moderator circuits, and turbine to fuel) are also set using the ntpThermo input deck.
The main advantage of integrating ntpThermo is its ability to calculate the local temperatures and stresses in order to identify potential failures, although the latter was not investigated in this paper.In addition, the ability to perform multi-channel analysis rather than just focusing on a single hot channel model enables the identification of failures in channels with different material properties (e.g., different loading of fuel in the matrix) as opposed to defining the radial power peaking as the limiting case.In other words, cooler channels may be prone to failure if their thermophysical or thermomechanical properties are less favorable when compared to hotter channels.Finally, the correct distribution of the mass flow rate or the required orifice pattern can only be obtained when multi-channel simulations are performed, making this integration highly desirable.The main advantage of integrating ntpThermo is its ability to calculate the local temperatures and stresses in order to identify potential failures, although the latter was not investigated in this paper.In addition, the ability to perform multi-channel analysis rather than just focusing on a single hot channel model enables the identification of failures in channels with different material properties (e.g., different loading of fuel in the matrix) as opposed to defining the radial power peaking as the limiting case.In other words, cooler channels may be prone to failure if their thermophysical or thermomechanical properties are less favorable when compared to hotter channels.Finally, the correct distribution of the mass flow rate or the required orifice pattern can only be obtained when multi-channel simulations are performed, making this integration highly desirable.

Comparison against Large Nerva Published Results
This section focuses on benchmarking the capabilities implemented in this research against the publicly available NTP models developed by the NASA Glenn Research Center (GRC).A specific comparison [7] relies on the Large Nuclear Engine for Rocket Vehicle Application (NERVA) system [30] that utilizes fuel elements made from the ZrC-graphite composite.The components include two RL-60 hydrogen turbopumps, 564 fuel elements, 241 tie tubes (moderating elements), shown in Figure 8, with a length of 132.1 cm, and a reactor power rating of 555 MW, produced in the core, nozzle, and reflector.The specifications of the fuel element include a 1.905 cm flat-to-flat distance containing 19 hydrogen coolant channels with a diameter of 0.257 cm.The fuel element itself is comprised of (U, Zr)C-graphite in a hexagonal pattern with an average coolant channel cladding of 100 µm of Zirconium Carbide.Meanwhile, the tie tube elements are modeled in two flow regions, namely the return and supply lines.The supply line is centered with the hydrogen coolant flow, which is surrounded by Inconel-718 and the Zirconium hydride moderator.The return line surrounds the supply line with another hydrogen flow channel encompassed by another layer of Inconel-718 clad with Zirconium carbide, surrounded by graphite filler.Furthermore, the nozzle and reflector are a hydrogen coolant channel with a 14.7 cm thick Beryllium exterior.

Comparison against Large Nerva Published Results
This section focuses on benchmarking the capabilities implemented in research against the publicly available NTP models developed by the NASA Glenn Research Center (GRC).A specific comparison [7] relies on the Large Nuclear Engine for Rocket Vehicle Application (NERVA) system [30] that utilizes fuel elements made from the ZrC-graphite composite.The components include two RL-60 hydrogen turbopumps, 564 fuel elements, 241 tie tubes (moderating elements), shown in Figure 8, with a length of 132.1 cm, and a reactor power rating of 555 MW, produced in the core, nozzle, and reflector.The specifications of the fuel element include a 1.905 cm flat-to-flat distance containing 19 hydrogen coolant channels with a diameter of 0.257 cm.The fuel element itself is comprised of (U, Zr)Cgraphite in a hexagonal pattern with an average coolant channel cladding of 100 µm of Zirconium Carbide.Meanwhile, the tie tube elements are modeled in two flow regions, namely the return and supply lines.The supply line is centered with the hydrogen coolant flow, which is surrounded by Inconel-718 and the Zirconium hydride moderator.The return line surrounds the supply line with another hydrogen flow channel encompassed by another layer of Inconel-718 clad with Zirconium carbide, surrounded by graphite filler.Furthermore, the nozzle and reflector are a hydrogen coolant channel with a 14.7 cm thick Beryllium exterior.Previous analysis on the system has documented the energy deposition breakdown as well as the pressure, temperature, and mass flow rates at each state point.The models include a conjoined flow for the moderating element (as supply and return lines) as well as the nozzle and reflector (in a conjoined flow line).Furthermore, the piping between components such as the fuel and pump lines are modeled to account for potential pressure losses that can occur in the transportation between components.A direct comparison between the reference state point values and the model values is provided through utilizing a constant mass flow split between the reflector/nozzle line and the moderating element line as 7.53 kg/s and 5.15 kg/s, respectively.By matching the desired chamber conditions of 2794 K, 6.89 MPa, and 12.68 kg/s along with a pump pressure ratio matching of 0.22 Previous analysis on the system has documented the energy deposition breakdown as well as the pressure, temperature, and mass flow rates at each state point.The models include a conjoined flow for the moderating element (as supply and return lines) as well as the nozzle and reflector (in a conjoined flow line).Furthermore, the piping between components such as the fuel and pump lines are modeled to account for potential pressure losses that can occur in the transportation between components.A direct comparison between the reference state point values and the model values is provided through utilizing a constant mass flow split between the reflector/nozzle line and the moderating element line as 7.53 kg/s and 5.15 kg/s, respectively.By matching the desired chamber conditions of 2794 K, 6.89 MPa, and 12.68 kg/s along with a pump pressure ratio matching of 0.22 MPa to 15.65 MPa, 21 K, and 12.68 kg/s, a comparison can be made.The percentage differences are found to be in relatively good agreement (Table 2), given that most of the differences are observed in temperatures.Additionally, the pressure and mass flow rates are in exceptional agreement, with the largest discrepancy residing in the turbine.The discrepancy results from the difference in methodology for the turbopump calculation.The reference utilizes a slightly different efficiency.Since the efficiency depends on the converged mass flow rate and required pressure drop for the operation, the turbine properties would not be the same as the reference and should rather match closer to the realities of the intricate nature of the turbopumps.Another discrepancy can be noted at the conjunction between the reflector line and moderator return line.The flow rates sum from 5.15 kg/s and 7.53 kg/s to 12.68 kg/s, which is expected.However, the pressure from the reflector is said to be 13.61MPa and the moderator is 14.5 MPa, while the conjunction has a pressure of 13.61 MPa.The pressure at the conjunction must be identical in both circuits.The solution calculated by the ntpSystem does consider the pressure inconsistency and accurately predicts the boundary conditions at the turbine inlet.In spite of the inconsistency reported by the NPSS, the largest percentage difference found in the relevant state points reported in Table 2 is 6.3%, which, for the purposes of the system code, is considered to be satisfactory.Additional comparisons specific to the turbopump components are presented in Table 3, which shows key factors such as the volumetric flow rate and impeller diameter.

System Sensitivity Studies
Sensitivity studies have been conducted in terms of perturbing several inputs that are vital in mission analysis to improve the holistic understanding of the engine, identify trade-offs, and possibly even converge on optimal design parameters.Some crucial design parameters ascertained include the convergence scheme, compression ratio, specific impulse, thrust, power share, and power shape; further descriptions of each terminology are provided in their respective sections.The main objectives utilized are the specific impulse and thrust because their inputs alter the chamber conditions, which propagate into many other components; because of this, physical feasibility must be ensured given that design limitations, such as excessive temperatures or insufficient enthalpy, may result.Maintaining consistency across all comparisons is imperative because differing multi-variable operations yield misleading and uncharacteristic results; therefore, in all scenarios investigated, all variables are held constant unless stated otherwise.One example of this includes the power shape, in which all perturbations consist of constant and uniform power generation within each flow channel, thus preventing analytical misdirection.Finally, in all the sensitivity studies, the heat fluxes are provided from external neutronic analyses and are fixed.In other words, the system-level simulations are performed in a standalone manner, but future research will couple the neutronic feedback with the system-level simulations.

Convergence Scheme
The ntpSystem has the capability for the user to specify fluid flow convergence methods.The fixed flow rate is based on the reference values to the moderating element and a nozzle/reflector of 5.15 and 7.53 kg/s; in this scenario, the mass flow rate is unable to change during the simulation (denoted here as no convergence).An alternative solution technique, denoted here as the mass convergence, allows the flow to be split naturally by ensuring the pressure drop across the parallel channels is identical.It must be pointed out that the latter approach is the correct one when different T/H channels share the same plenum, and thus the solution scheme must ensure identical pressure drops.In essence, using these different methods will generally result in a different mass flow rate split and may alter the local T/H parameters (e.g., surface/wall and propellant temperatures).This section focuses on demonstrating the variations in the results when a proper mass convergence scheme is not used.A reference case (presented in Section 3) is chosen using the performance parameters of 900 s Isp, a thrust of 25 klbf, and a pumping pressure ratio (PPR) of 70.Then, each parameter is perturbed separately from the reference point (by about 20%), and the two solution schemes, i.e., cases with (w) and without (w/o) mass convergence schemes, are executed and compared.Table 4 displays the pressure and mass flow rate within the different components (e.g., nozzle) of the system.It can be noticed that the results can be less sensitive for some specific cases (e.g., the reference point), but the solution method does have an impact on the flow (e.g., variation of Isp) within the system.There is no way to determine this sensitivity in advance but to use the mass flow rate convergence scheme.The comparison of convergence schemes is continued in Figure 9, utilizing the axial propellant bulk temperatures within the moderating element circuit (i.e., supply and return).This is to illustrate the local effects within the T/H channels and emphasize the need of the entire system solution as a means to determine the safety margins in various components.Although in some of the perturbed cases (e.g., variation of thrust) the deviations appear rather minor, the non-convergence scheme produces non-physical results as the pressure differential across the channels must be identical.

Compression Pressure
The first action the flow encounters is the transition from the cryogenic hydrogen tank through the pump.The pump increases the pressure sufficiently to provide adequate head to enforce the flow rate needed to drive the system.The output pressure of the pump is enforced by the user input for the pump pressure ratio, which is the ratio of the pressure after the coolant leaves the pump to the pressure in the cryogenic tank.The observational analysis of the system output is enabled through the adjustment of this pressure ratio.
Figure 10a shows that the system operates at near-constant turbine efficiencies with the increase in the pump pressure ratio and a direct decrease in the pump efficiency.Although there is only a slight decrease in the pump efficiency, this decrease would be extrapolated as the pump pressure ratios further increase.As expected, the shaft work increases (Figure 10b) with the pump pressure ratio, but this creates a decrease in the pump efficiency, given as the volumetric flow, and hence the mass flow rate admitted through

Compression Pressure
The first action the flow encounters is the transition from the cryogenic hydrogen tank through the pump.The pump increases the pressure sufficiently to provide adequate head to enforce the flow rate needed to drive the system.The output pressure of the pump is enforced by the user input for the pump pressure ratio, which is the ratio of the pressure after the coolant leaves the pump to the pressure in the cryogenic tank.The observational analysis of the system output is enabled through the adjustment of this pressure ratio.
Figure 10a shows that the system operates at near-constant turbine efficiencies with the increase in the pump pressure ratio and a direct decrease in the pump efficiency.Although there is only a slight decrease in the pump efficiency, this decrease would be extrapolated as the pump pressure ratios further increase.As expected, the shaft work increases (Figure 10b) with the pump pressure ratio, but this creates a decrease in the pump efficiency, given as the volumetric flow, and hence the mass flow rate admitted through the turbine must be adjusted (as seen in Figure 10c).The phenomenon being described here is derived from the fact that these new conditions are changing the specific speed and diameter, thus leading to a slightly different efficiency (as shown in Figure 4).A mass flow rate increase is seen in the turbine with a constant mass flow rate in the pump.The turbine efficiency remains constant, so, based on Equation ( 20), the volumetric flow rate increases, and therefore the mass flow rate must also increase.Meanwhile, the total mass flow rate within the system remains constant, thus the increases in the turbine mass flow are balanced by a decrease in the bypass mass flow, explain their opposing trendlines (Figure 10c).
The moderating elements and the nozzle/reflector are situated between the pump and turbine in terms of the system flow; therefore, the pressure increase in the pump also translates to a pressure increase in the moderating elements and nozzle/reflector circuit, as seen in Figure 11.The temperature (Figure 11a) in the components is inverted in comparison to the mass flow rate (Figure 11b) because a higher mass flow rate leads to lower enthalpy (and hence temperature) variation.As expected, the pressure that these components experience (Figure 11c) must increase because the outlet pump pressure increases.
rate increase is seen in the turbine with a constant mass flow rate in the pump.The turbine efficiency remains constant, so, based on Equation ( 20), the volumetric flow rate increases, and therefore the mass flow rate must also increase.Meanwhile, the total mass flow rate within the system remains constant, thus the increases in the turbine mass flow are balanced by a decrease in the bypass mass flow, explain their opposing trendlines (Figure 10c).The moderating elements and the nozzle/reflector are situated between the pump and turbine in terms of the system flow; therefore, the pressure increase in the pump also translates to a pressure increase in the moderating elements and nozzle/reflector circuit, as seen in Figure 11.The temperature (Figure 11a) in the components is inverted in comparison to the mass flow rate (Figure 11b) because a higher mass flow rate leads to lower enthalpy (and hence temperature) variation.As expected, the pressure that these components experience (Figure 11c) must increase because the outlet pump pressure increases.The moderating elements and the nozzle/reflector are situated between the pump and turbine in terms of the system flow; therefore, the pressure increase in the pump also translates to a pressure increase in the moderating elements and nozzle/reflector circuit, as seen in Figure 11.The temperature (Figure 11a) in the components is inverted in comparison to the mass flow rate (Figure 11b) because a higher mass flow rate leads to lower enthalpy (and hence temperature) variation.As expected, the pressure that these components experience (Figure 11c) must increase because the outlet pump pressure increases.The main observation here is that simulating this system is largely dependent on the pressure ratio, which might affect the flow characteristic within the analyzed system.In reality, however, the system should operate at the lowest possible pressure ratio value as it will require reduced shaft work and satisfy the needed chamber conditions, thus also satisfying the specific impulse and thrust.

Specific Impulse
A common practice is that the maximum achievable specific impulse is preferential given the increased utilization of the propellant; however, it is evident in Figure 12 that the chamber outlet temperature increases quadratically with the specific impulse.The lat- The main observation here is that simulating this system is largely dependent on the pressure ratio, which might affect the flow characteristic within the analyzed system.In reality, however, the system should operate at the lowest possible pressure ratio value as it will require reduced shaft work and satisfy the needed chamber conditions, thus also satisfying the specific impulse and thrust.

Specific Impulse
A common practice is that the maximum achievable specific impulse is preferential given the increased utilization of the propellant; however, it is evident in Figure 12 that the chamber outlet temperature increases quadratically with the specific impulse.The latter is expected as the specific impulse is known to be proportional to the square root of the chamber temperature [21].Therefore, increases in the specific impulse result in even larger increases in the temperature.One key consideration in designing an NTP core is the thermo-mechanical safety margin of the solid components (e.g., nuclear fuel solid elements).A high propellant temperature indicates that these solid components will exhibit even higher temperatures and potentially prohibitive stresses, which may require limiting the specific impulse.
The main observation here is that simulating this system is largely dependent on the pressure ratio, which might affect the flow characteristic within the analyzed system.In reality, however, the system should operate at the lowest possible pressure ratio value as it will require reduced shaft work and satisfy the needed chamber conditions, thus also satisfying the specific impulse and thrust.

Specific Impulse
A common practice is that the maximum achievable specific impulse is preferential given the increased utilization of the propellant; however, it is evident in Figure 12 that the chamber outlet temperature increases quadratically with the specific impulse.The latter is expected as the specific impulse is known to be proportional to the square root of the chamber temperature [21].Therefore, increases in the specific impulse result in even larger increases in the temperature.One key consideration in designing an NTP core is the thermo-mechanical safety margin of the solid components (e.g., nuclear fuel solid elements).A high propellant temperature indicates that these solid components will exhibit even higher temperatures and potentially prohibitive stresses, which may require limiting the specific impulse.In addition, the DRACO program does not necessarily target a specific impulse of 900 s as it only aims to demonstrate the feasibility of the NTP technology in space and collect valuable data [4] (e.g., operation at high temperatures).Designs having different specific impulse values would determine the boundary conditions, mass flow rate splits, and general flow characteristics.
Increasing the specific impulse requires increasing the power of the system (Figure 13a), but, at the same time, the required shaft work is reduced (Figure 13b) as a smaller mass flow rate is required to be admitted (Figure 14a) through the entire system.The lower work requirements to drive the turbopump component also results in a different breakdown of the mass flow rates through the turbine (Figure 14a), with less flow in the turbine itself.The mass flow rate drop is correlated with the pressure differential decrease mainly In addition, the DRACO program does not necessarily target a specific impulse of 900 s as it only aims to demonstrate the feasibility of the NTP technology in space and collect valuable data [4] (e.g., operation at high temperatures).Designs having different specific impulse values would determine the boundary conditions, mass flow rate splits, and general flow characteristics.
Increasing the specific impulse requires increasing the power of the system (Figure 13a), but, at the same time, the required shaft work is reduced (Figure 13b) as a smaller mass flow rate is required to be admitted (Figure 14a) through the entire system.The lower work requirements to drive the turbopump component also results in a different breakdown of the mass flow rates through the turbine (Figure 14a), with less flow in the turbine itself.The mass flow rate drop is correlated with the pressure differential decrease mainly in the reflector and moderator component, as seen in Figure 15.This figure also shows that the pressure drop in the fuel elements is not greatly affected.Finally, Figure 16a shows the axial variation of the solid fuel surface temperature within the fuel elements, which follows the propellant temperature (Figure 16b).As the power of the overall system was increased and the system's mass flow rate decreased, a larger temperature gain is seen for higher specific impulse values.in the reflector and moderator component, as seen in Figure 15.This figure also shows that the pressure drop in the fuel elements is not greatly affected.Finally, Figure 16a shows the axial variation of the solid fuel surface temperature within the fuel elements, which follows the propellant temperature (Figure 16b).As the power of the overall system was increased and the system's mass flow rate decreased, a larger temperature gain is seen for higher specific impulse values.

Thrust
The thrust in an NTP system is controlled in tandem with the specific impulse.Generally speaking, the higher the specific impulse, the less propellant is required to achieve a certain thrust.However, if the Isp is kept constant, achieving a higher thrust requires a higher mass flow rate to be admitted through the system, which in turn leads to increasing the total power of the system (Figure 17a).The higher mass flow rate requirements (Figure 18) also lead to an increased shaft power, as shown in Figure 17b.To satisfy the required shaft work, a higher mass flow rate must be passed via the turbine as well as the turbine bypass control valve (TBCV).A direct effect of feeding higher flow rates is the resulting pressure drop in the different flow circuits within the expander cycle (Figure 19).

Thrust
The thrust in an NTP system is controlled in tandem with the specific impulse.Generally speaking, the higher the specific impulse, the less propellant is required to achieve a certain thrust.However, if the Isp is kept constant, achieving a higher thrust requires a higher mass flow rate to be admitted through the system, which in turn leads to increasing the total power of the system (Figure 17a).The higher mass flow rate requirements (Figure 18) also lead to an increased shaft power, as shown in Figure 17b.To satisfy the required shaft work, a higher mass flow rate must be passed via the turbine as well as the turbine bypass control valve (TBCV).A direct effect of feeding higher flow rates is the resulting pressure drop in the different flow circuits within the expander cycle (Figure 19).
higher mass flow rate to be admitted through the system, which in turn leads to increasing the total power of the system (Figure 17a).The higher mass flow rate requirements (Figure 18) also lead to an increased shaft power, as shown in Figure 17b.To satisfy the required shaft work, a higher mass flow rate must be passed via the turbine as well as the turbine bypass control valve (TBCV).A direct effect of feeding higher flow rates is the resulting pressure drop in the different flow circuits within the expander cycle (Figure 19).higher mass flow rate to be admitted through the system, which in turn leads to increasing the total power of the system (Figure 17a).The higher mass flow rate requirements (Figure 18) also lead to an increased shaft power, as shown in Figure 17b.To satisfy the required shaft work, a higher mass flow rate must be passed via the turbine as well as the turbine bypass control valve (TBCV).A direct effect of feeding higher flow rates is the resulting pressure drop in the different flow circuits within the expander cycle (Figure 19).Finally, Figure 20 shows the axial variation of the solid fuel surface and propellant temperatures within the fuel elements.This figure implies that the power-to-mass flow rate ratio is kept roughly constant when modeling designs with differing thrust values; as a result, the boundary conditions may be fixed.It is observed from Figure 20b  The sensitivity studies previously discussed provide insight into the steady-state so- Finally, Figure 20 shows the axial variation of the solid fuel surface and propellant temperatures within the fuel elements.This figure implies that the power-to-mass flow rate ratio is kept roughly constant when modeling designs with differing thrust values; as a result, the boundary conditions may be fixed.It is observed from Figure 20b that the propellant temperature at the inlet of the fuel channel experiences a difference of about 0.48%.
The sensitivity studies previously discussed provide insight into the steady-state solution of the system with respect to the pump pressure ratio, specific impulse, and thrust.The analysis of each study individually provides a thermomechanical idea of the property development within the system, as each input acts as a hinge on the physical constraints.However, shifting to a design perspective requires an input value for each of the three properties themselves.Therefore, a relationship between all three inputs must be known for an optimized system.As discussed in Section 4.2, the minimum value for the pump pressure ratio that will enable the system's operation is utilized as the boundary condition as it is representative of the lowest power requirement for the turbopump.A convergence scheme was utilized to minimize the pump pressure ratio that provided a steady-state solution in the system for analysis.Evidently, Figure 21 shows that the pump pressure ratio and specific impulse lack a correlative relationship, with a slight increase in the lower specific impulse, whereas the thrust displays a clear direct relationship.The result is unsurprising as the sensitivity study analysis explains the characterization of the direct relationship between the thrust and chamber pressure and, henceforth, the minimum operable pump pressure ratio, while a direct non-linear correlation is present between the specific impulse and chamber temperature, and is therefore relatively unrelated to the minimal compression.Finally, Figure 20 shows the axial variation of the solid fuel surface and propellant temperatures within the fuel elements.This figure implies that the power-to-mass flow rate ratio is kept roughly constant when modeling designs with differing thrust values; as a result, the boundary conditions may be fixed.It is observed from Figure 20b that the propellant temperature at the inlet of the fuel channel experiences a difference of about 0.48%.The sensitivity studies previously discussed provide insight into the steady-state solution of the system with respect to the pump pressure ratio, specific impulse, and thrust.The analysis of each study individually provides a thermomechanical idea of the property development within the system, as each input acts as a hinge on the physical constraints.However, shifting to a design perspective requires an input value for each of the three properties themselves.Therefore, a relationship between all three inputs must be known for an optimized system.As discussed in Section 4.2, the minimum value for the pump pressure ratio that will enable the system's operation is utilized as the boundary condition as it is representative of the lowest power requirement for the turbopump.A convergence scheme was utilized to minimize the pump pressure ratio that provided a steady-state solution in the system for analysis.Evidently, Figure 21 shows that the pump pressure ratio and specific impulse lack a correlative relationship, with a slight increase in the lower specific impulse, whereas the thrust displays a clear direct relationship.The result is unsurprising as the sensitivity study analysis explains the characterization of the direct relationship between the thrust and chamber pressure and, henceforth, the minimum operable pump pressure ratio, while a direct non-linear correlation is present between the specific impulse and chamber temperature, and is therefore relatively unrelated to the minimal compression.

Power Share
The system package allows the power fractions for each component to be user-defined.As a result, the package can be used to analyze the systematic effects resulting from power deposition alterations.In the case presented in Table 5, a comparison is completed through altering the power share between the fuel elements and moderating elements (supply and return combined).The power share of the fuel is altered from 85% to 93% (in the reference case reported in Section 3, the power share was 88.4%).Although this change is arbitrary, it aims to represent cores having considerably higher fuel-to-moderator ratios (i.e., less moderated).The result portrays the expected increases in temperature within the moderating elements and the reduced inlet boundary conditions to the fuel elements.The temperatures, pressures, and mass flow rates reported in Table 5 maintain the tank and chamber conditions that are identical to those reported in Section 3. Additionally, the power deposition alters the minimum pump pressure ratio required to meet the necessary chamber conditions, where the 85% fuel deposition requires a compression ratio of 51, while the 93% fuel deposition lowers this to 46.The crux of this investigation is that the changes in power deposition are affected by neutronic considerations, which in turn are dependent on the ratio of moderating-to-fuel elements and the general core design.Therefore, the system level analysis and neutronic analysis cannot be isolated, but must rather be treated in tandem with each other.

Power Share
The system package allows the power fractions for each component to be user-defined.As a result, the package can be used to analyze the systematic effects resulting from power deposition alterations.In the case presented in Table 5, a comparison is completed through altering the power share between the fuel elements and moderating elements (supply and return combined).The power share of the fuel is altered from 85% to 93% (in the reference case reported in Section 3, the power share was 88.4%).Although this change is arbitrary, it aims to represent cores having considerably higher fuel-to-moderator ratios (i.e., less moderated).The result portrays the expected increases in temperature within the moderating elements and the reduced inlet boundary conditions to the fuel elements.The temperatures, pressures, and mass flow rates reported in Table 5 maintain the tank and chamber conditions that are identical to those reported in Section 3. Additionally, the power deposition alters the minimum pump pressure ratio required to meet the necessary chamber conditions, where the 85% fuel deposition requires a compression ratio of 51, while the 93% fuel deposition lowers this to 46.The crux of this investigation is that the changes in power deposition are affected by neutronic considerations, which in turn are dependent on the ratio of moderating-to-fuel elements and the general core design.Therefore, the system level analysis and neutronic analysis cannot be isolated, but must rather be treated in tandem with each other.

Multi-Channel Analysis
This section describes the impact of multi-channel analysis in predicting realistic boundary conditions and flow characteristics.The multi-channel capability can also be used to predict safety margins, but this is not discussed here.The case presented here relies on the large NERVA case from Section 3, with the exception of using a supercell adopted from ref. [31] rather than a single heated channel.In the large NERVA case, there are 564 fuel elements with a total propellant flow of 12.68 kg/s.Here, there are six channels which are representative of the 564 elements.An arbitrary radial power peaking was assigned with a maximum value of 1.2.The axial power distribution was obtained from the neutronic simulations conducted in ref. [31] and is shown in Figure 22.The mass flow rate is scaled per channel to preserve the total value of 12.68 kg/s.Two simulations are performed here, with the first having uniform power distribution across all six channels (i.e., this case is equivalent to performing a single T/H channel analysis).The second simulation is performed for the power distributions presented in Figure 22.This simulation also leverages the built-in capability [28] of converging on the orifice dimensions so that the exit propellant velocity is uniform.Without these orifices, the outlet temperatures are not uniformly distributed (Figure 23) as larger portions of the mass flow rates are diverted away (Figure 23) from the elements, producing more power to satisfy uniform pressure losses across all channels.The results for the uniform and non-uniform power distributions are presented in Figure 24, showing that the orifices are indeed producing a uniform outlet temperature, but are also introducing additional form pressure losses.The uniform outlet temperatures are achieved via introducing orifices with smaller entry diameters in channels producing less power.The greater flow resistance naturally pushes some of the flow to elements producing more power, thus flattening the fuel temperature distribution, as seen in Figure 25.
An additional comparison is provided in Table 6 where the turbine flow characteristics are presented for both the uniform and non-uniform test cases.As seen from the table, there is a slight difference in the mass flow rate through the turbine.The temperatures and pressures for all other components are very similar and thus are not reported in the table.The main point of this section is to demonstrate the fact that pressure drops might be strongly affected if multi-channel analysis is included but the impact on the system is case-dependent.Yet, without the inclusion of multi-channel modeling, the resulting simulations might be inconsistent with the physical system.tions are presented in Figure 24, showing that the orifices are indeed producing a uniform outlet temperature, but are also introducing additional form pressure losses.The uniform outlet temperatures are achieved via introducing orifices with smaller entry diameters in channels producing less power.The greater flow resistance naturally pushes some of the flow to elements producing more power, thus flattening the fuel temperature distribution, as seen in Figure 25.tions are presented in Figure 24, showing that the orifices are indeed producing a uniform outlet temperature, but are also introducing additional form pressure losses.The uniform outlet temperatures are achieved via introducing orifices with smaller entry diameters in channels producing less power.The greater flow resistance naturally pushes some of the flow to elements producing more power, thus flattening the fuel temperature distribution, as seen in Figure 25.An additional comparison is provided in Table 6 where the turbine flow characteristics are presented for both the uniform and non-uniform test cases.As seen from the table, there is a slight difference in the mass flow rate through the turbine.The temperatures and pressures for all other components are very similar and thus are not reported in the table.The main point of this section is to demonstrate the fact that pressure drops might be strongly affected if multi-channel analysis is included but the impact on the system is case-dependent.Yet, without the inclusion of multi-channel modeling, the resulting simulations might be inconsistent with the physical system.An additional comparison is provided in Table 6 where the turbine flow characteristics are presented for both the uniform and non-uniform test cases.As seen from the table, there is a slight difference in the mass flow rate through the turbine.The temperatures and pressures for all other components are very similar and thus are not reported in the table.The main point of this section is to demonstrate the fact that pressure drops might be strongly affected if multi-channel analysis is included but the impact on the system is case-dependent.Yet, without the inclusion of multi-channel modeling, the resulting simulations might be inconsistent with the physical system.

Summary and Conclusions
Nuclear thermal propulsion has been identified to be a beneficial and crucial technology in the advancement of cis-lunar and interplanetary space travel, given its capability of high thrust and specific impulse yields.Pushes to advance this technology are essential as the upcoming DRACO project is targeting demonstration by 2027 [4].Modeling and simulations are complementary to experimental efforts, but are necessary to mature this technology.Among the computational needs is the ability to model the engine as a whole.Previous works have demonstrated 0/1-dimensional analysis with turbomachinery implementation [7] or transient neutronic implementation, including temperature reactivity coefficients, delayed neutron fractions, and precursor decay constants [8].Additional efforts include MOOSE's fluid flow and heat transfer analysis with minimized computational cost through high-fidelity reduced-order simulations [16].However, achieving a comprehensive system simulation requires the coupling of turbomachinery with multi-dimensional T/H solutions, hence the development of the current iteration of the ntpSystem.
This paper illustrates the key methodologies necessary for the development of an accurate steady-state solution while eliminating or clarifying the assumptions utilized by previous works.Verification and analysis are accomplished through comparisons with the NPSS Large NERVA model in conjunction with sensitivity studies of the convergence method, compression ratio, specific impulse, thrust, power share, and power shape.The purpose of these sensitivity studies is to unfold the behavior of the system for a given set of assumptions and simplifications.The main driving question is how a variation in certain objectives (e.g., required thrust) or constraint changes the operational conditions within the expander cycle (e.g., mass admitted through the turbine or pressure drop across various components).In practice, many reported results disregard or do not justify various design considerations when performing decoupled core-level simulations.For example, modern NTP cores are envisioned to operate using HALEU fuel and require the careful optimization and placement of moderating elements to control the spatial power and achieve the desired critical operational time.As was shown here, the number of moderating elements will change the flow characteristics throughout the system and will induce a unique set of boundary conditions that must be applied to the core bearing these moderating elements.In other words, even if decoupled core simulations are performed, the boundary conditions must be adjusted to capture the specific core loading.
This paper indicates that the multi-channel approach together with the proper handling of the turbomachinery is vital in full core steady-state simulations.The results portray an accurate comparison and optimal design conclusions, such as evidence for the minimization of the pump pressure ratio and the maximization of the thrust and specific impulse within temperature constraints to yield an optimal rocket design.Further results also insist on the utilization of accurate convergences and power distributions as neutronic, and system analyses have been proven to require iterative development and cannot be treated in isolation.
As a result, the ntpSystem is shown to encapsulate an accurate, modular, multi-physics, multi-channel, multi-dimensional steady-state system solution of NTP technology that can be invoked for the advancement of space travel.Future studies will use the steady-state solution to initialize transient simulations of reactor startup, operation, and shut-down scenarios in order to provide a comprehensive NTP simulation.

Figure 2 .
Figure 2. Flowchart of the solution sequence.

Figure 2 .
Figure 2. Flowchart of the solution sequence.
3.1), nozzle, reflector, piping, and moderator elements.The operational conditions at SP7a/b and SP8a/b are known, but the mass flow rate split between the turbine ( .m t ) and the bypass .m bypass valve is not known at this stage.(b)

Figure 3 .
Figure 3. Correlation of cavitation with pump suction specific speed.

Figure 3 .
Figure 3. Correlation of cavitation with pump suction specific speed.
17,  x FOR PEER REVIEW 9 of 27 contours, and the outcome is a set of values that represent the specific speed and diameter as well as the pump efficiency ( ).

Figure 4 .
Figure 4. Specific speed and diameter pump diagram.

Figure 4 .
Figure 4. Specific speed and diameter pump diagram.

Figure 5 .
Figure 5. Specific speed and diameter turbine diagram.

Figure 6 .
Figure 6.Illustration of a modular setup using ntpThermo.

Figure 6 .
Figure 6.Illustration of a modular setup using ntpThermo.

Figure 7 .
Figure 7. Schematic description of the main component modeling in ntpThermo.(a) Reflector and moderator circuits; (b) fuel channels.

Figure 7 .
Figure 7. Schematic description of the main component modeling in ntpThermo.(a) Reflector and moderator circuits; (b) fuel channels.

Figure 8 .
Figure 8. Illustration of the fuel (left) and moderator (right) elements.

Figure 8 .
Figure 8. Illustration of the fuel (left) and moderator (right) elements.

Figure 10 .
Figure 10.Turbopump operational characteristics as a function of the pump pressure ratio.(a) Thermodynamic efficiencies; (b) shaft work, MW; (c) mass flow rate through the turbine.

Figure 10 .
Figure 10.Turbopump operational characteristics as a function of the pump pressure ratio.(a) Thermodynamic efficiencies; (b) shaft work, MW; (c) mass flow rate through the turbine.

Figure 10 .
Figure 10.Turbopump operational characteristics as a function of the pump pressure ratio.(a) Thermodynamic efficiencies; (b) shaft work, MW; (c) mass flow rate through the turbine.

Figure 11 .
Figure 11.Flow characteristics in the system for different pump pressure ratios.(a) Mass flow rates; (b) outlet temperature; (c) outlet pressure.

Figure 11 .
Figure 11.Flow characteristics in the system for different pump pressure ratios.(a) Mass flow rates; (b) outlet temperature; (c) outlet pressure.

Figure 13 .
Figure 13.System total power (a) and shaft work (b) for varying specific impulse values.Figure 13.System total power (a) and shaft work (b) for varying specific impulse values.

Figure 13 .
Figure 13.System total power (a) and shaft work (b) for varying specific impulse values.Figure 13.System total power (a) and shaft work (b) for varying specific impulse values.

Figure 13 .Figure 14 .
Figure 13.System total power (a) and shaft work (b) for varying specific impulse values.

Figure 15 .
Figure 15.Pressure drops for varying specific impulse values.

Figure 15 .
Figure 15.Pressure drops for varying specific impulse values.Figure 15.Pressure drops for varying specific impulse values.

Figure 15 .Figure 16 .
Figure 15.Pressure drops for varying specific impulse values.Figure 15.Pressure drops for varying specific impulse values.Energies 2024, 17, x FOR PEER REVIEW 20 of 27

Figure 16 .
Figure 16.Axial temperature variation for different specific impulse values.(a) Solid fuel surface temperature; (b) propellant temperature.

Figure 17 .
Figure 17.System power (a) and shaft work (b) for varying thrust values.

Figure 18 .
Figure 18.System power (left) and shaft work (right) for varying thrust values.

Figure 17 .
Figure 17.System power (a) and shaft work (b) for varying thrust values.

Figure 17 .
Figure 17.System power (a) and shaft work (b) for varying thrust values.

Figure 18 .
Figure 18.System power (left) and shaft work (right) for varying thrust values.Figure 18. System power (left) and shaft work (right) for varying thrust values.

Figure 18 . 27 Figure 19 .
Figure 18.System power (left) and shaft work (right) for varying thrust values.Figure 18. System power (left) and shaft work (right) for varying thrust values.

Figure 20 .
Figure 20.Axial temperature variation for different thrust values.(a) Solid fuel surface temperature; (b) propellant temperature.

Figure 19 .
Figure 19.Pressure drops for varying thrust values.

Figure 19 .
Figure 19.Pressure drops for varying thrust values.

Figure 20 .
Figure 20.Axial temperature variation for different thrust values.(a) Solid fuel surface temperature; (b) propellant temperature.

Figure 21 .
Figure 21.Specific impulse, thrust, and minimum pump pressure ratio comparison.

Figure 22 .
Figure 22.Multi-channel axial power distributions (left) and radial power peaking factors (right) with channels labeled in black.

Figure 23 .
Figure 23.Propellant exit (left) and max fuel (center) temperatures in K, and mass flow rates in kg/s (right) with no orifices.

Figure 22 .
Figure 22.Multi-channel axial power distributions (left) and radial power peaking factors (right) with channels labeled in black.

Figure 22 .
Figure 22.Multi-channel axial power distributions (left) and radial power peaking factors (right) with channels labeled in black.

Figure 23 .
Figure 23.Propellant exit (left) and max fuel (center) temperatures in K, and mass flow rates in kg/s (right) with no orifices.

Figure 23 .
Figure 23.Propellant exit (left) and max fuel (center) temperatures in K, and mass flow rates in kg/s (right) with no orifices.Energies 2024, 17, x FOR PEER REVIEW 24 of 27

Figure 25 .
Figure 25.Mass flow rates (left) in kg/s, orifice diameter in cm (center), and max fuel temperature in K (right).

Figure 25 .
Figure 25.Mass flow rates (left) in kg/s, orifice diameter in cm (center), and max fuel temperature in K (right).

Figure 25 .
Figure 25.Mass flow rates (left) in kg/s, orifice diameter in cm (center), and max fuel temperature in K (right).

Table 1 .
Material stress data.

Table 4 .
Comparison of flow characteristics with (w) and without (w/o) a convergence scheme.

Table 5 .
Comparison of flow characteristics for different power share values.

Table 6 .
Flow characteristics through a single turbine with and without orifices.Results are pre-

Table 6 .
Flow characteristics through a single turbine with and without orifices.Results are pre-

Table 6 .
Flow characteristics through a single turbine with and without orifices.Results are presented for single-and multi-T/H channels (non-uniform power).
* Flow from the turbine and bypass valve combine.