Thermal Hydraulics Design and Analysis Methodology for a Solid-Core Nuclear Thermal Rocket Engine Thrust Chamber

1. Introduction

Nuclear thermal propulsion [1] can carry far larger payloads and reduce travel time for astronauts going to Mars than is now possible with chemical propulsion. The most feasible concept, extensively tested during the Rover/NERVA (Nuclear Engine for Rocket Vehicle Application) era, is the solid-core concept [2]. It features a solid-core nuclear reactor consisting of hundreds of heat generating prismatic flow elements. Each flow element contains tens of tubular flow channels through which the working fluid, hydrogen, acquires energy and expands in a high expansion nozzle to generate thrust. Approximately twenty nuclear thermal engines with different sizes and designs were tested during that era. Using hydrogen as propellant, a solid-core design typically delivers specific impulses (ISP) on the order of 850 to 1000 s, about twice that of a chemical rocket such as the Space Shuttle Main Engine (SSME).

With the announcement of the Vision for Space Exploration on January 14, 2004, NASA Marshall Space Flight Center started a two-year solid-core nuclear rocket development effort in 2006. The tasks included, but not limited to, nuclear systems development, design methodology development, and materials development. In 2011, with the retirement of Space Shuttle fleets, and NASA’s shifting focus to further out places such as Mars and asteroids, nuclear thermal propulsion is likely to garner substantial interest again. It is therefore timely to discuss the design methodology development effort from 2006 to 2007, entitled “Multiphysics Thrust Chamber Modeling”, which developed an advanced thermal hydraulics computational methodology and studied a solid-core, nuclear thermal engine designed in the Rover/NREVA era.

One of the impacts made by this thermal hydraulics design methodology is the consideration of chemical reacting flow that addresses the effect of hydrogen decomposition and recombination. The advantage of using hydrogen as a propellant is well known in the chemical rocket due to its low molecular weight. However, molecular hydrogen decomposes to atomic hydrogen during high-temperature heating in the thermal nuclear reactor. Since atomic hydrogen has half the weight of that of molecular hydrogen, it was speculated by some that the total thrust could be doubled if all of the hydrogen is dissociated at very high temperatures, therefore leaning to the high power density reactor design. In actuality, the hydrogen conversion is often not uniform across the solid-core since the reactor temperature depends on the hydrogen flow distribution and the actual power profile generated by the nuclear fuel. In addition, the hydrogen atoms recombine in the nozzle since temperature decreases rapidly during the gas expansion, thereby negating the thrust gain. To the best of our knowledge, however, the effect of hydrogen decomposition in a thermal nuclear thrust chamber on the thrust performance has never been studied.

On the other hand, it is always desirable to decrease the reactor weight while one of the ideas is to reduce the reactor size, which increases the power density. One of the impacts of operating at the combination of high temperature and high power density is a phenomenon known as the mid-section corrosion [3], as reported during the legacy engine tests. It is the cracking of the coating layer deposited on the inner wall of the flow channel, coupled with an excessive mass loss of the material near the mid-section of a flow element. The purpose of the coating layer was to isolate the carbonaceous compound in the flow element matrix from the attack by hydrogen. The causes of mid-section corrosion were speculated as a mismatch in the thermal expansion of flow element and its coating material, high flow element web internal temperature gradients, and change of solid thermal property due to irradiation [3, 4]. Those are all possibilities related to the materials. We, however, wanted to trace the cause from a thermal-hydraulics design view point. That is, with the long, narrow flow channel design that is used to heat the hydrogen, the possibility of flow choking in the channel was never studied. One of the efforts in this task was therefore to investigate the possibility of chocked flow occurring in the long, narrow flow channel, and its implication on heat accumulation in the flow element and mid-section corrosion.

The objective of this study was to bridge the development of current design methods with those developed in the Rover/NERVA era, thereby providing a solid base for future development. This is because during the Rover/NERVA era, there was a wealth of lessons learned from those legacy engine tests. All those lessons learned culminated in the final design of the Small Engine [5], but was never built nor tested since it was designed near the end of that era. The Small Engine was therefore a ‘paper engine’ that bears the best features of lessons learned. By simulating and comparing the computed environments with those of the Small Engine design analysis and available test information from the legacy tests, the lessons learned during Rover/NERVA era may be revitalized and the effect of some important design features may be validated. This Chapter therefore reviews the thermal hydraulics computational methodology developed to help the design of a materials tester [4, 6], and the results reported in the design analyses of the Small Engine [7, 8]. Figure 1 shows the computational model [7] of the Small Engine that shows the surfaces of the hydrogen inlet duct, a pressure vessel that houses the solid-core reactor, and an exhaust nozzle that provides the thrust.

Figure 1.

Computational Model of the Small Engine.

A two-pronged approach was conducted to address the aforementioned efforts: a detailed modeling of a single, powered flow element that addresses possible, additional cause of the mid-section corrosion [8], and a global modeling approach that computes the entire thermal flowfield of the Small Engine thrust chamber [7]. The latter links the thrust performance with the effects of power profiles, chemical reactions, and overall heat transfer efficiency. The global approach solves the entire thrust chamber with a detailed modeling, except the thousands of flow channels in the hundreds of flow elements were lumped together as a porous media, for computational efficiency. The heat transfer between other supporting solid components and the working fluid is solved with the conjugate transfer methodology, which was developed in [4, 6]. Theoretical and neutronics provided power profiles were used in lieu of the coupled neutronics modeling. The computational methodology, the results of the simulations of a single flow element and that of the entire thrust chamber, are presented and discussed herein.

2. Computational heat transfer methodology

2.1. Fluid dynamics

The computational methodology was based on a multi-dimensional, finite-volume, viscous, chemically reacting, unstructured grid, and pressure-based, computational fluid dynamics formulation [9]. Time-varying transport equations of continuity, species continuity, momentum, total enthalpy, turbulent kinetic energy, and turbulent kinetic energy dissipation rate were solved using a time-marching sub-iteration scheme and are written as:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\right)=0$$E1

$$\frac{\partial \rho \text{}{\alpha }_i}{\partial \text{}t}+\frac{\partial }{\partial x_j}\left(\rho u_j\text{}{\alpha }_i\right)=\frac{\partial }{\partial x_j}\left[\left(\rho D+\frac{{\mu }_t}{{\sigma }_{\alpha }}\right)\frac{\partial {\alpha }_i}{\partial x_j}\right]+{\omega }_i$$E2

$$\frac{\partial \rho u_i}{\partial \text{}\text{}t}+\frac{\partial }{\partial x_j}\left(\rho \text{}{u}_j{u}_i\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}$$E3

$$\begin{array}{l}\frac{\partial \rho H}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_{j}H\right)=\frac{\partial p}{\partial t}+\frac{\partial }{\partial x_j}\left[\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\frac{\partial H}{\partial x_j}\right]+\frac{\partial }{\partial x_j}\left[\left(\left(\mu +{\mu }_t\right)-\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\right)\frac{\partial \left(V^2/2\right)}{\partial x_j}\right]\\ +\frac{\partial }{\partial x_j}\left[\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\left(u_k\frac{\partial u_j}{\partial x_k}-\frac{2}{3}{u}_j\frac{\partial u_k}{\partial x_k}\right)\right]+Q_r\end{array}$$E4

$\frac{\partial pk}{\partial t}+\frac{\partial }{\partial x_j}(p{u}_j\left[\left(\mu \frac{{\mu }_t}{{\delta }_k}\right)\frac{\partial k}{\partial x_j}\right]+p\left(\Pi -\epsilon \right)$E5

$\frac{\partial p\epsilon }{\partial t}+\frac{\partial }{\partial x_j}\left(p{u}_{j\epsilon }\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\delta }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+p\frac{\epsilon }{k}\left(C_j\Pi -C_2\epsilon +C_3{\Pi }^2/\epsilon \right)$E6

A predictor and corrector solution algorithm was employed to provide coupling of the governing equations. A second-order central-difference scheme was employed to discretize the diffusion fluxes and source terms. For the convective terms, a second-order upwind total variation diminishing difference scheme was used. To enhance the temporal accuracy, a second-order backward difference scheme was employed to discretize the temporal terms. Details of the numerical algorithm can be found in [8-13].

An extended k-ε turbulence model [14] was used to describe the turbulent flow and turbulent heat transfer. A modified wall function approach was employed to provide wall boundary layer solutions that are less sensitive to the near-wall grid spacing. Consequently, the model has combined the advantages of both the integrated-to-the-wall approach and the conventional law-of-the-wall approach by incorporating a complete velocity profile and a universal temperature profile [10].

2.2. Heat transfer in solids

A solid heat conduction equation was solved with the gas-side heat flux distributions as its boundary conditions. The solid heat conduction equation can be written as:

$\frac{\partial p{C}_p{T}_s}{\partial t}=\frac{\partial }{\partial x_j}\left(K\frac{\partial T_s}{\partial x_i}\right)+Q_p+Q_s$E7

The present conjugate heat transfer model [4] solves the heat conduction equation for the solid blocks separately from the fluid equations. The interface temperature between gas and solid, which is stored at interior boundary points, is calculated using the heat flux continuity condition. For solution stability and consistency, the gas and solid interface boundary temperature is updated using the transient heat conduction equation (7).

2.3. Flow and heat transfer in porous media

A two-temperature porosity model was formulated with separate thermal conductivities for the flow and the solid parts. The heat transfer between the flow and solid was modeled by using the empirical correlation of the heat transfer coefficient for circular pipes as a function of flow Reynolds numbers. Empirical multipliers for both the heat transfer and drag loss were determined by comparing solutions of flow passing through a porous flow element with those of a Small Engine 19-channel flow element using detailed conjugate heat transfer modeling [8]. The only affected fluid governing equations are Navier-Stokes and energy equations and can be rewritten as:

$$\frac{\partial \rho u_i}{\partial \text{}\text{}t}+\frac{\partial }{\partial \text{}{x}_j}\left(\rho \text{}{u}_j{u}_i\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}-\frac{L}{\beta }$$E8

$$\begin{array}{l}\frac{\partial \rho H}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_{j}H\right)=\frac{\partial p}{\partial t}+\frac{\partial }{\partial x_j}\left[\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\frac{\partial H}{\partial x_j}\right]+\frac{\partial }{\partial x_j}\left[\left(\left(\mu +{\mu }_t\right)-\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\right)\frac{\partial \left(V^2/2\right)}{\partial x_j}\right]\\ +\frac{\partial }{\partial x_j}\left[\left(\frac{K}{C_p}+\frac{{\mu }_t}{{\sigma }_H}\right)\left(u_k\frac{\partial u_j}{\partial x_k}-\frac{2}{3}{u}_j\frac{\partial u_k}{\partial x_k}\right)\right]+Q_r+\frac{1}{\beta }\text{}{Q}_s\end{array}$$E9

For the solid heat conduction in porous media,

$$\frac{\partial {\rho }_s{C}_{ps}{T}_s}{\partial t}=\frac{\partial }{\partial x_j}\left(K_s\frac{\partial T_s}{\partial x_i}\right)+\frac{1}{1-\beta }\left(Q_p-Q_s\right)$$E10

For the Small Engine 19-channel flow element heat-exchanger configuration, drag loss for flow in circular pipes can be used as a point of departure. That is,

$$L=\frac{1}{2}\rho f_L{c}_f\left|U\right|u_i/d$$E11

where $$c_f=0.0791R{e}^{-0.25}$$ is the Blasius formula for turbulent pipe flow [15]. Typical Reynolds numbers in a flow channel range from 10,000 to 40,000.

For the heat exchange source term,

$$Q_s=\frac{1}{2}\rho f_q\frac{c_f}{p_r^{2/3}}\left|U\right|C_p\left(T_s-T\right)/d$$E12

For the purpose of this study, the conjugate heat transfer module for solids was benchmarked with the analysis of a cylindrical specimen heated by an impinging hot hydrogen jet [4]. The computed solid temperature profiles agreed well with those of a standard solid heat transfer code SINDA [16]. The methodology for flow through porous media was verified through a particle-bed nuclear flow element [17] and the Space Shuttle Main Engine (SSME) main injector assembly [18]. The numerical and physical models for predicting the thrust and wall heat transfer were benchmarked with an analysis of the SSME thrust chamber flowfield, in which the computed axial-thrust performance, flow features, and wall heat fluxes compared well with those of the test data [12].

3. Small engine

The goals of this study were achieved by computing the thermal hydraulics of a flow element and the entire thrust chamber of an engine designed near the end of the Rover/NERVA era – the Small Engine. The thrust chamber of Small Engine composes of an inlet plenum, the solid-core nuclear reactor or heat exchanger, and an exhaust nozzle, as shown in Fig. 1. There are 564 flow elements and 241 support elements or tie-tubes designed for the thermal nuclear reactor. The flow element is shaped like a hexagonal prism, with a length of about 890 mm and a width of about 19 mm from flat to flat [5, 8]. The prismatic flow element contains 19 tubular coolant channels. Three coolant channel diameters were designed for the Small Engine. Each flow element is held in position by three tie-tubes and the corresponding hot-end support system (not modeled). General geometry and operating conditions were obtained from [5], while certain specific operating conditions and nozzle geometry were calculated and provided by the Systems Engineering group.

3.1. Power profiles

For the purpose of this study, theoretical and neutronics calculation provided power profiles were imposed onto the solid-core domain in lieu of the coupled neutronics calculations, for computational efficiency. Combinations of two axial and three radial power profiles, as shown in Figs. 2 and 3, were used to show the effect of which on the heat transfer and thrust performance. Figure 2 shows a Cosine profile and a clipped Cosine profile for the power distribution in the axial direction of the solid-core reactor, while Fig. 3 shows a Cosine profile, a flattened (Cosine) profile, and a flat profile for the power distributions in the radial direction of the nuclear heat exchanger. Given an example for the thrust chamber computations, three combinations of these power profiles were assumed. The first combination uses the shape of the Cosine curve (shown in Figs. 2 and 3) in both the axial and radial directions. That combined power distribution resembles the thermal flux distribution in bare reactors [19], and is simply named as the Cosine-Cosine power profile. By definition, the combined Cosine-Cosine power profile peaks at the middle of the core and drops to zero at the core boundary due to escaping neutrons. The second combination was prescribed by a neutronics calculation with varied Uranium loading. It features the clipped Cosine profile (shown in Fig. 2) in the axial direction and the flattened (Cosine) profile (shown in Fig. 3) in the radial direction, and is dubbed as the clipped Cosine-flattened power profile. The varied fuel loading flattens the (Cosine) power profile in the radial direction, but the power rises near the boundary to show the effect of the reflector, as shown in Fig. 3. The idea of flattening the radial profile is such that the flow in the channels is heated more uniformly, thereby improving the heat transfer efficiency. It is envisioned that the clipped Cosine-flattened power profile is probably the closest power profile to that intended for the Small Engine.

It can be seen that if a radially flattened power profile improves the heat transfer efficiency, then a theoretically flat, radial power profile should reach even higher heat transfer efficiency. We therefore proposed a flat power profile design for the radial direction. A theoretically flat radial power distribution may be achieved with a combination of varied fuel loading and working fluid flow distribution. The third combination therefore employs the clipped Cosine curve (shown in Fig. 2) in the axial direction and a flat curve (shown in Fig. 3) in the radial direction, and is called the clipped Cosine-flat power profile.

Figure 2.

Power profiles used in the axial direction

Figure 3.

Power profiles employed in the radial direction.

3.2. Thermal properties and kinetics

High-temperature real-gas thermodynamic properties were obtained from [20]. These properties were generated for temperatures up to 20,000 K. The peak gas temperature computed in this study did not exceed 10,000 K, hence is well within the applicable range. A 2-species, 2-reaction chemical kinetics mechanism was used to describe the finite-rate hydrogen dissociation and recombination reactions, as shown in Table 1. The first hydrogen recombination reaction is abridged from a large set of kinetics mechanism developed for kerosene combustion [21], while an irreversible, second reaction [22] is added to describe the hydrogen decomposition. The kinetics of the first hydrogen recombination step have been benchmarked through many kinetic mechanism studies, as described in Ref. [21], while the kinetic rates of the second hydrogen recombination reaction were measured [22]. Note the first reaction is a reversible reaction.

Reaction^a	A	B	E/R	M	Ref.
M + H + H = H2 + M	5.0E15	0	0	H2, H	21
M + H2 → H + H + M	8.8E14	0	48300	H2	22

Table 1.

Hydrogen reaction kinetics mechanism

^aM is third-body collision partner and forward rate constant K_f = AT^Bexp(-E/RT).

The importance of finite-rate chemistry was demonstrated in [7], where thermal hydraulics analyses were conducted for the Small Engine with and without finite-rate chemistry. One of the results show that, when the Cosine-Cosine power profile was imposed on the solid-core, the maximum solid temperature calculated with the finite-rate chemistry case was 5369 K, while that for the frozen chemistry case was much higher at 9366 K. This is because the hydrogen decomposition is endothermic. The frozen chemistry freezes molecular hydrogen throughout the thrust chamber and does not allow the hydrogen decomposition to occur, hence an artificially high maximum solid temperature was calculated. It can be seen that without finite-rate chemistry involved, the computed thermal environment and thrust performance are not physical [7].

The solid-core flow element material is assumed to be the (U, Zr)C-graphite composite A, which was tested as flow element material in a legacy reactor [23]. Properties of thermal conductivity, density, and heat capacity as a function of temperature were obtained for (U, Zr)C-graphite composite A from Ref. [23]. Those properties of Beryllium [19] were used for the reflector and slat in the thrust chamber computation. Slat acts as a buffer between the solid-core and the reflector. The thermal properties for the coating layer deposited on the inner wall of flow channels were also obtained from [23].

4. Small engine single flow element modeling

From the material properties point of view, the main cause of mid-section corrosion was speculated as a mismatch in the thermal expansion of flow element and its coating material [3, 4]. The solution is therefore to improve the material properties through materials development [4, 6]. In addition to the materials development, however, we feel further study through the thermal hydraulics analysis of the fundamental reactor design is also important. That is, reduction of the reactor size often started with reducing the diameter of the flow channels, which results in higher aspect ratio, or longer flow channels. According to the Rayleigh line theory, flow with continuous heat addition in a long tube could choke. When that happens, any further heat addition can only serve to reduce the mass flow rate in the tube or, in other words, to jump to another Rayleigh line of lower mass flow [24]. This phenomenon could cause unintended mass flow mal-distribution in the solid-core reactor, resulting in uneven and high local thermal load in the flow element matrix, thereby cracking the coating material. The goal is therefore to compute if choking could occur in one of the flow channels in the Small Engine. To achieve this goal, the worst case was pursued. That is, the Cosine-Cosine power profile which generates peak thermal load in the center of the reactor was assumed. As a result, the flow element located at the center of the solid-core reactor was selected for the computation. In addition, among the three diameters considered for the flow channel of the Small Engine [5], the smallest flow channel diameter was selected, for its highest aspect ratio. These choices were made such that the computed hot hydrogen flow inside one of the flow channels had the best chance to choke.

4.1. Computational grid generation

As mentioned above, the flow element is shaped like a hexagonal prism, containing 19 tubular coolant channels for the propellant hydrogen to acquire heat from the neutronic reactions, as sketched in Fig. 4. For computational efficiency, a 60º pie-section of a single flow element was computed by taking advantage of the symmetrical nature of the prism geometry. A hybrid volume grid was generated by extrusion [25, 26] from a 2-D cross-sectional mesh, as shown in Fig. 5. Structured grid cells were used for the inner and outer layers of the flow channel, and for the outer edges of the prismatic flow element, while unstructured grid cells were used for the rest of the internal web. The structured grid cells at the inner layer of the flow channel was to resolve the turbulent boundary layer, while the general strategy of the hybrid mesh was to minimize numerical diffusion and the number of cells.

Figure 4.

Sketch of the flow element geometry.

Figure 5.

Cross-sectional computational grid layout for the flow channels and the solid web.

In this study, a coating layer, which has different thermal properties from those of the flow element, was also modeled for the flow channel. This was done to imitate the design of those flow elements described in the legacy engine test [23], in which the coating layer was added to protect the carbonaceous compound in the flow element from the alleged chemical attack by the heated, high temperature hydrogen. It is noted that, in the present simulation, the thermal properties of the flow element and the coating layer vary with temperatures. Lastly, in the computational model, an entrance region and an exit region were added to the upstream and downstream of the flow element, respectively, to simulate the environments above and below the solid-core. Total number of cells used was 4.5 million. The heat transfer between the fluid flow, coating layer and the solid fuel was simulated with the conjugate heat transfer model.

4.2. Results and discussion

In this study, various power generation levels were simulated, and it was found that flow choking occurred in some coolant channels as the power level reached about 80% of the maximum power level. Theoretically, once the flow choking occurs, further increase of heat addition would lead to the reduction of mass flow rate in the flow channel, amount to shifting from one Rayleigh line to another as described in Ref. [24]. This could cause mass flow maldistribution in the flow-element matrix, resulting in uneven thermal load in the solid-core. That is, the temperature of the internal web houses the flow channels starved with coolant hydrogen may rise unexpectedly, potentially leading to the cracking of the coating material. Unfortunately, further increase of the power led to unrealistic numerical solutions because of the boundary conditions employed (fixed mass flow rate at the inlet and mass conservation at the exit). This set of boundary conditions was used because only one flow element was simulated, and thus, mass flow reduction at higher power level was precluded. To allow for local mass flow reductions, at least 1/12 of a pie shaped cross-section of the solid-core has to be computed such that fixed mass flow boundary condition need not be used. However, that option was out of the scope of this study. Nevertheless, the current setup is still acceptable since our goal is to determine whether flow choking could occur in the flow channel.

Figure 6 shows the computed temperature distribution in the radial direction (across the solid-fluid interface) at the mid-section of the flow channel. The radial temperature profile reveals that substantial thermal gradients occur at the fluid-coating interface and the coating-solid fuel interface. This kind of high thermal gradient occurring at the coating-solid interfaces could be an issue already without flow choking in the flow channel. With flow choking and the coolant flow rate reduced, the thermal gradient at the coating-solid interface could be even higher since heat is not carried away at the design point, and eventually mid-section corrosion could develop.

Figure 7 shows the computed axial temperature and Mach number profiles along centerlines 1 and 2 (see locations of centerlines 1 and 2 in Fig. 5). The axial temperature profiles almost overlap for centerlines 1 and 2, so are the axial Mach number profiles. The peak Mach number exceeds unity near the end of the flow channel. It can also be seen that the temperatures of propellant hydrogen, rise steadily but drop about 200 K to 2300 K near the end of the fuel port. That is because near the end of fuel port, the flows become choked in both channels. The predicted maximum temperature of both flow channels is about 2700 K. The most significant result from Fig. 7 is that the flow choked near the end of the fuel port, indicating the mass flow rates in these two channels could be reduced, resulting in higher flow element web internal temperature and higher thermal gradient at the coating-solid fuel interfaces and eventually, the possible cracking of the coating layer. Under the operating conditions and assumptions made in this study, the possibility of chocked flow occurring in the flow channels is therefore demonstrated. For future study, three-dimensional numerical simulations of multiple flow elements with a fixed upstream total pressure boundary condition, and a downstream nozzle to remove the fixed exit mass flow boundary condition are necessary to include the inter-element effect. Although the chocked flow occurring in the flow channel is a potential design flaw, the impact of which can readily be avoided by shortening the length of the flow channel to, say, 0.85 m (Fig. 7), without lowering the maximum gas temperature.

Figure 6.

Radial temperature profiles at the mid-section of the flow element.

Figure 7.

Axial temperature and Mach number profiles along centerlines 1 and 2.

5. Small engine thrust chamber modeling

The rationale of analyzing the entire thrust chamber was to bridge the Small Engine design in the Rover/NERVA era to the current modern thermal hydraulic computational methodology. Specifically, to relate the assigned power profiles to that of the Small Engine by comparing the computed thrust performance with that of the design. Since there are 564 flow elements and 241 support elements or tie-tubes, and each flow element contains 19 tubular flow channels, it is computationally cost-prohibitive to solve detailed thermal hydraulics for each and every flow channel. The 19 flow channels in each flow element were therefore lumped as one porous flow channel for computational efficiency, and the effect of drag and heat transfer were modeled with flow and conjugate heat transfer through porous media, or Eq. (8) through Eq. (12). The medium coolant channel diameter, out of the three coolant channel diameters designed for the Small Engine [5], was selected to derive the porosity of the flow elements.

5.1. Computational grid generation

Hybrid computational grids were generated using a software package GRIDGEN V15.07 [27]. The layout of the flow elements and tie-tubes of the solid-core is such that the whole thrust chamber is symmetrical about a 30 deg. pie-shaped cross-section, hence only the 30 deg. cross-section was computed assuming symmetry. The strategy of using hybrid mesh was again for computational efficiency. Figure 8 shows a 240 deg. view of the computational grid and geometry layout of the Small Engine thrust chamber. The thrust chamber includes the hydrogen inlet duct, pressure vessel, and a nozzle with an expansion ratio of 100. Only tie-tubes in the pressure vessel are shown for clarity. Figure 9 shows a cross-sectional view of the solid-core computational grid, depicting flow elements, tie-tubes, slat, and reflector.

Figure 8.

Computational grid and geometry layout of the solid-core Small Engine thrust chamber.

Figure 9.

A cross-sectional view of the solid-core. For illustration, circles are imposed onto the flow elements to differentiate them from the tie-tubes.

No-slip boundary condition was applied to all solid walls, while supersonic outflow boundary condition was employed at the nozzle exit. A fixed total pressure and temperature condition was used at the inlet. Inlet hydrogen flow properties were obtained from a system model simulation. Since the minor thermal effects of tie-tubes were included in the system model simulation, an adiabatic wall boundary condition was used for tie-tube walls. The core surrounding components such as the slat and reflector were treated as heat conducting solids to provide accurate boundary conditions for the solid-core boundary. The heat conducted from the core through slat and reflector to the outer thrust chamber walls was dissipated to a far-field temperature assumed to be 310 K. Table 2 shows the run matrix with three combinations of the axial and radial power profiles: the Cosine-Cosine, clipped Cosine-flattened (Cosine), and clipped Cosine-flat profiles. A series of mesh studies using grid sizes of 1,903,278, 1,925,742, 4,681,751, and 7,460,255 points were performed on case 1. It was found that the computed average pressure drops across the solid-core were very similar and the differences among the computed thrust values were less than 2.5%. Based on those findings and the consideration for computational efficiency, the grid size of 1,925,742 points (or 2,408,198 computational cells) was selected for the rest of the computations.

Case	Axial power shape	Radial power shape
1	Cosine	Cosine
2	Clipped Cosine	Flattened
3	Clipped Cosine	Flat

Table 2.

Run matrix.

5.2. Results and discussion

5.2.1. Temperature contours and profiles in the thrust chamber

Figure 10 shows the computed thrust chamber temperature contours with three different power profiles: the Cosine-Cosine, clipped Cosine-flattened, and clipped Cosine-flat profiles. In the plotted temperature contours, those in the solid-core represent the solid temperatures in the flow elements; the portions surrounding the solid-core depict the temperatures in the slat and reflector, while the rest are the gas temperature contours. It can be seen that the solid-core temperature contours on the left (case 1) reflect the effect of Cosine-Cosine power distribution, since the peak temperature is located near the middle of the reactor, except it is shifted downstream of the geometrical center in the axial direction. The hydrogen gas temperature contours in the solid core take the same shape as those of the solid temperature, except lower. The heat transfer delay in the axial direction is a result of the energy balance between the cooling from the incoming cold hydrogen in the flow channels and the heating from the nuclear material in the web of the flow elements. It is apparent that the effect of cold hydrogen won the fight between the two counter-acting phenomena, and pushes the peak flow element temperature downstream that shows up as a delay in heat transfer.

Figure 10.

Solid and gas temperature contours on the symmetry plane. From left to right: case 1, case 2, and case 3.

Figure 11 shows the computed solid and gas axial temperature profiles on the symmetry plane, along the centerline of each flow element. Only the temperature profiles for the 1^st, 3^rd, 5^th, 7^th, and 9^th flow elementfrom the center tie-tube are plotted for clarity. The solid temperatures (T_s shown as dashed lines) lead those of the gas temperatures (T_g, solid lines), showing the effect of the heat transfer lag between the flow element matrix and hydrogen gas. The gas temperatures appear to peak near the core exit, while the solid temperatures appear to peak more upstream for the 1^st, 3^rd, and 5^th flow element, comparing to the peak temperature locations of the hydrogen. The lead decreases at the 7^th flow element, and disappears completely at the 9^th flow element. These axial temperature profiles reflect the effect of the Cosine-Cosine power profile that concentrates the power in the middle of the core, and drops to zero power at the core boundary. The peak solid temperature for the 9^th flow element is the lowest at around 700 K, as expected. The gas temperature at the core-exit is about 4600 K for the 1^st flow element, and again about 700 K for the 9^th flow element. The 3900 K core-exit gas temperature spread between the 1^st and 9^th flow elements indicates a very non-uniform temperature distribution at the core-exit, or poor heat transfer efficiency, apparently caused by the power-centric Cosine-Cosine power distribution. Note the temperature gradient of the solid temperature is very steep for the 1^st flow element occurring at a location between the core entrance and the peak temperature.

Figure 11.

Solid and gas axial temperature profiles for the 1^st, 3^rd, 5^th, 7^th, and 9^th flow elements on the symmetry plane for case 1.

As discussed above, it is apparent that the Cosine-Cosine power distribution causes not only peak solid temperatures that are too high, but also very non-uniform core-exit gas temperatures, indicating inefficient heat transfer. Fortunately, by measures such as varying the fuel loading, a flattened (Cosine) distribution in the radial direction could be achieved, as indicated in Fig. 3. Note that this profile is not completely flat, since the normalized power is about 1.1 at the center, drops to a minimum of about 0.9 at two-thirds of the radius outward, then rise to a maximum power of 1.6 at the boundary. But comparing it to the Cosine-Cosine power profile, it is relatively flat. In fact, by our estimation, this clipped Cosine–flattened power profile is the most representative of that intended for the Small Engine. The middle picture of Fig. 10 shows the computed solid and gas temperatures, for this clipped Cosine-flattened power profile. It can be seen that the solid and gas temperature contours, reflect the effect of the clipped Cosine-flattened power distribution. Comparing to the temperature contours in the left of Fig. 10, those in the middle are more uniform. The peak solid and gas temperatures in the middle of Fig. 10 at about 3149 and 3081 K., respectively, are much lower than those in the left of Fig. 10, or 5369 and 4596 K., respectively.

Figure 12 shows the computed solid and gas axial temperature profiles on the symmetry plane. Comparing to those of Fig. 11, it can be seen that other than the 9^th flow element, the solid temperature generally leads the gas temperature only slightly, since the clipped Cosine-flattened power profile is more uniform than the Cosine-Cosine profile. The peak core-exit gas temperature of the 1^st fuel element is about 3037 K., while that of the 9^th fuel element is about 2514 K. The range of the gas temperatures at the core-exit is hence only about 523 K, much less than that of 3900 K in case 1, indicating the core-exit gas temperature of case 2 is much more uniform than that of case 1. The flattened radial power profile is therefore a vast improvement over the Cosine radial power distribution.

As mentioned before, it is possible to make the flattened radial power profile in Fig. 3 completely flat, by further fine tuning the fuel loading and by shaping the radial hydrogen mass flow rate profile to match that of the shape of the radial power profile; hence, the proposed theoretical flat radial power profile. The hydrogen mass flow rate profile may be shaped with the installation of various sizes of orifices at the entrance of each flow element. The result of the computed solid and gas temperature contours of such a case, or case 3, using the clipped Cosine-flat power profile, is shown in the right of Fig. 10. It can be seen in terms of the solid-core radial temperature distribution and the uniformity of the gas temperature beneath the solid core, this result is the most uniform. Figure 13 shows the solid and gas axial temperature profiles for the 1^st, 3^rd, 5^th, 7^th, and 9^th flow elements in the solid-core for case 3. It can be seen that the temperature spread among the 1^st, 3^rd, 5^th and 7^th flow elements are much smaller than those of earlier cases, except those of the 9^th flow element due to the heat loss to the slat and reflector. From the right picture of Fig. 10 and from Fig. 13, it can be seen this flat radial power profile is a very good power profile in terms of general uniformity of the temperatures for both solid and gas temperatures inside the core and for gas exit temperatures.

Figure 12.

Solid and gas axial temperature profiles for the 1^st, 3^rd, 5^th, 7^th, and 9^th flow elements on the symmetry plane for case 2.

Figure 13.

Solid and gas axial temperatures for the 1^st, 3^rd, 5^th, 7^th, and 9^th flow elements on the symmetric plane for case 3.

5.2.2. Hydrogen atom mass fraction contours in the thrust chamber

Figure 14 shows the computed thrust chamber hydrogen mass fraction contours with three aforementioned power profiles: the Cosine-Cosine, clipped Cosine-flattened, and clipped Cosine-flat profiles. The hydrogen atom mass fractions are important because the amount of which is directly related to the eventual temperature contours. As a result, the hydrogen atom mass fraction contours look qualitatively similar to those of the solid-core temperature contours, with the peak conversion pushed to near the end of the core due to the same reason that pushed the peak solid temperature away from the center of the core. For the Cosine-Cosine power profile case, as shown in the left picture of Fig. 14, coolant hydrogen enters the thrust chamber at about 370 K, heats up to about 4800 K in the solid-core, then cools and expands into the diverging nozzle to generate the thrust. Molecular hydrogen decomposes to atomic hydrogen at about 2400 K; hence, most of the hydrogen atoms are formed while in the flow elements. Once the local temperature starts to cool, i.e. during the expansion in the nozzle, hydrogen atoms recombine to become molecular hydrogen. The peak hydrogen atom mass fraction in case 1 is 0.40, or 40% conversion from hydrogen molecule decomposition.

The middle picture in Fig. 14 shows the computed hydrogen atom mass fraction contours for case 2. It surely reflects the effect of the clipped Cosine-flattened power distribution shown in Figs. 2 and 3. Since the peak solid and gas temperatures of the clipped Cosine-flattened power profile are lower than those of the Cosine-Cosine power profile, as displayed in the previous section, the peak hydrogen atom mass fraction also drops from a high of 0.40 in the Cosine-Cosine power profile case to a lowly 0.02, or a 2% conversion in the clipped Cosine-flattened power profile case. Note that the atomic hydrogen contours in case 2 exhibit a more pronounced striation in the pressure vessel beneath the core than those in case 1. This is caused by the higher power profile near the core boundary, resulting in higher local temperatures, thereby higher hydrogen molecule conversion near the core boundary for case 2 than that for case 1.

Figure 14.

Hydrogen atom mass fraction contours on the symmetry plane. From left to right: case 1, case 2, and case 3.

Shown in the right figure of Fig. 14 is the computed hydrogen atom mass fraction contours for the clipped Cosine-flat power profile case, or case 3. It can be seen that the hydrogen atom mass fraction contours are the most uniform throughout the thrust chamber, among the three cases. The effect of the flat radial power profile is apparently the driver. As a result, the hydrogen atom conversion drops from 2% in case 2 to 1.5% in case 3.

5.2.3. Summary

The computed heat transfer and thrust performance design parameters for all three cases are compared with those available of the Small Engine and summarized in Table 3. The computed core pressure drops are shown in the second column. The values for all three cases are close to the design number, although that of the Cosine-Cosine power distribution case is slightly lower. The computed peak solid temperatures are shown in the third column, in which the highest peak solid temperature belongs to case 1, followed by case 2, with the lowest peak solid temperature comes from case 3, as expected.

The fourth column shows the average core-exit gas temperatures. These were obtained by averaging the temperatures at the core-exit for all nine elements on the symmetry plane. It can be seen that the average core-exit gas temperature for the Cosine-Cosine power distribution at 3334 K, is much higher than that of the Small Engine design temperature at 2750 K. On the other hand, the average core-exit gas temperatures for the clipped Cosine-flattened and clipped Cosine-flat power distribution cases are both quite close to that of the Small Engine at 2785 and 2782 K, respectively.

Case	ΔPcore,atm	Ts,max, K	Tg,core exit, K	ΔTg,core exit, K	ISP
1	8.9	5369	3334	3900	811
2	9.1	3149	2785	523	868
3	9.1	3066	2782	842	900
Small Engine	9.0	-	2750	-	860~875

Table 3.

A summary of the heat transfer and performance parameters.

The computed temperature spread, or the difference of the core-exit gas temperature between the first and the 9^th flow elements, is shown in the fifth column. Lower temperature spread represents more uniform temperature distribution amongst the flow elements, or better heat transfer efficiency. It can be seen that the clipped Cosine-flattened power profile produces best heat transfer efficiency, the clipped Cosine-flat power profile ranks the second, while the Cosine-Cosine power profile has the worst heat transfer efficiency. Note that if without the outlier of the temperature profile at the 9^th fuel element for case 3, its temperature spread from the rest is more uniform than that of case 2, as indicated in Figs. 12 and 13.

From the last column comes with the comparison of the specific impulses. The best thrust performance belongs to case 3 that is much higher than the design value, followed by case2 which is within the design range, while case 1 has the lowest ISP. Yet the maximum atomic hydrogen conversions are 40%, 2% and 1.5%, for cases 1, 2, and 3, respectively. These results demonstrate that the notion of higher temperature, more atomic hydrogen, therefore higher thrust is not valid. Rather, it is the most uniform radial power profile that produces the highest thrust.

In summary, the computed core pressure drop, core-exit gas temperature, and specific impulse, from the clipped Cosine-flattened power distribution agree well with those of the Small Engine. The result indicates the fluid, thermal, and hydrogen environments computed with the present computational thermal hydraulics methodology closely emulate what was intended for the original Small Engine design. In addition, examining the results of our proposed clipped Cosine-flat power distribution, not only the computed core-exit gas temperature and core pressure drop agree equally well with those of the Small Engine, the computed specific impulse is 25 to 40 seconds higher than those of the Small Engine, demonstrating the present computational methodology can be used to guide the future design.

6. Conclusion

Thermal hydraulics computational design analyses were performed to investigate the mid-section corrosion issue occurred during the legacy engine tests, and to predict the heat transfer efficiency and thrust performance for a virtual solid-core, nuclear thermal engine thrust chamber – the Small Engine. Multiphysics invoked include the turbulent flow and heat transfer, finite-rate chemistry, power generation, and conjugate heat transfer for solids and porous media. The results of a detailed single flow element modeling show that under the assumption of the worst design condition, the hot hydrogen flow could choke near the end of the narrow, long flow channels, which could lead to hydrogen flow reduction and local hot spots in the solid-core; hence, originating the mid-section corrosion. In addition, a unified Small Engine thrust chamber thermal flow field constituting those of the inlet plenum, the pressure vessel, and the nozzle, was analyzed with three power profiles. It was found that the computed core pressure drop, core-exit gas temperature, and specific impulse for the clipped Cosine-flattened power profile closely agreed with those of the Small Engine design values. It was also found that with our proposed clipped Cosine-flat power profile, not only the computed core-ext gas temperature and core pressure drop closely agreed with those of the Small Engine, the corresponding specific impulse was much higher than those of the original Small Engine numbers. Design lesson learned from this effort is that high hydrogen molecule conversion to hydrogen atoms neither improves the heat transfer efficiency, nor increases the thrust performance. Rather, it is the more uniform radial power profile that produces lower peak solid temperature, higher heat transfer efficiency, and higher thrust performance.

Nomenclature

A=pre-exponential rate constant

B=temperature power dependence

C₁,C₂,C₃,C_μ=turbulence modeling constants, 1.15, 1.9, 0.25, and 0.09

C_p=heat capacity

D=diffusivity

d=flow channel diameter

E=activation energy

f_L, f_q=empirical multipliers

H=total enthalpy

K=thermal conductivity

k=turbulent kinetic energy

L=drag loss due to porous media

P=pressure

p_r=Prandtl number

Q=volumetric heat source

R=radial distance

R=gas constant

Re=Reynolds number

T=temperature

t=time, s

U=flow speed

u_i=mean velocities in three directions

x=Cartesian coordinates

Z=axial distance

α=species mass fraction

β=porosity or void of fraction

ε=turbulent kinetic energy dissipation rate

μ=viscosity

μ_t=turbulent eddy viscosity (=ρC_μk²/ε)

Π=turbulent kinetic energy production

ρ=density

σ=turbulence modeling constants

τ=shear stress

ω=chemical species production rate

Subscripts and superscripts

cl=centerline

p=nuclear power source

r=radiation

s=surface, solid, or porous media

t=turbulent flow

Acknowledgments

This study was partially supported by a Nuclear Systems Office task entitled “Multiphysics Thrust Chamber Modeling” with funding from the Prometheus Power and Propulsion Office, NASA Headquarter. Ron Porter was the Marshall Space Flight Center Nuclear Systems Office Manager. Wayne Bordelon was the Nuclear Thermal Propulsion manager. Michael Houts was the Nuclear Research Manager. Steve Simpson and Karl Nelson provided the clipped Cosine-flattened power profile, nozzle geometry and system modeling results. Bill Emrich suggested the Cosine-Cosine power profile. Solid material thermal properties provided by Panda Binayak, Robert Hickman, and Bill Emrich are also acknowledged.