Nuclear Science and Techniques

NUCLEAR ENERGY SCIENCE AND ENGINEERING

DEM-CFD simulation of modular PB-FHR core with two-grid method

Feng-Rui Liu，

Xing-Wei Chen，

Zhong Li，

Na-Xiu Wang

Nuclear Science and Techniques

Vol.28, No.7

Article number 100

Published in print 01 Jul 2017

Available online 06 Jun 2017

DOI：10.1007/s41365-017-0246-3

6201

For designing and optimizing the reactor core of modular pebble-bed fluoride-salt-cooled high-temperature reactor (PB-FHR), it is of importance to simulate the coupled fluid and particle flow due to strong coolant-pebble interactions. Computational fluid dynamics (CFD) and discrete element method (DEM) coupling approach can be used to track particles individually, while it requires a fluid cell being greater than the pebble diameter. However, the large size of pebbles makes the fluid grid too coarse to capture the complicated flow pattern. To solve this problem, a two-grid approach is proposed to calculate inter-phase momentum transfer between pebbles and coolant without the constraint on the shape and size of fluid meshes. The solid velocity, fluid velocity, fluid pressure and void fraction are mapped between hexahedral coarse particle grid and fine fluid grid. Then the total inter-phase force can be calculated independently to speed up computation. To evaluate HYPERLINK "http://cn.bing.com/dict/search?q=suitability&FORM=BDVSP6&mkt=zh-cn" suitability of this two-grid approach, the pressure drop and minimum fluidization velocity of a fluidized bed were predicted, and movements of the pebbles in complex flow field were studied experimentally and numerically. The spouting fluid through a central inlet pipe of a scaled visible PB-FHR core facility was set up to provide the complex flow field. Water was chosen as liquid to simulate the molten salt coolant, and polypropylene balls were used to simulate the pebble fuels. Results show that the pebble flow pattern captured from experiment agrees well with the simulation from two-grid approach, hence the applicability of the two-grid approach for the later PB-FHR core design.

PB-FHRPebble flowDEM-CFDTwo-grid approach

1 Introduction

As a novel reactor concept, pebble bed fluoride salt-cooled high-temperature reactor (PB-FHR) combines the advantages of liquid fluoride salts, separate buffer salt pool and pebble fuel together [1]. Fluoride salts has low neutron cross section and high heat capacity, which means higher operating temperature and better safety margin [2]. A closed primary-salt loop immersed in a buffer salt pool is used to control the temperature of primary loop components. Hundreds of thousands of fuel pebbles recirculate through reactor core. Fluoride is driven to pass through the pore space of pebble bed to extract the fission heat for electricity generation [3]. A full understanding of the pebble fuel and coolant behavior in pebble bed reactor (PBR) is important to assess the reactor core performance and operation safety. This is especially prominent when one tries to optimize the PBR design such as online refueling, loss-of-coolant accident and earthquake effect on the core stability [4].

Models for studying solid-fluid flows have been investigated with a large number of numerical algorithms at different levels of detail. Different numerical methods can be mainly classified into direct numerical simulation (DNS) models, the Euler–Lagrange (EL) models and two-fluid models.

DNS models are applied to resolve all details of the flow, heat and species transport at the smallest relevant length scales. The exchange of conserved quantities of mass, momentum and energy between a single particle and the fluid is modeled via, instead of empirical correlations, boundary conditions of the particle surface [5]. It resolves the flow between particles to get accurate correlations for coarse-grained models. DNS methods proposed so far include the Lattice Boltzmann Method (LBM), the Immersed Boundary Method (IBM), etc. [6-8], but they were expensive even for the reactor core at lab scale, owing to the complexity of particle-flow interactions. However, with increasing computer capacity and rapid developments in computational science, DNS methods have received considerable attention to study complex multiphase flows.

The numerical prediction of engineering-scale multiphase flow equipment can be achieved with continuum models such as the Two-Fluid Model (TFM) [9]. Both the fluid and solid phases are described as fully interpenetrating continua. The interaction between the two phases is incorporated by drag force correlations, which depend on the local relative velocity of the two phases and the local solid volume fraction. The internal momentum transport in the solids phase is commonly based on the kinetic theory of granular flow (KTGF) [10]. The drawback of this method is that it does not adequately model the details of particle-particle and particle-gas interactions.

The EL models such as computational fluid dynamics coupled with the discrete element model (DEM-CFD) can combine the advantages of DNS and TFM. The general theoretical framework has been established for DEM-CFD models [11]. In this model, the fluid phase is described as a continuous phase, while the solid phase is described by the Newtonian equations for each individual particle. EL models can capture the discrete nature of the particle phase while maintain the computational efficiency by ignoring the details of fluid field at particle level. Currently, DEM-CFD can be used to simulate systems with millions of particle numbers. Thus, it is suitable to simulate the pebble flow in modular PB-FHR core with tens of thousands pebble fuels efficiently [12]. Using DEM-CFD model, Yanheng LI and Wei JI simulated complicated behavior of a PBR reactor core with helium gas-cooled and fluoride salt-cooled designs [13]. Noticeable changes were observed, such as higher pebble density in the cylindrical core region and more uniform vertical fluid speed profile due to the coupling effect. These suggested good balance between fidelity and efficiency of the coupling method compared with uncoupled methodology in analyzing pebble flow in PBRs. However, the work did not demonstrate the grid independence of the algorithm and topology of fluid mesh, especially the mesh in the conical top boundary where the grid density was difficult to be large enough.

Generally, DEM-CFD requires a much larger control volume than the particle size for proper treatment of the drag force, which restraint the coolant resolution. However, fine fluid mesh is usually needed when one tries to capture fluid details especially for the modular PB-FHR core with low bed-to-particle ratio. The meshes with cell size comparable to or smaller than the particle diameter are often required to resolve complicated flow fields with large velocity gradients, so the size restriction from solid and fluid phases often contradict with each other. On the other hand, too few particles in a cell result in large spatial and temporal fluctuations of the solid fraction field, hence unstable numerical integration.

For breaking up the restriction of fluid mesh size, S. Deb and S.K. Tafti developed a typical two-grid formulation for systems involving DEM–CFD coupling in a parallel processing framework [14]. A coarser particle grid was introduced to transfer fluid field variables, interphase transfer terms and void fractions between particle and fluid grid. For simulation of circulating fluidized bed, Falah Alobaid et al. determined the physical values of fluid and particle phases in two separated grids, which allowed the resolution of fluid grid to be independent of the particle size, with improved calculation accuracy [15,16]. In this formulation, the fluid flow equations were solved on a fine grid, whereas the discrete particle equations were solved on a coarse grid. Fluid field variables, interphase transfer terms and void fractions were transferred between two grids using proper mapping methods. The results showed that this model could predict the particles motion and the pressure gradients in the circulating fluidized bed accurately.

However, the fluid particle grid mapping in this method requires that the two types of orthogonal meshes should alignment each other. This strong geometrical restriction of the mesh shape is difficult to be implemented in non-orthogonal domain. In order to use unstructured CFD mesh, S. Rickelt et al. employed an additional orthogonal fine transfer grid to transfer data between solid and fluid phase [17]. Temperatures, velocities and fluid properties of each fluid cell were transferred to the DEM code through transfer grid. This model requires that the transfer grid should be much finer than the fluid mesh to perform interpolation, which may constraint the fluid grid resolution.

In this paper, we extend the two-grid method to make it suitable for the fluid domain with irregular shapes and change the communication pattern between fluid and solid phase for higher efficiency. In this method, the coarse grid is defined on both sides (DEM and CFD). A conservative interpolation scheme is employed between arbitrarily polyhedral fluid mesh and hexahedral coarse mesh. Data are transferred only at the start of each coupling time step, and the governing equations of DEM and CFD are solved simultaneously until next coupling step. This algorithm is firstly validated against the case of fluidized bed with simple flow pattern [18]. In order to test the feasibility of the code for complex flow characteristics in the PB-FHR core, the corresponding laboratory-scale experiments with polyethylene spheres are conducted at a specially designed jet flow environment. The pebble flow under the influence of spouted fluid is studied experimentally and numerically under a range of operating conditions. Good agreement is observed in the motion of particle group, which provides a concrete base for further development and study.

2 Numerical model

2.1. Particle motions

The discrete element method, a powerful numerical tool proposed by Cundall and Strack [19], has been extensively used in granular systems. The granular media is represented as an assemblage of discrete particles interacting with one another. According to Newton's second law of motion, the motion of spheres is obtained by solving the translational and rotational equations. The corresponding governing equations of particle i can be written as:

m_{i} \frac{d {\vec{v}}_{i}}{d t} = {\vec{f}}_{p, i} + \sum_{j = 1}^{N_{c}} ({\vec{f}}_{n, i j} + {\vec{f}}_{t, i j}) + m_{i} \vec{g}

(1)

I_{i} \frac{d {\vec{w}}_{i}}{d t} = \sum_{j = 1}^{N_{c}} ({\vec{r}}_{i} \times {\vec{f}}_{t, i j})

(2)

where, v_i, w_i, m_i, and I_iare the linear velocity vector, angular velocity vector, mass and moment of inertia of particle i respectively; f_n,ij f_t,ij are normal and tangential collision contact forces between particle i and neighboring particle j; and f_p,i is the drag on individual particle i. According to the Hertzian spring-dashpot model, they can be written as:

{\vec{f}}_{n, i j} = (- k_{n} δ_{n, i j}^{3 / 2} - η_{n} {\overset{⇀}{v}}_{n, i j} \cdot {\overset{⇀}{n}}_{i j}) {\overset{⇀}{n}}_{i j}

(3)

{\vec{f}}_{t, i j} = - \min (k_{t} δ_{t, i j} + η_{t} {\overset{⇀}{v}}_{t, i j} \cdot {\overset{⇀}{t}}_{i j}, μ | {\vec{F}}_{n, i j} |) {\overset{⇀}{t}}_{i j}

(4)

where, n_ij and t_ij are normal and tangential unit vectors, respectively, from particle i to particle j; δ_n,ij and δ_t,ij are the displacements of particle caused by normal and tangential forces, respectively; the normal stiffness coefficient k_n,ij, tangential stiffness coefficient k_t,ij, normal damping coefficient η_n,ij, and tangential damping coefficient η_t,ij are calculated by:

k_{n} = \frac{2 E}{3 (1 - σ^{2})} \sqrt{\frac{R_{i} R_{j}}{R_{i} + R_{j}}}

(5)

k_{t} = k_{n}

(6)

η_{n} = η_{t} = - \frac{2 \ln e}{\sqrt{π^{2} + \ln^{2} e}} \sqrt{m k_{n}} δ_{n, i j}^{1 / 4}

(7)

where, E, e and σ are Young’s modulus, restitution coefficient and Poisson coefficient, respectively.

2.2 Navier–Stokes equations for fluids

The fluid phase is modeled by the incompressible volume-averaged Navier–Stokes equations, which can be described as:

\frac{\partial ρ_{f} ε_{f}}{\partial t} + \nabla \cdot (ρ_{f} ε_{f} \vec{U}) = 0

(8)

[\frac{\partial ρ_{f} ε_{f} \vec{U}}{\partial t} + \nabla \cdot (ρ_{f} ε_{f} \vec{U} \vec{U})] = - ε_{f} \nabla p - \vec{F} + \nabla \cdot (ε_{f} τ) + ρ_{f} ε_{f} \vec{g}

(9)

where, U is the fluid velocity, p is the pressure, ρ_f is the density, τ is the viscous stress tensor, g is acceleration due to gravity, and ε_f is the local voidage of the fluid cell. The turbulent effect is represented by

\nabla \cdot (τ) = \frac{\partial}{\partial x_{j}} [(μ_{f} + μ_{t}) (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})]

(10)

where μ_f is fluid viscosity and μ_t is turbulence viscosity.

F, the momentum exchange rate per volume between fluid and particles, is given by:

\vec{F} = - \frac{1}{Δ V} \sum_{i = 1}^{N} (α_{i} {\vec{f}}_{p, i})

(11)

where, ΔV is the volume of fluid cell, α_i is the volume fraction of the particle i that falls into the specific fluid cell, and N is the number of particles within the cell.

For simplicity, we ignore the effect of solid phase on liquid turbulent kinetic energy and dissipation. The governing transport equations for turbulent energy k and turbulent dissipation rate σ are:

\frac{\partial ε_{f} ρ_{f} k}{\partial t} + \nabla \cdot (ε_{f} ρ_{f} k \vec{U}) = \nabla \cdot (ε_{f} (μ_{f} + \frac{μ_{t}}{σ_{k}}) \nabla k) + ε_{f} ρ_{f} G_{k} - ε_{f} ρ_{f} σ

(12)

\frac{\partial ε_{f} ρ_{f} σ}{\partial t} + \nabla \cdot (ε_{f} ρ_{f} σ \vec{U}) = \nabla \cdot (ε_{f} (μ_{f} + \frac{μ_{t}}{σ_{σ}}) \nabla σ) + ε_{f} \frac{σ}{k} (C_{1} ρ_{f} G_{k} + C_{2} ρ_{f} σ)

(13)

where, G_k is the generation of turbulence kinetic energy owing to the mean velocity gradients, C₁ and C₂ are model constants of 1.44–1.92, and σ_k= 1.0 and σ_σ = 1.3 are the turbulent Prandtl numbers.

2.3. Fluid–particle interactions

The main contributions to the interphase exchange terms include the drag force, pressure gradient force, virtual mass term, Saffman liftforce, Magnus lift force and the Basset history term [11]. For dense and slow-moving pebble flows in liquid environment, the lift force and history term become negligible and is therefore ignored in the simulations presented herein. f_p,i in Eq. (1) denotes the forces due to fluid-particle interactions, e.g., pressure gradient force, viscous stress tensor gradient force, drag force and virtual mass force:

{\vec{f}}_{p, i} = {\vec{F}}_{s, i} + {\vec{F}}_{d, i} + {\vec{F}}_{v, i}

(14)

The pressure gradient force and viscous stress tensor gradient force due to variations of the fluid stress tensor on particle i can be written as:

{\vec{F}}_{s, i} = - V {}_{i}\nabla p_{i} + V_{i} \nabla \cdot τ

(15)

The drag force on particle i can be written as

{\vec{F}}_{d, i} = [\frac{β}{1 - ε_{f}} (\vec{u} - {\vec{v}}_{i})] V_{i}

(16)

where, u is gas velocity, v_i is particle velocity, p is gas pressure, V_i is particle volume, and β is the drag force coefficient, which can be expressed by a combination of the Ergun equation for the dense granular regime and the Wen-Yu correlation for the dilute regime:

β = {\begin{cases} 150 \frac{{(1 - ε_{f})}^{2} μ}{ε_{f} d_{p}^{2}} + 1.75 \frac{ρ_{f} (1 - ε_{f})}{d_{p}} | \vec{u} - {\vec{v}}_{i} | ε_{f} \leq 0 .8 \\ \frac{3}{4} C_{D} \frac{ε_{f} ρ_{f} (1 - ε_{f})}{d_{p}} | \vec{u} - {\vec{v}}_{i} | ε_{f}^{- 2.65} ε_{f} >0 .8 \end{cases}

(17)

where d_p is particle diameter, and C_D is the drag coefficient for a single unhindered particle evaluated by

C_{D} = {\begin{cases} 24 (1 + 0.15 R e_{p}^{0.687}) / R e_{p} R e_{p} \leq 1000 \\ 0.44 R e_{p} >1000 \end{cases}

(18)

Re_p, particle Reynolds number, can be written as:

R e_{p} = \frac{ε_{f} ρ_{f} | \vec{u} - {\vec{v}}_{i} | d_{p}}{μ_{f}}

(19)

The virtual mass force due to the pressure gradient on particle i can be written as:

{\vec{F}}_{v, i} = C_{vm} V_{i} ρ_{f} (\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} - \frac{\partial {\vec{v}}_{i}}{\partial t} - {\vec{v}}_{i} \cdot \nabla {\vec{v}}_{i})

(20)

where C_vmis the virtual mass effect coefficient.

2.4. Two-grid approach

In this work, we used two grid types: the hexahedral particle mesh and the polyhedral fluid mesh. The particle mesh is large enough to contain a sufficient number of particles for the locally averaged N-S equations. The size and shape of fluid grid depend on just the fluid resolution and stability of discretization scheme. The whole algorithm is divided into two parts (CFD and DEM). A fluid grid and a particle grid are defined on CFD part. The particle grid is defined on DEM part.

In CFD part, fluid velocity and fluid pressure on the fluid grid (U and p) are calculated based on N-S equation and mapped to the coarse grid as U_c and p_c. The U_c and p_c are then copied to the coarse grid in DEM part.

In DEM part, void fraction and solid velocity on the coarse grid (ε_c and Us_c) are evaluated based on particle centroid method or divided particle volume method and copied to the coarse grid in CFD part. The ε_c and Us_c are then mapped to the fine grid as ε and Us.

Once the data transfer has been finished, both the DEM part and CFD part maintain the four quantities (U_c, P_c, ε_c, and Us_c) defined on coarse grid respectively, so both sides can evaluate the interphase momentum without synchronization. The fluid-particle interaction force obtained will be used to update the movement of particles directly in DEM part. However, it should be mapped back to fine grid in CFD part since the N-S equations are solved on fluid fine mesh.

Compared with the traditional DEM-CFD scheme, data exchange only occurs at the start of coupling procedure in the algorithm we proposed, then the governing equation of fluid and solid can be solved independently until next coupling time step. Both parts do not wait for each other to exchange fluid-particle momentum any longer. The coupling procedures are outlined in Fig. 1.

Fig. 1

Coupling procedures.

2.4.1. Mapping scheme

In this article, mesh-to-mesh interpolation between hexahedron and polyhedron is based on the conservative interpolation [21]. The interpolation weight relies on the intersection between two grids. In order to calculate the intersection volume between a hexahedron and an arbitrary polyhedron, the lines and faces of two control volumes are traversed to evaluate the intersection points, then the convex hull can be constructed from these points and the cell vertices based on the Quickhull algorithm in CGAL [22]. A schematic illustration of the two-dimensional process is given in Fig. 2. The black circles represent particles. Suppose the quadrilateral A₀A₄B₆B₅ is the particle cell i and the three polyhedrons j₁, j₂and j₃ are fluid cells. Take the cell j₂ as example, the value of cell j₂ is given by:

Fig. 2

Coarse and fine 2D grid.

ϕ_{j 2} = ϕ_{i} \frac{S_{i, j 2}}{S_{i}}

(21)

where, S_i,_j2 is intersection area between cell j₂ and cell i. The intersection point of line of j₂ and cell i is A₃ and B₀. The source and target cell vertices that do not lie outside each other, is A₇, A₉ and A₄. The five vertices can construct a convex hull with the area of S_i and j₂. The weight of fine cell j₂ on the coarse cell i is the ratio between the j₂ and S_i. The value of cell i is given by:

ϕ_{i} = \frac{\sum_{k = 1}^{N} ϕ_{j} S_{i, j k}}{S_{i}}

(22)

where N is the total number of fluid cells that overlap with particle cell i.

2.4.2. Solid fraction

The ε_c and Us_c on DEM part are usually calculated based on divided particle volume method [20]. A particle's volume is divided into all cells it overlaps. For the cell k inside the computational domain, the solid volume fraction is calculated by summing up all the parts of particles that overlap with the cell k, which is given by:

ε_{s, k} = \frac{\sum_{i \in c e l l k} α_{i, k} V_{p, i}}{V_{k}}

(23)

where, V_p,_i is the volume of the particle i, V_k is the volume of the coarse cell k and α_i,k is the volume fraction of the particle i that falls into the cell k. The solid phase velocity Us_c in cell k is computed in a similar way.

For the particle cell that lies outside the fluid domain or occupies part of computational domain, solid volume fraction is calculated in a different way, as shown in Fig. 3. The dashed area represents the intersection part between the cylindrical fluid domain and one coarse cell. The black circles are particles next to the container wall. The solid volume fraction field is given by

Fig. 3

Coarse cell partially occupies the domain.

ε_{s, k} = {\begin{cases} \frac{\sum_{i \in cell k} α_{i, k} V_{p, i}}{V_{k}} / \frac{V_{cut}}{V_{k}} V_{cut} > 0 \\ 0 V_{cut} = 0 \end{cases}

(24)

where V_cut is the intersection volume between particle cell k and fluid domain.

3 Numerical simulations of fluidized bed

The proposed two-grid DEM-CFD method was firstly validated against the test case of C. Goniva et al. [18]. In their case, the bed-to-particle-diameter ratio was 27. When the fluid mesh was much larger than the particle diameter (4 times), there would be only seven grids in radial direction, which is too coarse to capture the flow characteristics. Using the two-grid approach, we used three types of fluid meshes in different resolutions to test the grid independence of algorithm. All the simulation parameters (geometry, boundary conditions, and particle properties) but the fluid grid were exactly the same as the case of Ref.[18]. The boundary condition of fixed value velocity inlet was applied on the bottom of the cylinder. Inlet velocities ranged from 0 to 0.02 m/s. The cylinder top was applied with a fixed pressure outlet. The side wall was applied with no slip condition. On the lower part of the cylinder, 10000 particles were randomly packed.

The pressure drop at different inlet velocities was predicted and compared with theoretical value. The conditions and parameters used in the simulations are listed in Table 1.

Simulation conditions and parameters

Items	Value
Domain radius and height (m)	0.0138, 0.553
Gas Density (kg/m³)	1
Gas viscosity(Pa·s)	1.5×10^-4
Gas superficial velocity (m/s)	0.002–0.02, in 0.002 steps
Fluid time step (s)	5×10^-4
Number of particles	10000
Particle diameter (m)	0.001
Particle density (kg/m³)	200
Normal spring-stiffness (N/m)	1×10⁴
Friction coefficient	0.3
Restitution coefficient	0.9
Solid time step (s)	1×10^-5

The three types of grids used in simulation are shown in Fig.4. The bold lines are coarse grid, and the fine lines are fluid grid. Their ratios of average fluid mesh size to particle diameter (L/d) are 2.5, 1.1 and 0.55, respectively. The length of coarse cell is 4 times larger than particle diameter.

Fig. 4

Top view of the three fluid meshes with L/d =2.5, 1.1 and 0.55, respectively.

The simulation results are shown in Fig. 5. The theoretical value of pressure drop for a fixed bed was evaluated by Ergun equation [23]. Below the minimum fluidization velocity of about 0.012 m/s, the pressure drops calculated at different fluid grids, and the theoretical value, agree well with each other. It means that, for small fluidized bed with simple flow pattern, the algorithm can quantify the particle-fluid momentum exchange based on the fluid grid with fine enough resolution.

Fig. 5

Pressure drop vs. inlet velocity at different L/d ratios.

4 Experimental procedures and the results

4.1. Experimental procedures

In the PB-FHR core, pebbles are removed from near the core top and fed back to the lower inlet plenum. Before being fully burnt, the pebbles go through the core several times in quit low, i.e. the slow continuous flow of pebbles through the reactor core. For reproducing the pebble transport hydrodynamics of the PREX facility, A. Griveau did scaled experiments at approximately 50% geometric scaling [24], using water and polyethylene spheres that matched the pebble-to-liquid density ratio, and the Reynolds and Froude numbers. However, the aim of this work is to validate the algorithm against complex pebble-coolant flow pattern instead of HYPERLINK "http://cn.bing.com/dict/search?q=addressing&FORM=BDVSP6&mkt=zh-cn" addressing an actual working condition in scaled facility, so we canceled the pebble recirculation passage, but kept the pebble-to-liquid density ratio and chose the jet flow to form complex flow environment.

Considering the limited space of laboratory, the polypropylene spheres were sized at Φ1.5 cm and the column-to-particle diameter ratio was 23.3. The scaled experimental facility is shown in Fig. 6. The water circulation system contained a water tank, an acrylic glass cylindrical container with hexahedral blanket, and pebbles. A pump connected to the water tank was used to control the flow rate. The water flew out of the cylindrical core from four outlets on the side wall of the top cover. The conic top had about 2000 small holes to facilitate the outflow of the water and keep the packing stable.

Fig. 6

Scaled PB-FHR facility.

A total of 18000 pebbles were poured into the cylindrical core, the conic top was closed and the pumping was started. The water flow rate was kept low enough so as not to interrupt the outline of packing. After the pebbles floated upward and the steady state was achieved, the flow rate was increased to produce different fluid flow patterns.

For quantitative investigation of the flow rate effect on particle movement, two indexes were chosen as indicators of degree of particle circulation. One is the average height of particles moving downward next to the wall, denoted as HT, while another is the average height of particles at the bottom of the packing, denoted as HB. Compared with other indexes, they are easy to be observed and measured. The definition of HT and HB are shown in Fig. 7.

Fig. 7

Schematics of the facility (in mm).

4.2. Model solution procedure

The symmetric section of the initial state of the facility is shown in Fig. 7. The computational domain was the region of dotted lines, in 800-mm height and 350-mm width. The conic top was of 100 mm height and the inlet pipe was of Φ50 mm.

The fluid equations were solved based on finite volume method. The semi-implicit method for pressure-linked equation (SIMPLE) algorithm was used to solve the fluid flow equations. The quadratic upwind interpolation of convective kinematics (QUICK) scheme and a central-differenced scheme was employed for the discretization of the convection term and the diffusion term respectively. The motion of fluid phase was solved with standard k-ε two equations turbulence model.

Proper boundary and initial conditions should be assigned for different boundaries of the domain to enclose the governing equations. The inlet boundary was the circular area at the bottom of the cylindrical container in radius of less than 0.05 m. The conic top was treated as outlet boundary. The remaining walls were wall boundary. Initially the particles were randomly packed in the upper part of domain in porosity of about 0.4. For the kinetic energy of turbulence and energy dissipation rate of fluid, the fixed value was assumed at the inlet and zero gradient boundary conditions at the other boundaries. A fixed value of inlet velocity, non-slip boundary condition and zero gradient boundary conditions were assigned for the velocity of inlet, wall and outlet boundary. Zero gradient boundary conditions were assigned for the pressure of inlet and wall boundary and the fixed value was applied for the outlet boundary.

To evaluate the flow rate effect on particle movement quantitatively, six inlet velocities in experiments were applied. The DEM and CFD time steps for satisfactory convergence were 1×10⁻⁵ s and 2×10^–3 s, respectively. The simulation conditions and parameters are summarized in Table 2.

Simulation conditions and parameters

Items	Values
Fluid Density (kg/m³)	1000
Fluid viscosity (Pa s)	1.0×10⁻³
Inlet velocity (m/s)	1.71, 1.95, 2.16, 2.33, 2.43, 2.59
Fluid time step (s)	2×10⁻³
Number of particles	18000
Particle diameter (m)	0.015
Particle density (kg/m³)	910.4
Normal spring-stiffness (N/m)	1×10⁴
Friction coefficient	0.2
Restitution coefficient	0.9
Solid time step (s)	1×10⁻⁵

All cases were simulated for 15 s to obtain statistically steady state of time-average values. Two symmetric hexahedral meshes of different resolutions were generated to analyze the grid independence of the algorithm. The size of coarse mesh was one-eighth the diameter of cylindrical container. The grids of the fluid and initial packing simulated are shown in Fig. 8. Square grids represent coarse mesh and red points are the center of particles of initial packing.

Fig. 8

Fluid mesh 1, mesh 2 and initial packing.

4.3. Simulation results

4.3.1. Particle circulation indexes

Unlike the gas-solid system, particles driven by net buoyance floated upward automatically until steady state of the facility was reached. When the fluid inlet velocity was less than 1.71m/s, particles next to the wall remained stationary. At higher fluid velocities, the drag force and pressure gradient could overcome the gravitational force acting on the particles. Pebble circulation developed gradually, which could be observed by movement of the pebbles next to the wall. In this case, particles in the central spout region were accelerated by the fluid to the upper part of the packing, split into the annular region and then flew downwards to the bottom of the packing to feed into the circulation back and forth.

Assuming that the container bottom is the zero height and Z-axis direction is positive, the HT and HB of particles at different inlet velocities are shown in Fig. 9. With increasing inlet velocities, HT moves upwards and HB moves downwards, which means more particles participate in pebble circulation. The overall trends of the HT and HB from experiments and simulations are very similar, though both show some fluctuations for these cases. This oscillation of the height of HB can be explained by relative errors of the method to locate the particles next to the wall.

Fig. 9

HT(a) and HB(b) as a function of inlet velocity.

When determining the HT and HB, we located the moving particles next to the wall at the front, back, left and right, and then averaged over them. For each direction, the deviation of the height of particle centre might exceed one particle diameter. The maximum deviation of HT and HB between experimental and simulation results did not exceed 2d (two times the diameter of particle), which means that the fluid-solid interactions can be evaluated quantitatively by the proposed two-grid DEM-CFD code. The deviation of the two meshes of different resolutions also falls into the acceptable error range, so the proposed algorithm allows fine enough fluid mesh to capture the details of complicated flow pattern.

4.3.2 Particle trajectories and fluid velocity profiles

Trajectories of tracer particles and detailed information of the movement behavior were obtained from the simulation results. Eight clusters of particles distributed uniformly in the radial direction were chosen as studying object. There were ten tracer particles next to the wall in each cluster. Fig. 10 shows the top view of the tracking particles.

Without loss of generality, the inlet velocity effect on the movement of tracer particles were investigated at 1.95, 2.33 and 2.59 m/s. The tracer particles were selected as follows: The chosen particles were those with their axial coordinate around HT at the last time step; the particle ID was determined, and the trajectories over the whole time were obtained. In this way, we could discover how particles reach the top area next to the wall eventually. In Fig.11, the particle trajectories at different operation conditions are shown on the left side, and the corresponding 3D velocity profiles of fluid on the right side. Different color helps to distinct particles in each cluster. It can be seen that fluid will split from the spout region to annulus region to form a vertex at large scale, which is accompanied by the formation of pebble circulation.

Fig. 10

Top view of sample particles.

Fig. 11

Particle trajectories and fluid velocity profile.

From Fig.11, the domain of vertex is larger than the region of pebble circulation since the movement of particle can only be activated by a large-enough fluid velocity field. Trajectories of some particles form an involute instead of a close-loop. Colored tracers next to the wall might move to the interior of the packing through a few cycles owing to radial mixing, but it was difficult to visualize this in the experiment. The affected area of vertex increased with the inlet superficial fluid velocity. At the inlet velocity of 1.95m/s, the particles just formed a half of loop; while at 2.33 and 2.59m/s, they circled at least one time and the circulation region expanded both upstream and downstream correspondingly. The amplitude of the upward circulation was greater than that of the downward circulation since the bottom of vertex is close to the bottom of container. This can be seen from the rate of change of HT and HB in Fig.9. The higher the inlet superficial velocity, the larger the domain of vertex and pebble circulation, agreeing well with experimental results.

5 Conclusion

DEM-CFD simulation based on two-grid approach can track the pebbles individually and have acceptable computational amount for tens of thousands of particles, which are two main advantages of the algorithm when one tries to simulate pebble flow and coolant flow in modular PB-FHR core with strong particle-fluid interactions. The proposed coupling method decouples the size restriction between fluid mesh and particle diameter, so the details of fluid field can be captured by the fine enough grid. With the proper mapping scheme between hexahedrons and polyhedrons, the shape of fluid grid is not limited to regular hexahedrons, which makes the algorithm suitable for irregular computational domain. Data exchange only occurs at the start of coupling procedure and the computation of both parts can overlap with each other, so the algorithm is more efficient than the traditional coupling scheme.

The complicated granular flow pattern in scaled facility can be predicted quantitatively. Actually, the flow field in reactor core is much smoother than the simulated flow field, which reduces the complexity of flow pattern significantly. It suggests that the proposed two-grid DEM-CFD algorithm can give correct granular dynamics for fluid-particle systems typically found in PB-FHR core.

References:

C.W. Forsberg, P.F. Peterson, P.S. Pickard et al.,

Molten-salt-cooled advanced high-temperature reactor for production of hydrogen and electricity

. Nucl. Technol. 144(3), 289-302 (2003). doi: 10.13182/nt03-1

Baidu Scholar

Google Scholar

J. Serp, M. Allibert, O. Beneš et al.,

The molten salt reactor (MSR) in generation IV: Overview and perspectives

. Prog. Nucl. Energ. 77, 308-319 (2014). doi: 10.1016/j.pnucene.2014.02.014

Baidu Scholar

Google Scholar

R. O. Scarlat, M. R. Laufer, E. D. Blandford et al.,

Design and licensing strategies for the fluoride-salt-cooled, high-temperature reactor (FHR) technology

. Prog. Nucl. Energ. 77, 406-420 (2014). doi: 10.1016/j.pnucene.2014.07.002

Baidu Scholar

Google Scholar

Z.J. Yang, J.L. Gou, J.Q. Shan et al.,

Analysis of SBLOCA on CPR1000 with a passive system

. Nucl. Sci. Tech. 28, (2016). doi: 10.1007/s41365-016-0154-y

Baidu Scholar

Google Scholar

S. Bu, J. Yang, M. Zhou et al.,

On contact point modifications for forced convective heat transfer analysis in a structured packed bed of spheres

. Nucl. Eng. Des. 270, 21-33 (2014). doi: 10.1016/j.nucengdes.2014.01.001

Baidu Scholar

Google Scholar

H. Zhang, Y. Tan, S. Shu et al.,

Numerical investigation on the role of discrete element method in combined LBM-IBM-DEM modeling

. Comput. Fluids. 94, 37-48 (2014). doi: 10.1016/j.compfluid.2014.01.032

Baidu Scholar

Google Scholar

H. Zhang, A. Yu, W. Zhong et al.,

A combined TLBM-IBM-DEM scheme for simulating isothermal particulate flow in fluid

. Int. J. Heat. Mass. Tran. 91, 178-189 (2015). doi: 10.1016/j.ijheatmasstransfer.2015.07.119

Baidu Scholar

Google Scholar

H. Kruggel-Emden, B. Kravets, M. Suryanarayana et al.,

Direct numerical simulation of coupled fluid flow and heat transfer for single particles and particle packings by a LBM-approach

. Powder. Technol. 294, 236-251 (2016). doi: 10.1016/j.powtec.2016.02.038

Baidu Scholar

Google Scholar

Z. Xia, Y. Fan, T. Wang et al.,

A TFM-KTGF jetting fluidized bed coal gasification model and its validations with data of a bench-scale gasifier

. Chem. Eng. Sci. 131, 12-21 (2015). doi: 10.1016/j.ces.2015.03.017

Baidu Scholar

Google Scholar

10.

L. Yu, J. Lu, X. Zhang et al.,

Numerical simulation of the bubbling fluidized bed coal gasification by the kinetic theory of granular flow (KTGF)

. Fuel. 86(5), 722-734 (2007). doi: 10.1016/j.fuel.2006.09.008

Baidu Scholar

Google Scholar

11.

H. Zhu, Z. Zhou, R. Yang et al.,

Discrete particle simulation of particulate systems: theoretical developments

. Chem. Eng. Sci. 62(13), 3378-3396 (2007). doi: 10.1016/j.ces.2006.12.089

Baidu Scholar

Google Scholar

12.

J. Vujić, R. M. Bergmann, R. Škoda et al.,

Small modular reactors: Simpler, safer, cheaper?

. Energy. 45(1), 288-295 (2012). doi: 10.1016/j.energy.2012.01.078

Baidu Scholar

Google Scholar

13.

Y. Li, W. Ji,

Effects of fluid-pebble interactions on mechanics in large-scale pebble-bed reactor cores

. Int. J. Multiphas. Flow. 73, 118-129 (2015). doi: 10.1016/j.ijmultiphaseflow.2015.03.004

Baidu Scholar

Google Scholar

14.

S. Deb, D. K. Tafti,

A novel two-grid formulation for fluid-particle systems using the discrete element method

. Powder. Technol. 246, 601-616 (2013). doi: 10.1016/j.powtec.2013.06.014

Baidu Scholar

Google Scholar

15.

F Alobaid, J Ströhle, B Epple et al.,

Extended CFD/DEM model for the simulation of circulating fluidized bed

. Adv. Powder. Technol. 24(1), 403-415 (2013). doi: 10.1016/j.apt.2012.09.003

Baidu Scholar

Google Scholar

16.

F Alobaid,

A particle-grid method for Euler-Lagrange approach

. Powder. Technol. 286, 342-360 (2015). doi: 10.1016/j.powtec.2015.08.019

Baidu Scholar

Google Scholar

17.

S Rickelt, F Sudbrock, S Wirtz et al.,

Coupled DEM/CFD simulation of heat transfer in a generic grate system agitated by bars

. Powder. Technol. 249, 360-372 (2013). doi: 10.1016/j.powtec.2013.08.043

Baidu Scholar

Google Scholar

18.

C. Goniva, C. Kloss, A. Hager et al.,

An open source CFD-DEM perspective

. Proceedings of OpenFOAM Workshop (2010).

Baidu Scholar

Google Scholar

19.

P. A. Cundall, O. D. Strack,

A discrete numerical model for granular assemblies

. Geotechnique. 29(1), 47-65 (1979). doi: 10.1680/geot.1979.29.1.47

Baidu Scholar

Google Scholar

20.

Z. Peng, E. Doroodchi, C. Luo et al.,

Influence of void fraction calculation on fidelity of CFD‐DEM simulation of gas‐solid bubbling fluidized beds

. Aiche. J. 60(6), 2000-2018 (2014). doi: 10.1002/aic.14421

Baidu Scholar

Google Scholar

21.

J. Su, Z. Gu, C. Chen et al.,

A two‐layer mesh method for discrete element simulation of gas‐particle systems with arbitrarily polyhedral mesh

. Int. J. Numer. Meth. Eng. 103(10), 759-780 (2015). doi: 10.1002/nme.4911

Baidu Scholar

Google Scholar

22.

A. Fabri, S. Pion,

CGAL: The computational geometry algorithms library

. in Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information systems. ACM (2009). doi: 10.1145/1653771.1653865

Baidu Scholar

Google Scholar

23.

S. Ergun,

Fluid flow through packed columns

. J. Ind. Eng. Chem. 41, 1179-1184 (1949). doi: 10.1021/ie50474a011

Baidu Scholar

Google Scholar

24.

A. Griveau.

Modeling and Transient Analysis for the Pebble Bed Advanced High Temperature Reactor~PB-AHTR, Master Thesis

, University of California, Berkeley (2007).

Baidu Scholar

Google Scholar