Analysis of the Flow Distribution in a Particle Bed Reactor for Nuclear Thermal Propulsion

Abstract

Nuclear thermal propulsion (NTP) is regarded as the preferred option for the upcoming crewed interstellar exploration due to its excellent performance compared to the current most advanced chemical propulsion systems. Over the past several decades, many novel concepts have been proposed, among which the particle bed reactor (PBR) is the most efficient, compact, and lightweight method. Its unique features, such as the extremely high power density and the radial flow path of coolant in the fuel region, introduce many challengeable issues to the thermal hydraulic design of PBR, with the flow distribution being representative. In this work, the flow distribution process within the core is analyzed based on the understanding of the axial pressure profile in a dummy PBR. A “flow shift” phenomenon leading to the hot spot in the core is introduced first, and three methods, i.e., decreasing the pressure drop within the hot gas channel, increasing the flow resistance on the cold frit or hot frit, and changing the flow pattern from “Z” to “U”, are proposed to reduce the “flow shift” and the consequent temperature mal-distribution. The pros and cons of using cold frit or hot frit to distribute the coolant are also discussed. Finally, by using three numerical examples, these analyses are demonstrated. The findings here may provide technical support for PBR design.

1. Introduction

Interplanetary exploration and deep space activities have attracted enormous attention over the past several decades. However, these projects are expensive even for developed countries or mature space technology organizations. Also, long-term space explorations impose many challenges and risks to astronauts, who suffer more cosmic ray exposure and bone loss. Therefore, a new, trustworthy technique is in great demand to reduce the mission duration, risk, and cost for future space exploration [1,2,3]. The nuclear thermal propulsion (NTP) system is a promising technology, which uses nuclear fission energy rather than the traditional chemical reactions to heat liquid hydrogen [4,5], thus its specific impulse shows significant improvement over chemical rockets due to the lower molecular weight of propellant. The specific impulse represents the efficiency of the rocket engines. More specifically, it is the change in momentum per unit mass for fuels, or rather how much more push accumulates as the fuel is used [6]. Increasing the specific impulse means less propellant is used for same mission or a heavier payload is achieved at the same propellant consumption [7]. Furthermore, NTP has a higher thrust (even up to 100 tons) and a longer lifetime than any other current in-space propulsion systems, which ensures a faster transmit procedure and longer traveling distance on space missions [8,9]. Because of its excellent transportation performance, NTP broadens the launch window of space missions (to months instead of days) [10], and is capable of returning to Earth at any time within three months of Earth departure burn, which increases the mission robustness and allows for potential abort scenarios. National Aeronautics and Space Administration (NASA) has anticipated that NTP could be the preferred option in outer space, especially for manned Mars exploration activities programmed in the Mars Design Reference Architecture 5.0 (DRA 5.0) in 2009 [11], as shown in Figure 1. In this sequence, the NTP system will be equipped on the cargo Mars transfer vehicles (MTVs) and crewed MTV. It will also be applied to the return journey of the crew to Earth. According to the architecture, the initial mass in low earth orbit (IMLEO) reduces by 33%, while the traveling time in space decreases by 40% when using NTP instead of a conventional chemical stage [10].

NTP has been studied for several decades. It was first conceived in 1940s and 1950s, and aggressively developed in the 1960s and the early 1970s. The principal research was carried out in the United States of America (USA) and the Union of Soviet Socialist Republics (USSR) [12]. In the USA, the most representative effort is the “Rover” project and the later Nuclear Engine for Rocket Engine Application (NERVA) [13]. During the projects, more than 20 different reactors were built, and numerous key technologies, including fuel fabrication and testing, were realized. The projects demonstrated that NTP could be a feasible and reliable tool to explore space, and the latest NERVA engine, the XE, was qualified for a human mission to Mars. Meanwhile, the USSR claimed to build the most developed prototype nuclear engine, RD-0410, which had a slightly better performance than NERVA. Regrettably, no nuclear thermal propulsion system has been launched into space yet. The NERVA project was terminated in 1973 due to non-technical reasons, but research still proceeded in some institutions. In the 1980s, a more efficient, compact, and lightweight design, namely the particle bed reactor (PBR), was proposed based on the knowledge derived from NERVA, which further enhanced exploration capabilities [14]. However, PBR is much more immature than NERVA in terms of the technology due to many inherent technical challenges.

A typical PBR mainly consists of a reactor core, a shield, a control system, a liquid hydrogen (LH₂) tank, a turbo-pump, and a nozzle, and it can be divided into two types according to the core layout, i.e., individual PBR (Figure 2a) and fuel elements-assembled PBR (Figure 2b). In the individual PBR, fuel pebbles (particles) are packed randomly in an annular region between two orifices called cold frit and hot frit, respectively. In the fuel elements-assembled PBR, the size of the fuel particles is reduced to that of a grain of sand and the thickness of the fuel region decreases as well; together with the moderator base, a hexagonal fuel element is formed, and several similar elements are assembled to form the core. When the propulsion system runs, cryogenic hydrogen is extracted from the liquid tank and pumped up to a higher pressure. After that, the pressurized hydrogen is preheated by cooling the regenerative channel of the nozzle, the reactor reflector, the control drums, and other important components. Subsequently, a portion of the preheated hydrogen expands in the turbo to drive the pump and the exhausted hydrogen flows into the reactor core. In the core, the coolant flows radially inward to the particle bed to cool the fuel particles, and then flows axially toward the thrust chamber through a hot gas channel. Finally, the high-temperature coolant is ejected from a convergent–divergent nozzle to generate thrust [15]. The partial flow path of the hydrogen is shown in Figure 2.

The particles increase the interfacial heat transfer area between the fuel and the coolant, which enhances the heat transfer capability in PBR. Thus, the power density of the reactor core is raised to an extremely high level, even more than 40 MW/L in the fuel elements, which requires an excellent thermal hydraulic design to avoid the hot spot. However, the flow distribution of coolant in the core is a challenging issue. In PBR, especially in the fuel elements-assembled PBR, the coolant is reasonably distributed among the fuel elements first, and then the coolant is distributed along the axis to match the heat-release profile in each fuel element. The first part does not present a problem, as the flow distribution occurs in the prism-type reactor and can be directly achieved by using orifices. The second part is more difficult, and there are only a few pieces of research that have reported the relative information.

Morley et al. [16] investigated the flow distribution in a critical pellet (pebble) bed reactor using an in-house code, NUTHAM-S, which is based on a two-dimensional thermal hydraulic model for calculating pressure, flow, and temperature. According to their findings, the hot frit contributed to control the amount of coolant flowing in any particular region of the reactor core, while the cold frit was inefficient in the distribution of flow through the core. Hence, the porosity profile on the hot frit was optimized to make the axial temperature in the region adjacent to the hot frit as uniform as possible. However, in another thermal hydraulic analysis performed by Benenati et al. [17], the conclusion was drawn that the cold frit porosity was an important factor which significantly affected the flow distribution in the core. Moreover, they also found that the hot gas channel diameter majorly affected the convective heat transfer process in the core. Tuddenham [18] developed a code to investigate the flow distribution and the thermal hydraulic characteristics in a fuel element of a packed bed nuclear reactor, where the cold frit was designated to control the coolant flow. They observed some significant axial flow in the fuel region, which led to a great variation in temperature along the hot gas channel. Ludwig et al. [14] presented a detailed introduction of three fuel elements-assembled PBRs, which consisted of 19, 37, and 61 fuel elements, respectively. In their work, the flow distribution in the fuel elements was mainly realized by the cold frit. It was also addressed that an optimized porosity profile on the cold frit was beneficial in terms of a better thermal hydraulic design, which was demonstrated through the work performed by Ji et al. [19], where the non-dominated sorting genetic algorithm with elitisms approach (NSGA-II) was adopted to optimize the porosity profile on the cold frit; the thermal performance was exactly improved after this optimization.

From the literature review above, it is noted that there was no consensus regarding the use of cold frit or hot frit to control coolant flow. In addition, the effects of other factors on the flow distribution have not been discussed yet. Therefore, in this paper, the typical axial pressure profile of a dummy PBR was analyzed, based on which, the flow distribution was studied and several approaches, including decreasing the pressure drop within the hot gas channel, increasing the pressure drop across the cold or hot frits, and changing the flow pattern, were proposed to improve the flow distribution and reduce temperature mal-distribution in the PBR core. Moreover, the pros and cons of using cold frit or hot frit for flow distribution were also analyzed. Following this, a computational model was developed and some numerical examples were applied to verify the analyses.

2. Flow Distribution in the PBR Core

2.1. Flow Shift Phenomenon

In an effort to avoid the tedious expression, the individual PBR core and the fuel elements in the assembled PBR core are collectively called the PBR core in the following analyses due to their similar geometry. The simplified view of the PBR core, including the inlet plenum, the cold frit, the particle bed, the hot frit, and the hot gas channel, is shown in Figure 3. From the figure, the coolant flows axially into the reactor core through an inlet plenum, and then turns its direction to flow radially across the cold frit, the particle bed, and the hot frit, and finally flows axially again out of the reactor from the hot gas channel. Therefore, reasonable distribution of the coolant along the axial direction in the PBR core is required to ensure efficient heat transfer. For convenience of analysis, a dummy reactor of uniform power, uniform porosity on both cold and hot frits, and uniform porosity in the fuel region was adopted to analyze the flow process; the conclusions could also be applied to a PBR of non-uniform profiles.

As is known to us all, flow is determined by the pressure difference. Hence, a pressure analysis is the base of the flow distribution in PBR. A schematic diagram of a typical axial pressure profile at several important radial locations, including the outer cold frit interface 1-1′, inner cold frit interface 2-2′, outer hot frit interface 3-3′, and inner hot frit interface 4-4′, is depicted in Figure 3. The pressure drop across the cold frit, the particle bed, and the hot frit is represented by ①, ②, and ③, respectively. The ideal pressure profile at the outer hot interface is depicted by 3-3″, which represents an equal pressure drop through the particle bed along the axial direction and a uniform flow distribution within the particle bed. However, the pressure decreases a lot within the hot gas channel (4-4′) due to frictional resistance to high speed flow and acceleration [16], while increases slightly (1-1′) due to the conversion of dynamic pressure to static pressure in the inlet plenum [20]. Thus, the pressure difference on the left side of the PBR is less than that on the right side. The mal-distribution of the pressure leads to more coolant always flowing toward the right side, which is called the “flow shift” phenomenon. Consequently, the real pressure drop across the fuel bed majorly deviates from the ideal, and the hot spot occurs on the left side of the particle bed due to insufficient cooling there, which negatively affects the thermal hydraulic performance of the PBR.

2.2. Methods to Reduce Flow Shift

There are three methods which use geometry tailoring or flow pattern change to reduce the negative effect of flow shift on the thermal hydraulic characteristics. Detailed information regarding these methods is discussed below.

2.2.1. Decreasing the Pressure Drop within the Hot Gas Channel

According to the previous introduction, flow shift is mainly caused by significant pressure variation within the hot gas channel, thus, the most direct approach is to decrease the pressure drop there. This can be achieved by increasing the diameter of the hot gas channel to decelerate the hydrogen flow, which was also addressed in the work by Benenati et al. [17]. This way, the pressure drop shows a smaller difference between the left side (p₄ − p₁) and the right side (p_4′ − p_1′), and the real pressure profile (3-3′) is closer to the ideal one (3-3″) when compared to the original one. The pressure drop in the hot gas channel can be further decreased to $\Delta p_{hot}^{\prime }$ by reducing the length of the reactor core, e.g., from L to L′. In this case, the difference in pressure drop between the left side and the right side reduces further and the pressure profile (3-3′) approaches the ideal pressure profile (3-3‴) much more, which contributes to the uniform flow distribution and even temperature fields in the reactor. These changes are all shown in Figure 4. However, increasing the diameter of the hot gas channel or decreasing the length of the PBR leads to increased size and weight, meaning the system is more expensive and difficult to control. The size of rockets in particular may limit the effects of this approach in some cases.

2.2.2. Increasing the Pressure Drop across the Cold or Hot Frits

As introduced previously, increasing the diameter of the hot gas channel or decreasing the length of the reactor are beneficial in the reduction of flow shift, but this method is usually limited by other physical requirements. Increasing the pressure drop across the cold or hot frits is another indirect approach and is not limited by these requirements significantly. This can be easily implemented by tailoring the pore size or porosity, and is effective in PBR design. Through this approach, the percentage of the pressure drop within the hot gas channel accounting for the total pressure drop reduces, and the relative discrepancy between the pressure drop on the left side and the right side reduces as well. The PBR of increased pressure drop across a hot frit is called a hot frit-dominated design, while the PBR of increased pressure drop across a cold frit is called a cold frit-dominated design. The choice of a hot frit-dominated design or a cold frit-dominated design is a controversial issue, as illustrated in Section 1, and should be analyzed in detail. In this work, the discussion is carried out from the viewpoint of power density levels in PBR core.

For a PBR of low power density, the fuel particles are usually of large size. Here, a cold frit-dominated design and a hot frit-dominated design are both feasible to distribute the coolant more uniformly within the reactor core. Specifically, if the design is cold frit-dominated, the flow is distributed directly before the cold frit, but some flow redistribution occurs in the fuel region due to the large particle size; even axial flow cannot introduce much additional pressure drop in this case. The flow redistribution may cause a slightly negative effect on the thermal performance of the core [18]. If the design is hot frit-dominated, the coolant flow still adjusts itself before the hot frit, which positively affects the thermal performance of the PBR. Figure 5 shows a schematic diagram of the pressure profile along the axial direction at several radial positions for a hot frit-dominated design and a cold frit-dominated design under such conditions. For both designs, the real pressure profile on the inner hot frit interface (3-3′) approaches the ideal (3-3″), which is marked as a dashed line in the figure.

For a PBR of high power density, small fuel particles (~500 μm) are required to increase the interfacial area and ensure efficient heat transfer from fuel to the coolant. Hence, the pore size of the cold and hot frits must be reduced as well. Under this condition, if cold frit-dominated is the selected method, flow redistribution will be restricted in the fuel region, otherwise a large additional pressure drop is introduced into the core. However, if a hot-dominated design is adopted, a severe hot spot occurs, which can be explained by the Ergun equation [21]:

$$\frac{\Delta p}{L}=\frac{150\mu }{d_p^2}\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3}u+\frac{1.75\rho }{d_p}\frac{\left(1-\epsilon \right)}{{\epsilon }^3}u\left|u\right|,$$

(1)

where d_p is the pore size, ε is the porosity on the frits, μ is the viscosity, ρ is the density, and u is the superficial velocity of the fluid. The first term on the right hand in Equation (1) is the viscous effect, while the second term is the inertia effect. If the pore size on the hot frit is reduced to ~100 μm and the hydrogen is raised to 2800 K, the pressure drop due to viscous effect is much more important than the inertia effect, thus, the flow is viscous-dominated. In this case, less coolant flows toward the hot region due to the increased viscosity of the hydrogen gas with increasing temperature, which is a flow instability phenomenon. This flow instability leads to a severe hot spot and damages the PBR core.

2.2.3. Changing Flow Pattern

In addition to geometry tailoring, flow shift can be also reduced by changing the flow pattern. In the previous introduction, the inlet and the outlet were set on different sides, and coolant flowed co-currently in the inlet plenum and the hot gas channel. In this section, the flow pattern involves the inlet and outlet set on the same side. The coolant flows counter-currently in the inlet plenum and the hot gas channel. The former is called the “Z” flow pattern, while the latter is called the “U” flow pattern. Figure 6 shows a schematic diagram of the pressure profile along the axial direction at several radial positions in under “U” flow pattern conditions. In the figure, the difference between the pressure drop on the left side and that on the right side decreases compared to the “Z” flow pattern, which was shown in the Figure 3. Therefore, the mal-distribution of the coolant is reduced, and the peak temperature moves to the right side of the PBR core. However, the effect of this approach on flow shift reduction is limited to the slight pressure variation in the inlet plenum. For a design of a significant pressure increase in the inlet plenum, changing the flow pattern from “Z” to “U” may be more effective.

3. Numerical Demonstrations

3.1. Numerical Approach

To verify the analyses in the previous section, some numerical examples based on the computational fluid dynamics (CFD) are presented in this section. The investigated domain consisted of an inlet plenum, a cold frit, a particle bed, a hot frit, and a hot gas channel, similar to the simplified view of a PBR as shown in Figure 3. The wall connecting the inlet plenum and radial reflector was set as adiabatic. The cold frit, the particle bed, and the hot frit were regarded as porous media. The flow within the inlet plenum and hot gas channel were both turbulent flows. Generally, the thermal hydraulic characteristics in PBR were investigated by solving the compressible steady Navier–Stokes equation numerically, and the averaged governing equations for conservation of mass and momentum were written as [22]:

$$\nabla ·\left(\rho \mathit{U}\right)=0,$$

(2)

$$\nabla ·\left(\rho \mathit{U}\mathit{U}\right)=-\nabla p+\nabla ·\left({\mu }_{eff}\nabla \mathit{U}\right)+\mathit{S},$$

(3)

where S is the source term of momentum equation and should be considered carefully for the porous media flow, U is the mean superficial velocity of the fluid, and μ_eff is effective viscosity, which is determined by the realizable k-ε model for the flow in the inlet plenum and the hot gas channel, while remaining equal to the molecular viscosity due to turbulence suppression in the particle bed and the cold and hot frits, according to the porous media model. In addition, the non-equilibrium thermal model was adopted to evaluate the heat transfer in the particle bed, and the energy equations for the fluid and solid matrixes were expressed as [23]:

$$\nabla ·\left(\rho \mathit{U}{H}_f\right)=\nabla ·\left(\alpha {\lambda }_f\nabla T_f\right)+h_{sf}{A}_{sf}\left(T_s-T_f\right),$$

(4)

$$0=\nabla ·\left(\left(1-\alpha \right){\lambda }_s\nabla T_s\right)+Q+h_{sf}{A}_{sf}\left(T_f-T_s\right),$$

(5)

where H_f is the enthalpy of the fluid, α is the porosity of the particle bed equaling to 0.38, λ_f is the thermal conductivity of the fluid, λ_s is the thermal conductivity of the solid, T_f is the temperature of fluid, T_s is the temperature of solid, Q is the power density in the fuel region, h_sf is the heat transfer coefficient between coolant and solid, and A_sf is the interfacial area density. A_sf was determined by:

$$A_{sf}=\frac{6\left(1-\alpha \right)}{d_p},$$

(6)

where d_p is the diameter of fuel particle. According to the Safety Standard of the Nuclear Safety Standard Commission (KTA) 3102 “Reactor Core Design for High-Temperature Gas-Cooled Reactors”, the source term S for the momentum equation was calculated as [24]:

$$S=-f·\frac{1-\alpha }{{\alpha }^3}·\frac{1}{d_p}·\frac{1}{2\rho }·{\left(\frac{\dot{m}}{A}\right)}^2,$$

(7)

where $\dot{m}$ is the mass flow rate, A is flow area, and f is the coefficient evaluating pressure loss in the particle bed,

$$f=\frac{320}{\left(\frac{Re}{1-\alpha }\right)}+\frac{6}{{\left(\frac{Re}{1-\alpha }\right)}^{0.1}},$$

(8)

and Re/(1 − α) is the equivalent Reynolds number. Similarly, the heat transfer coefficient was determined by the following relationships [25]:

$$h_{sf}=\frac{Nu·{\lambda }_f}{d_p},$$

(9)

$$Nu=1.27\frac{P{r}^{1/3}}{{\alpha }^{1.18}}R{e}^{0.36}+0.033\frac{P{r}^{1/2}}{{\alpha }^{1.07}}R{e}^{0.86}.$$

(10)

For the porous media flow across the cold and hot frit, the pressure loss was determined by the Ergun model instead, as Equation (1) shows.

To describe the variation of the thermo-physical and transport properties of the coolant, the real gas Peng–Robinson equation was applied to relate the pressure, the specific volume, and the temperature of hydrogen, through which the density could be directly determined [26]. By combining the thermodynamic differential relationships and properties of ideal hydrogen, the specific heat capacity of real gas hydrogen could also be calculated [27]. The molecular viscosity and the thermal conductivity could be both determined by semi-empirical correlations, which included the dilute gas contribution and the excess contribution, due to the effects of two-body interactions and many-body collisions of molecules [28,29].

All numerical simulations were performed using ANSYS FLUENT v19.2. Due to the simplified 2D geometry, the grid systems of quadrilateral cells were constructed easily using the computational domain. These cells were all orthogonal without skewness. In the simulation, the SIMPLEC algorithm was employed for pressure–velocity coupling. The second order upwind scheme was used to discretize the governing equations of continuity, momentum, and energy, while the first order upwind scheme was used for discretization of governing equations with regard to turbulence. Furthermore, the calculation was considered to be converged only when the two criteria were met simultaneously, i.e., (1) when the residual for solving energy equation was less than 10⁻⁸, while the residuals for other solutions were no greater than 10⁻⁵, and (2) when the mass-averaged temperature at the outlet, the area-averaged pressure at the inlet, and the maximum temperature within the computational domain were all nearly unchanged with iterations.

3.2. Verification and Validation of the Numerical Approach

Before the demonstrations, the numerical approach was verified and validated. According to the previous introduction, the conceptual model mainly consisted of three parts, i.e., the Reynolds Averaged Navier-Stokes (RANS) turbulence model, the hydrogen model, and the porous media model. The hydrogen model and the porous media model were embedded into the Computational Fluid Dynamics (CFD) solver through user-defined functions (UDFs). The partial validation of the numerical approach is shown in Figure 7.

Firstly, the hydrogen model was directly validated by the database from National Institute of Standards and Technology (NIST), and it was found that there was only some deviation near the pseudo-critical temperature [30]. In the thermal hydraulic design of PBR, the minimum temperature was not lower than 150 K, hence the hydrogen model was considered feasible. Detailed information regarding this can be found in [15,30].

Subsequently, the RANS model, accurately the k-ε model, was validated by the experimental data obtained by Hendricks et al. [31,32]. These experiments investigated the convection heat transfer of hydrogen flow in a straight tube, where various high heat fluxes were imposed on the tube wall and the properties of hydrogen varied violently due to strong heating. In these simulations, the NIST data of hydrogen were used to exclude possible discrepancies resulting from the user-defined hydrogen model. Two typical experimental conditions, i.e., the trans-critical flow Run 64-706 and the fully supercritical flow Run 103-315 were selected for validation. It is known that the most difficult thing in modeling supercritical flow is the simulation of flow near a pseudo-critical temperature through a numerical approach. In both Run 64-706 and Run 103-315, the reduced inlet temperature _Tr,in, i.e., T_in/T_c, was less than 2, which means that the hydrogen deviated the ideal gas significantly. Hence, it was aggressively inferred that the numerical model was sufficient to reproduce real flow if the derived results agreed well with the measured data in experiments Run 64-706 and Run 103-315. Figure 8 shows the comparison of the computed wall temperature and experimental data, and it can be seen that the calculated results agreed well with the experimental data regarding both tendency and peak temperature. As for the absolute results, the maximum errors were −27.12% and 21.03% for Run 64-706 and Run 103-315, respectively, corresponding to the first measured point. In the other points, the calculated results agreed to within ±10%, which was fairly good regarding the modeling of the convection heat transfer to supercritical flow. Detailed information regarding this is included in Appendix A.

Regarding the porous media model, no experimental or direct numerical simulation investigations have been performed to study hydrogen flow through the heated particle bed. The most well-accepted models are semi-empirical relations in the KTA Safety Standard [24,25], which are widely used in the design and evaluation of high-temperature gas-cooled reactors. In this work, it was assumed that the KTA standard was still viable, and the major objective was to verify the implementation of these models into the CFD solver. The verification was performed by comparing the results of rated High Temperature gas-cooled Reactor-Pebble Module (HTR-PM) [33] from CFD and THERMIX [34], which was a dedicated code for analyzing the thermal performance of the high-temperature gas-cooled reactor. Comparisons of the predicted results using these two codes are listed in

Table 1. Comparisons of the results from CFD and THERMIX.

Parameters	CFD	THERMIX	Discrepancy
Mass flow rate, kg/s	90.49	96.2 × 94% *	-
Inlet temperature of helium, K	523	523	-
Outlet temperature of helium, K	1055	1059	(−4.0) 0.5%
Pressure drop through core, kPa	80.0	76.7	(+3.3) 4.3%
Maximum temperature of helium, K	1210	1170	(+40) 4.5%
Maximum temperature of fuel surface, K	1215	1178	(+37) 4.1%

* 96.2 kg/s is the total mass flow rate at the inlet of core, in which 94% of coolant flows into the core, while the other 6% is by-pass or leaks.

. Through these strong agreements, the implementation of the porous media model was verified.

3.3. Numerical Examlpes and Discussions

3.3.1. Example A: Effects of the Flow Pattern

A particle bed reactor, which was similar to the concept of the critical reactor proposed at the University of New Mexico’s Institute for Space Nuclear Power Studies (ISNPS), was selected [16]. The reactor was a hydrogen-cooled fast reactor. The fuel particles, which were 6 mm in diameter, were packed randomly in an annular region between the hot frit and the cold frit. The pore sizes of the cold and hot frits were 4 mm. The fuel bed, the cold frit, and the hot frit were all treated as porous media of uniform porosity, and the averaged porosities were 0.38, 0.1, and 0.2 respectively. The main geometrical parameters are listed in

Table 2. Main geometrical parameters of the PBR.

Parameters	Value
Radius of the core R, mm	507
Length of the core L, mm	1200
Radius of the hot gas channel Rh, mm	100
Thickness of the hot frit δh, mm	3
Thickness of the fuel bed δb, mm	281
Thickness of the cold frit δc, mm	2
Thickness of the inlet plenum δin, mm	10
Thickness of radial reflector δr, mm	111

.

The thermal hydraulic characteristics in the reactor core were assumed indifferently along the circumference, thus the computational model was simplified to a two-dimensional axis-symmetric model. The heat deposition was assumed to be uniform in the core and the total thermal power was 1000 MW, thus, the mean power density was 1.938 MW/L. The inlet temperature of hydrogen was 200 K, while the mass flow rate of hydrogen at the inlet boundary was 21 kg/s to ensure the exhaust gas approached 2800 K. The operating pressure for the whole reactor was 8.0 MPa. Before the demonstration, a grid independence check was performed on three grids which consisted of 33,120, 13,2480, and 480,800 cells, respectively. These cells were uniformly spaced along the axial direction and non-uniformly spaced along the radial direction. For the inlet plenum and the hot gas channel in particular, the height of the cell adjacent to the wall or interface needed to be small enough for resolution of the boundary layer. The computed results showed that grid 2, with 132,480 cells, was the best regarding both the accuracy and the computational cost; thus, this grid was selected for the following investigations.

Figure 9 shows the temperature field in the PBR of two different flow patterns, i.e., the “Z” flow pattern and the “U” flow pattern. Form the figure, the maximum temperature within the fuel bed changed from the bottom left corner to the bottom right corner when the flow pattern changed from “Z” to “U”. The maximum temperature also significantly decreased. The temperature field was more uniform as the flow pattern changed.

Table 3. Thermal hydraulic characteristics of the two PBRs using different flow patterns.

Flow Pattern	Δp (kPa)	Tout (K)	Tmax (K)
“Z”	872.4	2772.2	7121.3
“U”	867.4	2775.6	4833.0

lists the total pressure drop in the PBR, the outlet temperature of coolant and the maximum temperature in the fuel bed. The hot spot was found to be reduced without any significant negative effects regarding the pressure drop and the outlet temperature when the flow pattern changed from “Z” to “U”. That is to say, the coolant was distributed more uniformly and the thermal hydraulic performance of the “U” flow pattern was better than that of the “Z” flow pattern, which coincided with the analyses in Section 2.2.3.

3.3.2. Example B: Effects of the Core Height

Using the previous PBR as a baseline, the effects of the core height were investigated by reducing the core height to 100 cm and 80 cm. The remaining geometry parameters and the power density were unchanged. The mass flow rates at the inlet boundary decreased to 17.5 kg/s and 14 kg/s, respectively, to keep the ratio of total power output to the mass flow rate constant. In addition, to characterize the influence of reducing core height on the flow distribution in the fuel bed, parameter β was defined as the ratio of the maximum pressure drop to the minimum pressure drop through the particle bed, i.e.,

$$\beta =\frac{\Delta p_{b,\max }}{\Delta p_{b,min}}.$$

(11)

This parameter was directly related to the mal-distribution of coolant. The greater the value of the parameter was, the more severe the mal-distribution of coolant was.

Table 4. Thermal hydraulic characteristics in three PBRs of various core heights.

Heights	Δp (kPa)	Tout (K)	Tmax (K)	β
L = 80 cm	551.8	2799.1	4161.2	2.3345
L = 100 cm	693.0	2786.2	5119.6	3.5258
L = 120 cm	865.7	2773.9	6940.4	5.4813

lists the thermal hydraulic characteristics in three PBRs of various core heights. As shown in the table, the pressure drop through the whole PBR reduced as the core height decreased. This was mainly due to the significant reduction of the pressure drop within the hot gas channel, as shown in Figure 10. The table also indicates that β was reduced, meaning the coolant was distributed more uniformly, which decreased the maximum temperature in the fuel bed from 6940.4 K to 4161.2 K when L decreased from 120 cm to 80 cm. Therefore, reducing the core height contributed to the improvement of the thermal hydraulic performance in PBR, which was in accordance with the analyses discussed in Section 2.2.1.

3.3.3. Example C: Effects of Cold and Hot Frit Resistance

Based on the previous baseline design, a cold frit-dominated PBR was achieved by increasing the resistance of the cold frit and decreasing the resistance of the hot frit. In this design, the pore size of the cold frit was reduced to 10 μm and the porosity of the cold frit was increased to 0.2. Using this PBR to establish a computational model, the pressure profile, streamline, and temperature field were derived, which are shown in Figure 11a and Figure 12a, respectively. Similarly, through the reduction of the pore size of the hot frit to 400 μm, a hot frit-dominated design was derived as well. By performing the CFD calculation, the pressure profile, streamline, and temperature field were obtained; these are shown in Figure 11b and Figure 12b. The maximum temperatures of both the cold frit-dominated PBR and the hot frit-dominated PBR were reduced compared to the previous baseline design from Figure 12 and Figure 9, meaning that the coolant was distributed more uniformly than before. The results confirm the analyses regarding the PBR of the low power density in Section 2.2.2.

Although these two evolutions both show improvement in the temperature profiles compared to the baseline design in the previous section, the hot frit-dominated design was slightly superior to cold frit-dominated design due to a lower peak temperature and more uniform temperature distribution. Figure 11 displays why this was the case. In the cold frit-dominated design, the coolant was redistributed in the fuel bed, especially in the region near the cold frit, where the flow rate was low. This redistribution did not introduce much of an additional pressure drop to the PBR due to the large fuel particle size, but did lead to a significant temperature variation along the hot gas channel. This result was in accordance with the phenomenon observed by Tuddenham [18] and Morley et al. [16]. In the hot-dominated design, a lower magnitude of flow redistribution was observed and the flow was almost fully radial in the fuel region. Logically, the coolant in the hot frit-dominated system was distributed more uniformly than in the cold frit-dominated system.

In addition to the PBR above, another design of higher power density was also selected to demonstrate the effect of increasing resistance of cold or hot frits. This design included a fuel element in an assembled PBR, which was similar to the design of Ludwig et al. in Brookhaven National Laboratory [13]. The main geometrical parameters are listed in

Table 5. Main geometrical parameters of the PBR.

Parameters	Value
Length of the fuel element L, mm	639.20
Radius of the hot gas channel Rh, mm	21.29
Thickness of the hot frit δh, mm	3.00
Thickness of the fuel bed δb, mm	6.15
Thickness of the cold frit δc, mm	2.00
Thickness of the inlet plenum δin, mm	3.00
Fuel particle diameter dp, mm	0.5
Fuel bed porosity, -	0.38

. Besides these parameters, the power density was raised to 30 MW/L and the mass flow rate was 0.396 kg/s to ensure high outlet temperatures approaching 3000 K.

Firstly, the pore sizes of the cold and hot frits of this fuel element were designated as 0.002 mm and 0.1 mm, respectively, and the porosities of the cold and hot frits were 0.4 and 0.5, respectively. In this way, the fuel element was a cold frit-dominated system. Using the CFD calculation, the temperature field, the pressure profile, and the streamline of this cold frit-dominated fuel element were determined, and these values are shown in Figure 13. A seen in the figure, the temperature was rather uniform, and the maximum temperature was 3163.8 K. The total pressure drop was 1189.2 kPa, which was calculated from the CFD results. According to the streamline, the flow redistribution was much weaker than that of the previous cold frit-dominated PBR due to smaller fuel particle sizes; any axial flow led to a significant increase in pressure drop throughout the whole fuel element. However, when the fuel element was designed as a hot frit-dominated one, the CFD calculation diverged, meaning that flow instability occurred. Under this condition, the maximum temperature exceeded 8000 K (8000 K is the pre-defined temperature limit in the CFD calculation), and the PBR core was damaged. Therefore, only the cold frit-dominated design was feasible to distribute the coolant in the fuel element of high power density and small particles, thereby demonstrating the analyses discussed in Section 2.2.2. and accounting for why the cold frit were adopted to control coolant flow in the work by Ludwig et al. [14] and Ji et al. [19].

4. Conclusions and Further Considerations

The particle bed reactor (PBR) is the most efficient, compact, and lightweight method among all current nuclear thermal propulsion (NTP) concepts, but the thermal hydraulic design is a challenging issue, where the flow distribution of coolant is an important component. Based on the understanding of the axial pressure profile in the PBR core, this paper thoroughly analyzed the flow process of coolant in a dummy design. The analyses revealed that more coolant always flows toward one side in this dummy PBR, which was called a “flow shift” phenomenon. This phenomenon was mainly caused by a significant pressure drop within the hot gas channel and led to a hot spot in the core. The stronger the flow shift was, the higher the maximum temperature in the core was. Three approaches, i.e., decreasing the pressure drop in the hot gas channel directly, increasing the resistance of cold or hot frits to reduce the percentage of pressure drop within the hot gas channel accounting for the total pressure drop, and changing the flow pattern from a “Z” to a “U” pattern were proposed to reduce flow shift. In particular, the pros and cons of using a cold frit-dominated design or a hot frit-dominated design were discussed. To restrict flow redistribution and avoid flow instability, the use cold frit is recommended in PBRs of high power density, while the use hot frit is recommended in PBRs of relatively low density to distribute the coolant. All analyses were demonstrated by numerical simulations.

Overall, this work is beneficial to understanding of the flow distribution in PBRs. However, these qualitative analyses are only part of the puzzle. In future work, more studies emphasizing the quantitative relationship of these geometry parameters, power density, and flow variables should be performed. In addition, almost all of work presented here was based on uniform heat release, hence, further analyses should be carried out based on the non-uniform power profile.