# Pore Scale Thermal Hydraulics Investigations of Molten Salt Cooled Pebble Bed High Temperature Reactor with BCC and FCC Configurations

Authors: Shixiong Song; Xiangzhou Cai; Yafen Liu; Quan Wei; Wei Guo
Journal: Science and Technology of Nuclear Installations
 (2014)

---

## Abstract

The present paper systematically investigated pore scale thermal hydraulics
characteristics of molten salt cooled high temperature pebble bed reactor. By using computational fluid dynamics (CFD) methods and employing simplified body center cubic (BCC) and face center cubic (FCC) model, pressure
drop and local mean Nusselt number are calculated. The simulation result shows that the high Prandtl number molten salt in packed bed has unique fluid-dynamics and thermodynamic properties. There are divergences between
CFD results and empirical correlations’ predictions of pressure drop and local Nusselt numbers. Local pebble surface temperature distributions in several default conditions are investigated. Thermal removal capacities of molten salt are confirmed in the case of nominal condition; the pebble surface temperature under the condition of local power distortion shows the tolerance of pebble in extreme neutron dose exposure. The numerical experiments of local pebble insufficient cooling indicate that in the molten salt cooled pebble bed reactor, the pebble surface temperature is not very sensitive to loss of partial coolant. The methods and results of this paper would be useful for optimum designs and safety analysis of molten salt cooled pebble bed reactors.


---

## Body

## 1. Introduction

The pebble bed advanced high-temperature reactor (PB-AHTR) is one class of the next generation advanced high temperature reactor nuclear plant (NGNP) which was primarily developed at UC-Berkeley [1–3] in collaboration with Oak Ridge National Laboratory and other national labs. The PB-AHTR is different from conventional gas cooled high temperature pebble bed reactors (HTR) because the primary coolant of the reactor uses molten fluoride salt rather than high pressure helium. With benefits of the high heat capacity, excellent heat transport characteristics, and neutron slowing down characteristic of fluorine-lithium-beryllium (FLiBe), the PB-AHTR could allow the increase operational power density, the improving of thermal electricity efficiency, the enhancing of passive safety, and the reducing of core dimensions [2]. The nuclear system could operate at about atmospheric pressure due to low vapour pressure liquid salt, which could significantly reduce the construction cost. More recently, the Shanghai Institute of Applied Physics (SINAP) launched five-year government funded research and development fluoride-salt-cooled high-temperature reactor (FHR) project and constructed an experiment facility to research the thermal hydraulics characteristic of the molten salt loop [4, 5]. The PB-AHTR type reactor may stand a good chance of becoming one of the next generation reactors with inherent safety [6].Based on the similarity of HTR and PB-AHTR, thermal hydraulics analysis of PB-AHTR core can be analogous to HTR, but systematically the study of thermal-hydraulics characteristics is necessary, such as heat removal capacity of fluid salt and overall pressure drop of random [1, 3] or ordered pebble bed [4, 5]. The molten salt has been recommended as a blanket coolant in the design of fusion reactors [7, 8], the high Prandtl number (Pr) fluid has unique fluid-dynamics and thermodynamic properties compared to high pressure helium and water. As indicated by Takeuchi et al. [9] the low thermal conductivity and high viscosity of FLiBe mean that heat transport from heated wall to flow center is dominated by turbulence heat transfer rather than thermal diffusion. Experimental and numerical heat removal capacity of high Pr number fluid for sphere-packed pipe were studied by Watanabe et al. and Shimizu et al. [8, 10]; the results show that intense turbulence near the wall is efficient to mix very thin thermal boundary layer of high Pr fluid. However, there is also lack of experimental and numerical study of ordered pebble bed cooled by FLiBe; the thermal-hydraulics characteristics such as Nusselt number (Nu) and pressure drop need more careful inspection.According to previous study, the thermal-hydraulics characteristics of ordered configuration pebble bed may be quite different from random pebble bed. The Ergun’s equation [11] and Wakao’s relation [12] over-predict the value of pressure drop and Nu number in the case of ordered pebble bed [13–15]. It is known that the thermal-hydraulics model in the standardized HTR code THERMIX [16, 17] and TINTE [18] is macroscopic porous media model, which calculate local block average temperature of pebble bed and coolant with local Nu number rather than pore scale realistic temperature distribution. The pressure drop and local Nu number recommended by different researchers show significantly discrepancies at high Reynolds number (Re); empirical correlations validation has crucial importance in PB-AHTR thermal-hydraulics studies and design optimization.In the past decades, with the development of computer technology, researchers were able to investigate pore scale flow and heat transport character of ordered and random configuration packed bed with the help of computational fluid dynamics (CFD). Initiated by Dalman et al. [19] and Nakayama et al. [20], numerical computational model to investigate fluid flow and heat transfer of packed beds was established and macroscopic correlations were derived. Logtenberg and Dixon [21] investigated wall heat transfer coefficient, Nu number, and the radial effective thermal conductivity ratio in 3D geometry of pseudorandomly stacked spheres in a cylindrical channel. Calis et al. [22] and Romkes et al. [23] investigated the flow and heat transfer characteristics in composite structured packed beds of spheres with CFD and experimental methods, and the CFD results show good agreement with experimental data and they also demonstrate that the effects of packing forms on macroscopic flow and heat transfer characteristics were remarkable.In the field of thermal-hydraulics studies of pebble bed reactor, in order to optimize design and analysis safe operation, microscopic level coolant flow and heat transfer were investigated by several researchers. The CFD method has been performed to investigate pore scale flow and heat transport. The large Eddy simulation (LES) and Reynolds stress model (RSM) were firstly introduced to analysis flow through pebble bed core by Yesilyurt and Hassan [24]. By comparing experimental data, LES, and a large number of Reynolds averaged Navier-Stokes (RANS) based models to investigate flow and temperature patterns, Lee et al. concluded that LES was better in predicting turbulence-induced heat transfer [25]. Hassan investigated velocity, pressure, and pebble surface temperature distributions [26] by using a Smagorinsky-based LES model and interpebble contact points. More recently, quasi-direct numerical calculation was first introduced to optimize a pebble bed configuration by Shams et al. [27, 28] and a broad range of qualitative and quantitative data throughout the computational domain were generated. Wu et al. [29] investigated the advantages and disadvantages of realistic packed pebbles and macroscopic porous media approach for closely packed pebbles in CFD simulation, averaged pressure gradient, and temperature rise across the pebbles using realistic model are similar to porous media model. Ferng and Lin investigating effects of BCC and FCC arrangements on flow and heat transfer characteristics in HTR, the result reveals that KTA recommend correlations developed from the fully-developed flow condition are appropriate for system codes such as THERMIX and TINTE [30].As stated above, these studies demonstrate CFD method is feasible and reliable for modeling flow and heat transfer phenomena of pore scale packed beds. Pore scale temperature distribution, local averageNu number, and overall pressure drop of ordered and random pebble bed nuclear reactor can be predicted by using of CFD method. In this present paper, pore scale thermal-hydraulics characteristic of PB-AHTR/FHR was studied, pressure drop and local average Nu number under forced convection flow are investigated. Temperature distribution of pebbles under normal operation, local power distortion, and condition of insufficient cooling were derived and calculated.The remainder of this paper is organized as follows. A detailed description of the geometry model, mesh, and CFD model are described in Section2. The pressure drop and forced convection heat transfer coefficient (Nu number) are calculated in Section 3. Temperature distribution of pebbles under normal operation, local power distortion, and local pebbles insufficient cooling are demonstrated in Section 4. Finally, we will give some concluding remarks and future work in Section 5.


## 2. Geometry Model and Computational Method

### 2.1. Geometry Model Description

In the conventional gas cooled high temperature pebble bed reactor, the reactor core is deposited by random distribution pebbles. The possible Local configurations of pebbles are simple cubic (SC), body-centered cubic (BCC), face-centered cubic (FCC), and hexagonal close-packed (HCP). The random core is convenient for engineering construction and refueling, but makes great challenge for accurate thermal-hydraulics prediction. The ordered configuration core was employed in [5]. In this present work microscopic thermal-hydraulics characteristics of BCC and FCC configuration of local PB-AHTR core are studied. Following Nijemeisland and Dixon [31], to avoid poor quality mesh cells at the points of contact between pebbles, the “near-miss” model is introduced by left small gaps (1% of pebble diameter) between surfaces. As can be seen in Figure 1 the computational domain consists of inlet block, packed cells, and outlet block. Both of the BCC and FCC models consist of seven unit of cells, the extended inlet and outlet domain guarantee fully developed flow and eliminate reserved flow; model dimensions can be seen in Table 1.Table 1
Geometry parameters of BCC and FCC model.
Parameters
L
1 (mm)
L
2 (mm)
L
3 (mm)
H (mm)
d
p (mm)
Porosity (ϵ)
BCC
80.0
275
145
3.499
30
0.3398
FCC
65.0
330
165
4.285
30
0.2813BCC (b) and FCC (a) packed cells with extend inlet and outlet blocks.
(a)
(b)


### 2.2. Computational Model, Theory, and Mesh Generation

#### 2.2.1. Turbulence Model

In this present work, theRe number ranges from 10 to 5000. There are three flow patterns in the range of the Re number, laminar flow regime, laminar-turbulent transition regime and turbulence flow regime, respectively. In the case of laminar-turbulent transition flow regime, there are no standard guidelines about applicability of CFD turbulence models. Experimental study has found that a transition from laminar to turbulence flow in a spheres packed bed occurs over the range from 110 to 150 for the Re number, fully turbulent was found at Re>300. In this paper, when Re locates at [10,200), the laminar model is used. For the transition flow and turbulence flow regime ([200,5000]), pressure drop and Nu number are calculated by using the renormalization group (RNG) k-ɛ model, which is applicable to the small scale eddies. A brief introduction of the turbulent models will be presented, which mainly focuses on the key feature of the models. Researchers could refer to a more comprehensive discussion [32]. The governing equations of Reynolds averaged Navier-Stokes method in Cartesian coordinates for steady-state incompressible flow are described as follows.For continuity,(1)∂∂xi(ρui)=0.
For momentum,
(2)∂∂xj(ρuiuj)=-∂p∂xi+∂∂xj[μ(∂ui∂xj+∂uj∂xi-23δij∂ul∂xl)]-∂∂xi(-ρuj′uj′¯).
For energy equations, laminar flow is
(3)∂∂xi(ui(ρE+P))=∂∂xi(k∂T∂xi)+Sh,
and turbulence flow is
(4)∂∂xi(ui(ρE+P))=∂∂xi[(k+cpμtPrt)∂T∂xi+ui(τij)]+Sh.
The Reynolds averaged values and the effects of the turbulence are represented by the Reynolds stress, which can be modeled by the Boussinesq hypothesis in order to close the constitutional equations,
(5)-∂∂xi(ρuj′uj′¯)=μt[∂ui∂xj+∂uj∂xi]-23[ρkt+μt∂uk∂xk]δij.kt and ɛ represent the turbulence kinetic energy and rate of turbulence dissipation, which can be obtained from transport equations,
(6)∂∂xi(ρktui)=∂∂xj(αkμeff∂kt∂xj)+μtS2-ρɛ,∂∂xi(ρɛui)=∂∂xj(αɛμeff∂ɛ∂xj)+C1ɛɛktμtS2-C2ɛρɛ2kt-Rɛ,
where the quantities αk=αɛ are the inverse Pr number for k and ɛ, which can be derived from
(7)|α-1.3929α0-1.3929|0.6321|α+2.3929α0+2.3929|0.3679=μmolμeff,
where α0=1. The effective viscosity can be obtained by integrating equation as follows:
(8)d(ρ2kɛμ)=1.72ν^ν^3-1+Cν,
and in the high-Reynolds number limit, the above equation gives
(9)μt=ρCμkt2ɛ,
where ν^=μeff/μ. Equation (8) allow the obtaining of an accurate description of the effective turbulent transport varying with the effective Reynolds number, which is suitable for capturing near wall effect in such complex geometry with a broadly ranged Re number. And compared to the standard k-ɛ model, the additional term Rɛ makes the model more responsive to the effects of rapid strain and streamline curvature. Constants of the RNG k-ɛ model are listed in Table 2.Table 2
Constants of RNGk-ε model.
C
1
ε
C
2
ε
C
ν
C
μ
1.42
1.68
100
0.0845The particle Reynolds number (Red), hydraulics Reynolds number of pebble bed (Reh), and pore scale hydraulics diameter (dh) are defined as follows:
(10)Red=dpρVDμ,Reh=dhρVDϵμ,dh=4ϵ1-ϵ(VpAp)cell,
where VD is the Darcy velocity, Vp is the pebble volume, Ap is the area of pebble surface in one cell, and ϵ is the porosity of packed cells ϵ=Vfluid/Vtotal.


#### 2.2.2.  Wall Treatment

It is of importance to notice that all of the standard turbulence models are only valid for turbulent flow sufficiently faraway from walls where the shear stress terms are small and turbulence viscosity is isotropic. In the near wall region, the viscous boundary layer dominates the flow behavior and there exists large gradient of momentum and temperature, which needs adequate treatment. In the case of pebble bed flow zone, we can expect that both the wall of pebble surface and configurations of pebbles will greatly affect the flow behavior and the wall treatment method should be carefully selected.In the computational domain of pebble bed, the presence of complex distribution solid surface strongly affects turbulent flows. The dissipation of turbulent kinetic energyɛ is much larger near a solid surface than near main flow domain. Large gradient of mean velocity appears in the viscous sublayer, which has a great influence on overall pressure drop. In order to accurately model the viscous sublayer, the enhanced wall treatment (EWT) [32] is employed when the k-ɛ model is activated. The EWT for the ɛ-equation is a two-layer model. In the near-wall viscous sublayer Rey=ρyk/μ the viscosity-affect is completely resolved, in the fully turbulent region (Rey>Rey*, Rey*=200), the k-ɛ model is employed. The EWT is applicable for capturing boundary layer effect of high viscosity fluid, especially in the case of moderate Reh where the viscous contribution can be analogous to inertial effect in pebble bed model.


#### 2.2.3.  Mesh and Boundary Condition

The mesh generation in such complex geometries needs more careful treatment, and in the laminar flow regime mesh convergence can be achieved by decreasing face mesh area and increasing body mesh density. For turbulent flow with EWT the near-wall grid should locate in the viscous sublayer, and the dimensionless wall distancey+ should preferably be equal to 1, at least under 5. We should note that in pebble bed domain the local Re number not only changes along with different calculational case corresponding to Reh numbers but also has large deviations for different local velocities near the pebble surfaces. In order to control y+ in reasonable value and avoid poor mesh quality, face mesh numbers of one pebble are up to 2.5×104 with five prism layers, although grid convergence can be achieved in a lower grid density. As shown in Figure 2 structured hexahedral meshes are used in the computational domain of inlet and outlet block in order to save mesh numbers. The mesh numbers can be seen in Table 3.Table 3
Mesh number for BCC and FCC model.
Pebbles configuration
BCC
FCC
Mesh elements number
4562145
8547926FCC mesh with part of inlet block, outlet block, and one cell.
(a)
(b)
(c)Constant inlet temperature is used for eachReh number; the outflow condition is used at the outlet. All the surfaces of pebbles use the nonslip wall condition. Fixed heat flux (q′′′) thermal conditions at pebble wall are assigned, if Re≤100, q′′′=104[W/m2]; when 100<Re≤500, q′′′=105[W/m2]; otherwise, q′′′=5×105[W/m2]. In the cases of AHTR under normal operation and other conditions, q′′′ is reference to the design of [2]. Symmetric conditions are used for the cell side boundaries.


## 2.1. Geometry Model Description

In the conventional gas cooled high temperature pebble bed reactor, the reactor core is deposited by random distribution pebbles. The possible Local configurations of pebbles are simple cubic (SC), body-centered cubic (BCC), face-centered cubic (FCC), and hexagonal close-packed (HCP). The random core is convenient for engineering construction and refueling, but makes great challenge for accurate thermal-hydraulics prediction. The ordered configuration core was employed in [5]. In this present work microscopic thermal-hydraulics characteristics of BCC and FCC configuration of local PB-AHTR core are studied. Following Nijemeisland and Dixon [31], to avoid poor quality mesh cells at the points of contact between pebbles, the “near-miss” model is introduced by left small gaps (1% of pebble diameter) between surfaces. As can be seen in Figure 1 the computational domain consists of inlet block, packed cells, and outlet block. Both of the BCC and FCC models consist of seven unit of cells, the extended inlet and outlet domain guarantee fully developed flow and eliminate reserved flow; model dimensions can be seen in Table 1.Table 1
Geometry parameters of BCC and FCC model.
Parameters
L
1 (mm)
L
2 (mm)
L
3 (mm)
H (mm)
d
p (mm)
Porosity (ϵ)
BCC
80.0
275
145
3.499
30
0.3398
FCC
65.0
330
165
4.285
30
0.2813BCC (b) and FCC (a) packed cells with extend inlet and outlet blocks.
(a)
(b)


## 2.2. Computational Model, Theory, and Mesh Generation

### 2.2.1. Turbulence Model

In this present work, theRe number ranges from 10 to 5000. There are three flow patterns in the range of the Re number, laminar flow regime, laminar-turbulent transition regime and turbulence flow regime, respectively. In the case of laminar-turbulent transition flow regime, there are no standard guidelines about applicability of CFD turbulence models. Experimental study has found that a transition from laminar to turbulence flow in a spheres packed bed occurs over the range from 110 to 150 for the Re number, fully turbulent was found at Re>300. In this paper, when Re locates at [10,200), the laminar model is used. For the transition flow and turbulence flow regime ([200,5000]), pressure drop and Nu number are calculated by using the renormalization group (RNG) k-ɛ model, which is applicable to the small scale eddies. A brief introduction of the turbulent models will be presented, which mainly focuses on the key feature of the models. Researchers could refer to a more comprehensive discussion [32]. The governing equations of Reynolds averaged Navier-Stokes method in Cartesian coordinates for steady-state incompressible flow are described as follows.For continuity,(1)∂∂xi(ρui)=0.
For momentum,
(2)∂∂xj(ρuiuj)=-∂p∂xi+∂∂xj[μ(∂ui∂xj+∂uj∂xi-23δij∂ul∂xl)]-∂∂xi(-ρuj′uj′¯).
For energy equations, laminar flow is
(3)∂∂xi(ui(ρE+P))=∂∂xi(k∂T∂xi)+Sh,
and turbulence flow is
(4)∂∂xi(ui(ρE+P))=∂∂xi[(k+cpμtPrt)∂T∂xi+ui(τij)]+Sh.
The Reynolds averaged values and the effects of the turbulence are represented by the Reynolds stress, which can be modeled by the Boussinesq hypothesis in order to close the constitutional equations,
(5)-∂∂xi(ρuj′uj′¯)=μt[∂ui∂xj+∂uj∂xi]-23[ρkt+μt∂uk∂xk]δij.kt and ɛ represent the turbulence kinetic energy and rate of turbulence dissipation, which can be obtained from transport equations,
(6)∂∂xi(ρktui)=∂∂xj(αkμeff∂kt∂xj)+μtS2-ρɛ,∂∂xi(ρɛui)=∂∂xj(αɛμeff∂ɛ∂xj)+C1ɛɛktμtS2-C2ɛρɛ2kt-Rɛ,
where the quantities αk=αɛ are the inverse Pr number for k and ɛ, which can be derived from
(7)|α-1.3929α0-1.3929|0.6321|α+2.3929α0+2.3929|0.3679=μmolμeff,
where α0=1. The effective viscosity can be obtained by integrating equation as follows:
(8)d(ρ2kɛμ)=1.72ν^ν^3-1+Cν,
and in the high-Reynolds number limit, the above equation gives
(9)μt=ρCμkt2ɛ,
where ν^=μeff/μ. Equation (8) allow the obtaining of an accurate description of the effective turbulent transport varying with the effective Reynolds number, which is suitable for capturing near wall effect in such complex geometry with a broadly ranged Re number. And compared to the standard k-ɛ model, the additional term Rɛ makes the model more responsive to the effects of rapid strain and streamline curvature. Constants of the RNG k-ɛ model are listed in Table 2.Table 2
Constants of RNGk-ε model.
C
1
ε
C
2
ε
C
ν
C
μ
1.42
1.68
100
0.0845The particle Reynolds number (Red), hydraulics Reynolds number of pebble bed (Reh), and pore scale hydraulics diameter (dh) are defined as follows:
(10)Red=dpρVDμ,Reh=dhρVDϵμ,dh=4ϵ1-ϵ(VpAp)cell,
where VD is the Darcy velocity, Vp is the pebble volume, Ap is the area of pebble surface in one cell, and ϵ is the porosity of packed cells ϵ=Vfluid/Vtotal.


### 2.2.2.  Wall Treatment

It is of importance to notice that all of the standard turbulence models are only valid for turbulent flow sufficiently faraway from walls where the shear stress terms are small and turbulence viscosity is isotropic. In the near wall region, the viscous boundary layer dominates the flow behavior and there exists large gradient of momentum and temperature, which needs adequate treatment. In the case of pebble bed flow zone, we can expect that both the wall of pebble surface and configurations of pebbles will greatly affect the flow behavior and the wall treatment method should be carefully selected.In the computational domain of pebble bed, the presence of complex distribution solid surface strongly affects turbulent flows. The dissipation of turbulent kinetic energyɛ is much larger near a solid surface than near main flow domain. Large gradient of mean velocity appears in the viscous sublayer, which has a great influence on overall pressure drop. In order to accurately model the viscous sublayer, the enhanced wall treatment (EWT) [32] is employed when the k-ɛ model is activated. The EWT for the ɛ-equation is a two-layer model. In the near-wall viscous sublayer Rey=ρyk/μ the viscosity-affect is completely resolved, in the fully turbulent region (Rey>Rey*, Rey*=200), the k-ɛ model is employed. The EWT is applicable for capturing boundary layer effect of high viscosity fluid, especially in the case of moderate Reh where the viscous contribution can be analogous to inertial effect in pebble bed model.


### 2.2.3.  Mesh and Boundary Condition

The mesh generation in such complex geometries needs more careful treatment, and in the laminar flow regime mesh convergence can be achieved by decreasing face mesh area and increasing body mesh density. For turbulent flow with EWT the near-wall grid should locate in the viscous sublayer, and the dimensionless wall distancey+ should preferably be equal to 1, at least under 5. We should note that in pebble bed domain the local Re number not only changes along with different calculational case corresponding to Reh numbers but also has large deviations for different local velocities near the pebble surfaces. In order to control y+ in reasonable value and avoid poor mesh quality, face mesh numbers of one pebble are up to 2.5×104 with five prism layers, although grid convergence can be achieved in a lower grid density. As shown in Figure 2 structured hexahedral meshes are used in the computational domain of inlet and outlet block in order to save mesh numbers. The mesh numbers can be seen in Table 3.Table 3
Mesh number for BCC and FCC model.
Pebbles configuration
BCC
FCC
Mesh elements number
4562145
8547926FCC mesh with part of inlet block, outlet block, and one cell.
(a)
(b)
(c)Constant inlet temperature is used for eachReh number; the outflow condition is used at the outlet. All the surfaces of pebbles use the nonslip wall condition. Fixed heat flux (q′′′) thermal conditions at pebble wall are assigned, if Re≤100, q′′′=104[W/m2]; when 100<Re≤500, q′′′=105[W/m2]; otherwise, q′′′=5×105[W/m2]. In the cases of AHTR under normal operation and other conditions, q′′′ is reference to the design of [2]. Symmetric conditions are used for the cell side boundaries.


## 2.2.1. Turbulence Model

In this present work, theRe number ranges from 10 to 5000. There are three flow patterns in the range of the Re number, laminar flow regime, laminar-turbulent transition regime and turbulence flow regime, respectively. In the case of laminar-turbulent transition flow regime, there are no standard guidelines about applicability of CFD turbulence models. Experimental study has found that a transition from laminar to turbulence flow in a spheres packed bed occurs over the range from 110 to 150 for the Re number, fully turbulent was found at Re>300. In this paper, when Re locates at [10,200), the laminar model is used. For the transition flow and turbulence flow regime ([200,5000]), pressure drop and Nu number are calculated by using the renormalization group (RNG) k-ɛ model, which is applicable to the small scale eddies. A brief introduction of the turbulent models will be presented, which mainly focuses on the key feature of the models. Researchers could refer to a more comprehensive discussion [32]. The governing equations of Reynolds averaged Navier-Stokes method in Cartesian coordinates for steady-state incompressible flow are described as follows.For continuity,(1)∂∂xi(ρui)=0.
For momentum,
(2)∂∂xj(ρuiuj)=-∂p∂xi+∂∂xj[μ(∂ui∂xj+∂uj∂xi-23δij∂ul∂xl)]-∂∂xi(-ρuj′uj′¯).
For energy equations, laminar flow is
(3)∂∂xi(ui(ρE+P))=∂∂xi(k∂T∂xi)+Sh,
and turbulence flow is
(4)∂∂xi(ui(ρE+P))=∂∂xi[(k+cpμtPrt)∂T∂xi+ui(τij)]+Sh.
The Reynolds averaged values and the effects of the turbulence are represented by the Reynolds stress, which can be modeled by the Boussinesq hypothesis in order to close the constitutional equations,
(5)-∂∂xi(ρuj′uj′¯)=μt[∂ui∂xj+∂uj∂xi]-23[ρkt+μt∂uk∂xk]δij.kt and ɛ represent the turbulence kinetic energy and rate of turbulence dissipation, which can be obtained from transport equations,
(6)∂∂xi(ρktui)=∂∂xj(αkμeff∂kt∂xj)+μtS2-ρɛ,∂∂xi(ρɛui)=∂∂xj(αɛμeff∂ɛ∂xj)+C1ɛɛktμtS2-C2ɛρɛ2kt-Rɛ,
where the quantities αk=αɛ are the inverse Pr number for k and ɛ, which can be derived from
(7)|α-1.3929α0-1.3929|0.6321|α+2.3929α0+2.3929|0.3679=μmolμeff,
where α0=1. The effective viscosity can be obtained by integrating equation as follows:
(8)d(ρ2kɛμ)=1.72ν^ν^3-1+Cν,
and in the high-Reynolds number limit, the above equation gives
(9)μt=ρCμkt2ɛ,
where ν^=μeff/μ. Equation (8) allow the obtaining of an accurate description of the effective turbulent transport varying with the effective Reynolds number, which is suitable for capturing near wall effect in such complex geometry with a broadly ranged Re number. And compared to the standard k-ɛ model, the additional term Rɛ makes the model more responsive to the effects of rapid strain and streamline curvature. Constants of the RNG k-ɛ model are listed in Table 2.Table 2
Constants of RNGk-ε model.
C
1
ε
C
2
ε
C
ν
C
μ
1.42
1.68
100
0.0845The particle Reynolds number (Red), hydraulics Reynolds number of pebble bed (Reh), and pore scale hydraulics diameter (dh) are defined as follows:
(10)Red=dpρVDμ,Reh=dhρVDϵμ,dh=4ϵ1-ϵ(VpAp)cell,
where VD is the Darcy velocity, Vp is the pebble volume, Ap is the area of pebble surface in one cell, and ϵ is the porosity of packed cells ϵ=Vfluid/Vtotal.


## 2.2.2.  Wall Treatment

It is of importance to notice that all of the standard turbulence models are only valid for turbulent flow sufficiently faraway from walls where the shear stress terms are small and turbulence viscosity is isotropic. In the near wall region, the viscous boundary layer dominates the flow behavior and there exists large gradient of momentum and temperature, which needs adequate treatment. In the case of pebble bed flow zone, we can expect that both the wall of pebble surface and configurations of pebbles will greatly affect the flow behavior and the wall treatment method should be carefully selected.In the computational domain of pebble bed, the presence of complex distribution solid surface strongly affects turbulent flows. The dissipation of turbulent kinetic energyɛ is much larger near a solid surface than near main flow domain. Large gradient of mean velocity appears in the viscous sublayer, which has a great influence on overall pressure drop. In order to accurately model the viscous sublayer, the enhanced wall treatment (EWT) [32] is employed when the k-ɛ model is activated. The EWT for the ɛ-equation is a two-layer model. In the near-wall viscous sublayer Rey=ρyk/μ the viscosity-affect is completely resolved, in the fully turbulent region (Rey>Rey*, Rey*=200), the k-ɛ model is employed. The EWT is applicable for capturing boundary layer effect of high viscosity fluid, especially in the case of moderate Reh where the viscous contribution can be analogous to inertial effect in pebble bed model.


## 2.2.3.  Mesh and Boundary Condition

The mesh generation in such complex geometries needs more careful treatment, and in the laminar flow regime mesh convergence can be achieved by decreasing face mesh area and increasing body mesh density. For turbulent flow with EWT the near-wall grid should locate in the viscous sublayer, and the dimensionless wall distancey+ should preferably be equal to 1, at least under 5. We should note that in pebble bed domain the local Re number not only changes along with different calculational case corresponding to Reh numbers but also has large deviations for different local velocities near the pebble surfaces. In order to control y+ in reasonable value and avoid poor mesh quality, face mesh numbers of one pebble are up to 2.5×104 with five prism layers, although grid convergence can be achieved in a lower grid density. As shown in Figure 2 structured hexahedral meshes are used in the computational domain of inlet and outlet block in order to save mesh numbers. The mesh numbers can be seen in Table 3.Table 3
Mesh number for BCC and FCC model.
Pebbles configuration
BCC
FCC
Mesh elements number
4562145
8547926FCC mesh with part of inlet block, outlet block, and one cell.
(a)
(b)
(c)Constant inlet temperature is used for eachReh number; the outflow condition is used at the outlet. All the surfaces of pebbles use the nonslip wall condition. Fixed heat flux (q′′′) thermal conditions at pebble wall are assigned, if Re≤100, q′′′=104[W/m2]; when 100<Re≤500, q′′′=105[W/m2]; otherwise, q′′′=5×105[W/m2]. In the cases of AHTR under normal operation and other conditions, q′′′ is reference to the design of [2]. Symmetric conditions are used for the cell side boundaries.


## 3. Pressure Drop and Nu Number Simulations of BCC and FCC Model

### 3.1. Basic Definition of Thermal-Hydraulic Correlations

The most widely used empirical correlations to predict pressure drop and local meanNu number are Ergun’s equation and Wakao’s correlation. Ergun’s equation was first proposed by Ergun [11],
(11)(ΔPΔL)Ergun=f·12(VDϵ)21dh,f=(c1Reh+c2),
where f is the friction factor, VD is Darcy velocity in porous medium, and c1=133 and c2=1.33 are constants in Ergun’s original work. Wakao’s type correlations were first introduced by Wakao et al. [12],
(12)NuWakao=hwf·dpkf=a1+a2Pr1/3·(ϵdpdh)n,
where hwf is the heat transfer coefficient of pebble wall to fluid, dp is particle diameter, a1, a2, and n are model constants with a1=2.0, a2=1.1, and n=0.6. The constant a1=2.0 corresponds to the heat transfer coefficient by solving the heat transfer equations for a single pebble immersed in a stagnant liquid with uniform temperature far from the pebble; the other term of (12) is similar to the correlations for smooth tube. In thermal hydraulics studies of gas cooled pebble bed reactors, correlations recommended by German regulatory commission (KTA correlations) [33, 34] were mostly used, which can be described as follows:
(13)(ΔPΔL)KTA=f·(1-ϵ)ϵ3·ρ2·VD2,f=320Red/(1-ϵ)+6(Red/(1-ϵ))0.1,NuKTA=1.27Pr1/3ϵ1.18Red0.36+0.033Pr1/2ϵ1.07Red0.86.
Pressure drop and Nu number predicted by correlations of Ergun, Wakao, and KTA show significant divergences especially the Reh number up to 103, which suggests careful validations in thermal hydraulics analysis of ordered and random molten salt cooled pebble bed reactors. In this present work, normalized pressure drop and local mean Nu number for each cell of BCC and FCC model can be defined as follows [14]:
(14)ΔpΔz=1N∑cellNicellNj‍1(zout-zin)Ain∫‍∫Ain‍(pout-pin)dA,Nu=dpkf1N∑cellNicellNj‍ρfcp∫‍∫Ain‍u(Tf,out-Tf,in)dAAp(T¯p-T¯f),
where u, Tf,in, Tf,out are flow velocity and mass flow temperature at (in,out) cross section of each cells, respectively. The fluid average temperature T¯f and pebble wall average temperature T¯wall in each cell can be defined as
(15)T¯f=12∫‍∫Ain‍u(Tf,in+Tf,out)dA∫‍∫Ain‍udA,T¯p=1N∑1N‍Twall.As shown in Figure1, the calculation domain consists of seven periodically packed cells in the flow direction. One important question is whether this simplified model can depict overall effects of “integrated bed” such as pressure drop and average Nu number. To illustrate this problem, Figure 3(a) shows the pressure drop of BCC and FCC model in the case of Reh=1000 of each cell. As can be seen in Figure 3(a), in the flow direction after the flow through two cells the pressure gradient is almost constants, this result indicating the flow appears to be periodically formed. This phenomenon was also observed in [22] by using numerical and experimental methods.Pressure drop and Nu number of BCC and FCC model.
(a)
(b)The local mean Nu number of each cell ofBCC and FCC model in the case of Reh=1000 is shown in Figure 3(b). In Figure 3(b), it can be found that the entrance effect of BCC model causing higher heat transfer capability of fluid flow, the cell averaged Nu number has significant change in Cell 1–3. The flow stream lines of BCC and FCC in the case of Reh=1000 are shown in Figure 4 for the sake of intuitive observation. Figure 4 shows that when the entrance effect decays along with the flow passes downstream through two cells, the flow is periodically developed, and thus the local Nu number tending towards stability in Cell 3–6 (Figure 3(b)). At the last cell of BCC model, because of a hollow located at the cell center, the turbulent mainstream can not be suppressed on the near-wall region which causes a decline of Nu number. In the FCC packed cells the entrance effect is observed at the very beginning of Cell 1; when flow passed through Cell 2, the entrance effect decays and the Nu number is almost equal in Cell 3–6.Streamline distributions in BCC and FCC packed cells.
(a)
(b)It is interesting to notice that, experimental results demonstrate that heat transfer from pebbles positioned in the entrance layers of a bed is lower than the average [35]. In [35] Achenbach indicated that the lower Nusselt numbers in first layer of packed bed are due to the fact that the velocity and turbulence levels of the incoming flow are smaller than the superficial velocity, and we can expect that the smaller velocity and turbulence levels should be owed to lower porosity in the entrance region. Unlike the random bed the ordered packed bed has a uniform porosity, flows acceleration, and strong turbulent mixing which induce higher heat transfer coefficient. This phenomenon has also been observed by several researches using experimental and numerical methods in the geometry of sphere and cylinder packed bed [36–38]. In Figure 4, it also can be found that the local flow velocity of the FCC cells is higher than BCC model. This observation supports the results of Figure 3(b); under the same Reh number, the Nu number of FCC cells is higher than BCC cells.In the next section, the normalized pressure drop and local average Nu number are defined by (14) with Ni=4 and Nj=5.


### 3.2.  Pressure Drop and Forced Convection Heat Transfer Investigations

In this section, density, dynamic viscosity, heat conductivity, andPr number are assumed as reference constant values when Reh>100. Thermal-physical properties are calculated by (16) at Tin=928.15K, which is the average operation temperature of PB-AHTR steady state [2]. Thermal properties are given in Table 4. Density and dynamic viscosity of FliBe [39] are as follows:
(16)ρ=(2.28-4.88*10-4*(Tin-273.15))*103[kg/m3],μ=0.116*103*exp⁡(3755Tin)[kg/ms].Table 4
Thermal-physical properties of FliBe atTin=928.15 K.
Densityρ [kg/m2]
1960.0
Dynamic viscosityμ  [kg/m s]
0.00663
Thermal conductivityk [W/(m K)]
1.1
Specific heatCv [J//(kg K)]
2415.0The pressure drop (Δp/Δz) of BCC and FCC model proportioned to Reh number is shown in Figure 5. For the case of BCC model, Figure 5(a) shows that Δp/Δz predicted by Ergun’s equation is higher than KTA’s correlation. When Reh<1000, Δp/Δz calculated by laminar model and RNG k-ɛ model show good agreement with Yang et al.’s experiment data fitting. As Reh increases to relatively high numbers ([2000,5000]), the inertial effect dominates the pressure drop and Δp/Δz tends to be underestimated compared to the Yang et al.’s experimental data. While in the case of FCC model (Figure 5(b)), pressure drop predictions by k-ɛ model show good agreement with Ergun and KTA’s correlations when Reh<1000. When Reh>1000, and Δp/Δz of FCC model shows deviation with the empirical correlations and tends to Yang et al.’s experimental data fitting. One interesting phenomenon can be found in Figure 5, when Reh tends to relatively high values, pressure drop predicted by CFD model shows divergence from empirical correlations. We should note that the empirical correlations fitted from experimental data almost with gas and water as working medium. The viscosity of molten salt is much higher than that of gas and water; in the case of relative high Reh number, the pressure drop is dominated by inertial effect and the viscous force can be ignored. Compared to water and gas, high viscosity of molten salt will result in less vorticities in the packed bed, which lead to lower pressure drop under the same Reh number.Pressure drop of BCC and FCC model.
(a)
(b)Pebble to molten salt heat transfer coefficients (Nu number) is defined by (14) and (15). While the experimental data of Yang et al. [15] is obtained in the condition of fixed wall temperature. In this present work, we discard their data and compare the CFD results with Wakao’s equation and KTA’s recommended correlation. The Nu numbers proportioned to Reh number of BCC and FCC model are shown in Figure 6.Nu number of BCC and FCC model.
(a)
(b)In Figure6, the Nu number predicted by Wakao’s equation and KTA’s correlation shows significant difference and the NuKTA is about twice over NuWakao under the same Reh. In the laminar flow zone of BCC model (Figure 6(a)), the Nu numbers of the present work are lower than Wakao’s correlation predictions. While in the laminar zone of FCC model, the points of present Nu are located at the region between Wakao and KTA’s predictions. In addition, in the case of low Re (Re<100) number, the natural convection of flow and heat transportation induced by density gradient (buoyancy effect) can not be ignored. In this present work, the buoyancy effect is included by activating the Boussinesq model. The Nu number in the laminar zone shows lower values compared to Wakao (BCC) and KTA’s (FCC) predictions. As discussed in Section 1, the molten salt has low conductivity, high viscosity, and high Pr number, which means a thin thermal boundary fluid. The thermal diffusivity (α=k/ρcp) of molten salt is weak than metal, water, and even pressured gas. With the laminar flow pattern, heat transportation is dominated by thermal diffusion and natural convection, and low volume expansion ratio and high Pr number result in lower Nu number than other low Pr number fluid.With the increasing ofReh number, the flow pattern comes into laminar-turbulent (100<Reh<300) and fully turbulent (Reh>300) flow regime. The Nu numbers of the BCC model show good agreement with that of the Wakao’s equation. And for the case of FCC model the Nu number is shown to be higher than Wakao’s equation and lower than the KTA’s correlation. In the turbulent flow regime (Re≥103) heat transport from heated wall to flow mainstream center is dominated by turbulence heat transfer rather than thermal conduction, and the diffusive heat transport is proportional to the turbulence production in the area close to the rim of the viscous layer and depends on the velocity difference between rim and flow core. One conclusion can be derived which is that heat transfer capability of molten salt is significantly enhanced by turbulent mixing effect. In Figure 6 it is found that fluid heat transfer coefficients of FCC model are much higher than those of BCC model. The first reason is that under the same Reh number, the local velocities of FCC are higher than those of BCC model. And the other reason is that there exist more vorticities of fluid flow in BCC model than in FCC model (Figure 4), which induced the degeneration of heat transfer capability.Results shown in Figures5 and 6 indicate that the pressure drop and Nu number of packed bed are strongly affected by the pebble configuration and flow pattern. Moreover, pressure drop and Nu number in random pebble bed with molten salt are not investigated in this present work, which need further studies.


## 3.1. Basic Definition of Thermal-Hydraulic Correlations

The most widely used empirical correlations to predict pressure drop and local meanNu number are Ergun’s equation and Wakao’s correlation. Ergun’s equation was first proposed by Ergun [11],
(11)(ΔPΔL)Ergun=f·12(VDϵ)21dh,f=(c1Reh+c2),
where f is the friction factor, VD is Darcy velocity in porous medium, and c1=133 and c2=1.33 are constants in Ergun’s original work. Wakao’s type correlations were first introduced by Wakao et al. [12],
(12)NuWakao=hwf·dpkf=a1+a2Pr1/3·(ϵdpdh)n,
where hwf is the heat transfer coefficient of pebble wall to fluid, dp is particle diameter, a1, a2, and n are model constants with a1=2.0, a2=1.1, and n=0.6. The constant a1=2.0 corresponds to the heat transfer coefficient by solving the heat transfer equations for a single pebble immersed in a stagnant liquid with uniform temperature far from the pebble; the other term of (12) is similar to the correlations for smooth tube. In thermal hydraulics studies of gas cooled pebble bed reactors, correlations recommended by German regulatory commission (KTA correlations) [33, 34] were mostly used, which can be described as follows:
(13)(ΔPΔL)KTA=f·(1-ϵ)ϵ3·ρ2·VD2,f=320Red/(1-ϵ)+6(Red/(1-ϵ))0.1,NuKTA=1.27Pr1/3ϵ1.18Red0.36+0.033Pr1/2ϵ1.07Red0.86.
Pressure drop and Nu number predicted by correlations of Ergun, Wakao, and KTA show significant divergences especially the Reh number up to 103, which suggests careful validations in thermal hydraulics analysis of ordered and random molten salt cooled pebble bed reactors. In this present work, normalized pressure drop and local mean Nu number for each cell of BCC and FCC model can be defined as follows [14]:
(14)ΔpΔz=1N∑cellNicellNj‍1(zout-zin)Ain∫‍∫Ain‍(pout-pin)dA,Nu=dpkf1N∑cellNicellNj‍ρfcp∫‍∫Ain‍u(Tf,out-Tf,in)dAAp(T¯p-T¯f),
where u, Tf,in, Tf,out are flow velocity and mass flow temperature at (in,out) cross section of each cells, respectively. The fluid average temperature T¯f and pebble wall average temperature T¯wall in each cell can be defined as
(15)T¯f=12∫‍∫Ain‍u(Tf,in+Tf,out)dA∫‍∫Ain‍udA,T¯p=1N∑1N‍Twall.As shown in Figure1, the calculation domain consists of seven periodically packed cells in the flow direction. One important question is whether this simplified model can depict overall effects of “integrated bed” such as pressure drop and average Nu number. To illustrate this problem, Figure 3(a) shows the pressure drop of BCC and FCC model in the case of Reh=1000 of each cell. As can be seen in Figure 3(a), in the flow direction after the flow through two cells the pressure gradient is almost constants, this result indicating the flow appears to be periodically formed. This phenomenon was also observed in [22] by using numerical and experimental methods.Pressure drop and Nu number of BCC and FCC model.
(a)
(b)The local mean Nu number of each cell ofBCC and FCC model in the case of Reh=1000 is shown in Figure 3(b). In Figure 3(b), it can be found that the entrance effect of BCC model causing higher heat transfer capability of fluid flow, the cell averaged Nu number has significant change in Cell 1–3. The flow stream lines of BCC and FCC in the case of Reh=1000 are shown in Figure 4 for the sake of intuitive observation. Figure 4 shows that when the entrance effect decays along with the flow passes downstream through two cells, the flow is periodically developed, and thus the local Nu number tending towards stability in Cell 3–6 (Figure 3(b)). At the last cell of BCC model, because of a hollow located at the cell center, the turbulent mainstream can not be suppressed on the near-wall region which causes a decline of Nu number. In the FCC packed cells the entrance effect is observed at the very beginning of Cell 1; when flow passed through Cell 2, the entrance effect decays and the Nu number is almost equal in Cell 3–6.Streamline distributions in BCC and FCC packed cells.
(a)
(b)It is interesting to notice that, experimental results demonstrate that heat transfer from pebbles positioned in the entrance layers of a bed is lower than the average [35]. In [35] Achenbach indicated that the lower Nusselt numbers in first layer of packed bed are due to the fact that the velocity and turbulence levels of the incoming flow are smaller than the superficial velocity, and we can expect that the smaller velocity and turbulence levels should be owed to lower porosity in the entrance region. Unlike the random bed the ordered packed bed has a uniform porosity, flows acceleration, and strong turbulent mixing which induce higher heat transfer coefficient. This phenomenon has also been observed by several researches using experimental and numerical methods in the geometry of sphere and cylinder packed bed [36–38]. In Figure 4, it also can be found that the local flow velocity of the FCC cells is higher than BCC model. This observation supports the results of Figure 3(b); under the same Reh number, the Nu number of FCC cells is higher than BCC cells.In the next section, the normalized pressure drop and local average Nu number are defined by (14) with Ni=4 and Nj=5.


## 3.2.  Pressure Drop and Forced Convection Heat Transfer Investigations

In this section, density, dynamic viscosity, heat conductivity, andPr number are assumed as reference constant values when Reh>100. Thermal-physical properties are calculated by (16) at Tin=928.15K, which is the average operation temperature of PB-AHTR steady state [2]. Thermal properties are given in Table 4. Density and dynamic viscosity of FliBe [39] are as follows:
(16)ρ=(2.28-4.88*10-4*(Tin-273.15))*103[kg/m3],μ=0.116*103*exp⁡(3755Tin)[kg/ms].Table 4
Thermal-physical properties of FliBe atTin=928.15 K.
Densityρ [kg/m2]
1960.0
Dynamic viscosityμ  [kg/m s]
0.00663
Thermal conductivityk [W/(m K)]
1.1
Specific heatCv [J//(kg K)]
2415.0The pressure drop (Δp/Δz) of BCC and FCC model proportioned to Reh number is shown in Figure 5. For the case of BCC model, Figure 5(a) shows that Δp/Δz predicted by Ergun’s equation is higher than KTA’s correlation. When Reh<1000, Δp/Δz calculated by laminar model and RNG k-ɛ model show good agreement with Yang et al.’s experiment data fitting. As Reh increases to relatively high numbers ([2000,5000]), the inertial effect dominates the pressure drop and Δp/Δz tends to be underestimated compared to the Yang et al.’s experimental data. While in the case of FCC model (Figure 5(b)), pressure drop predictions by k-ɛ model show good agreement with Ergun and KTA’s correlations when Reh<1000. When Reh>1000, and Δp/Δz of FCC model shows deviation with the empirical correlations and tends to Yang et al.’s experimental data fitting. One interesting phenomenon can be found in Figure 5, when Reh tends to relatively high values, pressure drop predicted by CFD model shows divergence from empirical correlations. We should note that the empirical correlations fitted from experimental data almost with gas and water as working medium. The viscosity of molten salt is much higher than that of gas and water; in the case of relative high Reh number, the pressure drop is dominated by inertial effect and the viscous force can be ignored. Compared to water and gas, high viscosity of molten salt will result in less vorticities in the packed bed, which lead to lower pressure drop under the same Reh number.Pressure drop of BCC and FCC model.
(a)
(b)Pebble to molten salt heat transfer coefficients (Nu number) is defined by (14) and (15). While the experimental data of Yang et al. [15] is obtained in the condition of fixed wall temperature. In this present work, we discard their data and compare the CFD results with Wakao’s equation and KTA’s recommended correlation. The Nu numbers proportioned to Reh number of BCC and FCC model are shown in Figure 6.Nu number of BCC and FCC model.
(a)
(b)In Figure6, the Nu number predicted by Wakao’s equation and KTA’s correlation shows significant difference and the NuKTA is about twice over NuWakao under the same Reh. In the laminar flow zone of BCC model (Figure 6(a)), the Nu numbers of the present work are lower than Wakao’s correlation predictions. While in the laminar zone of FCC model, the points of present Nu are located at the region between Wakao and KTA’s predictions. In addition, in the case of low Re (Re<100) number, the natural convection of flow and heat transportation induced by density gradient (buoyancy effect) can not be ignored. In this present work, the buoyancy effect is included by activating the Boussinesq model. The Nu number in the laminar zone shows lower values compared to Wakao (BCC) and KTA’s (FCC) predictions. As discussed in Section 1, the molten salt has low conductivity, high viscosity, and high Pr number, which means a thin thermal boundary fluid. The thermal diffusivity (α=k/ρcp) of molten salt is weak than metal, water, and even pressured gas. With the laminar flow pattern, heat transportation is dominated by thermal diffusion and natural convection, and low volume expansion ratio and high Pr number result in lower Nu number than other low Pr number fluid.With the increasing ofReh number, the flow pattern comes into laminar-turbulent (100<Reh<300) and fully turbulent (Reh>300) flow regime. The Nu numbers of the BCC model show good agreement with that of the Wakao’s equation. And for the case of FCC model the Nu number is shown to be higher than Wakao’s equation and lower than the KTA’s correlation. In the turbulent flow regime (Re≥103) heat transport from heated wall to flow mainstream center is dominated by turbulence heat transfer rather than thermal conduction, and the diffusive heat transport is proportional to the turbulence production in the area close to the rim of the viscous layer and depends on the velocity difference between rim and flow core. One conclusion can be derived which is that heat transfer capability of molten salt is significantly enhanced by turbulent mixing effect. In Figure 6 it is found that fluid heat transfer coefficients of FCC model are much higher than those of BCC model. The first reason is that under the same Reh number, the local velocities of FCC are higher than those of BCC model. And the other reason is that there exist more vorticities of fluid flow in BCC model than in FCC model (Figure 4), which induced the degeneration of heat transfer capability.Results shown in Figures5 and 6 indicate that the pressure drop and Nu number of packed bed are strongly affected by the pebble configuration and flow pattern. Moreover, pressure drop and Nu number in random pebble bed with molten salt are not investigated in this present work, which need further studies.


## 4. Pebble Wall Temperature under Normal Operation Condition, Local Power Distortion, and Local Pebbles Insufficient Cooling

The basic design of pebble bed AHTR/FHR fuel element named “Pebbles” is inherited from the HTR, which is made of pyrolytic graphite and contains thousands of microfuel particles called “TRISIO” particles [40]. According to [40] the safety thresholds for maintaining the integrity of pebbles are the maximum fuel temperature below 1250° for steady state and 1600° for loss of coolant accident scenarios. In the well-established standard thermal-hydraulics code THERMIX/TINTE, pebble surface temperature was evaluated by macroscopic porous media method, but there was inability to accurately depict intuitionistic temperature fields distribution of pebbles’ surfaces and coolant round them. As indicated by Moormann [41], occurrence of possible local hot spots in the pebble bed possibly affects the pebbles integrity. The CFD methods could be a reliable tool for simulations of local hot spots. As has been discussed in Section 3, pebbles configuration strongly affects the heat transportation capability and the Nu number of FCC configuration is much higher than that of BCC under the same Reh. And we can speculate that in the random core, due to random configurations of pebbles and disordered flow streamlines, there may exist a certain number of vorticities like BCC model, which is a disadvantage for heat removal.In this section, in order to confirm cooling capacity of molten salt under nominal operation and other default conditions, we investigate the temperature distribution of pebble surface of module AHTR designed by UC-Brecly [1, 2]. Three default cases are defined, the nominal operation, local power distortion, and local pebbles insufficient cooling, respectively. In these cases, it is assumed that the calculation domain is located at middle of the reactor core where highest neutron dose exists. The local power distortion conditions are considered by promoting the heat flux of pebbles’ surfaces in Cells 4 and 5 of BCC and FCC model. The local pebbles insufficient cooling condition is defined by decreasing the inlet flow velocity. Density and dynamic viscosity of molten salt are defined by (16); other thermal-physical properties is given in Table 4. RNG-k-ɛ turbulent model is activated with enhanced wall treatment. Detailed descriptions of these cases are stated as follows.Case
1 (nominal operation). Inlet velocity is 0.315m/s; heat flux of pebble surfaces of all cells is 831650W/m2.Case
2a (local power distortion). Inlet velocity is 0.315m/s; heat flux of pebble surfaces of cells 4 and 5 is 1663300W/m2; heat flux at pebble surfaces of other cells is 831650W/m2.Case
2b (local power distortion). Inlet velocity is 0.315m/s; heat flux at pebble surfaces of cells 4 and 5 is 4158250W/m2; heat flux at pebble surfaces of other cells is 831650W/m2.Case
2c (local power distortion). Inlet velocity is 0.315m/s; heat flux at pebble surfaces of cells 4 and 5 is 6653200W/m2; heat flux at pebble surfaces of other cells is 83165W/m2.Case
3a (local pebbles insufficient cooling). Inlet velocity is 0.252m/s; heat flux of pebble surfaces of all cells is 831650W/m2.Case
3b (local pebbles insufficient cooling). Inlet velocity is 0.1575m/s; heat flux of pebble surfaces of all cells is 831650W/m2.Case
3c (local pebbles insufficient cooling). Inlet velocity is 0.063m/s; heat flux of pebble surfaces of all cells is 831650W/m2.Figures8–12 show the temperature distribution of pebble located at center of Cell 5 for BCC model and inlet face of Cell 5 for FCC model (Figure 7).Figure 7
Typical pebble.Pebble temperature distribution ofBCC (a) and FCC (b) model of Case 1 (view from inlet (left), outlet (middle)) and temperature distribution histogram.
(a)
(b)Pebble temperature distribution ofBCC model of Cases 2a, 2b, and 2c (view from inlet (left), outlet (middle)) and temperature distribution histogram.
(a)
(b)
(c)Pebble temperature distribution ofBCC model of Cases 3a, 3b, and 3c (view from inlet (left), outlet (middle)) and temperature distribution histogram.
(a)
(b)
(c)Pebble temperature distribution ofFCC model of Cases 2a, 2b, and 2c, (view from inlet (left), outlet (middle)) and temperature distribution histogram.
(a)
(b)
(c)Pebble temperature distribution ofFCC model of Cases 3a, 3b, and 3c, (view from inlet (left), outlet (middle)) and temperature distribution histogram.
(a)
(b)
(c)Figure8 shows the temperature distribution of pebble surface of BCC model (a) and FCC model (b) of Case 1. As can be seen from Figure 8(a), temperature of four arc regions at the pebble bottom exceeds 1076 K, and the average temperature of whole pebble surface (Tpave) is 1013 K. The histogram in Figure 8(a) shows that the interval of temperature distribution is about 40 K, which indicates a quite uniform temperature distribution of pebble surface. Figure 9 shows the temperature distribution of pebble surface of BCC model of Cases 2a, 2b, and 2c. Tpave of pebble surface in Cases 2a, 2b, and 2c are 1083, 1269, and 1436 K, respectively. In Case 2a, the surface temperature is mainly distributed in the interval of 1000–1120 K (Figure 9(b)), which is under the limit of pebble temperature safety threshold. As the heat fluxes promote times of nominal case, the Tpave increase significantly. The default conditions of Cases 2b and 2c were unlikely to occur in steady state of AHTR, studies of these cases may provide useful references to some extreme conditions. When the heat flux is promoting to five times of nominal cases, the temperature of pebble surface about 10%  As is over 1450 K, which may induce a certain number of TRISIO fuel temperature exceed the temperature safety threshold. In Case 2c, the temperature of pebble surface about 25%  As is exceed 1500 K (Figure 9(c)), which is unacceptable for steady-state operation.One major cause of local hot spots of HTR is insufficient cooling of pebble caused by highly local flow resistance [41]. Detailed studies of 3D local thermal-hydraulics characteristics in the condition of local pebble insufficient cooling can provide useful suggestions to optimize of liquid salt pebble bed design and safety analysis.Figure10 shows the temperature distribution of pebble surface of Cases 3a, 3b, and 3c in BCC model. In Case 3a, compared to Case 1, the rise of Tpave is 13.4 K, and rise of maximum pebble surface temperature (Tpmax⁡) is 40 K, which is induced by decreasing 20% of coolant. In Case 3b, rise of Tpave is 47.5 K, and rise of Tpmax⁡ is 103.7 K, which is induced by decreasing 50% of coolant. In Case 3c, rise of Tpave is 52.5 K, and rise of Tpmax⁡ is 351.4 K; meanwhile, the inlet coolant is only 20% of nominal case. The results of Case 3 indicate that temperature of pebble surface is not very sensitive to loss of local coolant in molten salt cooled reactor. Moreover, it is worth noting that the density of liquid salt is higher than pebble fuel, which means densification phenomenon of pebbles in pebble bed FHR is much less than gas cooled pebble bed reactor.Figures11 and 12 show the temperature distribution of pebble surface of FCC model of Cases 2 and 3. It can be seen that temperature rise induced by power distortion and local pebbles insufficient cooling is lower than corresponding case of BCC model. As has been discussed in Section 3, the Nu number of FCC model is higher than that of BCC model under the same Reh number, and the heat removal capacity of FCC model is superior to BCC model, which results in a lower surface temperature of FCC model under the same condition.


## 5. Conclusion

Pore scale thermal-hydraulics studies of molten salt cooled pebble bed high temperature reactor withBCC and FCC model have been performed. Pressure drop and Nusselt number were calculated and compared to empirical correlations. Pebble surface temperature was evaluated in several default cases. The obtained results are stated as follows:(i)
The overall effects of packed bed such as pressure drop and meanNu number can be reasonably derived by simplified BCC and FCC model.(ii)
Thermal-hydraulics characteristics of packed bed are strongly affected by the spheres configuration and flow pattern. Simulation results show that the mostly used empirical correlations of fraction factor andNu number should be carefully verified in ordered or random pebble bed.(iii)
Individual pebble surface temperature has been calculated by using of RNG-k-ɛ turbulent model and enhanced wall treatment. The result of nominal case shows the well enough cooling capability of molten salt. Temperature distribution in the case of local power distortion shows that the pebble’s temperature may exceed safety threshold in case of fuel element long exposure in five times neutron dose of nominal operation condition.(iv)
Numerical experiments about local pebbles insufficient cooling show that temperature of pebble surface is not very sensitive to loss of partial local coolant in molten salt cooled reactor. We can expect that in molten salt cooled pebble bed reactor the local hot spots should not be a serious problem.(v)
The thermal-hydraulics characteristic of random packed bed is not investigated in this present work, and the applicability of empirical thermal-hydraulics correlations needs further studies.



---

*Source: 589895-2014-06-23.xml*