Improved delayed detached eddy simulation of a randomly stacked nuclear pebble bed

Improved delayed detached eddy simulation of a randomly stacked nuclear pebble bed

Computers and Fluids 122 (2015) 12–25 Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/compf...

7MB Sizes 0 Downloads 2 Views

Computers and Fluids 122 (2015) 12–25

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Improved delayed detached eddy simulation of a randomly stacked nuclear pebble bed A. Shams a,∗, F. Roelofs a, E.M.J. Komen a, E. Baglietto b a b

Nuclear Research and Consultancy Group (NRG), P.O. Box 25, 1755 ZG Petten, The Netherlands Massachusetts Institute of Technology (MIT), 77 Massachusetts Ave, Cambridge, MA 02139, United States

a r t i c l e

i n f o

Article history: Received 24 October 2014 Revised 15 June 2015 Accepted 14 August 2015 Available online 22 August 2015 Keywords: High temperature reactor Pebble bed Hot-spots IDDES

a b s t r a c t The high temperature reactor (HTR) design concept exhibits excellent safety features due to the low power density and the large amount of graphite present in the core which gives a large thermal inertia in the event of an accident such as loss of coolant. However, the possible appearance of hot spots in a pebble bed core design may affect the integrity of the pebbles and the fuel. This has drawn the attention of several scientists to understand this highly three-dimensional complex phenomenon. To obtain accurate predictions based on techniques such as DNS and LES, for a realistic (even for a small size) pebble bed flow, is still computationally too expensive and not foreseeable in the near future. On the other hand, the prediction capabilities of turbulence modelling approaches such as RANS and hybrid RANS-LES methods for such complex flow regime have not yet been rigorously evaluated. In this paper, numerical simulations of a limited sized randomly stacked bed of spherical pebbles are performed by using improved delayed detached eddy simulation (IDDES). The pebble bed configuration analysed consists of approximately 30 pebbles, which are randomly stacked to represent the core of an HTR. The obtained results are compared (qualitatively and quantitatively) with the available reference LES for validation purposes. The results show that the selected IDDES based on the SST k-ω model is able to reproduce the overall flow topology. The mean flow and thermal fields are found to be in good agreement with the LES. However, the second order statistics have shown significant disagreement. Nevertheless, these results are explicitly discussed in detail to understand the flow and thermal fields appearing in this complex flow configuration. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction High temperature reactors (HTR) are being considered for deployment around the world. An HTR uses helium gas as a coolant, while the moderator function is accomplished by graphite. The fuel is embedded in a graphite moderator and can sustain extremely high temperatures [7,8], basically preventing the fuel from melting. An important feature associated with the HTRs, is the high coolant temperature that can be achieved and therefore the high efficiency of the thermodynamic cycle. The high coolant outlet temperature also leads to optimal suitability for coupling the reactor to an industrial process requiring heat. The core can be designed using a graphite pebble bed. Some experimental and demonstration reactors have been operated over the world using this design [3] which have shown safe and efficient operation, however questions have been raised about the occurrence of potential local hot spots in the pebble bed, possibly affecting the pebble integrity [15]. Heat transfer around a curved surface shows a



Corresponding author. Tel.: +31 224568152. E-mail address: [email protected], [email protected] (A. Shams).

http://dx.doi.org/10.1016/j.compfluid.2015.08.015 0045-7930/© 2015 Elsevier Ltd. All rights reserved.

complex behaviour, which can be affected by both laminar and turbulent flow characteristics and by the effect of flow curvature. The narrow flow passages through the gaps between the pebbles can have concave and convex configuration which will cause augmentation or suppression of turbulence level [2]. In addition, pressure gradients strongly affect the boundary layers. Transition from laminar to turbulent, wake, flow separation and its respective reattachment make this flow configuration very complex, therefore a detailed evaluation of the pebble bed flow physics needs to be performed. Computationally, thermal-hydraulic aspects of a pebble bed configuration can be predicted either by a porous medium or by realistic geometry-resolving approaches. In case of a porous medium approach, an averaged concept of porosity is applied to simulate the closely packed geometry [8], whereas, in a realistic approach case, every pebble is accurately and separately modelled in a simulation, using direct numerical simulation (DNS), large eddy simulation (LES), Reynolds averaged Navier–Stokes (RANS) calculations or hybrid RANS/LES methods. However, applying a realistic approach to a reactor scale simulation is still computationally expensive and not foreseeable at an industrial scale. On the other hand porous media models are computationally efficient and represent a practical

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

13

Fig. 1. CFD triangle (left) and different pebble bed distributions (right) [21].

engineering approach for the simulation of large pebble cores. However porous methods could improve their predictive capabilities by incorporating complex behaviour from higher order simulation methods. To fulfil this gap, experimental studies or validated computational fluid dynamics (CFD) approaches could play a significant role in providing the necessary closure to the coarser porous solution. While a limited number of experimental studies have been performed, the available data are still insufficient in order to support rigorous validation of turbulence modelling approaches [1,12,13]. This limitation is directly related to the difficulty of extracting data from a complex random distribution of pebbles. For such scenarios, high fidelity CFD simulations are possible and/or available, and can be very helpful in quantifying the capabilities and limitations of the available turbulence closures. A large number of modelling attempts ([4–6,10,11]; [14,17,29]) have clearly evidenced the limitations of the adopted RANS methods. This limited success is not surprising, and is related to the intrinsic limitations of the RANS closures, as explained by Wu et al. [29]. Higher order LES based calculations have also been performed to investigate FCC and BCC type pebble bed configurations (see [2] and [11]). Unfortunately, no detailed validation on a reference database was provided in order to verify the accuracy of LES predictions for such complex flow physics. In 2013, Shams et al. [20] have proposed a stepwise validation procedure involving different types of pebble bed distributions, as shown in Fig. 1. This validation procedure consists of four steps, starting from an idealized pebble bed, via a limited sized random pebble bed, a limited sized random pebble bed with wall effects and finally the reactor scale pebble bed. The purpose of the method is to leverage the first three steps to validate a realistic numerical approach. Afterwards, this approach should be applied to derive parameters to incorporate in a porous medium representation in order to enable modelling of the full scale reactor pebble bed. As a step 1 of this validation procedure, an idealized face cubic centred pebble bed has been studied extensively. To generate the reference database quasi-direct numerical simulations (q-DNS) for this idealized pebble bed were performed, see Shams et al. [18,19]. The qDNS was further utilized to validate LES, hybrid RANS/LES and URANS methods for such complex geometries, see Shams et al. [20,21,22]. This extensive validation of step 1 has shown that among the used methods, LES is found to be in excellent agreement with the q-DNS with a maximum difference (for the first and second order statistics) of less than 6%. Hence, it can be said that for complex geometries where performing q-DNS calculations is not possible, LES can be used as a reference to validate the low order turbulence models. Lessons learned from step 1 have been carefully applied by Shams et al. [23] for step 2 to generate a reference database for a limited sized random

pebble bed. As a result, a high quality database for flow and thermal fields is obtained by using LES, for details see Shams et al. [23]. In the present study, the prediction capabilities of hybrid RANS–LES methods for this limited sized random pebble bed are evaluated. It is worth mentioning that in the validation procedure for step 1, it was found that among the tested hybrid methods, IDDES based in SST kω model performed the best. Hence, the same formulation of IDDES (i.e. IDDES SST k-ω model) was selected for the presented study. Details regarding the flow configuration and the considered numerical methods are given in Section 2. The obtained results of flow and thermal fields and their comparison with the reference LES are discussed and documented in Section 3. Finally, the summary and conclusions are presented. 2. Flow configuration and numerical strategies 2.1. Computational domain Following the reference LES case [23], the computational domain for the selected IDDES simulation is kept the same. It consists of approximately 30 pebbles with an average porosity (ε = volume of voids/total volume) level of about 0.4. These pebbles are randomly stacked in a rectangular domain of x = 0.177 m, y = 0.354 and z = 0.177 m, which is shown in Fig. 2. This random distribution of pebbles was obtained from [16]. The diameter of the pebbles is 6 cm. In the original geometry configuration, the pebbles were clustered in a cubic domain of 0.177 m per side. However, for the reference LES case, the inlet and outlet boundary conditions needed to be imposed. Hence, the domain was further extended in the y-direction, i.e. 0.0885 m for both inlet and outlet sections. As a result, the appearing pebbles on these sides are not truncated, as can be seen in Fig. 2 (Right). Additionally, in the original geometry the pebbles are in point contacts with each other. Whereas, to model such point contacts is problematic for meshing and consequently may induce numerical errors to the solution. Hence, in order to avoid difficulties, these point contacts are converted into small area contacts by scaling the pebbles by a factor of 1.034. This gives a maximum radial overlap of ∼1 mm between two pebbles and the corresponding contact area is 0.0019 cm2 . 2.2. Flow/simulation parameters Helium is considered as a working fluid with an imposed mass flow rate of 0.01607 kg/s at the inlet section. This is equivalent to the Reynolds number (based on the pebble diameter and the maximum velocity appearing in the computational domain) of ∼9753. The inlet

14

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 2. Schematic of the computational domain (left) isometric view (right) top view [23].

Table 1 Flow parameters. Properties

Values

Units

Density Viscosity Thermal conductivity Specific heat Mass flow rate Inlet temperature Heat flux

5.36 3.69 ∗ 10−5 0.3047 5441.6 0.01607 773 8317

kg/m3 N.s/m2 W/m.K J/kg.K kg/s K W/m2

temperature of helium is kept 773 K. Moreover, constant properties corresponding to this temperature are used and are given in Table 1. Following the reference LES case, to sustain the turbulent at the inlet section, a synthetic eddy method (SEM) [9] is used and the turbulence level is kept to 25%. In addition, periodic boundary conditions are imposed for mass flow rate and temperature on all four sides (i.e. in x- and z-directions) of the computational domain. A simplified assumption of constant heat flux of 8317 W/m2 is applied on the pebble surfaces. All the flow parameters used for this simulation are summarized in Table 1.

Fig. 3. Polyhedral mesh across the computational domain in (left) front view, (right) zoom near the pebbles and on a single pebble.

2.4. Turbulence modelling

2.3. Mesh generation STAR-CCM+ [28] is used to generate a fully polyhedral mesh for the selected computational domain. It is worth mentioning that the adopted meshing strategy is similar to the reference LES case. Fig. 3 displays the generated polyhedral mesh with an off-set layer near the pebble walls. This mesh consists of around 7 million grid cells with a dimensionless mesh size smaller than 1 in wall normal, around 10 in azimuthal and 12 in the cross-sectional directions. An extensive preceding comparative study of hybrid methods and q-DNS of an idealized pebble bed (see Shams et al. [20,21]) has shown that even a nondimensional mesh size of 10 and 23 in azimuthal and cross-sectional directions, respectively, is good enough to resolve the LES scales for such flow configurations. The results obtained corresponding to that mentioned mesh were found to be in good agreement with the reference q-DNS solution. Moreover, at the outlet section the mesh is coarsened to avoid the reverse flow effects.

In this paper, a hybrid (U)-RANS/LES turbulence modelling approach is selected to perform the numerical simulations. More specifically, a detached eddy simulation (DES) type hybrid method is considered. A DES formulation uses a blend of unsteady (U)RANS and LES, as it attempts to treat near-wall regions in a (U)RANS-like manner, and treat the rest of the flow in a LES-like manner. The RANS eddy viscosity transport equation contains a turbulence destruction term that is a function of the wall distance (d). In the subsequent DES approach developed by Spalart et al. [26], the wall distance was replaced by a length scale dependent on the grid size () which modifies the RANS model into a LES–SGS model. In this DES model the length scale switches between wall distances (RANS) and grid size (LES). For the current flow configuration, an improved formulation of DES method is selected, i.e. improved delayed DES (IDDES). The selected IDDES formulation is used in a combination with the SST k-ω model. This combination of IDDES with the SST k-ω model has been extensively tested and validated for a pebble bed flow configuration (for details see [21]) and has shown very good agreement with the reference

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

15

q-DNS database. The detailed description of the IDDES SST k-ω model is given below:

2.4.1. IDDES formulation The improved delayed detached eddy simulation (IDDES) of Shur et al. [25] is a further derivative of the delayed DES (DDES), for details see [21]. The IDDES formulation combines the DDES with the wall modelled LES (WMLES), and combines the applied RANS and LES length scales with blending functions. An essential element of IDDES is a new definition of the subgrid length-scale that includes explicit wall-distance dependence, which involves only the grid spacing [25]. This definition of the subgrid scale is given as:

 = min {max [Cw y, Cw max , wn ], max }

(1)

where, ‘Cw ’ is a constant, ‘y’ is the distance from the wall, ‘max ’ is the maximum local grid-spacing (= max (x , y , z )) and ‘wn ’ is the grid step in the wall normal direction. This sub-grid filter  is smaller than the conventional sub-grid filters in the near wall region. Far away from the wall, it is equal to the maximum local grid spacing. It thereby reduces the modelled stress in the near wall region, while boosting the stress far away from the wall. Principally, IDDES includes two branches, i.e. DDES and WMLES which are described below.

2.4.1.1. DDES branch of IDDES. This branch of the model is only active when the inflow conditions do not have any sort of turbulent content. The characteristic turbulent length scale lDDES in DDES is defined as:

lDDES = lRANS − fd max {0, lRANS − lLES }

(3)

Here, fd = 1 − F1 , is the delaying function for SST k-ω model, where F1 is shielding function of the SST model to protect the ‘RANS mode’ in the boundary layer from early switching to ‘LES mode’. In the length scale lDDES , lRANS is the RANS length scale, and for SST kω model, it is defined as lRANS = k1/2 /(cμ ω). Whereas, lLES is the LES length scale and is defined (via the subgrid length-scale ) as:

lLES = CDES 

(4)

Here, CDES is the constant of DES [24] whereas represents a lowReynolds number correction (see [27]).

2.4.1.2. WMLES branch of IDDES. This branch of IDDES is intended to be active only when the inflow conditions are unsteady and impose some turbulent content. The characteristic length scale of WMLES is defined as:

lWMLES = fB (1 + fe )lRANS + (1 − fB )lLES

(5)

Here, fB is a blending function RANS and LES, whereas, fe is an elevating function which is necessary in order to compensate for excessive reduction of the turbulent stresses near the wall region, for details see [25].

2.4.1.3. Blending of DDES and WMLES branches. The length scale for IDDES can now be calculated as:





d˜hyb = f˜d (1 + fe )lRANS + 1 − f˜d lLES

(6)

Where f˜d is the blending function and is given as f˜d = max{(1 − fdt ), fB }, and here fdt is a shielding, for details see [25]. For numerical simulations which contains inflow turbulent content: fdt is close to 0; f˜d is equal to fB , and thus d˜hyb = lWMLES . Otherwise, fe becomes 0 and d˜ = lDDES . hyb

Fig. 4. Isometric view of (left) the predicted flow features (at Re = 9753) in the computational domain (middle) the selected z-plane along (right) with the pebbles connected to this plane [23].

2.5. Numerical schemes The selection of numerical schemes for the hybrid methods, such as the selected IDDES based on SST k-ω model, poses a dilemma. The LES type calculations usually require a low dissipative numerical scheme as the sub-grid scale model provides the dissipation of the small turbulent scales. Hence, central differencing type schemes are frequently used to facilitate this requirement as upwind differencing schemes typically add too much dissipation. Such a central numerical scheme was used for the reference LES of a random pebble bed case, for details see [23]. On the other hand, the use of these central differencing schemes to a RANS calculation can lead to excessive numerical oscillations. To overcome this problem, a blend of central and upwind differencing schemes is recommended and has been frequently used for the hybrid methods. Hence, for the present study a second order hybrid upwind/central scheme for space discretization is selected to perform the calculations. Similar to the IDDES formulation, where a blending function is used between RANS and LES, in this hybrid scheme also a blending function is used. The formulation of the blending function is designed to deliver upwind differencing in areas of unidirectional flow or relatively coarse grids. Whereas, the central differencing is used for finer grids with higher vorticity and lower strain. As such, in areas where LES is possible and desirable, low-dissipative central differencing is achieved, and more stable upwind differencing is used elsewhere, for more details see [28]. In addition, an implicit three time level second order scheme is used for temporal discretization. The present discretization is performed by using collocated and a Rhie-and-Chow type pressure-velocity coupling combined with a SIMPLE-type algorithm. 3. Results and discussions The simulations are performed for a sufficiently long integration time in order to reach the statistical convergence for mean and RMS fields. The integration time is non-dimensionalized as τ + = Tsim × Uinlet(ave) /LD , where LD is the length of the computational domain, Tsim is the total simulation time and Uinlet(ave) is average inlet velocity. Simulations are performed for approximately τ + = 8.5, with a time step of ∼1.0 × 10−4 τ + , and 5 sub-iterations per time step. This time step corresponds to the average CFL of less than 0.5. The simulation has run on 120 processors (2.66 GHz) for ∼1.2 months. It is worth reminding that on a similar computational power, the reference LES case took 6 months to obtain a statistically converged solution. This suggests that the current IDDES computation is 5 times faster than the reference LES.

16

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 5. Iso-contours of time average velocity magnitude (U+ = U/Uinlet (ave) ) across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

Fig. 7. Iso-contours of time average v-velocity (v+ = v/Uinlet (ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

Fig. 6. Iso-contours of time average u-velocity (u+ = u/Uinlet (ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

Fig. 8. Iso-contours of time average w-velocity (w+ = w/Uinlet(ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

This section reports the detailed flow and thermal field analyses (on qualitative and quantitative basis) obtained from the IDDES SST k-ω and are extensively compared with the reference LES database [23]. These results are sub-divided into two main sections individually presenting the flow and thermal field analyses in the following sections.

pebbles. The flow from top to bottom is depicted by the streamlines, showing the principle direction. The pressure loss coefficient across this computational domain is 2.82 in case of IDDES application which is slightly over-predicted compared to the reference LES, which results in a pressure loss coefficient of 2.73. A mid-section across the computational domain (in the zdirection) is selected, as shown in Fig. 4. It is noticeable that this selected z-plane contains a number of pebbles on either side and for clarity these pebbles have been assigned numbers. The resulting complex flow is highlighted by the mean velocity distribution in Fig. 5. The iso-contours of mean velocity magnitude display a wide range of low and high velocity streaks including different flow complexities like stagnation regions and high velocity jets appearing in the narrow inter-pebble gaps. The overall flow topology is very well predicted by IDDES SST k-ω model and is in good agreement with the

3.1. Flow field analyses The reference LES results have shown that the selected flow domain exhibits a complex flow, hence, they pose a challenge for Hybrid and (U)RANS turbulence modelling approaches to predict it correctly. The complexities of the flow are evident from the LES solution shown in Fig. 4, which displays an isometric view and is highlighted by the instantaneous temperature distribution on the

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 9. Iso-contours of time average RMS of u-velocity (u+ rms = urms /Uinlet(ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

17

Fig. 11. Iso-contours of time average RMS of w-velocity (w+ rms = wrms /Uinlet(ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

Fig. 10. Iso-contours of time average RMS of v-velocity (v+ rms = vrms /Uinlet(ave) ) component across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

LES results. The prediction capabilities of the used hybrid method are further tested for a mid-section in the y-direction, i.e. xz-plane, see Fig. 5. In comparison with LES, IDDES shows good agreement in predicting the highly three-dimensional flow with strong rotational and streamline curvatures. This is mainly because of the improved variant of DES (particularly the WMLES branch of IDDES), which triggers the turbulent contents to nicely impose the boundary conditions for the LES part. Moreover, the explicit comparison of all the velocity components is given in Figs. 6–8. It is worth mentioning that U, u, v and w represent the velocity magnitude, x-, y- and z-velocity components, respectively. The u-velocity component exhibits strong reverse flow regions as shown in Fig. 6 and a similar behaviour is also noticeable for the

Fig. 12. Iso-contours of time average velocity field with the selected regions and the lines [23].

w-velocity component in Fig. 8. This clearly gives an indication of a strong cross and rotational flow with highly three-dimensional features. This strong cross and rotational flow is quite nicely captured by IDDES SST k-ω model and is in good agreement with the LES. Moreover, the appearance of shear flow in the centre of the computational domain is elaborated by the principle (v-) velocity component, see Fig. 7. Once again, IDDES displays a very good prediction in comparison with the reference LES. This highlights the importance of a blending function for a hybrid method, which is an improved one as

18

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 13. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 1, as in Fig. 12 (S+ = distance/X).

Fig. 14. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 3, as in Fig. 12 (S+ = distance/X).

Fig. 15. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 4, as in Fig. 12 (S+ = distance/X).

compared to the classical DDES method. Hence, the separated shear flow is nicely captured and shows the robustness of the selected IDDES SST k-ω model. The RMS of all the velocity components is compared with the LES and is given in Figs. 9–11. By looking at the RMS of u- and w- velocity components, it is noticeable that overall topology is very well predicted by the IDDES. However, it displays an under-prediction for this strong cross and rotational flow. On the other hand, the RMS of the principle velocity component highlights the shear-dominated regions (locations marked by “s” in Fig. 10). The flow characteristics of the separated shear layers are nicely reproduced by the IDDES, however, again an under-prediction is clearly noticeable. This is

understandable, since in the near wall region the IDDES formulation is switched into its URANS mode. To understand the combined influence of these velocity components on the overall flow-field, this z-plane is further sub-divided into three regions, i.e. A, B and C (similar to the reference LES) as shown in Fig. 12. The flow physics appearing in these respective regions are discussed in the following sections. 3.1.1. Region A In region A, the contribution of each velocity component is investigated at three different locations, i.e. lines 1, 3 and 4, as shown in Figs. 13–15. Time averaged velocity profiles for line 1 (as shown in

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

19

Fig. 16. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 2, as in Fig. 12 (S+ = distance/X).

Fig. 17. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 7, as in Fig. 12 (S+ = distance/X).

Fig. 13) clearly indicate the v-velocity as a dominant velocity component in this accelerated flow zone. It is worth reminding that the negative sign of v-velocity indicates that the flow is in negative ydirection, rather than the reverse flow. However, a significant contribution by u- and w-velocity components can also be observed. Where u- and w-velocities being constantly positive and negative, respectively, indicating the resultant flow direction is into the paper. The overall flow is accelerated to around 5.5 times higher than the average inlet velocity. The used IDDES model excellently captures this highly-accelerated stream. Fig. 14 displays the evolution of velocity profiles for line 3, running from pebble 5 to 4 and intersecting the two recirculation zones. Starting from S+ = 0.625 till 0.675, with the negative w-velocity component indicating flow going into the paper with low velocity magnitude. This is followed by high velocity streak with the resulting direction outward of the paper. By approaching pebbles 4 and 7, the magnitude of z-velocity component increases and remains positive, indicating a resultant flow outward of the paper. This complex phenomenon is well predicted by the IDDES SST k-ω model and shows a good agreement with the reference LES. Moreover, for line 4, the flow regime adjacent to pebble 5 is directed towards the paper with w-velocity being the dominant velocity component. This is followed by a high velocity stream (directed outward of the paper) with a velocity magnitude higher than the other zones in region A. This results in a slight mismatch (especially in the centre region) in the velocities predicted by the IDDES SST k-ω model. The rest of the regions are in very good agreement with the LES. Furthermore, by looking at the RMS profiles for these lines, it can be noticed that for lines 3 and 4 velocity fluctuations are higher than for line 1. Hence, the RMS for lines 3 and 4 are under-predicted by

the used IDDES model. Whereas, the RMS levels for line 1 are overpredicted. This is understandable, as line 1 exhibits an accelerating region, which is a resultant of the three-dimensional converging flow over the pebble, hence, requires a more sophisticated RANS model (used in the RANS mode of IDDES formulation) to correctly predict this anisotropic turbulent behaviour in the near wall region. Since, the current formulation of IDDES is based on an isotropic RANS model, which eventually over-predicts the RMS of principle velocity component. Nonetheless, these RMS profiles depict highly turbulent flow with strong three-dimensional effects and are not very well captured by the IDDES method. However, most importantly, the used IDDES model has preserved the trend of the turbulent fluctuations. 3.1.2. Region B The velocity profiles for 4 lines, i.e. 2, 5, 6 and 7, are extracted in region B and are given in Figs. 16–19. Line 2 is extracted at the accelerated zone as indicated by the strong dominance of v-velocity component and is nicely reproduced by the hybrid method. This predicted accelerated region is found to be in good agreement with the LES. Line 7 also displays a similar behaviour with strong dominance by the y-velocity component along with the w-velocity being negative and flow seems to be converging in between pebbles 10 and 11. The overall flow direction is going into the paper. The flow topology is well preserved by the IDDES SST k-ω model, however, with slight mismatch for all the velocity components. For line 5, a high velocity stream is guided by the strong v-velocity component, which is further decreased along the line shows the existence of a strong cross-flow. Line 5 displays an alternating occurrence of positive and negative u- and w-velocities, simultaneously, which

20

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 18. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 5, as in Fig. 12 (S+ = distance/X).

Fig. 19. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 6, as in Fig. 12 (S+ = distance/X).

indicates the spiralling in and out pair of vortices. This strong crossflow is further elaborated by the velocity profiles extracted for line 6. The IDDES SST k-ω model is able to reproduce the existence of strong cross-flow as indicated for line 5, however, with a slight mismatch in comparison with the LES. Whereas, for line 6 this mismatch increases and highlights the discrepancy of hybrid methods to accurately represent such strong cross-flow regions. Furthermore, the RMS profiles predicted by the IDDES SST k-ω model show significant under-prediction for the lines 5, 6 and 7 in region B. Whereas, line 2 is found to be in fair agreement with the LES. Overall, region B exhibits strong three-dimensional turbulence flow with significant contribution by u- and w- velocity components to maintain strong cross-flow streams. The flow topology is well represented by the IDDES SST k-ω model, however, with significant mismatch in the mean and RMS profiles of velocities. It is important to realize that the current flow configuration exhibit a highly three dimensional flow with several localized stagnation regions, hence, the use of an isotropic model (in the RANS mode of this hybrid formulation) makes the flow prediction more sensitive and results in a significant disagreement with respect to the reference LES. 3.1.3. Region C Two lines are extracted in region C and are shown in Figs. 20 and 21. Line 8 shows a high velocity magnitude dominated by the v-velocity component, which is followed by a large low-velocity region, where the contribution of all velocity components is of similar strength. The alternating positive and negative appearance of velocity components traces the existence of highly three-dimensional cross-flows. The appearance of a strong cross-flow in region C is very well captured by the IDDES SST k-ω and is found to be in good

Table 2 Average quantitative results (flow field) of IDDES w.r.t the reference LES. Property

Average % difference

Pressure drop Mean velocity Velocity RMS

3.3 3 15

Table 3 Predicted Nusselt number over pebbles. Pebble number

Reference LES

IDDES SST k-ω

% Difference

2 3 5 11 13 14

608.2 607.4 614.8 616.2 614.7 615.3

604.8 604.1 613.5 615.7 614.1 614.8

0.55 0.54 0.21 0.08 0.09 0.08

agreement with the reference LES. In addition, the velocity profiles for line 9 depict less turbulent flow governed by a dominant uvelocity which are excellently represented by the hybrid method. Moreover, the RMS profiles predicted by the IDDES SST k-ω model are in good agreement with LES, with only a slight under-prediction. In summary, the average percentage difference of the IDDES’s prediction with respect to the LES is given in Table 2. Although, the RMS levels throughout the computational domain are highly over-predicted by the IDDES SST k-ω model, the overall pressure drop and the mean velocity level are in very good agreement with the reference LES.

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

21

Fig. 20. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 8, as in Fig. 12 (S+ = distance/X).

Fig. 21. Evolution of (left) time average velocities and (right) the RMS of principle velocity (v-) component along Line 9, as in Fig. 12 (S+ = distance/X).

Fig. 22. Iso-contours of time average temperature distribution (T+ = (T – Tin )/Tconv ) across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

3.2. Thermal field analyses In the previous section, detailed flow-field analyses have elaborated the complexity in this small flow domain. In this section, the effect of the flow over the heat transfer phenomenon of these pebbles is further studied. It is worth mentioning that the simulations are performed (similar to the reference LES) by imposed heat flux rather

Fig. 23. Iso-contours of time average RMS of temperature (T+ rms = Trms /Tconv ) across (top) z- and (bottom) y-planes for the (left) LES and (right) IDDES SST k-ω.

than solving as a conjugate heat transfer case. Hence, as a result there appear peaks in the temperature prediction where these pebbles are in contact with each other. It is important to realize that such hotspots can also appears for a conjugate heat transfer, however, these would not be as pronounced as are observed here. Nevertheless, this results in high Nusselt numbers (Nu) over the pebble surfaces. The heat transfer in terms of average Nusselt number of a few pebbles is

22

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 24. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 1, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 25. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 2, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 26. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 3, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

computed and is given in Table 3. The results indicate that the averaged Nu number predicted by IDDES SST k-ω is in very good agreement with the reference LES with a maximum average difference well below 1%. The resulting mean temperature of the helium coolant predicted by IDDES SST k-ω model for the z- and y- planes and their respective RMS fields are given in Figs. 22 and 23. In comparison with the LES results, the IDDES SST k-ω model shows a very good agreement and is capable of capturing the complex heat transfer phenomena appearing in this complex flow regime.

By looking at the thermal fluctuations (T+ rms = Trms /Tconv. ), it is clearly noticeable that throughout the computational domain the hybrid method highly over-predicts the RMS level. This is further shown by the extracted profiles of mean temperature, its RMS and the turbulent heat flux for the principle velocity component (i.e. v-velocity) along the selected 9 lines, see Figs. 24–32. Starting with line 1, which appears in Region A, it can be noticed that the flow stream passing this line exhibits a very low temperature field except near the wall region. Similarly, subdued thermal fluctuations are observed as can be seen in Fig. 24 along with respective

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

23

Fig. 27. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 4, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 28. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 5, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 29. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 6, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

turbulent heat flux. It is clearly noticeable that the profile of mean temperature is very well predicted by the IDDES SST k-ω model. However, the RMS levels are highly over-predicted. On the other hand, the turbulence heat flux (THF) of the principal velocity component is not very well predicted by the hybrid method. In line 2, the passing fluid seems to be heated-up in the bulk and also exhibits an increase in the level of thermal fluctuations and the THF. The mean temperature

and its RMS predicted by the IDDES SST k-ω model are in good agreement with the reference database. Whereas, the turbulence heat flux shows a significant mismatch. Line 3 exhibits high velocities followed by a big recirculation zone, which as a consequence affects the temperature distribution and the predicted by the hybrid method is found to be in very good agreement. However, the RMS levels are significantly over-predicted by the

24

A. Shams et al. / Computers and Fluids 122 (2015) 12–25

Fig. 30. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 7, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 31. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 8, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

Fig. 32. Evolution of (left) time average temperature and its RMS and (right) turbulent heat flux of principle velocity (v-) component along Line 9, as in Fig. 12 (S+ = distance/X, THF = v’ T’ /Uinlet(ave) Tconv. ).

IDDES SST k-ω model. In addition for line 3, the resulting THF displays enhanced thermal activities in the centre and are nicely represented by the hybrid method with a slight mismatch. The overall topology of the mean temperature for lines 4, 5, 6, 7, 8 and 9 is nicely reproduced by the IDDES SST k-ω model, however, with a slight mismatch in the centre regions. The RMS level is always over-predicted by the hybrid method. On the other hand, a slight improvement in the predicted of THF is observed for all the aforementioned lines, however, with significant over- and under-predictions.

Table 4 lists the average percentage difference of the IDDES SST k-ω predictions with respect to the reference LES. The used hybrid method shows an excellent prediction of average Nu number in comparison with the LES. In addition, it displays only 5% difference with the reference database, which shows the prediction capability of the IDDES SST k-ω model for such a complex flow regime. However, the thermal fluctuations show significant discrepancies and rise up to 15% and 25% for the level of RMS and THF, respectively. It is worth mentioning that for such flow configurations,

A. Shams et al. / Computers and Fluids 122 (2015) 12–25 Table 4 Average quantitative results (thermal field) of IDDES w.r.t the reference LES. Property

Average % difference

Nu number Mean temperature Temperature RMS Turbulent heat flux

0.28 5 15 25

the main focus is on the prediction of the mean temperature throughout the computational domain. 4. Conclusions A three-dimensional numerical simulation of a limited sized randomly stacked spherical shaped pebble bed has been performed using the IDDES SST k-ω model. Extensive flow-field and heat transfer analyses are performed to understand this complex flow regime and are compared with the reference LES database. The predicted flow-field is complex and influenced by strong three-dimensional cross and rotational flows. Flow over the pebbles undergoes strong attenuation and enhancement of the turbulence and leads to flow separation and its subsequent reattachment. An extensive comparison of the flow-field, both qualitative and quantitative, is performed. This shows that the overall flow topology is nicely reproduced by the used hybrid method. The average percentage difference for the overall pressure drop and the mean velocities is found to be in the range of 3%. Whereas, a significant mismatch in the RMS of the velocities is observed and is found to in the order of 15%. In addition, similar comparisons for the thermal fields are performed. The average Nusselt number over the selected pebbles and the mean temperature throughout the computational field are very well predicted by the used IDDES SST k-ω model. The average percentage difference for these two properties is found in the range of 0.3% and 5%, respectively. On the other hand, the IDDES SST k-ω model has shown significant mismatch for the thermal fluctuations with an average difference of the order of 15% and 25% for the RMS and turbulence heat flux, respectively. Although, the IDDES SST k-ω model has shown some mismatch for the fluctuating part of velocity and thermal field, the mean and average field are very well predicted. This shows the significance of the use of this hybrid method for more complex geometries or larger flow domains, where the use of q-DNS and LES is challenging in terms of computer power requirements and meshing. Acknowledgements The work described in this paper is funded by the Dutch Ministry of Economic Affairs and the FP7 EC Collaborative Project THINS No. 249337. The limited sized periodic random pebble bed used in paper was provided by Dimitrios Pavlidis and Danny Lathouwers of the Technical University of Delft and is acknowledged here. References [1] Dominguez OEE, Perez CE, Hassan YA. Measurements of flow modification by particle deposition inside a packed bed using time-resolved PIV. In: Proceedings of 4th international meeting on high temperature reactor technology HTR 2008, September 28 – October 1., Washington DC, USA; 2008. [2] Hassan YA. Large eddy simulation in pebble bed gas cooled core reactors. Nucl Eng Des 2008;238:530–7.

25

[3] IAEA. 2001. Current status and future development of modular high temperature gas cooled reactor technology, IAEA-TECDOC-1198. [4] In WK, Chun TH, Lee WJ, Chang JH. CFD analysis of the fuel temperature in high temperature gas-cooled reactors. In: Transactions of the Korean nuclear society autumn meeting, Busan, Korea; 2005 October. [5] In WK, Lee SW, Lim HS, Lee WJ. Three-dimensional analysis of the hot-spot fuel temperature in pebble bed and prismatic modular reactors. In: Transactions of the Korean Nuclear Society Spring Meeting, Chuncheon, Korea; 2006 May. [6] In WK, Lee WJ, Hassan YA. CFD simulation of a coolant flow and a heat transfer in a pebble bed reactor. In: Proceedings of the 4th International Topical Meeting on High Temperature Reactor Technology HTR2008 September 28-October 1, 2008, Washington, DC USA; 2008. [7] Janse van Rensburg JJ, Viljoen C, van Staden MP. CDF modeling of high temperature gas cooled reactors. In: Proceedings of ICAPP’06, Reno, NV, USA; 2006. p. 6304. [8] Janse van Rensburg JJ, Kleingeld M. An integral CFD approach for the thermal simulation of the PBMR reactor unit. Nucl Eng Des 2011;241(8):3130–41 August. [9] Jarrin N, Benhamadouche S, Laurence D, Prosser R. A synthetic-eddy method for generating inflow conditions for large-eddy simulations. Int J Heat Fluid Flow 2006;27(4):585–93. [10] Lee JJ, Kang SK, Yoon SJ, Park GC. Assessment of turbulence models in CFD code and its application to pebble bed reactor. In: Proceedings of the Fourth International Conference on Heat Transfer, Fluid Mechanics and Thermodynamics, Cairo, Egypt; 2005 September. [11] Lee JJ, Yoon SJ, Park GC, Lee WJ. Comparison of LES and RANS in turbulenceinduced heat transfer and application to PBMR analysis. J Nucl Sci Technol 2007;44:7. [12] Lee JY, Lee SY. Flow visualization of pebble bed HTGR. In: Proceedings of the 4th International Meeting on High Temperature Reactor Technology HTR2008 September 28-October 1, 2008, Washington DC, USA; 2008. [13] Lee JY, Lee SY. Measurement of flow and thermal field in the pebble bed core of the VHTRG. In: Proceedings of the 5th International Meeting on High Temperature Reactor Technology HTR2010 October 18-20, 2010, Prague, Czech Republic; 2010. [14] Matzner HD. PBMR project status and the way ahead. In: Proceedings of the Second International Topical Meeting on High Temperature Reactor Technology, Beijing, China; 2004 September. [15] Moormann R., A safety re-evaluation of the AVR pebble bed reactor operation and its consequences for future HTR concepts, Jül-4275, Jülich, Forschungszentrum, Germany. [16] Pavlidis D, Lathouwers D. Realistic packed bed generation using small numbers of spheres. Nucl Eng Des 2013;263:172–8. [17] PBMR (Pty) Ltd., 2005, Data and boundary conditions to be used in VSOP, TINTE and MCNP PBMR 400MW (th) Reactor Models, 022028, Rev.01. [18] Shams A, Roelofs F, Komen EMJ, Baglietto E. Quasi-direct numerical simulation of a pebble bed configuration. Part I: Flow (velocity) field analysis. Nucl Eng Des 2013;263:473–89. [19] Shams A, Roelofs F, Komen EMJ, Baglietto E. Quasi-direct numerical simulation of a pebble bed configuration. Part II: Temperature field analysis. Nucl Eng Des 2013;263:490–9. [20] Shams A, Roelofs F, Komen EMJ, Baglietto E. Large eddy simulation of a nuclear pebble bed configuration. Nucl Eng Des 2013;261:10–19 August 2013. [21] Shams A, Roelofs F, Komen EMJ, Baglietto E. Numerical simulations of a pebble bed configuration using hybrid (RANS–LES) methods. Nucl Eng Des 2013;261:201–11 August 2013. [22] Shams, A., Roelofs, F., Komen, E.M.J., Baglietto, E., 2013. Towards numerical simulation of a nuclear pebble bed: guidelines for turbulence model selection, Paper no. 277, NURETH-15, Pisa, Italy. [23] Shams A, Roelofs F, Komen EMJ, Baglietto E. Large eddy simulation of a randomly stacked nuclear pebble bed. Comput Fluids 2014;96:302–21 June 2014. [24] Shur ML, Strelets MK, Spalart PR, Travin AK. Detached-eddy simulation of an airfoil at high angle of attack. Engineering turbulence modelling and measurements, 4. Elsevier (Conference Paper); 1999. p. 669–78 http://www.researchgate.net/publication/236888801. [25] Shur ML, Spalart PR, Strelets MK, Travin AK. A hybrid RANS-LES approach with delayed-des and wall-modeled les capabilities. Int J Heat Fluid Flow 2008;29:1638–49. [26] Spalart PR, Jou W-H, Strelets MK, Allmaras SR. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Advances in DNS/LES. Liu C, Liu Z, editors. Ruston, Louisiana: Greyden Press; 1997. [27] Spalart PR, Deck S, Shur ML, Squires KD, Strelets M, Travin A. A new version pf detached-eddy simulation, resistant to ambiguous grid densities. Theoret Comput Fluid Dyn 2006;20(3):181–95. [28] STAR-CCM+. User manual. London: CD Adapco; 2013. [29] Wu CY, Ferng YM, Chieng CC, Liu CC. Investigating the advantages and disadvantages of realistic approach and porous approach for closely packed pebbles in CFD simulation. Nucl Eng Des 2010;238:530–7.