The Influence of the Packing Factor on the Fuel Temperature Hot Spots in a Particle-Bed GCFR

Abstract

In the recent past the so-called GCFR has been again a subject of study by the international scientific community. This type of reactors, although still in a preliminary stage of development, is a very interesting perspective because combines the positive characteristics common to all the fast reactors with those of the reactors cooled by helium. Up to now, almost all the analyses on the GCFR thermodynamic aspects have been performed starting from a “global” point of view: generally the core has been modelled as a porous medium and only the global parameters have been taken into account. The local effects have been included in adhoc corrective peak factors. The analyses carried out in the present research will be devoted to the characterization of the local effects, on a microscopic scale. In order to have reliable “global” nuclear and thermal-fluid-dynamic data, the performed analyses will be based on simulations previously performed using the RELAP5-3D code, assuming as input parameters the ETDR core ones. For each considered case, the variation ranges of the evaluated parameters have been estimated on the basis of the “best” and the “worst” cases. To summarize the obtained results, even in transient conditions, the variations of the considered input parameters are less significant for the local output values if compared to those due to the assumed packing factor. As a consequence, in a more general core calculation, the obtained local temperature (and velocity) values will have to be corrected by a proper factor that would have to take into account the results of this research.

1. Introduction

Thermal-Fluid-Dynamics (TFD) is fundamental for the analysis of the Nuclear Power Plants (NPPs). A nuclear reactor is (virtually) capable to produce all the wanted power provided that it is possible to remove it from the core. This characteristic, already known since the first studies on nuclear energy, is sometimes indicated as the “First Reactor Theorem” [1]. Therefore a complete analysis of the reactor energy balance is always mandatory [2].

Among the six Gen IV reactor projects, there are two with a thermal neutronic flux (High-Temperature Gas Reactor or HTGR, Supercritical Water Reactor or SCWR), three with a fast flux (Gas-Cooled Fast Reactor or GCFR, Sodium Fast Reactor or SFR, Lead Fast Reactor or LFR), and one with epithermal flux (Molten Salt Reactor or MSR).

Generation IV reactor projects have to answer to new challenges that nuclear energy must face in this new century. The selected systems have to meet Gen IV criteria [3], namely∶

(1)Sustainability, (2)economics,(3)proliferation-resistance (4)safety and reliability.

The capability of removing the heat produced in the core (both in nominal and accidental conditions) is a key point of criterion 4, particularly for GCFR. In fact this reactor has a high power density (up to 100MW/m³) similar to an LWR, but in this case the power has to be removed by a noncondensable gas (helium) instead of water.

Up to now, almost all the analyses on the GCFR thermodynamic aspects have been performed starting from a “global” point of view: generally the core has been modelled as a porous medium and only the global parameters have been taken into account. The local effects have been included in adhoc corrective peak factors.

On the contrary, the analyses carried out in the present research will be devoted to the characterization of the local effects, on a microscopic scale (by using the FLUENT CFD code [4]). The first step will be to setup a proper “microcell” that is used to evaluate the effects induced by the variation of some parameters (packing factors, mean temperature level, etc.) on the local parameters (temperature, coolant velocity, etc.). Then some parametric analyses will be performed in order to find which are the most significant one with reference to the temperature distribution and velocity field inside the microcell. The final purpose is to supply further data in order to develop a more realistic and accurate way to derive the already mentioned corrective peak factors.

In order to have reliable “global” nuclear and TFD data on a GCFR core, the performed analyses will be based on simulation previously performed [5] using RELAP5-3D code [6], assuming input parameters from the ETDR core [5] (Please note that the Experimental Test Demonstration Reactor (ETDR) is now known as ALLEGRO). For each considered case, the variation range of the evaluated parameters has been estimated on the basis of the “best” and the “worst” cases.

2. Preliminary Considerations on the GCFR Energy Balance

The gas-cooled nuclear reactors are usually characterized by a (relatively) low power density. On the contrary in the GCFR power density ranges from 25 to 100MW/m³. It is clear that the more demanding heat transfer conditions makes accurate TFD calculations necessary. As an example, the hot spots characterization issue could be relevant already at the nominal operating conditions.

On the basis of a wide bibliographical research, it can be stated that some studies exist on the gas-cooled nuclear reactors TFD (e.g., [7, 8]) and on the particle-bed configuration. However all of them are substantially different from this paper, where the following has been taken into account∶

(i)very small particle dimensions (diameter < 1 mm),(ii)internal heat production,(iii)high power density.

Particularly studies on the following were found [2]∶

(i)particle-bed configurations without internal heat production,(ii)CFD applied to gas-cooled nuclear reactors but almost all of them are referred to block type cores and not to pebble (or particle) bed cores (so, among the other features, the stochastical particle arrangement is missing),(iii)CFD-Neutronic coupled studies on pebble (or particle) bed gas-cooled nuclear reactors; however the bed is not modelled in detail, mainly due to the great number of elements and it is usually simulated by the porous medium approximation: as a consequence, it is not possible to calculate the maximum local temperature and the real temperature gradient inside the fuel elements.

3. TFD Simulations

The considered GCFR core was made of annular fuel elements (FEs) containing beds of coated particles (CPs), stochastically packed. Particularly the considered FE is shown in Figure 1. The CPs are double-sized in order to limit the pressure drops inside the core. As a first simplification, in this study the elementary cells taken into account are made only by the smaller CPs (diameter ∼1 mm). The CPs are TRISO-like but the ratio of kernel to external layer dimensions is quite large with respect to “classical” TRISO (i.e., typical CPs for HTR). The coolant (helium) flows through the CPs.

Figure 1 

Particle-bed GCFR fuel element [9].

So, it is necessary to use a model, specifically designed for these conditions.

As a first step, by search of the available scientific literature, suitable physical models were selected. Particularly the models able to evaluate the helium pressure drops through the core and the heat transfer between CPs and coolant were researched.

Regarding the pressure drops through cylindrical beds (some modifications would be needed for taking into account also radial flow through an annular bed), it has been found that they can be evaluated (at least in the porous medium approximation) by Kugeler and Schulten model [7]. The latter model has been designed in order to describe the pressure drops in a reactor similar (at least from this point of view) to the considered GCFR. Particularly the pressure drops Δ𝑝 could be evaluated by using the following formula: Δ𝑝=𝜓1−𝜀𝜀³𝐻𝑑1(𝑚̇2𝜌𝐴)², (1) where 𝜓=320+6(Re/(1−𝜀))(Re/(1−𝜀))^0.1,𝜀=packingfactor. (2) Alternatively other models could be used, as the formula by Ergun [10]: 𝐹=(Δ𝑃Δ𝑧)=150(1−𝜀)²𝜀³𝜇𝑣𝐷²_𝑝+1.75(1−𝜀)𝜀³𝜌|𝑣|𝑣𝐷_𝑝. (3) Looking at the heat transfer, there have been considered some models.

Among the others are Kugeler’s [7] and Shirai’s [11] formulae: Nu_𝑝=1.27Pr^0.33𝜀^1.18Re_𝑝^0.36+0.033Pr^0.5𝜀^1.07Re_𝑝^0.86,Nu_𝑝=2.0𝜀+0.75𝜀Pr^1/3(Re_𝑝)^1/2. (4)(5)

4. Models

All the models shown in the previous paragraphs have been designed considering the whole core as a continuum (usually modelled by porous medium approximation). Therefore they are useful to evaluate the core’s global behaviour, both in terms of pressure drops and temperature distribution. However those models cannot deal with the local (microscale) phenomena, like the hot spots. While in the thermal HTRs these phenomena are not so relevant, they could assume a key role in the fast GCFR, where the power density is much higher.

The considered problem (1-phase coolant flow in complex geometry with the aim to obtain very accurate data) suggests that CFD codes should be used. So, the problem was simulated using FLUENT code [4].

In order to define preliminary fluid-dynamic conditions (laminar or turbulent regime) of the problem, it was necessary to estimate the Re value for the considered cases. It is important to highlight that this is a different formulation of Re by that required by (1)—(5). Starting from already published data (see as an example [5] or [12]), it was chosen the formula for Re in pebble-bed proposed in [12]: 𝐺Re=^′𝑒𝐷^′_𝑒𝜇, (6) where 𝐺^′_𝑒=𝑚̇^′𝐴^′_flows=𝑚̇^′(𝜋(𝐷^′_𝑒,𝐷/4))_′𝑒=4(𝑉^′−𝑉^′_𝑠)𝐴^′_𝑠,𝑚̇^′=(𝑚̇_𝑡𝐴)(^′core𝐴_core), (7)

where 𝑉^′ is the total core volume, 𝑉^′_𝑠 is the CPs volume, 𝐴^′_𝑠 is the CPs surface, 𝑚̇_𝑡 is the coolant flow, 𝐴^′_core is the CPs cross section, 𝐴_core is the whole core cross section.

As a result, for all the considered cases (adopting the above mentioned formula) Re well above 10000 have been obtained, so they are characterized by a fully turbulent regime.

5. Methods

5.1. Minimum Cell Dimensions

In order to model the turbulence correctly, the three “classical” approaches have been considered:

(i)Direct Numerical Simulation (DNS),(ii)Reynolds Averaged Navier Stokes (RANS) equations [13],(iii)Large Eddy Simulation (LES) [12, 14].

For the calculations a segregated numerical method has been chosen and the turbulence has been treated by an RANS approach, using an RNG 𝑘-𝜀, model. This latter, in comparison with the simple 𝑘-𝜀 model could be extended also to the laminar region near the wall without adopting the so-called wall functions [15]. RNG 𝑘-𝜀 model is based on a well-known mathematical theory as described in [16].

As a first step, in order to define the minimum cell dimensions for a meaningful model, the influence of the number of CPs number in the considered area has been investigated: particularly many different cell sizes (ranging from 4 to 32 CPs) have been considered. Due to the limited geometrical dimensions, as an approximation, a constant and uniform heat production inside the CPs has been imposed equivalent for the whole core to 40MW/m³ [9] (this relatively “low” value was assumed taking into account that the maximum power density for a particle-bed GCFR is limited if compared to a GCFR with a different fuel form, as indicated also in [13, 17]). The previous assumption could be plainly justified if we look at the (considered) CP dimensions (diameter ∼1 mm): those dimensions are so limited that a single CP cannot significantly modify the incident fast neutron flux; as a consequence the reaction rates remain substantially unchanged.

For boundary conditions, the He inlet velocity has been set at 30 m/s (opposite to gravity direction) while the coolant outlet pressure has been set at 69 bar [9]. In addition symmetry (cyclic) boundary conditions along 𝑥 and 𝑦 directions. have been imposed(This hypothesis was based on the assumption to consider the cell well inside the core and not near the physical boundaries.).

As a result we have obtained that no significant differences have been found for a CP numbers inside the cell ranging from 8 to 32. So the minimum basic cell is composed by 8 CPs.

5.2. Packing Factor Influence

As a second parameter to be investigated, the packing factor (PF) has been considered. So, in order to identify the influence of this parameter, two basic cells (with the highest and the lowest packing factors) have been considered. The first one (Simple Cubic, SC with a PF = 0.52) is reported in Figure 2, while the second (Face-Centred Cubic, FCC with a PF = 0.74) in Figure 3. For comparison, the experimental average PF in a typical PBMR core is ∼0.61 [18].

Figure 2 

Basic GCFR core cell (Simple Cubic, SC).

Figure 3 

Basic GCFR core cell (Face-Centred Cubic, FCC).

For this comparison the boundary conditions have been assumed, as anticipated, on the basis of simulations previously performed[5] using RELAP5-3D code[6], assuming the input parameters from the ETDR core[5]. Particularly some relevant points have been selected from the results obtained for the steady state and the transient RELAP calculations, distinguishing between the average and the hot channels.

So, some of the global values supplied by the RELAP simulations (i.e., the pressure level, the average temperature at the inlet of the cell) have been assumed as input parameters for the FLUENT local calculations while the others have been evaluated starting from the same global values. For example, looking at Figure 5, in order to evaluate the mean outlet temperature for a cell located at the point “new2” z-level (𝑧^new2), a linear temperature trend between the points “new2” (𝑇=𝑇^RELAP_new2) and “new3” (𝑇=𝑇^RELAP_new3) has been assumed; then, for the cell located in “new2”, the outlet temperature (i.e., the 𝑇 value at 𝑧=𝑧^cell_new2_out) has been calculated with the following formula:

𝑇^cell_new2_out=𝑇^RELAP_new2+(𝑇^RELAP_new3−𝑇^RELAP_new2)⋅𝑧^cell_new2_out−𝑧^new2𝑧^new3−𝑧^new2. (8)

For each point two different FLUENT simulations (one for SC and the other for FCC cells) have been performed.

6. Calculations

As already anticipated, all the CFD calculations were performed using as boundary conditions the results obtained by previous RELAP5-3D simulations on the ETDR plant [5, 19], both for Steady State and Transient conditions. So, in the following subparagraphs, all the general temperature trends reported in the figures were obtained in those RELAP simulations [5, 19] and were used as basis for the CFD simulations.

6.1. Steady State—Average Channel

The temperature axial trend in the average channel for the steady state condition is reported in Figure 4.

Figure 4 

Temperature axial profile—Steady State—Average Channel.

Figure 5 

Temperature axial profile (Steady State—Average Channel) with the points selected for the CFD simulations identified.

In Figure 5 the same profile is reported with the points selected for the CFD simulations evidenced by a circle and identified by an abbreviation.

The main results obtained for the “new1” point are reported in the following Figures 6, 7, 8 and 9.

Figure 6 

Two views of the wall temperature distribution (in Kelvin degrees) for SC cell “new1” (Steady State—Average Channel).

Figure 7 

Two views of the wall temperature distribution (in Kelvin degrees) for FCC cell “new1” (Steady State—Average Channel).

Figure 8 

Bar chart showing (normalized) number of points versus velocity for SC cell “new1” (Steady State—Average Channel).

Figure 9 

Bar chart showing (normalized) number of points versus velocity for FCC cell “new1” (Steady State—Average Channel).

In order to evaluate (at least qualitatively) the velocity field and the amount of stagnation points (where the velocity magnitude is zero or just above zero), some velocity magnitude bar charts are reported in Figures 8 and 9. Please note that the normalization was automatically performed by the code.

The results obtained for “new2” and “new3” points show that temperature profiles and velocity fields are quite similar from a qualitative point of view.

6.2. Steady State—Hot Channel

The temperature axial trend in the hot channel for the steady-state condition is reported in Figure 10.

Figure 10 

Temperature axial profile—Steady State—Hot Channel.

In Figure 11 the same profile is reported with the points selected for the CFD simulations evidenced by a circle and identified by an abbreviation.

Figure 11 

Temperature axial profile (Steady State—Hot Channel) with the points selected for the CFD simulations identified.

The obtained results for the “hot1,” “hot2,” and “hot3” points show that temperature profiles and velocity fields are quite similar to those shown in Figures 6, 7, 8 and 9 (at least from a qualitative point of view).

6.3. Transient (LOFA)—Average Channel

In order to extend the range of potential cases to be studied by FLUENT code, there has been simulated [5, 19] (by RELAP5-3D) also a transient, namely, a Loss Of Flow Accident (LOFA). This particular transient was chosen for two reasons.

(i)We had all the needed data to perform this calculation [20].(ii)It is a quite serious transient for the considered plant.

The temperature axial trends in the average channel for the transient condition are reported in Figure 12.

Figure 12 

Temperature axial profile—LOFA—Average Channel.

In Figure 13 the same profiles are reported with the points selected for the CFD simulations evidenced by a circle and identified by an abbreviation.

Figure 13 

Temperature axial profile (LOFA—Average Channel) with the points selected for the CFD simulations identified.

Again the results obtained for all the selected points show temperature profiles and velocity fields quite similar, at least from a qualitative point of view. So, in order to avoid to show too many similar figures, it will be shown in the next Figures 14, 15, 16 and 17 some results related only to “trs1” case.

Figure 14 

Two views of the wall temperature distribution (in Kelvin degrees) for SC cell “trs1” (LOFA—Average Channel).

Figure 15 

Two views of the wall temperature distribution (in Kelvin degrees) for FCC cell “trs1” (LOFA—Average Channel).

Figure 16 

Bar chart showing (normalized) number of points versus velocity for SC cell “trs1” (LOFA—Average Channel).

Figure 17 

Bar chart showing (normalized) number of points versus velocity for FCC cell “trs1” (LOFA—Average Channel).

Again, in order to evaluate (at least qualitatively) the velocity field and the amount of stagnation points (where the velocity magnitude is zero or just above zero), some velocity magnitude bar charts are reported in Figures 16 and 17. Please note that the normalization was automatically performed by the code.

6.4. Transient (LOFA)—Hot Channel

The last considered case is the hot channel for the LOFA transient. The related temperature axial trends are reported in Figure 18.

Figure 18 

Temperature axial profile—LOFA—Hot Channel.

In Figure 19 the same profiles are reported with the points selected for the CFD simulations evidenced by a circle and identified by an abbreviation.

Figure 19 

Temperature axial profile (LOFA—Hot Channel) with the points selected for the CFD simulations identified.

Once again the results obtained for all the selected points shows temperature profiles and velocity fields that are quite similar, (at least from a qualitative point of view) to those already shown.

7. Discussion

From the obtained results there appears the great influence of the PF on the local parameters. Even if we look at the variations of other “macroscopic” parameters (average helium temperature in the cell, mean coolant temperature difference between inlet and outlet, etc.) in reasonable conditions (i.e., as determined by RELAP simulations), it appears that the PF remains the key factor. To underline this consideration it is useful to have a look to Tables 1 and 2, where for each considered point, the following has been reported:

Table 1 

Parameters for SC cell.

Table 2 

Parameters for FCC cell.

(i)maximum temperature difference between two points inside the cell (Δ𝑇^max),(ii)absolute value of the mean coolant temperature difference between inlet and outlet (|Δ𝑇_out-in|),(iii)average helium inlet temperature (𝑇_in).

Looking at the tables it could be noted that, even for the wide variations of the considered “global” parameters, Δ𝑇^max ranges for the SC cell between 9.46^∘C and 10.92^∘C while for the FCC cell it ranges between 14.78^∘C and 15.34^∘C.

From the obtained results it can be also noted the presence of some stagnation points, where the velocity magnitude is zero or just above zero. As evidenced in some figures (look as an example at Figures 8 and 9), due to the geometrical configuration, the number of stagnation points is higher for FCC cell than for SC one.

To summarize, even in transient conditions, the variations of the considered parameters (average helium temperature in the cell, mean coolant temperature difference between inlet and outlet, pressure drops, etc.) are less significant for the local parameters when compared with those due to the assumed PF. As a consequence, in a more general core calculation (usually modelled by a porous medium), the obtained temperature (and velocity) values will have to be corrected by a factor that takes into account in the results of this research. In fact it is important to highlight that the results have been obtained for a cell with height ℎ≤2 mm, while a real core will have a total height >2000 mm. Although the corrective peak factors could not be calculated in a simple additive way, the global values (and relative differences) will be obviously much higher than those calculated here.

8. Conclusions

The issue related to the cooling by helium of a high power density GCFR core is complex and, up to now, not extensively studied. Also for this reason this research seems to be quite original and innovative, although the obtained results will have to be further validated by code-to-experiment comparisons, too.

From a deep analysis of the analyzed problem, the influence of the PF on the local temperatures distributions and velocity fields looks quite evident. Particularly looking at the obtained results it could be noted that, even for wide variations of the considered “global” parameters, Δ𝑇^max ranges for SC cell between 9.46^∘C and 10.92^∘C while for FCC cell it ranges between 14.78^∘C and 15.34^∘C. So, even in transient condition, the variations of the considered parameters (average helium temperature in the cell, mean coolant temperature difference between inlet and outlet, pressure drops, etc.) are less significant for the local parameters when compared with those due to the assumed PF.

In a real case the CP distribution will be stochastical. Therefore it will be necessary, during the design stage, to assume (in a conservative way) the higher hot spot values, at least until it will be possible to have a reliable way to shrink the PF variation range. At least at the moment, in a more general core calculation (usually modelled by a porous medium), the obtained temperature (and velocity) values will have to be corrected by a factor that would have to take into account the results of the present research.

Finally, as another possible future development, it is important to highlight that, in order to have a reliable and accurate TFD evaluation of an innovative gas-cooled reactors, it will be mandatory to quantify adequately the neutronic source term. So it will be important to perform coupled iterative neutronic-TFD calculations, where the power generation will be the output for the neutronics and a boundary condition for TFD, while, the temperature distribution will be the output for TFD and a boundary condition for the neutronics. That could be made either by multi-physics coupled codes or by coupling separate neutronic and TFD codes which are run in parallel.

Acknowledgments

The work presented in this paper was partly funded by the European Union Sixth Framework Program, under contracts GCFR, PuMA and RAPHAEL. First of all the authors like to thank Engineer Castelliti of SCK•CEN, Dr. Kuijper and Dr. van Heek both of NRG, Professor Kloosterman of TUD, Dr. von Lensa of FZJ and Dr. Mitchell and Dr. Stainsby both of AMEC-NNC for their support. Finally they want to thank Professor Forasassi of DIMNP for his very appreciated collaboration and Dr. Bufalino of SORIT for his precious suggestions and help.