- Email: [email protected]

Contents lists available at ScienceDirect

Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes

Numerical simulation of turbulent ﬂow in microscopic pore scale of pebble bed by large-eddy simulation S. Ebara a,∗ , T. Yokomine b , A. Shimizu b , H. Hashizume a a b

Department of Quantum Science and Energy Engineering, School of Engineering, Tohoku University, 6-6-01-2, Aramaki-aza aoba Aoba-ku, Sendai, Miyagi 980-8579, Japan Division of Energy Engineering and Science, Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, 6-1 Kasuga-kouen, Kasuga, Fukuoka 816-8580, Japan

a r t i c l e

i n f o

Article history: Available online 15 June 2010 Keywords: Pebble bed Large-eddy simulation Finite element method

a b s t r a c t Performance of a fusion reactor using pebble bed in its blanket depends on averaged properties such as packing density of the bed. On the other hand, crucial phenomena that affect safety issues, such as blockage of ﬂow area and the outset of a heat spot caused by thermal creep of constituent particles, depend on local properties of the bed. Conventional researches for pebble bed have been done in terms of average operation or coarse graining of the bed and can never capture the above local phenomena. In this study, the pebble bed is not coarse-grained but is treated such that small scale phenomena are treated as they are. Large-eddy simulation is used for simulation of turbulent ﬂow in microscopic pore scale of pebble bed in the frame of the ﬁnite element method. The bed is modeled as one with regular conﬁguration of spherical particles and ﬂuid motion in the pore is simulated in detail. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Currently, in order to develop breeding blanket systems for fusion DEMO plants, the test blanket program has been planned and advanced as one of the engineering test programs in ITER [1]. A water-cooling system with a solid breeder is one of the major candidates of Test Blanket Module concepts [2]. In this solid breeder blanket, solid pebble beds have characteristics of generating volumetric heating and tritium by incident neutron ﬂux and they are separated by membrane panels equipped with cooling channels. Generally speaking, performance of a fusion reactor using pebble bed in its blanket depends on averaged properties such as packing density of the bed. On the other hand, crucial phenomena that affect safety issues, e.g., stress concentration in the pebble bed, blockage of ﬂow area and the outset of a heat spot caused by thermal creep of constituent particles, depend on local properties of the bed such as relative position of the particles in the bed. Conventional researches for pebble bed have been done in terms of average operation or coarse graining of the bed and can never capture the above local phenomena because individual constitutive particle of the bed or ﬂuid occupying space among them is not treated as separated phase from the other at all but the bed is considered as a continuum and space everywhere in the bed is a homogeneous, mixed phase, e.g., every volume prescribed in the bed is composed of ﬂuid and solid phase in the portion of ε to

∗ Corresponding author. Tel.: +81 22 795 7905; fax: +81 22 795 7906. E-mail address: [email protected] (S. Ebara). 0920-3796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fusengdes.2010.04.065

1 − ε, where ε is porosity of the bed. In this study, the pebble bed is not coarse-grained but is treated such that each phase is separated from the other and small scale phenomena are treated as they are. In the pebble bed of the fusion blanket, many phenomena such as heat transfer, tritium transport, ﬂuid ﬂow, thermal stress, and so on will occur simultaneously and closely correlated. Firstly, in this paper, ﬂuid ﬂow is focused on and large-eddy simulation is used for simulation of turbulent ﬂow in microscopic pore scale of pebble bed. 2. Computational model and numerical method Pebble bed is assumed to consist of spherical particles with the same diameter, dp , arranged regularly in a space and have bodycentered cubic (bcc) structure. The bcc structure with a lattice constant of a can have various packing ratios according to √ dp , and its maximum packing ratio of 0.680 appears when dp = 3/2a. In this simulation, this packing condition is used. In addition, a part of the bed with unbounded region is computed. More speciﬁcally, the periodic boundary is imposed on every computational boundary in the simulation for the reason of numerical convenience. The ﬂow in the bed is driven by applied average pressure gradient in the ﬂow direction, which is given by Ergun’s expression [3] (Eq. (1)) in the simulation, where f is non-dimensional pressure gradient, ε is porosity of the bed, Re is Reynolds number based on microscopic (pore-averaged) ﬂuid velocity and pebble diameter, respectively. f =

150 Re

1 − ε 2 ε

+ 1.75

1−ε ε

(1)

S. Ebara et al. / Fusion Engineering and Design 85 (2010) 1638–1641

1639

Fig. 2. Comparison of mean velocity proﬁles with DNS data.

Fig. 1. Computational model of pebble bed.

Fig. 1 shows the computational model of the pebble bed. Assuming an incompressible ﬂow, governing equations are the following ﬁltered forms of continuity, momentum equations; ∂u¯ i =0 ∂xi ∂u¯ i ∂p¯ ∂ ∂u¯ + u¯ j i = − + ∂t ∂xj ∂xi ∂xj

(2)

2 Re

S ij − ij

+f

(3)

where all quantities are non-dimensionalized by using a particle diameter, dp , microscopic (pore-averaged) velocity, and p¯ denotes pressure deviation from the averaged pressure and f is added as a driving force of ﬂow. S¯ ij = 12 ((∂u¯ i /∂xj ) + (∂u¯ j /∂xi )) and Re denote the ﬁltered strain rate tensor and Reynolds number, respectively while ij is subgrid-scale stress tensor deﬁned as ij = ui uj − u¯ i u¯ j . We adopt the dynamic mixed model of Zang et al. [4] for ij as follows: ij −

1 1 m ı = −2SGS S ij + Lijm − ıij Lkk 3 ij kk 3

(4)

we solve the relaxation-transport equations for Lagrangian averaging, the SUPG stabilization is also used. Discretized equations are solved by the predictor-multicorrector algorithm (PMA) [8]. The ﬁnite element formulations are implemented in parallel using message passing interface (MPI) library. In parallel computations, the ﬁnite element mesh is partitioned into contiguous subdomains using METIS library [9]. The benchmark of a developed numerical code is ﬁrst done by comparing direct numerical simulation (DNS) database of a fully developed two-dimensional turbulent channel ﬂow between two parallel walls where Reynolds number based on a friction velocity, u , and a half channel width is 150 [10]. Mean velocity proﬁles are compared in Fig. 2 while root mean square (rms) values of velocity ﬂuctuations and Reynolds stress are shown in Figs. 3 and 4, respectively. Quantities of the vertical axes are normalized by using u , and the horizontal axes, y+ , are distance from the channel wall normalized by u and dynamic viscosity of ﬂuid. In Fig. 2, angle brackets indicates an average over time and u¯ indicates a velocity component in the ﬂow direction while in Fig. 3 urms , vrms and wrms correspond to intensities of velocity ﬂuctuation in the ﬂow direction, in the normal direction to the channel walls and in the parallel direction to the walls, respectively. In these ﬁgures, two kind of numerical data are also compared. Triangles (LDM) denote the data obtained from the Lagrangian-averaged stabilizing largeeddy simulation (LES) while circles (DMM) are those obtained from LES stabilized by spatial averaging in statistically uniform direction. We verify that all of statistical data in turbulent ﬂow obtained from

where Lijm is the modiﬁed Leonard term described as Lijm = u¯ i u¯ j −

¯ 2 |S| ¯ where u¯¯ i u¯¯ j , SGS , is an eddy viscosity deﬁned as SGS = (CS ) ¯ is grid-ﬁlter width. The model coefﬁcient of ¯ = |S| 2S¯ ij S¯ ij while

CS is obtained by use of the approach proposed by Germano et al. [5]. When it is dynamically determined, some kind of spatial averaging is introduced in order to prevent the numerical instability [5]. In the present simulation, the Lagrangian averaging procedure developed by Meneveau et al. [6] is used for that purpose. In this procedure, relaxation-transport equations are additionally to be solved. The detail is seen in the reference. In the present simulation, the stabilized ﬁnite element formulations based on SUPG (streamline-upwind/Petrov–Galerkin) and PSPG (pressure-stabilization/Petrov–Galerkin) [7] are developed and applied to solve the momentum equation (Eq. (3)) coupled with the continuity equation (Eq. (2)). Equation systems obtained after spatial discretization of the weak form of governing equations with stabilizing terms are solved by using the T1 formulation of Ref. [7]. The elements used are eight-node three-dimensional solid ones with continuous trilinear velocity and pressure. When

Fig. 3. Comparison of turbulent intensities in the ﬂow direction, urms ; wall-normal direction, vrms and spanwise direction, wrms .

1640

S. Ebara et al. / Fusion Engineering and Design 85 (2010) 1638–1641

Fig. 4. Comparison of proﬁles of Reynolds stress.

the present LES code show a moderate agreement with the reference data and it is found that LDM shows better agreement with DNS data than DMM.

Fig. 6. Velocity vector in the disjunct channel.

3. Results and discussion Fluid ﬂows in the pebble bed, circumventing pebbles as blockage. So, ﬂow path alternates between bifurcation and conﬂuence. In the body-centered packing used in the simulation, there are some paths which run straight through from entrance to exit. Figs. 5 and 6 show the velocity vector on the cross-section parallel to the ﬂow direction in the case of Re = 100. In Fig. 5, there are two “open” paths where ﬂuid goes straight without circumvention. The width of the paths repeatedly becomes wide and narrow, but ﬂow resistance is relatively small and the paths become the bypass of ﬂow in the bed. On the other hand, there is not any bypass in Fig. 6. Flow velocity is only partially allowed to become large. When a pore is large in the pebble bed but does not connect the bypass of ﬂow, the pore cannot contribute the ﬂow in the bed. It is found from the ﬁgure that there are many regions which have almost nothing to do with the ﬂow in the bed. These results, however, are considered to be limited for the cases without heat and mass transfer from the peb-

bles. If the pebbles produce heat by nuclear heating, stagnant ﬂow region can become heat spot. Then high temperature due to this heat spot will make compressible ﬂuid like He expand and affect ﬂow ﬁeld. Besides ﬂow ﬁeld, the high temperature can affect the structure of the pebble bed through thermal creep of pebbles and consequently ﬂow paths deform and then ﬂow ﬁeld also changes. So, thermal issue is very important for the soundness of the pebble bed and this should be incorporated into the ﬂow calculation in the bed. 4. Conclusion Turbulent ﬂow in microscopic pore scale of pebble bed was simulated by large-eddy simulation. For complex geometry in the bed, the ﬁnite element method was adopted and Lagrangian averaging is used to stabilize the dynamic procedure in LES. In this study, the pebble bed is assumed to have bcc structure and ﬂuid in the bed tends to ﬂow in uninterrupted paths with less resistance. Even if there is a large pore in the bed, it is not useful for ﬂowing ﬂuid when it is not connecting the “bypass”. This means that it is important for ﬂuid ﬂow in the bed not how large ﬂuid is occupying in the bed but how ﬂuid’s path is connecting in the bed. After all, averaged quantities of the bed such as porosity cannot explain the local phenomena and the conﬁguration or structure of the bed is very important to analyze small scale phenomena in the pebble bed. Acknowledgement This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Young Scientists (B), 20760577, 2008. References

Fig. 5. Velocity vector in bypass channels in the bed.

[1] V.A. Chuyanov, ITER Test Blanket Working Group, ITER Test Blanket Working Group activities: a summary recommendations and conclusions, Fusion Eng. Des. 61–62 (2002) 273. [2] M. Enoeda, T. Hatano, K. Tshuchiya, K. Ochiai, Y. Kawamura, K. Hayashi, T. Nishitani, M. Nishi, M. Akiba, Development of solid breeder blanket at JAERI, Fusion Sci. Technol. 47 (2005) 1060. [3] S. Ergun, Fluid ﬂow through packed column, Chem. Eng. Progr. 48 (1952) 89–94. [4] Y. Zang, R.L. Street, J.R. Koseff, A dynamic mixed subgrid-scale model and its application to turbulent recirculating ﬂows, Phys. Fluids A 5 (1993) 3186–3196.

S. Ebara et al. / Fusion Engineering and Design 85 (2010) 1638–1641 [5] D.K. Lilly, A proposed modiﬁcation of the Germano subgrid-scale closure method, Phys. Fluids A 4 (1992) 633–635; M. Germano, U. Piomelli, P. Moin, H. Carbot, A dynamic subgred-scale eddy viscosity model, Phys. Fluids A 3 (1991) 1760–1765. [6] C. Meneveau, T. Lund, W. Cabot, A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech. 319 (1996) 353–385. [7] T.E. Tezduyar, S. Mittal, S.E. Ray, R. Shih, Incompressible ﬂow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements, Comput. Methods Appl. Mech. Eng. 95 (1992) 221–242.

1641

[8] A.N. Brooks, Streamline upwind/Petrov–Galerkin formulations for convection dominated ﬂows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 32 (1982) 199–259. [9] http://glaros.dtc.umn.edu/gkhome/views/metis. [10] N. Kasagi, Y. Tomita, A. Kuroda, Direct numerical simulation of passive scalar ﬁeld in a turbulent channel ﬂow, ASME J. Heat Transf. 114 (1992) 598–606.

Copyright © 2018 KUNDOC.COM. All rights reserved.