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

DEM-CFD simulation of modular PB-FHR core with two-grid method 1 2 first-author 1 1 1 corresp corresponding author:wangnaxiu@sinap.ac.cn corresponding author:wangnaxiu@sinap.ac.cn wangnaxiu@sinap.ac.cn Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China 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-FHR Pebble flow DEM-CFD Two-grid approach 1 Introduction 1 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 1 ]. Fluoride salts has low neutron cross section and high heat capacity, which means higher operating temperature and better safety margin [ 2 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 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 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 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 6 - 8 8 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 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 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 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 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 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 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 15 , 16 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 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 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 1 2.1. Particle motions 2 The discrete element method, a powerful numerical tool proposed by Cundall and Strack [ 19 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: <math display="block" id="M1"><mrow><msub><mi>m</mi><mi>i</mi></msub><mfrac><mrow><mi>d</mi><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>N</mi><mi>c</mi></msub></mrow></munderover><mrow><mrow><mo>(</mo><mrow><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>+</mo><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>t</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow></mstyle><mo>+</mo><msub><mi>m</mi><mi>i</mi></msub><mover accent="true"><mi>g</mi><mo stretchy="true">→</mo></mover></mrow></math> <math display="block" id="M2"><mrow><msub><mi>I</mi><mi>i</mi></msub><mfrac><mrow><mi>d</mi><msub><mrow><mover accent="true"><mi>w</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>N</mi><mi>c</mi></msub></mrow></munderover><mrow><mrow><mo>(</mo><mrow><msub><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub><mo>×</mo><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>t</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></math> where, v i , w i , m i , and I i are 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: <math display="block" id="M3"><mrow><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><mo>−</mo><msub><mi>k</mi><mi>n</mi></msub><msubsup><mi>δ</mi><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow><mrow><mn>3</mn><mo>/</mo><mn>2</mn></mrow></msubsup><mo>−</mo><msub><mi>η</mi><mi>n</mi></msub><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>⋅</mo><msub><mrow><mover accent="true"><mi>n</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow><msub><mrow><mover accent="true"><mi>n</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> <math display="block" id="M4"><mrow><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>t</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mo>−</mo><mi>min</mi><mrow><mo>(</mo><mrow><msub><mi>k</mi><mi>t</mi></msub><msub><mi>δ</mi><mrow><mi>t</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>+</mo><msub><mi>η</mi><mi>t</mi></msub><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>t</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub><mo>⋅</mo><msub><mrow><mover accent="true"><mi>t</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>,</mo><mi>μ</mi><mrow><mo>|</mo><mrow><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>|</mo></mrow></mrow><mo>)</mo></mrow><msub><mrow><mover accent="true"><mi>t</mi><mo stretchy="true">⇀</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> 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: <math display="block" id="M5"><mrow><msub><mi>k</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mi>E</mi></mrow><mrow><mn>3</mn><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msup><mi>σ</mi><mn>2</mn></msup><mo stretchy="false">)</mo></mrow></mfrac><msqrt><mrow><mfrac><mrow><msub><mi>R</mi><mi>i</mi></msub><msub><mi>R</mi><mi>j</mi></msub></mrow><mrow><msub><mi>R</mi><mi>i</mi></msub><mo>+</mo><msub><mi>R</mi><mi>j</mi></msub></mrow></mfrac></mrow></msqrt></mrow></math> <math display="block" id="M6"><mrow><msub><mi>k</mi><mi>t</mi></msub><mo>=</mo><msub><mi>k</mi><mi>n</mi></msub></mrow></math> <math display="block" id="M7"><mrow><msub><mi>η</mi><mi>n</mi></msub><mo>=</mo><msub><mi>η</mi><mi>t</mi></msub><mo>=</mo><mo>−</mo><mfrac><mrow><mn>2</mn><mi>ln</mi><mi>e</mi></mrow><mrow><msqrt><mrow><msup><mi>π</mi><mn>2</mn></msup><mo>+</mo><msup><mrow><mi>ln</mi></mrow><mn>2</mn></msup><mi>e</mi></mrow></msqrt></mrow></mfrac><msqrt><mrow><mi>m</mi><msub><mi>k</mi><mi>n</mi></msub></mrow></msqrt><msubsup><mi>δ</mi><mrow><mi>n</mi><mo>,</mo><mi>i</mi><mi>j</mi></mrow><mrow><mn>1</mn><mo>/</mo><mn>4</mn></mrow></msubsup></mrow></math> where, E , e and σ are Young’s modulus, restitution coefficient and Poisson coefficient, respectively. 2.2 Navier–Stokes equations for fluids 2 The fluid phase is modeled by the incompressible volume-averaged Navier–Stokes equations, which can be described as: <math display="block" id="M8"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>ε</mi><mi>f</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>ε</mi><mi>f</mi></msub><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn></mrow></math> <math display="block" id="M9"><mrow><mrow><mo>[</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>ε</mi><mi>f</mi></msub><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>ε</mi><mi>f</mi></msub><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>=</mo><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub><mo>∇</mo><mi>p</mi><mo>−</mo><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ε</mi><mi>f</mi></msub><mi>τ</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>ε</mi><mi>f</mi></msub><mover accent="true"><mi>g</mi><mo stretchy="true">→</mo></mover></mrow></math> 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 <math display="block" id="M10"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mi>τ</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>μ</mi><mi>f</mi></msub><mo>+</mo><msub><mi>μ</mi><mi>t</mi></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>j</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></math> where μ f is fluid viscosity and μ t is turbulence viscosity. F , the momentum exchange rate per volume between fluid and particles, is given by: <math display="block" id="M11"><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>V</mi></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mi>i</mi></msub><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow></mstyle></mrow></math> 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: <math display="block" id="M12"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>k</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>k</mi><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ε</mi><mi>f</mi></msub><mrow><mo>(</mo><mrow><msub><mi>μ</mi><mi>f</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>k</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>G</mi><mi>k</mi></msub><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>σ</mi></mrow></math> <math display="block" id="M13"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>σ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>σ</mi><mover accent="true"><mi>U</mi><mo stretchy="true">→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msub><mi>ε</mi><mi>f</mi></msub><mrow><mo>(</mo><mrow><msub><mi>μ</mi><mi>f</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>σ</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>σ</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>ε</mi><mi>f</mi></msub><mfrac><mi>σ</mi><mi>k</mi></mfrac><mrow><mo>(</mo><mrow><msub><mi>C</mi><mn>1</mn></msub><msub><mi>ρ</mi><mi>f</mi></msub><msub><mi>G</mi><mi>k</mi></msub><mo>+</mo><msub><mi>C</mi><mn>2</mn></msub><msub><mi>ρ</mi><mi>f</mi></msub><mi>σ</mi></mrow><mo>)</mo></mrow></mrow></math> where, G k is the generation of turbulence kinetic energy owing to the mean velocity gradients, C 1 and C 2 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 2 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 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: <math display="block" id="M14"><mrow><msub><mrow><mover accent="true"><mi>f</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>s</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>d</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>v</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> 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: <math display="block" id="M15"><mrow><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>s</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mo>−</mo><mi>V</mi><mmultiscripts><mo>∇</mo><mprescripts/><mi>i</mi><mi/></mmultiscripts><msub><mi>p</mi><mi>i</mi></msub><mtext>+</mtext><msub><mi>V</mi><mi>i</mi></msub><mo>∇</mo><mo>⋅</mo><mi>τ</mi></mrow></math> The drag force on particle i can be written as <math display="block" id="M16"><mrow><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mtext>d</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mrow><mo>[</mo><mrow><mfrac><mi>β</mi><mrow><mn>1</mn><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>−</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><msub><mi>V</mi><mi>i</mi></msub></mrow></math> 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: <math display="block" id="M17"><mrow><mi>β</mi><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mn>150</mn><mfrac><mrow><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mi>μ</mi></mrow><mrow><msub><mi>ε</mi><mi>f</mi></msub><msubsup><mi>d</mi><mi>p</mi><mn>2</mn></msubsup></mrow></mfrac><mo>+</mo><mn>1.75</mn><mfrac><mrow><msub><mi>ρ</mi><mi>f</mi></msub><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>d</mi><mi>p</mi></msub></mrow></mfrac><mrow><mo>|</mo><mrow><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>−</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mo>|</mo></mrow><mtext>            </mtext><msub><mi>ε</mi><mi>f</mi></msub><mo>≤</mo><mtext>0</mtext><mtext>.8</mtext></mtd></mtr><mtr><mtd><mfrac><mn>3</mn><mn>4</mn></mfrac><msub><mi>C</mi><mtext>D</mtext></msub><mfrac><mrow><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>ε</mi><mi>f</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>d</mi><mi>p</mi></msub></mrow></mfrac><mrow><mo>|</mo><mrow><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>−</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mo>|</mo></mrow><msub><mi>ε</mi><mi>f</mi></msub><msup><mrow/><mrow><mo>−</mo><mn>2.65</mn></mrow></msup><mtext>                         </mtext><msub><mi>ε</mi><mi>f</mi></msub><mtext>>0</mtext><mtext>.8</mtext></mtd></mtr></mtable></mrow></mrow></math> where d p is particle diameter, and C D is the drag coefficient for a single unhindered particle evaluated by <math display="block" id="M18"><mrow><msub><mi>C</mi><mtext>D</mtext></msub><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mn>24</mn><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mn>0.15</mn><mi>R</mi><msub><mi>e</mi><mi>p</mi></msub><msup><mrow/><mrow><mn>0.687</mn></mrow></msup><mo stretchy="false">)</mo><mo>/</mo><mi>R</mi><msub><mi>e</mi><mi>p</mi></msub><mtext>       </mtext><mi>R</mi><msub><mi>e</mi><mi>p</mi></msub><mo>≤</mo><mtext>1000</mtext></mtd></mtr><mtr><mtd><mn>0.44</mn><mtext>                                       </mtext><mi>R</mi><msub><mi>e</mi><mi>p</mi></msub><mtext>>1000</mtext></mtd></mtr></mtable></mrow></mrow></math> Re p , particle Reynolds number, can be written as: <math display="block" id="M19"><mrow><mi>R</mi><msub><mi>e</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>ε</mi><mi>f</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mrow><mo>|</mo><mrow><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>−</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mo>|</mo></mrow><msub><mi>d</mi><mi>p</mi></msub></mrow><mrow><msub><mi>μ</mi><mi>f</mi></msub></mrow></mfrac></mrow></math> The virtual mass force due to the pressure gradient on particle i can be written as: <math display="block" id="M20"><mrow><msub><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mi>v</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><msub><mi>C</mi><mrow><mtext>vm</mtext></mrow></msub><msub><mi>V</mi><mi>i</mi></msub><msub><mi>ρ</mi><mi>f</mi></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>⋅</mo><mo>∇</mo><mover accent="true"><mi>u</mi><mo stretchy="true">→</mo></mover><mo>−</mo><mfrac><mrow><mo>∂</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub><mo>⋅</mo><mo>∇</mo><msub><mrow><mover accent="true"><mi>v</mi><mo stretchy="true">→</mo></mover></mrow><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow></math> where C vm is the virtual mass effect coefficient. 2.4. Two-grid approach 2 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 . 2.4.1. Mapping scheme 3 In this article, mesh-to-mesh interpolation between hexahedron and polyhedron is based on the conservative interpolation [ 21 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 22 ]. A schematic illustration of the two-dimensional process is given in Fig. 2 Fig. 2 . The black circles represent particles. Suppose the quadrilateral A 0 A 4 B 6 B 5 is the particle cell i and the three polyhedrons j 1 , j 2 and j 3 are fluid cells. Take the cell j 2 as example, the value of cell j 2 is given by: <math display="block" id="M21"><mrow><msub><mi>ϕ</mi><mrow><mi>j</mi><mn>2</mn></mrow></msub><mo>=</mo><msub><mi>ϕ</mi><mi>i</mi></msub><mfrac><mrow><msub><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mn>2</mn></mrow></msub></mrow><mrow><msub><mi>S</mi><mi>i</mi></msub></mrow></mfrac></mrow></math> where, S i , j 2 is intersection area between cell j 2 and cell i . The intersection point of line of j 2 and cell i is A 3 and B 0 . The source and target cell vertices that do not lie outside each other, is A 7 , A 9 and A 4 . The five vertices can construct a convex hull with the area of S i and j 2 . The weight of fine cell j 2 on the coarse cell i is the ratio between the j 2 and S i . The value of cell i is given by: <math display="block" id="M22"><mrow><msub><mi>ϕ</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>ϕ</mi><mi>j</mi></msub><msub><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mi>k</mi></mrow></msub></mrow></mstyle></mrow><mrow><msub><mi>S</mi><mi>i</mi></msub></mrow></mfrac></mrow></math> where N is the total number of fluid cells that overlap with particle cell i . 2.4.2. Solid fraction 3 The ε c and Us c on DEM part are usually calculated based on divided particle volume method [ 20 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: <math display="block" id="M23"><mrow><msub><mi>ε</mi><mrow><mi>s</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>i</mi><mo>∈</mo><mi>c</mi><mi>e</mi><mi>l</mi><mi>l</mi><mtext> k</mtext></mrow></msub><mrow><msub><mi>α</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><msub><mi>V</mi><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mstyle></mrow><mrow><msub><mi>V</mi><mi>k</mi></msub></mrow></mfrac></mrow></math> 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 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 <math display="block" id="M24"><mrow><msub><mi>ε</mi><mrow><mi>s</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>i</mi><mo>∈</mo><mtext>cell </mtext><mi>k</mi></mrow></msub><mrow><msub><mi>α</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><msub><mi>V</mi><mrow><mi>p</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mstyle></mrow><mrow><msub><mi>V</mi><mi>k</mi></msub></mrow></mfrac><mo>/</mo><mfrac><mrow><msub><mi>V</mi><mrow><mtext>cut</mtext></mrow></msub></mrow><mrow><msub><mi>V</mi><mi>k</mi></msub></mrow></mfrac><mtext>              </mtext><msub><mi>V</mi><mrow><mtext>cut</mtext></mrow></msub><mo>></mo><mn>0</mn></mtd></mtr><mtr><mtd><mtext>0                                        </mtext><msub><mi>V</mi><mrow><mtext>cut</mtext></mrow></msub><mo>=</mo><mn>0</mn><mtext>        </mtext></mtd></mtr></mtable></mrow></mrow></math> where V cut is the intersection volume between particle cell k and fluid domain. 3 Numerical simulations of fluidized bed 1 The proposed two-grid DEM-CFD method was firstly validated against the test case of C. Goniva et al . [ 18 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 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 Table 1 . 1 1 Items 1 1 Value 1 1 Domain radius and height (m) 1 1 0.0138, 0.553 1 1 Gas Density (kg/m 3 ) 1 1 1 1 1 Gas viscosity(Pa·s) 1 1 1.5×10 -4 1 1 Gas superficial velocity (m/s) 1 1 0.002–0.02, in 0.002 steps 1 1 Fluid time step (s) 1 1 5×10 -4 1 1 Number of particles 1 1 10000 1 1 Particle diameter (m) 1 1 0.001 1 1 Particle density (kg/m 3 ) 1 1 200 1 1 Normal spring-stiffness (N/m) 1 1 1×10 4 1 1 Friction coefficient 1 1 0.3 1 1 Restitution coefficient 1 1 0.9 1 1 Solid time step (s) 1 1 1×10 -5 The three types of grids used in simulation are shown in Fig.4 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. The simulation results are shown in Fig. 5 Fig. 5 . The theoretical value of pressure drop for a fixed bed was evaluated by Ergun equation [ 23 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. 4 Experimental procedures and the results 1 4.1. Experimental procedures 2 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 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 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. 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 . 4.2. Model solution procedure 2 The symmetric section of the initial state of the facility is shown in Fig. 7 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 −5 s and 2×10 –3 s, respectively. The simulation conditions and parameters are summarized in Table 2 Table 2 . 1 1 Items 1 1 Values 1 1 Fluid Density (kg/m 3 ) 1 1 1000 1 1 Fluid viscosity (Pa s) 1 1 1.0×10 −3 1 1 Inlet velocity (m/s) 1 1 1.71, 1.95, 2.16, 2.33, 2.43, 2.59 1 1 Fluid time step (s) 1 1 2×10 −3 1 1 Number of particles 1 1 18000 1 1 Particle diameter (m) 1 1 0.015 1 1 Particle density (kg/m 3 ) 1 1 910.4 1 1 Normal spring-stiffness (N/m) 1 1 1×10 4 1 1 Friction coefficient 1 1 0.2 1 1 Restitution coefficient 1 1 0.9 1 1 Solid time step (s) 1 1 1×10 −5 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 Fig. 8 . Square grids represent coarse mesh and red points are the center of particles of initial packing. 4.3. Simulation results 2 4.3.1. Particle circulation indexes 3 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 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. 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 2 d (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 3 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 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 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. From Fig.11 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 Fig.9 . The higher the inlet superficial velocity, the larger the domain of vertex and pebble circulation, agreeing well with experimental results. 5 Conclusion 1 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 http://doi.org/10.13182/nt03-1 Molten-salt-cooled advanced high-temperature reactor for production of hydrogen and electricity 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 http://doi.org/10.1016/j.pnucene.2014.02.014 The molten salt reactor (MSR) in generation IV: Overview and perspectives 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 http://doi.org/10.1016/j.pnucene.2014.07.002 Design and licensing strategies for the fluoride-salt-cooled, high-temperature reactor (FHR) technology 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 http://doi.org/10.1007/s41365-016-0154-y Analysis of SBLOCA on CPR1000 with a passive system 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 http://doi.org/10.1016/j.nucengdes.2014.01.001 On contact point modifications for forced convective heat transfer analysis in a structured packed bed of spheres 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 http://doi.org/10.1016/j.compfluid.2014.01.032 Numerical investigation on the role of discrete element method in combined LBM-IBM-DEM modeling 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 http://doi.org/10.1016/j.ijheatmasstransfer.2015.07.119 A combined TLBM-IBM-DEM scheme for simulating isothermal particulate flow in fluid 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 http://doi.org/10.1016/j.powtec.2016.02.038 Direct numerical simulation of coupled fluid flow and heat transfer for single particles and particle packings by a LBM-approach 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 http://doi.org/10.1016/j.ces.2015.03.017 A TFM-KTGF jetting fluidized bed coal gasification model and its validations with data of a bench-scale gasifier 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 http://doi.org/10.1016/j.fuel.2006.09.008 Numerical simulation of the bubbling fluidized bed coal gasification by the kinetic theory of granular flow (KTGF) 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 http://doi.org/10.1016/j.ces.2006.12.089 Discrete particle simulation of particulate systems: theoretical developments 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 http://doi.org/10.1016/j.energy.2012.01.078 Small modular reactors: Simpler, safer, cheaper? 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 http://doi.org/10.1016/j.ijmultiphaseflow.2015.03.004 Effects of fluid-pebble interactions on mechanics in large-scale pebble-bed reactor cores 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 http://doi.org/10.1016/j.powtec.2013.06.014 A novel two-grid formulation for fluid-particle systems using the discrete element method 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 http://doi.org/10.1016/j.apt.2012.09.003 Extended CFD/DEM model for the simulation of circulating fluidized bed F Alobaid , A particle-grid method for Euler-Lagrange approach . Powder. Technol . 286 , 342 - 360 ( 2015 ). doi: 10.1016/j.powtec.2015.08.019 http://doi.org/10.1016/j.powtec.2015.08.019 A particle-grid method for Euler-Lagrange approach 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 http://doi.org/10.1016/j.powtec.2013.08.043 Coupled DEM/CFD simulation of heat transfer in a generic grate system agitated by bars C. Goniva , C. Kloss , A. Hager et al ., An open source CFD-DEM perspective . Proceedings of OpenFOAM Workshop ( 2010 ). An open source CFD-DEM perspective 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 http://doi.org/10.1680/geot.1979.29.1.47 A discrete numerical model for granular assemblies 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 http://doi.org/10.1002/aic.14421 Influence of void fraction calculation on fidelity of CFD‐DEM simulation of gas‐solid bubbling fluidized beds 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 http://doi.org/10.1002/nme.4911 A two‐layer mesh method for discrete element simulation of gas‐particle systems with arbitrarily polyhedral mesh 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 http://doi.org/10.1145/1653771.1653865 CGAL: The computational geometry algorithms library S. Ergun , Fluid flow through packed columns . J. Ind. Eng. Chem . 41 , 1179 - 1184 ( 1949 ). doi: 10.1021/ie50474a011  http://doi.org/10.1021/ie50474a011  Fluid flow through packed columns A. Griveau . Modeling and Transient Analysis for the Pebble Bed Advanced High Temperature Reactor~PB-AHTR, Master Thesis , University of California , Berkeley ( 2007 ). Modeling and Transient Analysis for the Pebble Bed Advanced High Temperature Reactor~PB-AHTR, Master Thesis 10.1007/s41365-017-0246-3 This work was supported by the "Strategic Priority Research Program" of the Chinese Academy of Sciences (No.XD02001002) 2016-10-31 2017-01-03 2017-01-08 2017-06-06 2017-07-01 2021-04-25 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0246-3 http://dx.doi.org/10.1007/s41365-017-0246-3 . 2017 Nuclear Science and Techniques 07 28