Effect of pebble size and bed dimension on the distribution of voidages in pebble bed reactor

Abstract

The HTR-10 built at Tsinghua University is an advanced pebble bed reactor because of its inherent safety and economic efficiency. It is fundamental to explore the voidage of the pebble bed. The existing experimental bed is limited in depth and contains mono-size pebbles. The effects of pebble size and bed dimension of voidage distribution are still not well known. In this work, the discrete element method is used to simulate the static packing of pebbles of three sizes in 2D and 3D pebble beds under the same load. The effects of bed dimension and pebble size on voidage distribution are analyzed. The results are useful for better understanding of the voidage distribution of pebble bed reactor and the effects of bed dimension and particle size as well as the wall effects.

Keywords

High temperature gas-cooled reactor pebble flow discrete element method voidage bed dimension particle size

High-temperature gas-cooled reactor (HTGR) is a type of nuclear reactor which has inherent safety features that produce minimal environmental impact for effective power generation. Thus, it is an excellent candidate of next generation of advanced nuclear reactors. Due to its many excellent characteristics, such as modularity, relatively low cost, short construction period, and other advantages, it has attracted more and more attention. As nuclear power technology to mature, the economic competitiveness of nuclear power has become increasingly strong. Nuclear power compared to other power generation is a relatively clean energy. For example, compared with the thermal power plant, the exhaust of nuclear power plants does not contain sulfur or other harmful pollutants, and the emission of carbon dioxide is only 4.5%.¹ With the increasingly serious haze and lack of water resources in the north of China, the nuclear energy has the advantage of environmental protection. At the same time, on account of the inherent safeties of nuclear reactors, the process of promoting nuclear power has attracted the increasingly attention of developers and the public.

The HTGR core is a spherical bed composed of a large number of spherical particles: fuel element coated ceramic and graphite balls as a moderator, through the accumulation of fuel pebbles and graphite pebbles to form the reactor core. The researchers of Tsinghua University have carried out experimental research and theoretical analysis in the whole process of the pebbles’ loading and unloading in the helium to achieve continuous operation without stopping, while the United States, Germany, South Africa, and other countries also do the relevant experiments and theoretical research on the status of the core.^2–5 However, many mechanical problems on the characteristics of HTR-10 core are still to be studied.⁶

Tsinghua University established the two-dimensional visualization of the pebble bed experimental apparatus in 2005,^7,8 and a number of scientific researchers are published via discrete element method (DEM) as the theoretical basis of the experiment and simulation to verify each other.⁹ Experiments and simulation results show that the current HTR-10 core flow is uniform, and there is no retention area in the loading and unloading process, which indicates that the structure of core is reasonable.¹⁰

In general, the study of the flow and heat transfer properties of the entire pebble bed has a significance of basic guiding, and the distribution of voidage can affect the pressure drop and heat transfer between the pebble bed and the reactor wall. Therefore, the voidage of the pebble bed in HTR has become the focus of many current research projects. Benenati and Brosilow’s studies reveal that there is oscillating state in the distribution of radial voids of the bed and the trend of the voidage has an asymptotic value finally,¹¹ which is consistent with that observed by Bomboniet in the experiment.¹² Although this oscillation does not affect the average parameters of the reactor, it may have a significant effect on local flow properties and thermodynamic reaction characteristics. For example, studies by Schertz and Bischoff show that the turbulence of voids causes the coolant flow in the gas-cooled reactor to tend to flow near the wall. Du Toit has extensive studies on distribution of voidage via simulation and experimental results.^14–18 Previous majority studies on the voidage are based on the voidage of the radial or axial distribution.^{11,12,14–22} The void distribution is determined by the filling characteristics of the core. So far, there have been many arguments about radial variations in near-wall voidage, such as Benenati and Brosilow, who argued that the radial voidage distribution of the near-wall surface is fluctuating.^11–18 In contrast, Vortmeyer and Schuster were considered that void velocities were nonvolatile and exponential correlation functions.^19–22 Due to the accumulation of structures, such as stacking can affect the internal and near-wall region of the specific value of the voidage by a simple cube, face-centered cube or body-centered cube, and the permutation structure may further affect the flow and thermo-physical properties in the bed, such as effective thermal conductivity.

In addition, it should be noted that the voidage is also affected by the movement of the core in the reactor. Therefore, the voidage during operation should be time-dependent. In other words, due to the fundamental difference between the slow moving status and the fixed status, the results obtained from the fixed bed cannot be used directly to estimate the voidage during dynamic loading of the bed. However, even in fixed-bed studies, there is still a complete description of the voidage characteristics in the radial and axial directions. Therefore, it is necessary and important to study and analyze the characteristics of porosity in detail in fixed spheroids, including the influencing effect and mechanism of porosity distribution.

DEM is widely used to simulate the movement of discrete spheres. Many researchers have studied the two-dimensional and three-dimensional models via the DEM method.^23–26 In this paper, the software which is named OpenFOAM and the soft sphere model which is one of the models in the DEM are used to simulate the static accumulation of the three kinds of core pebbles under the same load under two-dimensional and three-dimensional conditions, and the influence of the size and dimension of the core pebbles was studied quantitatively to improve the voidage study in the fixed bed of HTR-10.

Numerical model

The basic principle of DEM is to regard the element as a separate unit. The use of dynamic relaxation iterative method is to calculate the dynamics of each unit and the kinematic information. Then, the force and motion of each unit body in each time step is computed as well as the position of the unit body. Therefore, the whole motion of the macroscopic system can be obtained by analyzing the motion state of each individual unit. During simulation, the interaction between the particles is simplified as instantaneous equilibrium states.

According to Newton’s second law, the complete control equation for DEM particle motion is divided into two parts: translation and rotation. In this model, considering the interaction between each core pebble in the reactor bed and the gravity, the expressions are as follows

m_{i} \frac{{dv}_{i}}{dt} = \sum_{j} F_{ij}^{c} + F_{i}^{g}

(1)

I_{i} \frac{{dw}_{i}}{dt} = \sum_{j} M_{ij}

(2)

In the equation (1), m_i and v_i represent the mass and velocity of the pebble of the serial number i. $F_{ij}^{c}$ and $F_{i}^{g}$ are the contact force between the pebbles and the gravity. In the equation (2), I_i is the moment of inertia of pebbles. w_i and $M_{ij}$ are the angular velocity of the core and the moment of contact force. In the nonlinear viscoelastic model proposed by Tsuji, the formula for the normal force of the contact force between the units is

{\binom{}{F}}^{n} = {\binom{}{F}}_{et}^{n} + {\binom{}{F}}_{diss}^{n} = - k^{n} ξ^{3 / 2} - γ^{n} \overset{\cdot}{ξ} ξ^{1 / 4}

(3)

In the equation (3), ${\binom{}{F}}_{et}^{n}$ and ${\binom{}{F}}_{diss}^{n}$ represent the elastic repulsive force and viscous dissipative force, respectively in the process of collision and contact, ξ and $\overset{\cdot}{ξ}$ respectively represent the overlap item and instantaneous velocity in normal direction when the collisions among unit balls happen. $k^{n}$ represents the stiffness coefficient of the spring vibrator in the simulation of elastic repulsive force while $γ^{n}$ represent the damping coefficient of the dampers in the simulation of viscous dissipative force. This model is one of the most mature and widely used normal contact force models in DEM. The value of $k^{n}$ can be derived from the Hertzian contact theory²⁸

k^{n} = \frac{4}{3} (\frac{R_{1} R_{2}}{R_{1} + R_{2}})^{1 / 2} (\frac{1 - λ_{1}^{2}}{E_{1}} + \frac{1 - λ_{2}^{2}}{E_{2}})^{- 1}

(4)

In the equation (4), R represents the radius of the collision unit ball, E is the Young’s modulus of the unit ball; λ is the Poisson's ratio of the unit ball. The model is based on the assumption that the surface of the collision unit ball is smooth. The contact surface is ellipsoid, and the friction of the contact surface can be ignored and only considering the small deformation. The same as the normal force, the nonlinear viscoelastic model is used for tangential contact force in this paper.²⁹

{\binom{}{F}}^{t} = {\binom{}{F}}_{et}^{t} + {\binom{}{F}}_{diss}^{t} = - k^{t} β - γ^{t} \overset{\cdot}{β}

(5)

In the equation (5), β and $\overset{\cdot}{β}$ , respectively represent the overlap item and velocity in the tangential direction in the process of collisions among unit balls. Similarly, the calculation method of tangential stiffness coefficient is consistent with the normal stiffness coefficient. When the tangential force is greater than the dynamic friction, slither of the unit ball and the wall is considered.³⁰

{\binom{}{F}}^{t} = μ {\binom{}{F}}^{n}

(6)

where μ is the coefficient of friction between the particles and the wall. The DEM model used in this paper has been validated by previous comparison³¹ between experiments and simulation.

Simulation conditions and function definitions

Simulation conditions

The computational model of this paper has the same geometrical dimensions as the actual HTR-10 experimental system. As shown in Figure 1, the gravity direction is −y. In this paper, it simulated three kinds of pebbles, namely d_p = 0.06, 0.09, 0.12 (m). In order to ensure the stability of the calculation, the simulation time step must be less than the contact time of the pebbles collision, which will lead to very much calculation time. In order to reduce the computational cost of the simulation, the inherent stiffness coefficient of the material can be increased the collision contact time. Studies have shown that the use of this method does not change the trend of the results of the calculation.¹⁶ And previous studies have shown that the boundary and the size of the pebbles are the main factor affecting the voidage. The reduced contact force does not change the boundary conditions and the size of pebbles, and only reduces the internal force evenly. It can be reasonably assumed that the tendency of the relative distribution characteristics of the voidage is not affected by the reasonable reduction in the stiffness coefficient. According to the previous study, the stiffness coefficient used in this paper is about 1% of the stiffness coefficient of the graphite material used in the actual pebbles. The simulation parameters are shown in Table 1 and the simulation parameters of pebbles are shown in Table 2.

Table 1.

Simulation parameters.

Parameters	Value
Diameter of pebble bed, D_bed (m)	1.8
Height of the bed, H (m)	1.97
Base cone angle, α (°)	30
Diameter of discharge (m)	0.5
Height of discharge (m)	0.5
Restitution coefficient, e	0.9
Stiffness factor, kⁿ (N/m)	1 × 10⁴
Friction coefficient, μ	0.3
Poisson rate	0.3
Time step (s)	1 × 10⁻⁴
Total simulated time (s)	20

Table 2.

Simulation parameters of pebbles.

	d_p = 0.06 m	d_p = 0.09 m	d_p = 0.12 m
Pebble number N_p in three-dimensional bed	27,000	8000	3375
Pebble number N_p in two-dimensional bed	1260	560	315

In the initial stage of the simulation, in order to eliminate the effect of the injection mode on the voidage, the continuous injection time of each condition is 0–5 s through the random injection of the pebbles at the entrance of the top of the bed. The pebbles are injected at 1.5 s as shown in Figure 2. The injection is stop at t = 5 s, then the whole process of the pebbles, only by gravity and free falling to the bottom of the bed, is a natural and random accumulation process.

Function of voidage definitions

As aforementioned, voidage is mainly divided into radial voidage and axial voidage. When calculating the radial voidage, this article divide pebble bed into two parts to solve the upper cylinder and cone base, respectively. It should be pointed out that, in this paper, the two-dimensional model that we study is a quasi-two-dimensional simulation of ball diameter thickness in the Z direction. So the two-dimensional simulation and three-dimensional simulation can use the same computational formula of voidage, it can also ensure that the calculation results of different dimensions are comparability. The voidage calculation formula of the upper cylinder is shown as below

In the equation (7), $δ h = h_{2} - h_{1}$ is height of the cylinder, $r_{i}$ shows the division of pebble bed cylinder radius. In this paper, the core ball bed is divided into 360 equal parts along the radial direction, and the volume calculation formula of sphere segment is used to simulated V_j, as follows

V_{j} = π H_{p}^{2} (R_{p} - \frac{H_{p}}{3})

(8)

where R_p represents the element sphere radius, and H_p represents the height of sphere segment. Accordingly, radial voidage of the cone base can be calculated as equation (9)

Unlike $δ h = h_{2} - h_{1}$ in equation (7), $Δ h (r)$ in equation (9) is a variable changed with r. Similarly, the axial voidage divide y-axis into 900 equal parts, shown as formula (10)

ɛ (y) = 1 - \frac{1}{A^{2} (y) δ y} \sum_{j = 1}^{N} δ V_{j} |_{y_{i - 1} \leq y (δ V_{j}) \leq y_{i}}

(10)

where

A (y) = π r^{2} (y_{i}')

δ y = y_{i} - y_{i - 1}

, and

y'_{i}

is the middle height between y_i and

y_{i - 1}

Simulation results and analysis

The static accumulation of the core pebbles is a natural random accumulation state, and the states of accumulation of various conditions are shown in Figure 3. It can be seen from the figure that the vertical and horizontal pebbles are almost confluent and have little local exchange.

Radial voidage of upper cylinder

The radial voidage distribution of the upper cylinder of the bed is shown in Figure 4. It can be seen from the figure that the void ratio of the three-dimensional condition is oscillating in the near wall region. After a distance from the wall, the voidage will approach a set value that is consistent with previous studies.¹⁴ In contrast, the void rate in two-dimensional conditions will always be shocked. Void ratio is due to the reason that the spherical objects cannot be tightly accumulated among the pebbles and the gap between the pebble and the pebble is greatest.

With the increase of the distance from the wall of the pebbles, the effect of the single wall on the accumulation of the ball is weakened by mutual cancellation and the random distribution of the core around each layer in the radial direction will counteract the effect of the spherical shape on the porosity.

In the two-dimensional bed, the influence of the wall on the accumulation state of pebbles cannot be neglected except the center of the pebble bed, and the void ratio at the diameter of the particles from the wall is significantly larger than it in the remaining layers. This is also consistent with the results obtained in Figure 3.

Radial voidage of conical substrate

The radial voidage distribution of the conical base is shown in Figure 5. It can be seen from the figure that the radial voidage distribution of the conical base is similar to that of the cylindrical radial velocities, and the void ratio in the three-dimensional case approaches a fixed value, however, since the radial direction is no longer the normal direction of the conical base, the oscillation phenomenon near the tapered wall surface is not as periodic as the void ratio in the cylinder. In practice, it should vary similarly in the direction of the inner normal of the tapered bottom, while the void in the two-dimensional conditions are still has been oscillating.

Compared with Figures 4 and 5, it can be seen that the void ratio of the conical substrate should be smaller than the radial void of the cylinder, which is consistent with the actual pressure of the core pebble due to the contraction of the cone, it is because that the void ratio has a characteristic that decreases as the pressure increases. The boundary of the two-dimensional static stacking conical base amplifies the oscillatory effect of the porosity, and the oscillation period is still related to the pebble diameter.

Axial voidage of z axial

The three-dimensional z-axis axial voidage distribution is shown in Figure 6. It can be seen from the figure that the bottom of the tapered substrate subjected to the maximum pressure, which means the axial porosity in the connecting region between the tapered base and the unloaded outlet is minimized.

The comparison of the two-dimensional and three-dimensional z-axis axial void ratio distributions is shown in Figure 7. It is not difficult to find that the void fraction distribution above the tapered base has little correlation with the diameter of the pebble and is decreasing as the stacking height decreases. This is because the greater the squeeze is, the smaller the void rate is. At the same time the higher the height of the core pebble due to the reduction of the upper region, thus the squeeze of the gravity reduced.

While the axial voids above the circular cone bottom of the two-dimensional working conditions remain oscillating, the magnitude of the shock is much larger than that in the three-dimensional condition, and the oscillation period is related to the diameter of the core ball. It should be pointed out that under the three-dimensional conditions, the axial voidage of the small unloading cylinder is also shocked. The tendency is that it becomes smaller when the diameter of the core pebble decreases. It shows that with the side wall effect, there is also the bottom wall effect.

Conclusions

The soft sphere model based on DEM is applied to simulate the pebbles’ static stacking of different scales of 2D model and 3D model. Radial and axial voidages have been calculated. The conclusions are as follows:

Radial and axial void ratios oscillate in both 2D model and 3D model and the oscillation of it mainly distribute nearby the surface. While void ratio far away from the surface, void ratio can be treated as a constant. Attenuation of the 2D model’s void ratio is rather slow, furthermore, amplitudes and frequencies of void ratio under different conditions are related to the diameter of pebbles.

The void fraction of the upper part of the bed is dominated by the wall factor in the change of the near wall. Under the three-dimensional condition, the wall factors could cancel out each other. Near the center of the area, the wall factor can be ignored. Moreover, the size of the area is related to the size of the pebble diameter, however two-dimensional conditions can hardly ignore the wall factor.

The minimum voidage in the three-dimensional condition is close to the small cylinder at the base of the cone. The void rate in the two-dimensional condition is always in the oscillating state, and the oscillation amplitude of it is much larger than the three-dimensional condition.

The z-axis axial void ratio of the different pebble diameters under three-dimensional conditions has the same trend and similar distribution. The void fraction distribution above the tapered base has little correlation with the diameter of the pebble and is decreasing as the stacking height decreases. The voidages of the pebble at the small unloading cylinder are also dependent on the diameter of the core. The z-axis axial porosity of the different core diameters in the two-dimensional case also presents the same trend.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by the National Natural Science Foundations of China (Grant No. 51576211), the Science Fund for Creative Research Groups of National Natural Science Foundation of China (Grant No. 51321002), the National High Technology Research and Development Program of China (863) (2014AA052701), and the Foundation for the Author of National Excellent Doctoral Dissertation of PR China (FANEDD, Grant No. 201438)

Appendix 1

References

Crossland I. The economics of nuclear power. Nuclear Fuel Cycle Science and Engineering. Woodhead Publishing, 2012, pp.100–126.

Caram

Hong

. Random-walk approach to granular flows. Phys Rev Lett 1991; 67: 828–831.

Rycroft

Grest

Bazant

et al.

Analysis of granular flow in a pebble-bed nuclear reactor. Phys Rev E 2006; 74: 021306.

Kadak AC and Ballinger RG. MIT modular pebble bed reactor (MPBR). In: 2nd international topical meeting on high temperature reactor technology, Beijing, China, 22–24 September 2004.

Pebble Bed Modular Reactor Limited. PBMR safety analysis report, 001929-207/4, Rev B, South Africa: PBMR, 1999.

Choi

Kudrolli

Rosales

et al.

Diffusion and mixing in gravity-driven dense granular flows. Phys Rev Lett 2004; 92: 174301.

Yang

Jiang

. Experimental investigation on feasibility of two-region-designed pebble-bed high-temperature gas-cooled reactor. J Nucl Sci Technol 2009; 46: 374–381.

Jiang

Yang

Tang

et al.

Experimental and numerical validation of two-region-designed pebble bed reactor with dynamic core. Nucl Eng Des 2012; 246: 277–285.

Thornton

. On the relationship between the modulus of particulate media and surface energy of the constituent particles. J Phys D Appl Phys 1993; 26: 1587–1591.

10.

Yang

Liu

. DEM Simulation of pebble flow in HTR-10 core by phenomenological method. Atom Energy Sci Technol 2013; 12: 2231–2237.

11.

Benenati

Brosilow

. Void fraction distribution in beds of spheres. AIChE J 1962; 8: 359–361.

12.

Bomboni

Cerullo

Lomonaco

. Simplified models for pebble-bed HTR core burn-up calculations with Monteburns2.0. Ann Nucl Energy 2012; 40: 72–83.

13.

Schertz

Bischoff

. Thermal and material transport in non isothermal packed beds. AIChE J 1969; 15: 597–604.

14.

du Toit CG. The numerical determination of the variation in the porosity of the pebble-bed core. In: Proceedings of 1st international topical meeting on high temperature reactor technology (HTR2002), Petten, Netherlands, 22–24 April 2002.

15.

du Toit CG. Analysis of the radial variation in the porosity of annular packed beds. In: Proceedings of the 3rd international topical meeting on high temperature reactor technology (HTR2006), Johannesburg, South Africa, 1–4 October, 2006.

16.

du Toit

. Radial variation in porosity in annular packed beds. Nucl Eng Des 2008; 238: 3073–3079.

17.

du Toit CG and Rousseau PG. The flow and heat transfer in a packed bed high temperature gas cooled reactor. In: Proceedings of the international heat transfer conference (IHTC14), Washington, DC, 8–13 August, 2010.

18.

du Toit

Rousseau

Greyvenstein

et al.

A systems CFD model of a packed bed high temperature gas-cooled nuclear reactor. Int J Therm Sci 2006; 45: 70–85.

19.

Cheng

Hsu

. Fully-developed, forced convection flow through an annular packed-sphere bed with wall effects. Int J Heat Mass Transfer 1986; 29: 1843–1853.

20.

Hunt

Tien

. Non-Darcian flow, heat and mass transfer in catalytic packed-bed reactors. Chem Eng Sci 1990; 45: 55–63.

21.

Sodré

Parise

JAR

. Fluid flow pressure drop through an annular bed of spheres with wall effects. Exp Therm Fluid Sci 1998; 17: 265–275.

22.

Vortmeyer

Schuster

. Evaluation of steady flow profiles in rectangular and circular packed beds by a variational method. Chem Eng Sci 1983; 38: 1691–1699.

23.

Arntz

Den Otter

Briels

et al.

Granular mixing and segregation in a horizontal rotating drum: a simulation study on the impact of rotational speed and fill level. AIChE J 2008; 54: 3133–3146.

24.

Gui

Fan

. Numerical study of heat conduction of granular particles in rotating wavy drums. Int J Heat Mass Transfer 2015; 84: 740–751.

25.

Gui

Yang

et al.

Numerical simulation and analysis of particle mixing and conduction in wavy drums. Dry Technol 2016; 34: 91–104.

26.

Gui

Fan

Cen

. A macroscopic and microscopic study of particle mixing in a rotating tumbler. Chem Eng Sci 2010; 65: 3034–3041.

27.

Tsuji

Tanaka

Ishida

. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol 1992; 71: 239–250.

28.

Mishra

. A review of computer simulation of tumbling mills by the discrete element method: part I – contact mechanics. Int J Miner Process 2003; 71: 73–93.

29.

Mishra

. A review of computer simulation of tumbling mills by the discrete element method: Part II – practical applications. Int J Miner Process 2003; 71: 95–112.

30.

Weerasekara

Powell

Cleary

et al.

The contribution of DEM to the science of comminution. Powder Technol 2013; 248: 3–24.

31.

Gui N. DEM-LES/DNS coupling simulation of particle collision in complex two phase flow. Hangzhou, China: Zhejiang University, 2010.