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

International Journal of Multiphase Flow 73 (2015) 118–129

Contents lists available at ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

Effects of fluid–pebble interactions on mechanics in large-scale pebble-bed reactor cores Yanheng Li, Wei Ji ⇑ Department of Mechanical, Aerospace, and Nuclear Engineering, Rensselaer Polytechnic Institute, 110 8th Street, Troy, NY 12180-3590, USA

a r t i c l e

i n f o

Article history: Received 29 October 2014 Received in revised form 3 March 2015 Accepted 8 March 2015 Available online 24 March 2015 Keywords: Pebble bed reactors Granular flow Fluid–particle system Fluid–pebble interaction Discrete Element Method Computational Fluid Dynamics

a b s t r a c t In Pebble-Bed Modular Reactor (PBMR) systems, typical representatives of large-scale fluid–particle energy systems, the gaseous coolant has strong interactions with the fuel pebbles that slowly move downward. To accurately predict the coolant dynamics and fuel bed mechanical properties, the fluid– pebble interactions should be taken into consideration in the mechanical modeling of the core. However, previous research usually neglected these interactions in the pebble flow analysis, assuming a fixed bed that is not affected by the coolant. The effects of these interactions on the core behavior remain unknown and have never been fully addressed. To investigate the interaction effect on the pebble flow mechanics, a coupled multi-physics model is employed which accounts for the fluid–pebble interactions in the dynamics simulations of a large annular nuclear reactor system. Discrete Element Method is used to simulate the pebble motion and predict important granular properties such as pebble packing density, coordination number and contact stress distributions. Computational Fluid Dynamics is employed to predict the locally averaged coolant properties. The simulations of pebble flow and coolant flow are coupled through the calculation and exchange of fluid–pebble interaction forces between two simulation solvers at a fixed time interval. To observe the interaction effect, decoupled simulations are also performed as references for comparison, i.e. a single pebble flow simulation without interactions with coolant, followed by a single coolant flow simulation in a static fixed packed bed. The comparison indicates that the interactions have significant effects on the local mean pebble contact force and observable effects on the axial and radial distributions of pebble volume packing fraction, as well as the pebble flow velocity/pressure profiles. The interaction effects are especially noticeable in the regions near the wall and close to the core inlet. The change in the pebble packing distribution will consequently bring in discrepancies in the distributions of axial and radial fluid velocity and pressure. Therefore, the interactions between coolant and pebble cannot be negligible in the study of the thermal–hydraulic properties of large-scale fluid–particle energy systems, although some macroscopic properties of the pebble flow (such as the coordination number distribution) are barely changed by the fluid–pebble interactions. Ó 2015 Elsevier Ltd. All rights reserved.

Introduction Packed beds of spherical components can be found in many applications to the energy and biological/chemical industry. Examples of these systems include trickling catalytic bed reactor (Propp et al., 2000), moving bed drying (Barrozo et al., 1998) and very high temperature gas-cooled pebble bed reactors (Gougar and Davis, 2006). In these systems, the particulate phase is densely packed, either fixed or slowly-moving. Strong fluid–particle interactions exist between the packed bed and the gas. Understanding these interactions is crucial to the interpretation of system behaviors. This is especially important when one studies the Pebble Bed ⇑ Corresponding author. E-mail address: [email protected] (W. Ji). http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.03.004 0301-9322/Ó 2015 Elsevier Ltd. All rights reserved.

Modular Reactor (PBMR400) system, which is a promising candidate for the next generation nuclear energy systems (Koster et al., 2003). PBMR400 design is a typical large-scale dense fluid– particle system: a total of about 440,000 tennis ball-sized fuel pebbles are densely packed and slowly circulating within the reactor core while high speed fluid coolant (mostly helium gas (Koster et al., 2003)) passes downward through the packed fuel bed. Therefore, a quasi-static, dense granular flow (pebble flow) and a continuum gas flow co-exist within the core. A clear understanding of the interaction between pebble flow and coolant flow, as well as its effect on the mechanical behaviors of both flows in such a largescale nuclear reactor system is crucial for the accurate prediction and analysis of the reactor thermal–hydraulic and neutronic performances.

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

In quasi-static dense granular flow like the pebble flow in the PBMR400, enduring contacts between particles are dominant over particle–particle collisions. The overall pebble flow moves very slowly with small shear rate throughout the core, which is similar to the property of solid (Forterre and Pouliquen, 2008). Based on this reason, previous research effort toward the modeling of pebble-bed reactors primarily used decoupled approaches (du Toit, 2008; du Toit et al., 2006; Gao and Shi, 2002; Gougar, 2006; Lee and Gardner, 2013; Reitsma et al., 2006; Rycroft et al., 2006). The decoupled approach assumes that solid pebble flow and the resultant static arrangement of solid pebbles are not affected by the fluid flow and the fluid–pebble interaction effect on granular properties is not considered. Only pebble–pebble or pebble–wall interactions are taken into consideration (Gougar, 2006; Rycroft et al., 2006). In the decoupled approach, a pebble (particle) flow simulation is first performed to obtain a packing arrangement of pebbles. Then the snapshot of pebble positions or pebble packing statistics is used as input to the fluid flow simulations for the prediction of fluid dynamics and energy transfer. However in large-scale pebble bed reactors, the interaction force between the pebble and fluid is usually significant and cannot be ignored when simulating pebble flows. For example, in the PBMR400 system, the mass flow rate of helium coolant is as large as 192 kg/s (Reitsma et al., 2006). The helium will exert significant force on pebbles, which can be more than one pebble’s weight and will have direct impact on the contact force among pebbles. Even for a very dense, solid-like granular flow, there are still zones behaving like fluid or gas (Forterre and Pouliquen, 2008) which will have a higher tendency to be affected by an external force such as the fluid force. From the rheological perspective of granular flow, the local packing fraction of granular flow is closely related to the local pebble flow shear rate, which is again correlated with the ratio between the volume-averaged local shear stress and the volume-averaged local normal stress on the granular assembly (Forterre and Pouliquen, 2008; Jop et al., 2006). For the large scale PBMR design, the fluid force is large compared with the granular stress of the pebble flow and is capable of altering the pebble flow shear rate and packing fraction distribution significantly. To accurately predict the system behavior, a full understanding of the fluid–pebble interaction effect and its degree on both flows is of importance, as the fluid flow dynamics is very sensitive to pebble packing fraction distribution (Bai et al., 2009). Investigation of this fluid effect has been lacking in previous research using decoupled approaches. Our recent research (Li and Ji, 2013b) has investigated the interaction effects on pebble flow and coolant flow in a small-scale pebble bed reactor system and showed certain degree of the effect on pebble motion and fluid dynamics. It is still unknown for a large-scale system, such as PBMR400, how significant these interactions can be to the system dynamic and thermal behaviors. Therefore, we aim to study the fluid–pebble interaction effect in a large-scale PBR system (PBMR400) in this paper, by simulating pebble and coolant flow using decoupled and coupled approaches. The first key effort is to find a multi-physics computational model with appropriate fidelity that can identify the difference between decoupled and coupled simulations. Many approaches have been proposed to simulate the fluid flow passing through a densely packed bed. The most widely used model is the two-fluid model (TFM). In the TFM, both particle and fluid flow phases are modeled based on the continuum theory. Particles are homogenized spatially using continuum theory such as kinetic model, and the fluid–particle system is treated as an interpenetrating mixture of two fluids (Artoni et al., 2011; Ishii, 1975; Lun, 1991). The coupling between the particle phase and fluid phase is realized by introducing an interaction term, which is usually obtained based on an empirical correlation (Ishii, 1975). Due to the fidelity loss in the continuum approximation of

119

particle phase (especially for densely packed particles), the exact location and velocity of each particle cannot be tracked, and some important phenomena, such as the near-wall ordering that is usually observed in fluid–particle systems, cannot be precisely accounted for by the two-fluid model (Rycroft et al., 2006). Based on these reasons, the TFM has limitations for both coupled fluid– particle modeling and decoupled single phase particle flow modeling in dense particle phase situations. These limitations restrict the applicability of the TFM for the study of the interaction effect in large-scale fluid–particle systems. On the other hand, high fidelity models based on the Direct Numerical Simulations (DNS) for both phases have been proposed (Hu et al., 2001). DNS does not adopt any empirical turbulence models or correlations that are assumed in TFM. The particle motion, fluid motion around each particle, and the fluid–particle interaction are explicitly simulated and fully resolved at micro-scale with a high degree of accuracy. Therefore, the DNS is highly capable of facilitating the study of the interaction effect in a fluid–particle system via both decoupled approach and coupled approach. However, the computational cost of DNS is extremely high especially at the condition of high Reynolds number in the fluid phase, which prevents DNS from modeling largescale systems such as the PBMR400 design. Current DNS-based simulations for fluid–particle systems can treat on the order of 103–104 particles (Hu et al., 2001; Tsuji, 2007; Wang et al., 2010) which is less than the particle number on the order of 105 in the PBMR400 design. As a compromise between efficiency and accuracy, many moderate-fidelity models are proposed for the coupled modeling of fluid–particle systems (Tsuji, 2007). The typical example is the coupled Discrete Element Method (DEM)-Computational Fluid Dynamics (CFD) approach (Zhu et al., 2007). In the DEM–CFD approach, the DEM (Cundall and Strack, 1979) is used to simulate the particle flow at the micro-scale level, incorporating all the essential particle contact physics, which enables the DEM model to track each pebble’s location and velocity, and to account for aforementioned phenomena that two-fluid model cannot capture. For the fluid phase, the continuum CFD approach is usually used to compute the ensemble-averaged fluid quantities at the meso-scale (Drew and Lahey, 1993). The meso-scale approach utilizes various cell- and time- averaged fluid–particle interaction models, removing the computational cost of explicitly resolving the detail of the flow fluctuation around each particle. With a good balance between fidelity and efficiency, the DEM–CFD approach is the most acceptable and pragmatic method that is used in analyzing large fluid–particle systems (Kafui et al., 2002; Zhu et al., 2007). By comparing DEM–CFD results with the DEM only results, it is very convenient to study the effect of the fluid–particle interaction. As for the PBMR400 system studied in this work, its large-scale computational domain requires good modeling efficiency, while the individual pebble’s motion and contacts should be resolved in order to obtain precise pebble distribution statistics and perform accurate pebble material strength analyses. Based on above reasons, the DEM–CFD approach is a good candidate for the full core analysis of large-scale PBR systems. The approach has a better efficiency than the DNS and a better accuracy than the TFM. In the DEM–CFD approach, the following two factors need to be carefully considered in order to obtain good computational fidelity: the turbulence effect at the high Reynolds number in fluid phase and the drag force model between solid and fluid phases. For the CFD simulations, in order to account for the turbulence effect at high Reynolds number, two general approaches have been proposed and used in simulations of the fluid dynamics in large-scale packed beds. One is to use an effective viscosity in the cell- and time-averaged fluid equations to account for the turbulence stress. The effective viscosity is determined based on either experiment measurements (Bey and Eigenberger, 1997) or solutions to the turbulent kinetic energy and dissipation (k–e)

120

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

transport equations (Jiang et al., 2002; Nakayama and Kuwahara, 2008; Pedras and de Lemos, 2001). The other approach is to incorporate contributions from turbulence stress and the drag force into a single 1D pressure drop model (a friction factor) without fully solving the momentum conservation equation (Gao and Shi, 2002; Park et al., 2010). This approach is simple but less accurate than the first approach. We employ the first approach with an empirical turbulence viscosity correlation to account for the turbulence effect in our study. The drag force between fluid and pebble dominates the fluid– pebble interaction effect and plays a key role in the overall DEM– CFD modeling accuracy. Various correlations have been proposed to seek an accurate model of the effective drag force in a densely packed bed. However, most available correlations are valid only under a limited range of physical conditions, e.g. low to intermediate Reynolds number (<1000) and/or volume packing fraction (less than 0.5), which are less than those in the PBMR400 application. For classical experiment-based drag model correlations, the Wen/Yu model (Wen and Yu, 1966) is valid for low packing fractions (<0.2) (Gidaspow, 1994; Zhu et al., 2007) and is used primarily in fluidized bed applications without specific restrictions on the particle Reynolds number, although some literature (Kandhai et al., 2003) claim that it is valid for packing fractions up to 0.6, close to the pebble volume packing fraction range in a PBR. For higher packing fractions (>0.2), the Ergun type model was used (Zhu et al., 2007) and it is in a quadratic form f p  rp  c1 u þ c2 juju (Ergun and Orning, 1949), where fp is the fluid–particle force density, p is the fluid pressure and u is the relative velocity between fluid and particle. The linear term accounts for the viscous effect and the quadratic term accounts for the inertial effect of the flow. The original Ergun model is a 1D pressure drop correlation validated by experiments in a segment of a packed bed at the conditions of the Reynolds number less than 1000 and a geometry aspect ratio greater than 50 (Dw/d > 50), where Dw is the hydraulic diameter of the container and d is the particle diameter. For the fluid–particle systems at the higher Reynolds number and smaller geometry aspect ratio, there is prominent near-wall channeling effect (flow bypasses through high porosity channels). The original Ergun model does not account for this effect and results in an overestimation of the pressure drop and the drag force (Mehta and Hawley, 1969), although Ergun-type models have been widely used in the flow through a densely packed bed. The problem for using the original Ergun equation in the PBR application is that for the PBMR400 case, the pebble Reynolds number is about 104 and the annular reactor pressure vessel wall has an aspect ratio less than 50. The higher Reynolds number and smaller aspect ratio lead to a considerable wall effect that the classical Ergun model does not fully account for. Because of this reason, the German Nuclear Safety Standards Commission (KTA) developed an alternative pressure drop model for PBR applications (KTA, 1981), which is a modification of the Ergun equation that extends the applicability of the original Ergun model to the high Reynolds number range (up to Re  105). The KTA model uses the quadratic term in the Ergun equation and introduces a variable coefficient that is a function of Reynolds number to replace the constant coefficients in the original Ergun model. By this modification, the channeling effect that makes the original Ergun model inaccurate at high Reynolds number regime is fixed. The valid packing fraction range for this model is demonstrated from 0.58 to 0.64 (KTA, 1981), which is the typical range in a PBR core. Same as the Ergun model, the KTA model was originally developed for PBRs with a large Dw/d ratio and the wall effect due to a small geometry aspect ratio is still not fully accounted for. A recent work has shown that the 1D pressure drop predicted from the KTA model has observable errors compared with the experimental data in a Heat Transfer Test Facility (HTTF) with a finite size annular

geometry (du Toit and Rousseau, 2012; Rousseau and van Staden, 2008). However, for the modeling of fluid–particle systems in 2D and 3D, the wall effect can be automatically accounted for in a viscous stress term. This allows to use the KTA model as an acceptable drag force model in the multi-dimensional modeling of fluid–particle systems (du Toit and Rousseau, 2012; Park et al., 2010). As suggested by Bey and Eigenberger (1997) and du Toit and Rousseau (2010, 2012), for multi-dimensional cases, we can extend the applicability of an Ergun-type pressure loss model that was developed for large perspective ratios (such as the KTA model) from a 1D system to a 2D or 3D system as a drag force model. The extended drag model can be combined with the effective viscous stress term which takes the molecular and turbulence stress into account in the fluid momentum balance equation to fully account for the wall effect mentioned earlier. du Toit and Rousseau have compared simulation predictions based on the KTA drag model plus effective viscous stress with the HTTF experiment measurements and good agreement was observed (du Toit and Rousseau, 2012). Based on this comparison, we select the KTA model as our drag model, and adopt the effective viscous stress to account for the turbulence stress. The effective viscous stress term can be obtained via available correlations, such as Nakayama and Kuwahara’s correlation which is well studied and an acceptable model for the fluid–particle system simulations (Guo et al., 2006; Nakayama and Kuwahara, 1999, 2008). In PBR designs, several important pebble flow and fluid flow quantities are crucial and of great interest to designers and operators, such as the pebble flow velocity distribution, pebble contact force, pebble flow packing fraction distribution, fluid pressure drop and fluid velocity distribution. The major work presented in this paper is to investigate the effect of pebble–fluid interactions in an annular large-scale PBR system based on the PBMR400 design. This is done by comparing these important pebble flow and coolant flow characteristics obtained from the coupled simulations with those obtained from traditional decoupled approach where fluid calculation is based on single phase pebble flow snapshots. The investigation shows that for a slow moving densely packed granular system such as PBMR400, the extra force exerted from the gaseous coolant has effects on the local quantities such as local average packing fraction and local mean normal contact force, especially near the upper wall region. The coolant force, however, hardly affects the pebble coordination number distribution. For the fluid velocity distribution and pressure drop, the interaction will result in noticeable differences between simulations based on the coupled and decoupled approaches and cannot be neglected. Coupled pebble flow–coolant flow modeling The pebble–coolant interaction, especially the interphase drag force, brings in the coupling between the pebble flow and coolant flow. The coupled model in the DEM–CFD approach is realized by adding a term for the fluid–pebble interaction in each flow’s dynamic equations. The equations of motion for the pebble flow are based on the Newton’s law of motion. The contact forces between pebbles are determined by the Hertzian contact mechanics theory (Johnson, 1985). The time- and cell-averaged coolant flow equations are based on fluid mass continuity and momentum conservation (time- and cell-averaged Navier–Stokes equations). The inter-phase drag force is provided based on experimentally established correlations. Meanwhile, to study the coupling effect, we also perform decoupled simulations, which only account for the inter-phase drag in the fluid phase calculation and do not include the drag term in the particle phase modeling. For coupled approach, since two phases are solved at different scales (micro-scale for particle phase and meso-scale for fluid phase), it is quite inefficient to solve the combined pebble–coolant

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

equations at the same time step. In practice, in order to resolve the contact force which is based on micro-scale inter-pebble overlap (Li and Ji, 2013b), the particle phase equation is solved at a much smaller time step than that for solving fluid phase equations. This implies that over a certain amount of particle phase time steps, the fluid phase quantities are constant. After each fluid phase time step, the updated particle quantities are passed to the fluid solver to solve for the new fluid quantity distribution. By this means the overall simulation efficiency is greatly enhanced without loss of fidelity (Li and Ji, 2013a). The exchange frequency of inter-phase interaction is determined by the relative time evolution rate between two phases. For PBMR400 applications, since the particle phase moves much slower than the fluid phase, the inter-phase information is updated after 1000 DEM time steps to balance between accuracy and efficiency (Li and Ji, 2013a). The equations of motion for the pebbles are:

m

N X dv i ¼ Fi ¼ F ij þ W i þ mi g þ F f ;i ; dt j–i

ð1Þ

j¼1

Ji

dxi ¼ Ti; dt

ð2Þ

where vi is the velocity of the ith pebble, Fi is the net force on the ith pebble including Fij the contact force from the jth pebble, Wi the wall contact force, mig the gravitational force, and Ff,i the fluid force which will be discussed in detail later. The pebble–pebble contact force Fij, is composed of normal contact force Fn,ij which is along the direction of nij (from the center of pebble j to the center of pebble i) and tangential contact force Ft,ij which is along the plane perpendicular to nij. Ti = rinij  Ft,ij is the torque on ith pebble due to the tangential components of contact forces, where ri is the radius of ith pebble. Ji is the moment of inertia, and xi is the angular velocity. According to the viscoelastic contact mechanics model (Rycroft et al., 2006; Shäfer et al., 1996), the normal contact force Fn,ij between two pebbles i and j with radius ri, rj and velocity vi, vj, is given by 0:5 F n;ij ¼ SðrÞðkn d1:5 ij  cn dij vn;ij Þnij ;

ð3Þ

where kn and cn are material elastic and viscoelastic constants respectively, S(r) is a shape function dependent on the shape of the pebble, and kn characterizes the normal contact stiffness of the pebble. vn,ij is the relative normal speed defined as (vi  vj)nij, and dij is the overlap size. According to (Shäfer et al., 1996), 2E for mono-sized pebbles, kn ¼ 43 Eeff ¼ 3ð1 m2 Þ and SðrÞ ¼ r eff ¼ ffi qffiffiffiffiffiffiffi p ffiffiffiffiffiffiffi ffi ri rj ¼ r=2, where E is the Young’s modulus and m is the r þr i

j

Poisson ratio. As for the viscoelastic constant cn, it is closely related to the restitution coefficient of the pebbles and there is no rigorous expression. In this work since we use a faster-than-reality recirculation, the value of cn is set to be high at 5000 kg/(m s) in order to dissipate the energy of the newly inserted pebbles and to avoid excessive bouncing. For the tangential contact force, the stick–slip model is used combined with the Coulomb dry friction law:

F t;ij ¼

8 Rt < kt d0:5 ij

t0

v t;ij dt; kF t;ij k 6 ls kF n;ij k

: ld kF n;ij k

vt;ij vt;ij

; kF t;ij k > ls kF n;ij k

;

ð4Þ

where kt is the tangential elastic constant, ls is the static friction coefficient and ld is the kinetic friction coefficient. vt,ij is the relative tangential velocity defined as v i  vj  ½ðvi  v j Þ  nij nij þ

r i xi ðdij nij Þþr j xj ðdij nij Þ , r i þrj

where dij is the distance between the center

of sphere i and the center of sphere j. Since it is very difficult to

121

identify the true value of ls and ld under operational condition within the helium environment, in this work ls is assumed to be equal to ld for simplicity (denoted as l). The wall-pebble contact force Wi follows similar models as used in the pebble–pebble contact. The fluid–pebble interaction force Ff,i represents the total force that is exerted on the pebble from the fluid. In the case of dense and slowly moving granular flow, the major interaction force from the gas coolant is the drag force FD which is usually obtained through empirical correlations. For the pebble flow modeling, it is noted by previous pebble flow works (Rycroft et al., 2006) that a much smaller-than-realistic kn (or equivalently, normal stiffness) can be chosen to speed up the simulation without loss of accuracy about the pebble distribution statistics. However, for the large-scale PBMR400 modeling, since the contact force at the bottom of the core is very large, a closeto-reality material stiffness (which can be characterized by the material Young’s modulus that can be as large as 10 GPa for graphite) should be used so that the inter-pebble overlap at the bottom core region would not be unrealistically high. Since the critical time step of DEM simulation (Dtc) has the relation of pffiffiffiffiffiffiffiffiffiffiffiffi Dtc  m=kn (Cundall and Strack, 1979), it is determined that a DEM time step of Dt  106 s should be adopted which is one order smaller than the time step used for small-scale pebble-bed reactor applications such as the HTR10 system modeling (Li and Ji, 2013b). For most pebble-bed reactor system designs, the reactor core is either cylindrical or annular. Moreover, the coolant and pebble flow are injected by forces along the axial direction and are confined by the forces along the radial direction from the inner and outer walls (du Toit et al., 2006). Such a highly axisymmetric geometry and dominant flow driven forces in the axial/radial direction would lead to a much smaller spatial variation in flow quantities along the tangential direction than that along the axial and radial directions. Therefore, it is common to average the fluid quantities over the tangential direction to obtain simplified 2D fluid equations. If we denote fp as the total fluid–pebble interaction force density exerted on the fluid by all the pebbles that fall into a fluid cell in three dimension, and hf p ih as the tangentially-averaged R 2p fluid–pebble interaction force density defined by hf p ih ¼ 21p 0 f p dh h iT ¼ hf p ih;r ; hf p ih;z , then the cell- and time-averaged 2D fluid equations in the radial and axial dimensions can be written as:

@ðeqf Þ 1 @ðeqf rur Þ @ðeqf uz Þ þ þ ¼ 0; @t @r @z r   @ðeur Þ 1 @ðerur ur Þ @ðeuz ur Þ qf þ þ @t r @r @z " #   @p 1 @ @ur @ 2 ur ur þ 2  2 þ hf p ih;r ; r ¼ e þ eleff @r r @r @r @z r

qf

  @ðeuz Þ 1 @ðeruz ur Þ @ðeuz uz Þ þ þ @t r @r @z " #   @p 1 @ @uz @ 2 uz þ 2 þ hf p ih;z þ eqf g; r ¼ e þ eleff @r @z @z r @r

ð5Þ

ð6Þ

ð7Þ

where e is the average porosity in a fluid cell, ur and uz are the radial and axial components of the fluid velocity vector u, g is the magnitude of gravitational acceleration, qf is the average fluid density in a fluid cell, leff is the effective fluid viscosity which accounts for the effect of turbulence, and p is the average fluid pressure in a fluid cell. In order to account for pebbles residing within multiple fluid cells, fp is calculated based on the particle-source-in-cell method using the following equation (Parry and Millet, 2010):

122

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

f p;k ¼ 

X

ai;k F f ;i =V f ;k ;

ð8Þ

i2fluid cell k

where ai,k is the volume fraction of the ith pebble that falls into the kth fluid cell, and Vf,k is the volume of the kth fluid cell. Eq. (8) correlates fp with Ff,i, and its accuracy determines the overall computation fidelity of the coupling DEM–CFD model (Parry and Millet, 2010). For pebble–coolant modeling, the major fluid–pebble interaction can be expressed as (Kafui et al., 2002)

F f ;i  F B;i þ F D;i þ F VM;i ;

ð9Þ

where F B;i is the buoyancy force, FD,i is the drag force on the ith pebble which is related with the pebble–fluid interaction by f

F D ¼  ep V s , with Vs as the volume of a single pebble. FVM,i is the virtual mass force on the ith pebble that can be written as:

F VM;i ¼ C VM V s qf ;i

  @ui @vi þ ui  r ui   v i  rv i ; @t @t

ð10Þ

where CVM is the virtual mass effect coefficient and ui represents the fluid velocity around ith pebble. In this work, the Ergun-based KTA model is adopted to calculate the fluid–pebble force density, which can be written as (Gao and Shi, 2002; KTA, 1981):

fp ¼ 

Wð1  eÞjuj d

qf u;

ð11Þ

 Re 1  Re 0:1 where W ¼ 160 1 þ 3 1 , 100 6 Re/(1  e) 6 105, 0.36 6 e e e qu d

6 0.42. The Reynolds number is defined as Re ¼ fl 0 , and u0 is the f magnitude of local superficial velocity of the coolant. We can see that, in this model the constant coefficients in the original Ergun equation are replaced by a variable coefficient dependent on the Reynolds number. It should be noted that in principle, a relative velocity is needed to replace the velocity variable in Eq. (11) due to the motion of the pebbles. Because pebbles move much slower than the coolant, it is acceptable to use the coolant velocity in Eq. (11). The effective viscosity based on Nakayama and Kuwahara’s correlation (Nakayama and Kuwahara, 2008) is adopted in the simulations. This correlation has been used by other researchers for the large-scale fluid–particle systems and is demonstrated to be a reliable correlation for PBR-like system applications (Guo et al., 2006). The Nakayama and Kuwahara’s correlation is:

leff e ¼ 0:0233 Re; lf 1e

equilibrium state in a much less runtime. This accelerated recirculation is a common practice for large-scale granular flow modeling and will not change the in-core distribution of pebble packing fraction and normalized pebble flow velocity based on the previous research (Rycroft et al., 2006). The accelerated pebble flow inevitably results in a larger pebble flow shear rate near the outlet region and thus affects the packing fraction in the conic bottom. In this work, we focus on the coupling effect study in the bulk region (in-core region) of the annular core, which is barely affected by the pebble recirculation rate (Rycroft et al., 2006). The side view of the reactor is shown in Fig. 1, the core parameters are listed in Table 1, the physical parameters for the pebble are listed in Table 2, and the coolant properties are listed in Table 3, which are based on the realistic PBMR400 design (Reitsma et al., 2006). Since the thermal effect is not considered in this work, a constant temperature of T = 942 K is assumed in the core, which is the average PBMR400 coolant temperature based on an Idaho National Laboratory report (Gougar, 2006). In the simulations of pebble flow and coolant flow, an initial static packing of 440,000 pebbles are firstly generated within the core using the collective dynamics-based initialization algorithm (Li and Ji, 2012; Li and Ji, 2013c, 2014), with the initial average packing fraction of around 61% in the active core region (a core region that contains exactly all the fuel pebbles). Then the bottom outlet is opened to let the pebble fuel drop out. Once a pebble drains out of the core, it will re-enter the core from the top, forming a recirculation process. Meanwhile, helium gas coolant enters at high pressure from the top and leaves at the bottom, with the boundary conditions of constant velocity inlet and constant pressure outlet. Although the pebble recirculation is a dynamic process that exhibits stochastic nature, the pebble flow and coolant flow will reach a steady state or a dynamic equilibrium state after a long enough simulation runtime. The ‘‘dynamic equilibrium state’’ in this work is defined as a state where the maximum variation of

200 cm

370 cm

ð12Þ

where leff is the effective viscosity and lf is the molecular viscosity. The Reynolds number in the original correlation can reach up to 1E4 which is within the applicable range of PBMR400. Investigation of fluid–pebble interaction effect in the PBMR400 system In the PBMR400 design (Reitsma et al., 2006), the reactor core has an annular geometry with an outer diameter of Ro = 370 cm and an inner diameter of Ri = 200 cm, hence it has an crosssectional area of A = 7.6105 m2. The height of the active core is around 1100 cm in the original PBMR400 design which is equivalent to a total pebble number of around 440,000. In this work, the same geometry of PBMR400 core is adopted. Fuel pebbles are inserted at the top of the core. Since there is no readily available information regarding the outlet geometry and the pebble extraction mechanism, a 45° conic bottom is assumed in this work to let the pebble drain out of the core at an accelerated rate compared with the realistic one (180 days) in order to achieve the

H=1100 cm

Bulk region

=73cm z=0 Fig. 1. Geometry of PBMR400 system.

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129 Table 1 PBMR400 core design parameters (Reitsma et al., 2006). Thermal power

Inlet temperature (Ti)

Outlet temperature (To)

System pressure (hpi)

400 MWt

773 K

1173 K

9 MPa

123

steady state. The small stochastic fluctuations of the calculated results are smoothed out over 5 s. In order to study the significance of fluid–pebble interaction in a PBMR400 core, the simulation solutions based on the decoupled approach, which is adopted by other researchers, are also presented and compared with the results from the coupled approach. Interaction effect on pebble flow at steady state

Table 2 Pebble physical properties. Radius (r)

Mass (mp)

Young’s modulus (E)

Friction coefficient (l)

Poisson ratio (m)

3 cm

210 g

1E9 Pa

0.7

0.3

Table 3 Helium coolant physical properties at the temperature of 942 K. Density (q)

Viscosity(lf)

Specific heat (cp)

Mass flow rate _ (m)

Outlet pressure (po)

4.54 g/L

4.43E5 Pa s

5.19E3 J/kg K

192 kg/s

Constant

pebble packing fraction over consecutive time steps is less than a threshold (0.1%) within any fixed local fluid cell in the bulk core region and averaged over 5 s of time. The determination of running time for steady state and the time period for averaging is primarily physics-based and can be understood from the perspective of granular dynamics. Considering a small volume of granular material in the granular flow, two physical parameters can be defined in an average sense for such a small volume: volume-averaged shear rate c_ and volume-averaged normal stress (or granular pressure) ps (Forterre and Pouliquen, 2008). Related to these two parameters, two types of time parameters for the granular flow system can be defined. One is a microscopic time parameter d t micro ¼ pffiffiffiffiffiffiffiffi , characterizing the time needed for an individual parps =qs

ticle in a small volume assembly of particles to move a distance of one particle diameter d under the granular pressure ps, where qs is the mass density of the granular material. At the scale of tmicro, there are small fluctuations in granular parameters such as volume fraction within a small volume region. The other type of time parameter is the macroscopic time parameter t macro ¼ 1c_ , representing the time needed for a granular assembly (consisting of many particles) to undergo a unit deformation. For a given mesoscale representative volume element (RVE), tmacro can be regarded as the (upper bound) characteristic time needed to have a significant change in volume fraction and other granular quantity in RVE. For large scale PBR such as PBMR400, the packing is very dense and the overall shear rate of the pebble flow is very small, typically less than 102 s1, hence tmacro  100 s, while tmicro  102 s. In our simulations, we use tmacro to determine whether the system has reached the steady state (i.e., no significant change in volume fraction and other granular parameters). Therefore, the simulation time should be larger than tmacro. Once the pebble flow reaches the steady state, there will be only very small fluctuations characterized by a time scale of tmicro  102 s. If we further denote Tcyl as the average time for a pebble fuel to travel from the top to the bottom at the steady state, which is about 110 days of CPU time using a single CPU core in a Dell T7500 workstation, then the dynamic equilibrium state at the scale of fluid cells is reached within about one Tcyl from the initial configuration. The simulation stops at Tcyl and the final simulation solutions, such as pebble packing fraction distribution, contact stress distribution, and coolant pressure distribution, shown in this paper are the averaged solutions over a very short time period (5 s) at the end of the simulation time, when both the coolant flow and pebble flow reach the

A rough evaluation can be made to show the significance of the fluid drag force on pebbles. If we choose the coolant density at an average in-core condition (942 K and 9 MPa) and use the superficial coolant speed u0 = 5.56 m/s and the average bulk region porosity e0  0.40, according to established models (Eq. (9)) for fluid drag force exerted on packed beds, the helium drag force on a single pebble can be estimated as F D0  3:2 N, which is over 50% larger than the gravity force of a single pebble (which is 2.06 N). It justifies the necessity to study the pebble–coolant interaction effect in PBMR400 and other similar applications. Axial pebble packing fraction distribution comparison The axial distribution of pebble volume packing fraction (/r ðzÞ) is calculated by uniformly dividing the core into 6 cm-thick horizontal layers and summing up the volume of pebbles in each layer. The sphere cap volume formula (Lamarche and Leroy, 1990) is used in order to account for spheres that partially fall into a layer. Our simulations show that the coupled approach has a higher core height at the equilibrium state compared with the decoupled approach. However, the core region at the height above z = 980 cm is not fully occupied by pebbles, while for the decoupled simulation, the core region is fully occupied below z = 1095 cm. We can see from the axial packing fraction (which is the axial-packing fraction averaged over tangential and radial directions) distribution comparison (Fig. 2a) that the coupled and decoupled results are quite different from each other. Fig. 2b shows the difference in packing fraction between coupled and decoupled simulations. From Fig. 2a and b, we can see that the interaction effect on axial packing fraction is not noticeable in the bulk region. The decoupled result shows that /r ðzÞ is almost uniform across the bulk region, which is at 60.3%. As for the coupled approach, /r ðzÞ is near-uniform from z = 195 cm to z = 800 cm at a higher value around 62.1%, and when z > 800 cm /r ðzÞ gradually decreases as z increases. The reason is that at the upper core, the fluid–pebble interaction force is comparable with the pebble contact force. The radial distribution of the drag force near the boundary will cause significant shear stress on local pebble flow and will consequently reduce local pebble volume fraction. This shear-induced pebble flow dilation is dominant over the compression exerted from the coolant in regions close to the inlet. For lower z, the pebble contact force is much larger than the fluid–pebble interaction and the drag-induced shear stress is not significant compared with local pebble flow stress. Hence at lower height (still above the converging outlet region), pebble packing fraction is nearly uniform, which is the same as the decoupled case but at a slightly higher packing fraction due to the compression from the fluid. As a result of the loosely packed upper zone, the coupled approach results in a core height larger than the proposed 11 m by other researchers (Reitsma et al., 2006). It should be noted that the Young’s modulus is chosen on the order of 1E9 Pa for both models, which is one order smaller than the realistic value (1E10 Pa) but large enough such that the change in the pebble packing fraction from decoupled model to the coupled model is not due to the increase of the pebble overlap. For the decoupled situation, the total overlapped volume is about 60 cm3, which accounts for only 0.0001% of the total pebble volume. For coupled situation, this total overlapped volume will be increased to around 100 cm3 but still less

124

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

1

1

0.8

0.8 0.6

φz (r)

φr (z)

0.6 0.4 0.2 0

coupled decoupled 0

200

0.4 0.2

400

600

800

1000

coupled decoupled

0 100

z (cm)

120

160

180

(a) 150cm
10

1

5

0.8

0

0.6 0.4

-5

0.2 -10

140

r (cm)

φz (r)

Δφr (z) (x100)

(a) Comparison of axial distributions of packing fraction between coupled and decoupled models

0

200

400

600

800

coupled decoupled

1000 0 100

z (cm)

120

(b) The difference of axial distribution of packing fraction between coupled and decoupled models

Radial pebble packing fraction distribution comparison To obtain the radial distribution of packing fraction (/z ðrÞ), the core is uniformly divided into 160 annular zones. The packing profile is obtained by counting the accumulative volume of pebbles in each zone and averaging over the axial direction for a certain height using the sphere-cylinder intersection volume formula (Lamarche and Leroy, 1990). In Fig. 3, /z ðrÞ from both coupled and decoupled approach are given and compared at different heights. From the figure, we can see that the fluid–pebble interaction effect on /z ðrÞ is most prominent near the boundary. Throughout the bulk region, the effect is visible near the inner wall and becomes observable near the outer wall at the upper bulk region (z  800 cm). As discussed by other researchers (Forterre and Pouliquen, 2008; Jop et al., 2006), the local packing fraction of granular flow is closely related to the local pebble flow shear rate, which is again correlated with the ratio between the volume-averaged local shear stress and the volume-averaged local normal stress on the granular assembly (Li, 2014) (defined as granular friction coefficient). In other words, in order for the fluid force to change the pebble flow packing distribution, the fluidinduced shear stress must be significant compared with local volume-averaged assembly normal stress. It is obvious that the fluid-induced shear stress is most prominent near the wall, and the pebble flow normal stress is lower as the height goes higher. Therefore in the upper wall region the fluid effect is more noticeable. This difference in the packing fraction near the wall region will consequently bring in a difference in the fluid velocity and pressure profile, which will be discussed in the next section.

160

180

r (cm)

(b) 450cm
Fig. 2. Interaction effect on axial distribution of pebble volume packing fraction.

1 0.8

φz (r)

than 0.0002% of the total pebble volume. Therefore, the slightly increased overlap volume is negligible and is not responsible for the increase in the pebble packing fraction.

140

0.6 0.4 coupled

0.2

decoupled 0 100

120

140

160

180

r (cm)

(c) 750cm
Mean normal contact force distribution comparison The normal contact force is the dominant force within a largescale packed bed system. It is crucial to use the peak normal contact force (averaged over the local volume) as a criterion to select pebble material (graphite) with appropriate yield strengths, and the magnitude of the normal contact force is directly related with the graphite dust generation (Cogliati et al., 2011). The normal contact force between two identical spheres (Fn) and its consequent normal contact stress (rn) can be calculated as (Johnson, 1985):

8 2 > < Fn ¼ 3 > : rn ¼ p1

E ð1m2 Þ

h

pffiffiffiffiffiffiffiffi 1:5 r=2d i1=3 ; 2

6F n E 2 ð1m2 Þ r 2

ð13Þ

where parameters have the same physical meanings as those in Eq. (3) and thereafter. It should be noted that the viscoelastic component of the original normal contact force in Eq. (3) is much smaller

125

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

than the elastic component and is neglected for simplicity in Eq. (13). The normal contact stress should not be confused with the mean normal stress of the pebble assembly (Rycroft et al., 2013), which is the normal stress of the local pebble assembly homogenized within a finite volume and is usually much smaller than normal contact stress, although these two types of stresses are highly dependent on each other. Although simulations based on the DEM can accurately predict many granular flow physics such as near wall ordering and axial mixing, discrepancy exists between the predicted contact force by the DEM and the realistic measured normal contact force (Vidyapati et al., 2012). This discrepancy is mainly due to the soft sphere model adopted by the DEM. Because of this discrepancy, it is more meaningful to use the local mean normal contact force, instead of individual contact force between two contact pebbles, to analyze the criteria for pebble material design. The local mean normal contact force is defined as the average contact force of pebbles within a local region. Similar to the calculation of the packing fraction distribution, the PBMR400 core is divided into mesh cells in the 2D cylindrical coordinate system, and the mean normal contact force per pebble within each cell is calculated and the spatial distribution is shown in Fig. 4a for coupled approach and Fig. 4b for decoupled approach. The calculated normal contact force is normalized by pebble gravity (Fn⁄ = Fn/mg). Since the normal contact force is more significant at the lower core, we only give the distribution in the bottom core. From the figure we can see that, for the coupled simulation, the local normal contact force peaks near the converging outlet region close to the wall boundary, at a value of about 150 mg (300 N). For decoupled approach the Fn distribution is not sensitive to z and stays below 30 mg. This ‘‘high normal contact force’’ region for coupled approach is of high interest to the PBR system analyzers and is the region where material failure and dust generation likely occur. For the PBMR400 application, the calculated average normal contact stress in the ‘‘high

normal contact force’’ region from the coupled approach is within the order of 100 MPa, which reaches the yield strength of graphite material. According to (Cogliati et al., 2011), the local dust generation (wear) rate is approximately proportional to the mean local normal contact force. Hence this high normal contact force region is also a concern for pebble wear. For de-coupled approach, since the fluid-pebble force is not included in the pebble flow modeling, the local pebble wear rate cannot be obtained accurately, which makes the coupled approach indispensable in pebble wear analysis. Interaction effects on pebble coordination number distribution The coordination number of a packed bed (Nc) is a key parameter to measure the structure of a granular system (Egami and Waseda, 1984). There are various definitions for coordination number based on how the closest neighbors of a sphere are defined, and in this work it is referred to as the mechanical coordination number, which is the average number of contacted spheres with nonzero contact force (Song et al., 2008). The radial distribution of the coordination number is calculated in this study and is shown in Fig. 5. We compare the distribution of local coordination number obtained based on the decoupled (Fig. 5a) and coupled (Fig. 5b) simulations at different heights of the bulk region. From these two figures it can be seen that both models yield similar radial distribution of Nc, indicating that the local packing structure of the pebble flow is not significantly altered by the large fluid force exerted on the pebble flow despite the visible difference in packing fraction between two approaches. The reasons are (1) coordination number distribution is a physical property represented at larger scale compared with local packing fraction and is not that sensitive to local packing fraction distribution; (2) the fluid–pebble

6 5

Fn*

100

z=55cm z=65cm z=75cm z=100cm z=140cm

4

Nc

150

3 2 1

50

z=360cm z=480cm z=600cm z=720cm

0 100 0 100

z=240cm

120

140

160

180

r (cm) 120

140

160

180

(a) radial distribution of coordination number at different heights (decoupled)

r (cm)

(a) Spatial distribution of mean normal contact force in the bottom core (coupled)

6 5

150

Fn*

Nc

100

4

z=55cm z=65cm z=75cm z=100cm z=140cm

2

50

0 100

3

1 0 100 120

140

160

180

r (cm)

(b) Spatial distribution of mean normal contact force in the bottom core (decoupled) Fig. 4. Mean normal contact forces in the PBMR400 system.

z=240cm z=360cm z=480cm z=600cm z=720cm 120

140

160

180

r (cm)

(b) radial of coordination number at different heights (coupled) Fig. 5. Interaction effects on the probability distribution function of the coordination number.

126

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

interaction is still not significant enough to drastically change the pebble flow regime, and all the effect from the fluid–pebble interaction is mostly limited within the change in local properties of the pebble flow.

Interaction effect on coolant flow dynamics Total pressure drop comparison In order to obtain the pressure drop profile, the pressure values at a height z are averaged over all the radial fluid cells at the same height, denoted as hpðzÞir , and then the hpðzÞir  p0 profile is calculated and plotted in Fig. 7, where p0 is the constant outlet pressure and 0 < z < 1100 cm. From the comparison it can be concluded that the coupled approach and the decoupled approach produce similar total pressure drop across the core (2.35E5 Pa and 2.5E5 Pa, respectively), although they have different pressure gradient profile along the axial direction. The predicted average pressure loss within the adopted uniform-temperature core is slightly less than the

x 10

2

p-pout(Pa)

Interaction effects on pebble flow axial velocity distribution Finally we compare the radial distributions of pebble flow axial velocity between the coupled and decoupled models at different heights at the lower core, which are shown in Fig. 6. The local pebble flow axial velocity (vz) is normalized by the bulk region mean pebble flow axial velocity (the normalized velocity is denoted as vz⁄). Fig. 6 shows that in the bulk region both models yield similar z  0 throughout the inner core and vz distributions, where @v @r slightly smaller near the wall due to the wall friction. This behavior is very similar to the viscous fluid flow within a constrained pipe. The two models have a noticeable difference at the bottom core (near the cone–cylinder interface), where the result from coupled z approach has a larger radial gradient (smaller @v ), meaning a larger @r shear rate.

5

2.5

1.5 1 0.5 0

coupled decoupled

0

200

400

600

800

1000

z (cm) Fig. 7. Interaction effect on pressure drop.

reported value from other researchers, such as 2.6E5 Pa reported by Hossain et al. (2008) in which the thermal effect is considered. Coolant velocity distribution comparison The axial coolant velocity uz is the major component of the coolant velocity and plays an important role in the convective heat transfer within the core. Therefore, the accurate evaluation of radial and axial uz profiles is crucial for a good prediction of the in-core temperature distribution. In order to investigate the significance of fluid–pebble interaction effect on the uz, the radial profile at different heights in the lower core (z = 80 cm, 100 cm, 140 cm, 180 cm and 240 cm) are calculated and plotted in Fig. 8, where each data point represents a value of uz in a fluid cell. From Fig. 8, it can be seen that, as the fluid force tends to reduce the packing fraction near the outer boundary of the core, the coupled approach shows a more prominent channeling effect near the outer wall, where significant amount of fluid passes via the low 1800

3

vz*

2 1.5

1400 1200

uz (cm/s)

2.5

1600 z=75cm z=85cm z=100cm z=125cm z=180cm

1

800 600

z=80cm z=100cm z=140cm z=180cm z=240cm

400

0.5 0 100

1000

200 120

140

160

0 100

180

110

120

130

170

180

1600 1400

uz (cm/s)

vz*

z=75cm z=85cm z=100cm z=125cm z=180cm

1.5

1200 1000 800 600

1

400

0.5 0 100

160

1800

3

2

150

(a) u z radial profile for coupled modeling

(a) coupled results

2.5

140

r (cm)

r (cm)

200

120

140

160

180

r (cm)

(b) decoupled results Fig. 6. Interaction effect on radial distribution of pebble axial velocity.

0 100

z=80cm z=100cm z=140cm z=180cm z=240cm

110

120

130

140

150

160

170

180

r (cm)

(b) uz radial profile for decoupled modeling Fig. 8. Interaction effect on radial distribution of uz at different heights.

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

porosity radial channels. This channeling effect can be observed in both coupled approach and decoupled approach. However decoupled approach shows a balanced inner and outer channel while for coupled approach the high porosity near the outer-wall region caused by the fluid force makes the outer channel much more prominent than the inner channel. Since the axial component of the fluid flow is the primary factor that takes away the fission heat, this unique unbalanced channeling effect will cause a significant difference in thermal calculations between the two approaches, which calls for utilizing the coupled approach for a more accurate thermal modeling of the PBMR400 system.

127

For the interaction effect on the axial distribution of uz, we first _ as the coolant mass flow rate and A as the cross sectional define m area of the core. For the axial distribution of huz ðzÞir , from coolant _ mass conservation huz ðzÞir ¼ m=ðA  heðzÞir Þ, we can know that it is uniquely determined by heðzÞir , the radially averaged porosity at the height z, which has been shown in Fig. 2a. As for the radial coolant velocity (ur), due to the highly symmetric geometry of PBMR400, it is much smaller than the axial coolant velocity for both coupled and decoupled models, with the maximum magnitude of less than 1 m/s within the cylinder core. At this level of radial velocity, the difference from both approaches will not cause significant effects on the pebble flow or fluid flow properties.

1.8 10

1.6

20

1.4 1.2

(H-z)/dH

30 1 40 0.8 50

0.6

60

0.4 0.2

70 8

6

4

2

r/dR

(a) coupled model 1.8 1.6

10

1.4 20 1.2

(H-z)/dH

30 1 40

0.8

50

0.6

60

0.4 0.2

70 8

6

4

Drag force distribution comparison The distributions of the drag force exerted on the pebble for both the decoupled and coupled models are presented in Fig. 9. From the figure we can see that, the drag force in the PBMR400 case is very large, accounting for 150% of the pebble gravity in the bulk region, which significantly increases the inter-pebble contact force. The consequent issues, such as pebble jamming and pebble wear, due to the increased pebble contact force should be thoroughly investigated with the consideration of the larger fluid force onto the pebble flow. Compared with coupled approach, the decoupled model yields a more uniform drag force distribution along the axial direction, which is related to its uniform packing fraction distribution across the bulk region. Since the order of magnitude of the mean normal stress of the pebble flow can be approximated using the hydrostatic approach, we can find it on the order of 1E4  1E5 Pa. On the other hand, the turbulence shear stress which is determined primarily by the radial gradient of the axial fluid velocity, can be found several orders of magnitude smaller than the local mean volume-averaged normal stress of the pebble flow (0.01–0.1%). Therefore the turbulence shear stress, although much more prominent than smaller Reynolds number fluid–particle systems like the small scale HTR10, can hardly affect the pebble flow dynamics (especially in the lower core where the contact force is larger) and is not the primary factor responsible for the fluid-pebble interaction effect. Instead, the fluid–pebble interaction effect mostly comes from the radial gradient of the fluid drag force which is much more significant than the turbulence shear stress and more comparable with the mean normal stress of the pebble flow. From Fig. 9a we can see that the drag force shows a more prominent radial gradient in regions near the inner and outer walls throughout the core. A higher radial gradient value of the drag force would introduce a larger shear stress to the pebble flow and can appreciably change its rheological behavior if the pebble flow has a relatively small mean volumeaveraged normal stress. For the upper-core (z > 700 cm) region and the region very close to the outlet (73 cm < z < 129 cm), the local mean volume-averaged normal stress of the pebble flow is lower than that elsewhere. The fluid drag force mostly affects the pebble flow behavior in these regions, which explains the differences of the predictions shown in previous sections between the coupled and decoupled approaches.

2

r/dR

(b) decoupled model Fig. 9. Fluid drag force (FD) distribution within the cylinder core (normalized by pebble gravity mg). In both figures, the active core region (primarily cylindrical region) are uniformly divided into 7 and 74 meshes along r and z directions respectively. The r-axis denotes the edge numbers of the radial meshes starting from 1 (corresponding to inner wall at r = 100 cm) to 8 (corresponding to the outer wall at r = 185 cm). The z-axis denotes the edge numbers of the axial meshes starting from 1 (corresponding to z = 1109 cm) to 75 (corresponding to z = 73 cm).

Conclusion In this work, the pebble–coolant interaction effects on pebble motions and coolant dynamics are investigated in a large-scale PBMR400 system. Simple calculations have shown that the average pebble–coolant interaction force is about 50% larger than the pebble gravity force, with an even larger value at regions where the pebble volume packing fraction is higher than the average value. It is such a large interaction force (primarily the drag force of the

128

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129

coolant on the pebbles) that leads to some noticeable changes in the system behaviors. For the radially-averaged pebble packing fraction, the interaction causes a non-uniform distribution along the axial direction with significant differences near the inlet and outlet regions compared with the decoupled results. For the radial distribution of the axially-averaged pebble packing fraction, the coolant-induced shear stress on the pebble flow generates a more obvious difference in packing fraction for the coupled approach near the upper and outer wall region. This result demonstrates that a coupled approach is necessary in order to obtain more accurate thermal-hydraulic and neutronic predictions in the near wall region. The difference in the pebble packing fraction for different models consequently introduces discrepancies to the fluid calculations for both axial and radial fluid velocity and pressure distributions based on fluid mass and momentum conservation. The radial distribution of uz at different height is obtained and the result indicates that, compared with the balanced channeling effect from the decoupled approach, the fluid–pebble interaction brings in an unbalanced channeling effect which is more prominent near the outer wall. The difference between two models also exists in the prediction of the mean normal pebble contact force distribution. The coupled approach reveals a more noticeable high normal contact force area near the converging conic region, due to the strong fluid–pebble interaction force in the PBMR400 system. This high normal contact force region raises concerns about the graphite material strength and the dust generation issue in the PBMR400 design. Despite these noticeable effects in local pebble and fluid quantities, the interaction has little effect on the coordination number distribution, which indicates that the overall packing structure of the packed bed is not significantly changed by the strong fluid– pebble interaction. In addition, from the pebble velocity distribution, we can conclude that the fluid–pebble interaction does not have substantial impact on the overall movement pattern of the pebbles within the bulk region, although there is difference near the converging conic region between two models. Although a uniform-temperature core is adopted and the pebble and coolant temperature distributions are not considered in this work, the presented analyses still reveal many important phenomena in a large-scale PBR system and most of the results (especially the pebble flow results) are expected to hold true with the thermal effect considered based on the granular flow physics. Future work will address the thermal effect with quantitative analysis. Another future work involves with parallelizing the computational code to further enhance the efficiency of the DEM–CFD modeling for a large scale PBR. In the coupled DEM–CFD simulations, the total computation time consists of three parts: (1) particle dynamics evaluation, (2) porosity evaluation, and (3) fluid dynamics evaluation. Each of the parts can be parallelized. For particle dynamics evaluation, since multi-pebble contact exists ubiquitously in the core, shared memory based multi-threading (such as OpenMP) is preferable compared with message passing interface (MPI) based parallelization. For porosity evaluation, MPI is a better choice as there is no message passing between neighboring parallelization blocks. For the CFD calculation, both shared memory and MPI based parallel mechanisms are applicable. Our current implementation is to use OpenMP-based multi-threading in particle and fluid solvers. A speedup of 4–5 has been observed on a 12-core machine in the preliminary implementation and test.

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ijmultiphaseflow. 2015.03.004.

References Artoni, R., Santomaso, A., Canu, P., 2011. Coupling between countercurrent gas and solid flows in a moving granular bed: the role of shear bands at the walls. Int. J. Multiph. Flow 37, 1209–1218. Bai, H., Theuerkauf, J., Gillis, P.A., Witt, P.M., 2009. A coupled DEM and CFD simulation of flow field and pressure drop in fixed bed reactor with randomly packed catalyst particles. Ind. Eng. Chem. Res. 48, 4060–4074. Barrozo, M.A.S., Murata, V.V., Costa, S.M., 1998. The drying of soybean seeds in countercurrent and concurrent moving bed dryers. Drying Technol. 16, 2033– 2047. Bey, O., Eigenberger, G., 1997. Fluid flow through catalyst filled tubes. Chem. Eng. Sci. 52, 1365–1376. Cogliati, J.J., Ougouag, A.M., Ortensi, J., 2011. Survey of dust production in pebble bed reactor cores. Nucl. Eng. Des. 241, 2364–2369. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Drew, D.A., Lahey, R.T.J., 1993. Analytical Modeling of Multiphase Flow. Butterworth-Heinemann, Boston, MA. du Toit, C.G., 2008. Radial variation in porosity in annular packed beds. Nucl. Eng. Des. 238, 3073–3079. du Toit, C.G., Rousseau, P.G., 2010. The flow and heat transfer in a packed bed high temperature gas-cooled reactor. In: 2010 14th International Heat Transfer Conference. American Society of Mechanical Engineers, pp. 193–207. du Toit, C.G., Rousseau, P.G., 2012. Modeling the flow and heat transfer in a packed bed high temperature gas-cooled reactor in the context of a systems CFD approach. J. Heat Transfer 134, 031015-1–031015-12. du Toit, C.G., Rousseau, P.G., Greyvenstein, G.P., Landman, W.A., 2006. A systems CFD model of a packed bed high temperature gas-cooled nuclear reactor. Int. J. Therm. Sci. 45, 70–85. Egami, T., Waseda, Y., 1984. Atomic size effect on the formability of metallic glasses. J. Non-Cryst. Solids 64, 113–134. Ergun, S., Orning, A.A., 1949. Fluid flow through randomly packed columns and fluidized beds. Ind. Eng. Chem. 41, 1179–1184. Forterre, Y., Pouliquen, O., 2008. Flows of dense granular media. Annu. Rev. Fluid Mech. 40, 1–24. Gao, Z., Shi, L., 2002. Thermal hydraulic calculation of the HTR-10 for the initial and equilibrium core. Nucl. Eng. Des. 218, 51–64. Gidaspow, D., 1994. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions. Academic Press. Gougar, H.D., 2006. The Application of the PEBBED Code Suite to the PBMR-400 Coupled Code Benchmark-FY 2006 Annual Report. Gougar, H.D., Davis, C.B., 2006. Reactor Pressure Vessel Temperature Analysis for Prismatic and Pebble-bed VHTR Designs. INL/EXE-06-11057. Guo, B., Yu, A., Wright, B., Zulli, P., 2006. Simulation of turbulent flow in a packed bed. Chem. Eng. Technol. 29, 596–603. Hossain, K., Buck, M., Said, N.B., Bernnat, W., Lohnert, G., 2008. Development of a fast 3D thermal-hydraulic tool for design and safety studies for HTRS. Nucl. Eng. Des. 238, 2976–2984. Hu, H.H., Patankar, N.A., Zhu, M.Y., 2001. Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian–Eulerian technique. J. Comput. Phys. 169, 427–462. Ishii, M., 1975. Thermo-Fluid Dynamic Theory of Two-Phase Flow. Eyrolles, Paris, France. Jiang, Y., Khadilkar, M.R., Al-Dahhan, M.H., Dudukovic, M.P., 2002. CFD of multiphase flow in packed-bed reactors: I. k-Fluid modeling issues. AIChE J. 48, 701–715. Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press, Cambridge. Jop, P., Forterre, Y., Pouliquen, O., 2006. A constitutive law for dense granular flows. Nature 441, 727–730. Kafui, K.D., Thornton, C., Adams, M.J., 2002. Discrete particle-continuum fluid modeling of gas-solid fluidized beds. Chem. Eng. Sci. 57, 2395–2410. Kandhai, D., Derksen, J.J., Van den Akker, H.E.A., 2003. Interphase drag coefficients in gas–solid flows. AIChE J. 49, 1060–1065. Koster, A., Matzner, H.D., Nicholsi, D.R., 2003. PBMR design for the future. Nucl. Eng. Des. 222, 231–245. KTA, 1981. Reactor Core Design of High-Temperature Gas-Cooled Reactors Part 3: Loss of Pressure Through Friction in Pebble Bed Cores. KTA rule 3102.3, Kerntechischen Ausschuss. Lamarche, F., Leroy, C., 1990. Evaluation of the volume of intersection of a sphere with a cylinder by elliptic integrals. Comput. Phys. Commun. 59, 359–369. Lee, K.O., Gardner, R.P., 2013. Prediction of pebble motion in pebble-bed reactors using Monte Carlo molecular dynamics simulation. Nucl. Sci. Eng. 174, 264– 285. Li, Y., 2014. Multiscale Granular-Fluid Dynamics in Pebble Bed Reactor Cores, Ph.D. Thesis, Department of Mechanical, Aerospace, and Nuclear Engineering. Rensselaer Polytechnic Institute, Troy, New York. Li, Y., Ji, W., 2012. A collective dynamics-based method for initial pebble packing in pebble flow simulations. Nucl. Eng. Des. 250, 229–236. Li, Y., Ji, W., 2013a. Acceleration of coupled granular flow and fluid flow simulations in pebble bed energy systems. Nucl. Eng. Des. 258, 275–283. Li, Y., Ji, W., 2013b. Pebble flow and coolant flow analysis based on a fully coupled multi-physics model. Nucl. Sci. Eng. 173, 150–162. Li, Y., Ji, W., 2013c. Stability and convergence analysis of a dynamics-based collective method for random sphere packing. J. Comput. Phys. 250, 373– 387.

Y. Li, W. Ji / International Journal of Multiphase Flow 73 (2015) 118–129 Li, Y., Ji, W., 2014. Convergence analysis of a heuristic collective sphere packing algorithm. Int. J. Appl. Nonlinear Sci. 1, 136–155. Lun, C.K.K., 1991. Kinetic theory for granular flow of dense slightly inelastic, slightly rough spheres. J. Fluid Mech. 233, 539–559. Mehta, D., Hawley, M.C., 1969. Wall effect in packed columns. Ind. Eng. Chem. Proc. Des. Dev. 8, 280–282. Nakayama, A., Kuwahara, F., 1999. A macroscopic turbulence model for flow in a porous medium. J. Fluids Eng. 121, 427–433. Nakayama, A., Kuwahara, F., 2008. A general macroscopic turbulence model for flows in packed beds, channels, pipes, and rod bundles. J. Fluids Eng. 130, 101205-1–101205-7. Park, H., Knoll, D.A., Gaston, D.R., Martineau, R.C., 2010. Tightly coupled multiphysics algorithms for pebble bed reactors. Nucl. Sci. Eng. 166, 118–133. Parry, A.J., Millet, O., 2010. Modeling blockage of particles in conduit constrictions: dense granular-suspension flow. J. Fluids Eng. 132, 011302-1–011302-10. Pedras, M.H.J., de Lemos, M.J.S., 2001. Macroscopic turbulence modeling for incompressible flow through undeformable porous media. Int. J. Heat Mass Transf. 44, 1081–1093. Propp, R.M., Colella, P., Crutchfield, W.Y., Day, M.S., 2000. A numerical model for trickle bed reactors. J. Comput. Phys. 165, 311–333. Reitsma, F., Strydom, G., de Haas, J.B.M., Ivanov, K., Tyobeka, B., Mphahlele, R., Downar, T.J., Seker, V., Gougar, H.D., Da Cruz, D.F., 2006. The PBMR steady-state and coupled kinetics core thermal-hydraulics benchmark test problems. Nucl. Eng. Des. 236, 657–668.

129

Rousseau, P.G., van Staden, M., 2008. Introduction to the PBMR heat transfer test facility. Nucl. Eng. Des. 238, 3060–3072. Rycroft, C.H., Grest, G.S., Landry, J.W., Bazant, M.Z., 2006. Analysis of granular flow in A pebble bed nuclear reactor. Phys. Rev. E 74, 021306. Rycroft, C.H., Dehbi, A., Lind, T., Güntaye, S., 2013. Granular flow in pebble-bed nuclear reactors: Scaling, dust generation, and stress. Nucl. Eng. Des. 265, 69– 84. Shäfer, J., Dippel, S., Wolf, D.E., 1996. Force schemes in simulations of granular materials. J. Phys. I 6, 5–20. Song, C., Wang, P., Makse, H.A., 2008. A phase diagram for jammed matter. Nature 453, 629–632. Tsuji, Y., 2007. Multi-scale modeling of dense phase gas-particle flow. Chem. Eng. Sci. 62, 3410–3418. Vidyapati, V., Kheiripour Langroudi, M., Sun, J., Sundaresan, S., Tardos, G.I., Subramaniam, S., 2012. Experimental and computational studies of dense granular flow: transition from quasi-static to intermediate regime in a Couette shear device. Powder Technol. 220, 7–14. Wang, L., Zhou, G., Wang, X., Xiong, Q., Ge, W., 2010. Direct numerical simulation of particle-fluid systems by combining time-driven hard-sphere model and lattice Boltzmann method. Particuology 8, 379–382. Wen, C.Y., Yu, Y.H., 1966. A generalized method for predicting the minimum fluidization velocity. AIChE J. 12, 610–612. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2007. Discrete particle simulation of particulate systems: theoretical developments. Chem. Eng. Sci. 62, 3378–3396.