Nuclear Engineering and Technology 49 (2017) 1310e1317

Contents lists available at S

Nuclear Engineering and Technology

journal homepage: www.elsevier .com/locate/net

Original Article

Large eddy simulation of turbulent flow using the parallel computational fluid dynamics code GASFLOW-MPI

Han Zhang*, Yabing Li, Jianjun Xiao, Thomas Jordan Institute of Nuclear and Energy Technologies, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

a r t i c l e i n f o

Article history: Received 2 June 2017 Received in revised form 31 July 2017 Accepted 3 August 2017 Available online 30 August 2017

Keywords: GASFLOW-MPI large eddy simulation M&C2017 scalable linear solver

* Corresponding author. E-mail address: [email protected] (H.

http://dx.doi.org/10.1016/j.net.2017.08.003 1738-5733/© 2017 Korean Nuclear Society, Published licenses/by-nc-nd/4.0/).

a b s t r a c t

GASFLOW-MPI is awidely used scalable computational fluid dynamics numerical tool to simulate the fluid turbulence behavior, combustion dynamics, and other related thermalehydraulic phenomena in nuclear power plant containment. An efficient scalable linear solver for the large-scale pressure equation is one of the key issues to ensure the computational efficiency of GASFLOW-MPI. Several advanced Krylov subspace methods and scalable preconditioningmethods are compared and analyzed to improve the computational performance. With the help of the powerful computational capability, the large eddy simulation turbulent model is used to resolve more detailed turbulent behaviors. A backward-facing step flow is performed to study the free shear layer, the recirculation region, and the boundary layer, which is widespread in many scientific and engineering applications. Numerical results are comparedwith the experimental data in the literature and the direct numerical simulation results by GASFLOW-MPI. Both time-averaged velocity profile and turbulent intensity are well consistent with the experimental data and direct numerical simulation result. Furthermore, the frequency spectrum is presented and a e5/3 energy decay is observed for a wide range of frequencies, satisfying the turbulent energy spectrum theory. Parallel scaling tests are also implemented on the KIT/IKET cluster and a linear scaling is realized for GASFLOW-MPI. © 2017 Korean Nuclear Society, Published by Elsevier Korea LLC. This is an open access article under the

CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

Hydrogen risk evaluation in nuclear power plant (NPP) containment is an important issue for severe accident analysis, because hydrogen could be ignited anddamage the integrity of NPPs during a severe accident, such as the Fukushima accident in 2011 [1]. Accurate prediction of the hydrogen turbulent transport is the first, but crucial, step to determine hydrogen distribution. Extensive experimental research [2e4] and numerical simulation [5e7] were carried out to study the turbulent transport and other related phe- nomena in NPP containment. With the help of the powerful computational capability, computational fluid dynamics (CFD) be- comes a practical numerical tool to analyze the hydrogen behavior in NPP containment [8e10]. In order to reduce the computational time or enable high-fidelity predictionwithmore detailed geometry information for large-scale containment simulations, efficient scalable parallel computation, especially the high-performance scalable parallel linear solver, is one of the key points for a

Zhang).

by Elsevier Korea LLC. This is an

successful hydrogen distribution prediction. In this article, a large eddysimulation (LES) of thebackward-facing stepflow is performed by using the parallel CFD code GASFLOW-MPI to study the wall- bounded turbulent flow as well as to preliminarily evaluate the parallel performance of the GASFLOW-MPI.

GASFLOW-MPI is an advanced parallel CFD numerical tool that has been developed based on the message passing interface (MPI) library and domain decomposition technique. GASFLOW-MPI [11,12] has been developed, validated, and widely used to predict the complex thermalehydraulic behavior in NPP containments. In the past decades, GASFLOW was applied to simulate the hydrogen cloud distribution and evaluate the risk mitigation strategy for different types of nuclear plants, such as the EPR [13], the Inter- national Thermonuclear Experimental Reactor (ITER) [14], the German Konvoi-type PWR [15], the VVER [16], and the APR1400 [17], as well as the well-known open tests [18] and blind tests [19].

Obviously, an efficient scalable parallel linear solver for the large- scale symmetrical sparse equations derived from the pressure equation is one of the key issues to ensure the computational effi- ciency of GASFLOW-MPI. The preconditioned Krylov subspace iter- ative method is a type of efficient linear solver. The first Krylov

open access article under the CC BY-NC-ND license (http://creativecommons.org/

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e1317 1311

subspace method, the conjugate gradient (CG) method, was pro- posed by Hestenes and Stiefel [20]; it is used for solving the sym- metrical linear system. Then,many variantmethods for symmetrical linear equations were presented, such as the minimal residual (MINRES)method [21] and symmetric LQ (SYMMLQ)method [21]. At the same time, the Krylov subspace method also extended to the asymmetrical cases, such as the generalized minimal residual (GMRES) method [22]. Besides the choice of Krylov subspace methods, how to construct the preconditioning matrix is another important technology for the linear solver performance. In this pa- per, several advanced Krylov subspace methods and scalable pre- conditioningmethods are implemented and compared. The optimal combination of the Krylov subspace solver and preconditioning method is presented to improve computational performance.

In order to accurately predict the turbulent behavior, a suitable turbulent model should be used. Several turbulent models had been developed and validated in GASFLOW-MPI, including the algebraic turbulent model [23], the keepsilon turbulentmodel [23], and the LES turbulent model [24]. The LES turbulent model can resolve the large-scale turbulent fluctuations directly, wherein only the unresolved subgrid scale fluid motion should be modeled. Meanwhile, for the Reynolds-averaged NaviereStokes (RANS)- based turbulent models, such as the algebraic model and the keepsilon model, all turbulent fluctuations are unresolved. With the help of the powerful parallel computational capability, the LES turbulent model is used in this paper to resolve more detailed turbulence information. The standard Smagorinsky subgrid scale (SGS) model [25] is used to model the effects of unresolved small- scale fluidmotions in the LES turbulentmodel. The turbulent inflow boundary based on white noise is used to consider the turbulent information at the inlet. The LES turbulent model in GASFLOW-MPI had been validated by turbulent jet flows, which are the free shear turbulence [24]. In this study, we simulate a backward-facing step turbulent flow by the LES turbulent model to study the wall- bounded turbulent flows, a widespread turbulence phenomenon in scientific and engineering applications.

This paper is organized as follows. The physical model in GASFLOW-MPI is described in Section 2. The conservation equation, LES turbulent model, and numerical methods are discussed here. The parallel linear solver and parallel computing capability are discussed in Section 3. The turbulent simulation is presented and discussed in Section 4. The conclusions are presented in Section 5.

2. Numerical methodology in GASFLOW-MPI

2.1. Governing equations

GASFLOW-MPI is a powerful CFD numerical tool used to simu- late the complicated thermalehydraulic behavior in NPP contain- ments where a three-dimensional (3-D) transient compressible multicomponent NaviereStokes equation system is solved [11]. However, as only single-species isothermal gas flow is carried out in this paper, the radiation transfer model, combustion model, and mass/heat transfer model are therefore not considered in the following conservation equations, which include the volume equation, mass equation, momentum equations, and internal en- ergy equation [Eqs. (1e4)]. General thermodynamic equation of state, Eq. (5), and the general caloric equation of state, Eq. (6), are also used to close the governing equation system.

Volume equation

vV vt

¼ VV$ðb� uÞ (1)

Mass equation

vr vt

¼ V$½rðb� uÞ� (2)

Momentum equations

vðruÞ vt

¼ V$½ruðb� uÞ� � Vpþ V$sþ rg� V$~s (3)

Internal energy equation

vðrIÞ vt

¼ V$½rIðb� uÞ� � pV$u� V$q� V$~q (4)

General thermodynamic equation of state

p ¼ Zðr; TÞr R M

T (5)

General caloric equation of state

I ¼ Iðr; TÞ (6)

2.2. LES turbulent model

In this section, SGS turbulent models are introduced to model the unresolved terms ~s and ~q in the momentum equation and the energy equation, respectively.

SGS Reynolds stresses ~s could be expressed by Eq. (7) based on the Boussinesq hypothesis.

~sij ¼ �mt � 2Sij �

2 3 Skkdij

� (7)

Sij ¼ � vui=vxj þ vuj=vxi

� =2 (8)

In this paper, the standard Smagorinsky model is used to calculate the SGS turbulent viscosity, as shown in Eqs. (9e12).

mt ¼ rL2s jSj (9)

Ls ¼ CsD (10)

D ¼ V1=3 ¼ ðDxDyDzÞ1=3 (11)

jSj ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi 2SijSij

q (12)

In theory, the Smagorinsky constant Cs can be derived by assuming that the production and dissipation of subgrid-scale turbulent kinetic energy are in equilibrium [25]. In practice, it has been found that the best results could be obtained when the Smagorinsky constant is set as 0.1 [26].

For another unclosed term, SGS heat flux term ~q in the energy equation, it could be expressed by Eq. (13) based on the gradient hypothesis [27]. It has been found that for air-like forced flow in this paper, the simulation result is not sensitive to the value of turbulent Prandtl number, Prt [27]. According to the suggestions from the existing research [26], the turbulent Prandtl number Prt is set as 0.90 in this paper.

~qj ¼ �lt vT vxj

(13)

lt ¼ rmtPrt (14)

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e13171312

2.3. Implicit continuous-fluid Eulerian arbitraryeLagrangianeEularian algorithm

GASFLOW-MPI uses a robust all-speed numerical methodology, the implicit continuous-fluid Eulerian arbitraryeLagrangiane Eularian method, to solve the 3-D transient compressible Naviere Stokes equations in Cartesian or cylindrical coordinates [11]. The solutionofNaviereStokes equation is time-split into three stepsvia an operator-splitting technique [11]. At the first step and the third step, thediffusion term, convection term, and source termare treated in the explicit scheme to reduce the computational cost. At the second step, an implicit Lagrangian scheme is used to avoid the time step being limited by acoustic wave, as shown by Eqs. (15e18). By defining edp¼ pn¼ pII, the symmetrical elliptic pressure equation Eq. (19) can be derived from Eqs. (15e18) [11].

V II � V I Dt

¼ Vn$ h ðu$AÞII � ðu$AÞn

i (15)

ðrVÞII � ðrVÞI Dt

¼ 0 (16)

ðrVuÞII � ðrVuÞI Dt

¼ �Vn h V � pII � Pn

� i (17)

ðrVIÞII � ðrVIÞI Dt

¼ �VnpnV$ ðu$AÞII � ðu$AÞn

(18)

Dt2V$

” AVnVdp

ðrVÞI # � V

I

Vn � pI þ C� dp

¼ V I�pn � pI�

Vn � pI þ C� þ DtV$

h� ðu$AÞI � ðu$AÞn

�i ; (19)

where the coefficient, for an ideal gas, is

C ¼ pn �

R M,cvðTÞ

�A (20)

3. Pressure equation parallel solution

In order to provide more accurate and detailed solutions for containment simulations, the advanced parallel CFD software GASFLOW-MPI had been developed based on the original serial version GASFLOW and had been validated by the well-known test cases [18,19].

Obviously, an efficient scalable linear solver plays an important role for GASFLOW-MPI because a large-scale symmetrical sparse linear equation system derived from the elliptic pressure equation, Eq. (19), should be solved per time step. First, we evaluate the computational performance of the different Krylov subspace methods that are preconditioned by the block Jacobi method. Then, the effect of scalable preconditioning methods is analyzed and compared. Finally, the parallel scalability test is implemented on the KIT/IKET cluster.

A hypothetical 3-D H2 cloud in the airsteam mixture is calcu- lated in this section to evaluate the linear solver computational performance. The feature of this case is that the gas velocity is relatively small, whereas the pressure field drastically changes per time step. Therefore, a considerable number of iterations are required to achieve convergence. In this case, the total number of computational nodes is 200 � 200 � 200.

3.1. Computational performance of Krylov subspace methods

The preconditioned Krylov subspace iteration method is a family of advanced efficient linear solvers. The basic idea of Krylov subspace iteration method is that, at iteration step k, the method extracts a “good” approximate solution to the linear system A , x ¼ b from the Krylov subspace Kk by imposing the PetroveGalerkin condition [28,29].

Krylov subspace:

Kk ¼ KkðA; r0Þ ¼ span � r0;Ar0;…;A

k�1r0 � ; (21)

where x0 is the initial guess solution and r0 ¼ A , x0 e b is the initial residual.

PetroveGalerkin condition:

rk ¼ b� A$xk ⊥ LkðA; r0Þ; (22)

where Kk is the approximation space and Lk is the constraint space. Different versions of Krylov subspace methods stem from

different choices of constraint space Lk and approximation space Kk. In this paper, three Krylov subspace iteration methodsdCG method, MINRES method, and SYMMLQmethoddare used to solve the large-scale symmetrical sparse pressure equation. The choice of constraint space and approximation space for each Krylov subspace method is as follows.

CG method [20]:

Kk ¼ KkðA; r0Þ and � Lk ¼ KkðA; r0Þ (23) MINRES method [21]:

Kk ¼ KkðA; r0Þ and � Lk ¼ AKkðA; r0Þ (24) SYMMLQ method [21]:

�Kk ¼ ATKk � AT ; r0

� and Lk ¼ Kk

� AT ; r0

� (25)

The computational performance test is implemented on the KIT/ IKET server, which provides eight nodes in total and 16 cores on each node using the Intel Xeon Processor E5-2667 v2 CPU (Intel Corporation, Santa Clara, CA, USA) [18]. The convergence rate of three Krylovmethods preconditioned by the block Jacobi method is shown in Fig. 1. The result shows that, for the medium convergence criterion, the performance of all three Krylovmethods is almost the same. However, for the tight convergence criterion, the perfor- mance of CG takes priority over that of MINRES and SYMMLQ because of its strong stability. Because gas velocity is sensitive to pressure gradient, the CG iteration method is recommended to evaluate the pressure field with sufficient accuracy.

3.2. Parallel computing performance of preconditioner

Because there is only a limited amount of parallelism in stan- dard preconditioners such as symmetric successive overrelaxation method and incomplete LU factorization method, the alternative preconditioning techniques that are specially designed for parallel environments should be used. Comparedwith the additive Schwarz method preconditioning, less communication cost is required in the parallel block Jacobi preconditioning method. Thus, the parallel block Jacobi preconditioning method is used in this paper, where each processor has one block. Here, we focus on the effect of block solver on the total computational efficiency. Several block solvers, including the Cholesky factorization and the incomplete Cholesky factorization with fill level k [ICC(k); k ¼ 0, 1,…, 4] [30], are used to solve each individual block. On one hand, with the increase of the

0 100 200 300 400 500 600 10–25

10–20

10–15

10–10

10–5

100

No.of krylov iteration

N or

m al

iz ed

re si

du al

CG MINRES SYMMLQ

Fig. 1. Convergence rate of preconditioned Krylov methods. CG, conjugate gradient; MINRES, minimal residual; SYMMLQ, symmetric LQ.

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e1317 1313

fill level k, the iteration number of ICC(k) continues to decrease and is close to that of exact cholesky factorization, as shown in Fig. 2 and Table 1. On the other hand, the computational cost of block solver ICC(k) also increases with the increase of the parameter k. For our symmetrical sparse pressure equation system, the best balance is achieved at the point of ICC(2), as shown in Table 1. Therefore, the computational performance of CG solver preconditioned by ICC(2) block Jacobi method is the most efficient combinations in our test

0 20 40 60 80 100 10–20

10–15

10–10

10–5

100

No.of krylov iteration

N or

m al

iz ed

re si

du al

ICC(0) ICC(1) ICC(2) ICC(3) ICC(4) Cholesky

Fig. 2. Convergence rate of CG preconditioned by PBJ with ICC(k). CG, conjugate gradient; ICC(k), incomplete Cholesky factorization with fill level k; PBJ, parallel block Jacobi.

Table 1 Computational performance.

Averaged iteration #

Cost per iteration (sec)

Total cost (sec)

PBJ þ ICC(0) 909.1 0.0295 26.8 PBJ þ ICC(1) 694.5 0.0332 23.1 PBJ þ ICC(2) 588.7 0.0368 21.7 PBJ þ ICC(3) 507.1 0.0500 22.8 PBJ þ ICC(4) 453.4 0.0596 27.0 PBJ þ Cholesky 360.6 0.2260 81.5

ICC, incomplete Cholesky factorization; PBJ, parallel block Jacobi.

case. And the CG together with ICC(2) block Jacobi is the recom- mended solver for the elliptic pressure equation in the current GASFLOW-MPI. More efficient preconditioners are also in the plan, such as the algebraic multigrid preconditioner in Hypre library.

The parallel scalability test is performed in this section. To eliminate the effect of the limited memory bandwidth per core, the number of processes is increased from 16 to 128 at intervals of 16. Fig. 3 demonstrates the speedup with respect to one node (16 processes). The computational performance of 128 processes is about 7 times higher than that of 16 processes, which is a little less than the theoretical prediction value because of the additional communication cost. The numerical test shows that a linear scaling is achieved in our 3-D test problem, which indicates that the advanced CFD code GASFLOW-MPI could improve the computa- tional capability efficiently when increased numbers of parallel cores are used.

4. Experimental results

A turbulence benchmark problem, turbulent backward-facing step flow, is performed in this section to study the parallel scal- ability and to further validate the LES turbulent model. The backward-facing step flow is chosen because of its simple geometry but relatively complicated turbulent flow behaviors. Furthermore, extensive experimental research on this flow has been conducted and a large bibliographic database exists, allowing comparisons. In this paper, the backward-facing step flow at Re ¼ 5100 is simulated by LES turbulent model in GASFLOW-MPI, and the results are compared with the Jovic and Driver’s experiment [31] and the direct numerical simulation (DNS) results by GASFLOW-MPI. The LES turbulent model used in this paper had been validated by a wide range of Re number turbulent jet flows [24], where the tur- bulent free shear flow is the dominant phenomenon. Besides the free shear layer, wewill further validate the performance of the LES turbulent model in the recirculation region, the reattachment zone, and the boundary layer in this backward-facing step case, which is highly important for many scientific and engineering applications.

4.1. Problem description

The computational domain used for this simulation is shown in Fig. 4. The longitudinal length Lx downstream of the step is 20h, where h is the step height (h ¼ 9.8 mm). There is a channel

16 32 48 64 80 96 112 128 1

2

3

4

5

6

7

8

No.of processes

Sp ee

d- up

GASFLOW-MPI Practical linear scaling Theoretical linear scaling

Fig. 3. Speedup relative to one node (16 processes).

Fig. 4. Backward-facing step flow configuration.

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e13171314

ahead of the step whose length Li is equal to 10h ahead of the step. The computational domain in the vertical direction Ly and spanwise direction Lz is 6h and 4h, respectively.

The uniform mesh size is chosen in the streamwise and span- wise directions, and a nonuniformmesh distribution that is refined at the boundary layer near the wall and the free shear layer at the step is used in vertical direction, as shown in Fig. 5. Specifically, there are 384 computational meshes in the longitudinal direction and 32 meshes in the spanwise direction. We use 96 computational cells to resolve the vertical flow behavior where 35 computational cells are imposed within the step (y < h). The normalized mesh sizes in the three directions based on the inlet boundary layer shear velocity are Dxþ z 20, Dyminþ z 0.6, Dymaxþ z 62, and Dzþ z 30, respectively. And a total number of 1.2 million computational cells are used in this simulation.

In our research, the turbulent velocity boundary is used as the inlet boundary to consider the turbulent intensity at the inlet, as shown in Eq. (26). The terms u, v, and w represent velocity com- ponents in x, y, and z directions, respectively. In this case, we use the velocity profile proposed by Spalart [32] at Req ¼ 670 as the mean turbulent profile at the inlet umean(y,z), where the thickness of themomentum boundary layer q is equal to d99¼1.2h. The mean u outside the boundary layer is 7.72 m/s. The mean values for v and w are set equal to zero. In this paper, the white noise method is used to model the turbulent fluctuation u0(y,z), because it is easy to implement in a parallel environment, as shown in Eq. (27). This method had been used to study the backward-facing step flow in Ref. [33]. According to Yang’s suggestions [34], a channel ahead of the step is used in our research to develop a desired realistic

0 20 40 60 80 100 0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

Index no.of nodes

δ y /h

Free shear layer at the step

Boundary layer near the wall

Fig. 5. Mesh distribution in the wall-normal direction.

turbulence. A free-slip boundary condition v ¼ 0, vu/vy ¼ 0, vw/ vy ¼ 0 is used at the upper wall and spanwise boundaries, whereas the no-slip boundary conditions are applied at all other horizontal boundaries. A continuous boundary is used at the outflow boundary, which means that the outflow is completely developed and the velocities gradients vui/vxi are set equal to zero. The time step in our case is self-adjusted to limit the CFL number below 0.25.

uðy; zÞ ¼ umeanðy; zÞ þ u0ðy; zÞ (26)

u0ðy; zÞ ¼ ffiffiffiffiffiffiffi 3 2 k

r ,ε (27)

where k is the turbulent energy at the inlet and ε is the random number obeying standard normal distribution.

The viscous stress treatment near the wall is also an important issue to accurately simulate the profile within the boundary layer. In this paper, the grid spacing is fine enough to resolve the profile within the boundary layer. Therefore, a direct discretization of the viscous stress term rather than a wall function is used, as shown in Eqs. (28e29).

� sijAj

��� wall ¼ m

” vui vxj

þ vuj vxi

!# wall

Aj ¼ twAj; isj (28)

vv

vx

���� wall

¼ v Dx=2

; vw vx

���� wall

¼ w Dx=2

(29)

4.2. Numerical results

LES simulation over a backward-facing step is carried out and discussed in this section. First, the parallel scalability test is per- formed. Second, the characteristics of the unsteady flow are analyzed. Then, the important flow parameters, such as the reat- tachment length, mean longitudinal velocity, and turbulent in- tensity profiles, are calculated and compared with Jovic and Driver’s experiment [31] and the DNS results by GASFLOW-MPI. These parameters are averaged by 30 sampling points in time and 10 sampling points along the spanwise direction. Lastly, the fre- quency spectrum of velocity fluctuations is computed and compared with the theory frequency spectrum.

4.2.1. Parallel scalability analysis The parallel scalability test is further performed based on the

backward-facing step case. The CG solver preconditioned by ICC(2) block Jacobi method is used in this section. The numerical results show that the computational performance of 128 processes is 7 times higher than that of 16 processes, as shown in Fig. 6, which is similar to the hypothetical 3-D H2 cloud case described in Section 3. Fig. 6 demonstrates that the linear scaling is achieved. According to Eq. (30), the parallel efficiency of turbulent backward-facing step case is about 0.88, which indicates that high parallel performance is achieved for the typical turbulent flow using the CFD code GAS- FLOW-MPI.

E ¼ performance on N CPUs N � performance on one CPU (30)

4.2.2. Numerical result analysis The eddy behavior is resolved by the LES method, as shown in

Fig. 7. The oscillatory flow behavior is a phenomenon that is not

(A) (B) (C)

–5 0 5 10 0

10

20

30

40

50

60

Un (m/s)

Po si

tio n

(m m

)

X * =–0.33

LES DNS Exp

–5 0 5 10 0

10

20

30

40

50

60

Un (m/s) Po

si tio

n (m

m )

X * = 0.00

–5 0 5 10 0

10

20

30

40

50

60

Un (m/s)

Po si

tio n

(m m

)

X * = 0.66

Fig. 8. Mean longitudinal velocity profiles at three different streamwise positions. (A) Recirculation region. (B) Reattachment region. (C) Recovery region. LES, large eddy simulation; DNS, direct numerical simulation.

16 32 48 64 80 96 112 128 1

2

3

4

5

6

7

8

No.of processes

Sp ee

d- up

GASFLOW-MPI Practical linear scaling

Fig. 6. Parallel scalability test.

Fig. 7. Instantaneous velocity distribution. (A) Instantaneous velocity at 0.40 s. (B) Instantaneous velocity at 0.50 s.

(A) (B) (C)

0 1 2 0

10

20

30

40

50

60

UU (m2/s2)

Po si

tio n

(m m

)

X * =–0.33

LES DNS Exp

0 1 2 0

10

20

30

40

50

60

UU (m2/s2)

Po si

tio n

(m m

)

X * = 0.00

0 1 2 0

10

20

30

40

50

60

UU (m2/s2)

Po si

tio n

(m m

) X * = 0.66

Fig. 9. Square roots of Reynolds stresses components u02 at three different streamwise positions. Recirculation region. (B) Reattachment region. (C) Recovery region. LES, large eddy simulation; DNS, direct numerical simulation.

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e1317 1315

observed in the RANS model, because all velocity fluctuations are also averaged and only time-averaged results could be resolved in RANS. However, the reattachment location always varies over time in a certain region owing to the velocity fluctuation, as shown in Fig. 7. A similar phenomenon is observed and analyzed in previous investigations [31].

According to the methods in Ref. [35], the reattachment loca- tion, Xr, is determined where the time-averaged longitudinal ve- locity near wall equals to zero. The computed reattachment length is Xr ¼ 6.9h, which is consistent with the LES result Xr ¼ 6.8h in Refs. [36,37] and the experimental value (Xr ¼ (6 ± 0.15)h) [31].

Westphal and Johnston [38] concluded that the averaged ve- locity distribution is independent of the initial/boundary condi- tions and geometrical parameters under the normalized coordinate X* ¼ (x e Xr)/Xr. These results are further confirmed by Dubief and Delcayre [39] and Panjwani et al. [33]. In this paper, the normalized coordinates X* are used for comparison of our LES results with the experiment of Jovic and Driver [31] and the DNS results by GASFLOW-MPI. For the DNS, the mesh sizes in three directions are only half of those in the LES simulation. The comparisons are implemented between the computed LES results and the

experimental data of Jovic and Driver [31] for the nondimensional mean longitudinal velocity, as shown in Fig. 8. Three locations are used to make the comparison: (1) where the recirculation (X* ¼ e0.333) is located, (2) where the reattachment (X* ¼ 0.000) is located, and (3) where the recovery region (X* ¼ e0.666) is located. In detail, there is a slight backflow near the wall in the recirculation region. The mean longitudinal velocity changes to be positive at the height of about 5 mm, and then reaches the core region velocity at the height of about 15 mm. On the reattachment point, the velocity is equal to zero near the wall. In the recovery region, the velocity profile near the wall is the boundary layer. The computational

0.4 0.5 0.6 0.7 0.8 0.9 1 –4

–2

0

2

4

F lu

ct u

a tin

g v

e lo

ci ty

( m

/s )

Time (s)

10 0

10 1

10 2

10 3

10 4

10 –5

10 0

10 5

10 10

E (f

) (m

2 / s)

Frequency (/s)

Velocity signal

Frequency spectrum

(A)

(B)

Fig. 10. Frequency spectrum of velocity fluctuations. (A) Velocity signal. (B) Frequency spectrum.

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e13171316

results of LES are virtually coincident with that of DNS and agree well with the experimental data at all three different locations. This indicates that the LES method can achieve the same accuracy level as the DNS method for the averaged physical quantities.

The averaged turbulent intensity u02 computed with the current LES is compared with the experimental data and the DNS results, as shown in Fig. 9, where u0 is the velocity fluctuation in the x direction. The comparison is made at the same three locations. It can be observed that, for each measuring location, the turbulent intensity reaches the maximum value at the height slightly below the step where a large amount of turbulence is generated by the free shear mechanism. Then, the turbulent intensity begins to delay as the height increases owing to the turbulent dissipation effect. Generally speaking, the computational turbulence intensity predicts well with the experimental data and DNS results; however, it should be noted that the DNSmethod results are more accurate than those of the LES method. As shown in Fig. 9A, the LES result overestimates the tur- bulent intensity in the region close to the wall (y < h), whereas the DNS result almost overlaps with the experimental data. The analo- gous phenomenon has also been observed at the reattachment re- gion and recovery region. This indicates that the DNS method may have certain advantages for turbulent intensity prediction since all turbulent information is resolved directly in DNS calculation. Meanwhile, the unresolved subgrid scale motion has to be modeled using the SGS models in the LES model. However, we should note that the computation time for DNS is almost 16 times higher than that of the LES in this case. For the higher Re number turbulent flow, the computational cost of DNS will be prohibitively high.

The sign of the longitudinal velocity fluctuations during the period between 0.3 s and 1.0 s is plotted in Fig. 10. The distribution of velocity fluctuations is random, which is not captured by the RANS model. The frequency spectrum is analyzed, as shown in Fig. 10. In the corresponding energy spectrum, we observe a e5/3 energy decay for a wide range of frequencies (red line), followed by a steeper decay for f > 103, which is characteristic for the energy

cascade decay in turbulence. This also indicates that most turbulent energy is contained in lower frequencies, satisfying the turbulent energy spectrum theory.

5. Conclusions

An LES of a typical turbulent flow, backward-facing step flow, is performed using the parallel CFD code GASFLOW-MPI to study the wall-bounded turbulent flow. The choice of the parallel linear solver and the preconditioning technique for the large-scale pres- sure equation is discussed to pursue high computation efficiency. The numerical results show that the best computational perfor- mance is achieved when a CG solver preconditioned with ICC(2) block Jacobi is used. Furthermore, two parallel scaling tests, including a backward-facing step flow case and a hypothetical 3-D H2 cloud case, are also implemented on the KIT/IKET cluster. A linear scaling is realized for both test cases, which indicates that the computational capability of GASFLOW-MPI could be improved efficiently when increased numbers of parallel cores are used.

The numerical results of backward-facing step flow are analyzed and compared with the reference results. The reattachment length in the longitudinal direction calculated by GASFLOW-MPI is consistent with the other LES results and experimental data. The mean velocity profile and turbulence intensity agree well with the experimental data set and DNS results by GASFLOW-MPI. The DNS method possesses certain advantages for the turbulent intensity prediction in this case; however, the computation time for DNS is much higher than that of the LES. The frequency spectrum for the longitudinal velocity fluctuations u0 is analyzed. A e5/3 energy decay is observed for a wide range of frequencies, and it also in- dicates that most turbulent energy is contained in low frequencies. More advanced turbulent inlet models, such as the synthetic eddy method (SEM), have been implemented in the serial version of GASFLOW. The parallelization of the SEM is under development for GASFLOW-MPI.

H. Zhang et al. / Nuclear Engineering and Technology 49 (2017) 1310e1317 1317

Nomenclature

A outward normal fractional area vector Cs Smagorinsky constant I fluid internal energy Ls mixing length for subgrid scales M fluid molecular weight Prt turbulent Prandtl number R universal gas constant Sij rate-of-strain tensor jSj inner product of strain rate tensor T absolute fluid temperature V discretized fluid control volume Z fluid compressibility factor b control volume velocity surface vector g gravitational vector p pressure q heat flux vector ~q SGS heat flux term t time u fluid velocity vector D filter width Dt time step r fluid density lt turbulent conductivity coefficient mt SGS turbulent viscosity s viscous stress tensor ~s SGS Reynolds stresses term

Superscript I result at the first operator-splitting step II result at the second operator-splitting step

Conflicts of interest

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influ- ence our work.

References

[1] Y.S. Chin, P. Mathew, G. Glowa, R. Dickson, Z. Liang, B. Leitch, D. Duncan, A. Vasic, A. Bentaib, C. Journaeau, J. Malet, E. Studer, N. Meynet, P. Piluso, T. Gelain, N. Michielsen, S. Peillon, E. Porcheron, T. Albiol, B. Clement, M. Sonnenkalb, G. Weber, W. Klein-Hessling, S. Arndt, J. Yanez, A. Kotchourku, M. Kuznetsov, M. Sangiorgi, J. Fontanet, L.E. Herranz, C. Garcia de la Rua, A.E. Santiago, M. Andreani, J. Dreier, D. Paladino, R.Y. Lee, Organisation for Economic Co-Operation and Development, Nuclear Energy Agency-OECD/NEA, Committee on the Safety of Nuclear Installations-CSNI, Containment Code Validation Matrix, Le Seine Saint-Germain, 12 boulevard desIles, F-92130, Issyles-Moulineaux, France, 2014.

[2] S. Gupta, E. Schmidt, B. von Laufenberg, M. Freitag, G. Poss, F. Funke, G. Weber, THAI test facility for experimental research on hydrogen and fission product behaviour in light water reactor containments, Nucl. Eng. Des. 294 (2015) 183e201.

[3] J. Malet, E. Porcheron, P. Cornet, P. Brun, O. Norvez, B. Menet, L. Thause, ISP-47, International Standard Problem on Containment thermal-hydraulics, Step 1: TOSQAN-MISTRA, Phase A: airsteam mixtures TOSQAN experimental results, IRSN Rapport DPEA/SERAC/LPMAC/02e45, D�ecembre, 2002.

[4] J. Malet, E. Porcheron, P. Cornet, P. Brun, B. Menet, J. Vendel, ISP-47, Inter- national Standard Problem on Containment thermal-hydraulics, Step 1: TOSQAN-MISTRA, TOSQAN Phase B: air-steamehelium mixtures, Comparison Code-Experiments, IRSN Rapport DSU/SERAC/LEMAC/05e19, 2005.

[5] GOTHIC, Containment Analysis Package User Manual, Version 7.2a, Numerical Applications, Inc., 2006. Report NAI 8907e02, Revision 17.

[6] SNL, MELCOR Computer Code Manual, vol. 1: Primer and Users’ Guide, Version 2.1, Sandia National Laboratories, NUREG/CR-6119, vol. 1, Rev 3179, 2011.

[7] SNL, MELCOR Computer Code Manual, vol. 2: Reference Manual, Version 2.1, Sandia National Laboratories, NUREG/CR-6119, vol. 2, Rev 3194, 2011.

[8] S. Kudriakov, F. Dabbene, E. Studer, A. Beccantini, J.P. Magnaud, H. Paillere, A. Bentaib, A. Bleyer, J. Malet, E. Porcheron, C. Caroli, The TONUS CFD code for hydrogen risk analysis: physical models, numerical schemes and validation

matrix, Nucl. Eng. Des. 238 (2008) 551e565. [9] M. Houkema, N.B. Siccama, J.A. Lycklama Nijeholt, E.M.J. Komen, Validation of

the CFX-4 CFD Code for Containment Thermal-Hydraulics, Nucl. Eng. Des. 238 (2008) 590e599.

[10] H. Wilkening, L. Ammirabile, Simulation of helium release in the Battelle Model Containment facility using OpenFOAM, Nucl. Eng. Des. 265 (2013) 402e410.

[11] J. Xiao, J.R. Travis, T. Jordan, P. Royl, G.A. Necker, R. Redlinger, A. Svishchev, GASFLOW-III: A Computational Fluid Dynamics Code for Gases, Aerosols and Combustion, Volume 1: Theory and Computational Model, Karlsruhe Institute of Technology Report, 2014. Revision 3.5.

[12] J. Xiao, J.R. Travis, T. Jordan, P. Royl, G.A. Necker, R. Redlinger, A. Svishchev, GASFLOW-III: A Computational Fluid Dynamics Code for Gases, Aerosols and Combustion, Volume 2: User’s Manual, Karlsruhe Institute of Technology Report, 2014. Revision 3.5.

[13] H. Dimmelmeier, J. Eyink, M.A. Movahed, Computational validation of the EPR combustible gas control system, Nucl. Eng. Des. 249 (2012) 118e124.

[14] J. Xiao, J.R. Travis, W. Breitung, T. Jordan, Numerical analysis of hydrogen risk mitigation measures for support of ITER licensing, Fusion Eng. Des. 85 (2010) 205e214.

[15] P. Royl, H. Rochholz, W. Breitung, J.R. Travis, G. Necker, Analysis of steam and hydrogen distributions with PAR mitigation in NPP containments, Nucl. Eng. Des. 202 (2000) 231e248.

[16] P. Kostka, Z. Techy, J.J. Sienicki, Hydrogen mixing analyses for a VVER containment, in: 10th International Conference on Nuclear Engineering, Arlington, VA, USA, 2002. CD-ROM.

[17] J. Kim, U. Lee, S.-W. Hong, S.-B. Kim, H.-D. Kim, Spray effect on the behavior of hydrogen during severe accidents by a loss of coolant accidents in the APR1400 containment, Int. Commun. Heat Mass Transfer 33 (2006) 1207e1216.

[18] J. Xiao, J.R. Travis, P. Royl, G. Necker, A. Svishchev, T. Jordan, Three-dimen- sional all-speed CFD code for safety analysis of nuclear reactor containment: status of GASFLOW parallelization, model development, validation and application, Nucl. Eng. Des. 301 (2016) 290e310.

[19] P. Royl, J. Xiao, T. Jordan, Blind simulations of THAI test TH27 with GASFLOW- MPI for participation in the international benchmark conducted within the German THAI Program, Workshop Application of CFD/CFMD Codes (CFD4NRS_6) September 13e15, Cambridge, Ma, USA, 2016.

[20] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (1952) 409e436.

[21] C.C. Paige, M.A. Saunders, Solution of sparse indefinite systems of linear equations, SIAM J. Numer. Anal. 12 (1975) 617e629.

[22] Y. Saad,M.H.Schultz,GMRES:ageneralizedminimal residual algorithmforsolving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856e869.

[23] P. Royl, J.R. Travis, J. Kim, GASFLOW-III: A Computational Fluid Dynamics Code for Gases, Aerosols and Combustion, Volume 3: Assessment Manual, Karlsruhe Institute of Technology Report, 2014. Revision 1.0.

[24] H. Zhang, The development and validation of the large eddy simulation (LES) turbulent model in GASFLOW-MPI, GASFLOW Users’ meeting, Karlsruhe, Germany, 2016. CD-ROM.

[25] P. Sagaut, Large Eddy Simulation for Incompressible Flows: An Introduction, Scientific computation, Springer-Verlag, Berlin, 2001.

[26] Ansys Fluent, Theory Guide, Release 14.0, Ansys Inc, 2011. [27] G. Grotzbach, M. Worner, Direct numerical and large eddy simulations in

nuclear applications, Int. J. Heat Fluid Flow 20 (1999) 222e240. [28] Y. Saad, A. Henk Van der Vorst, Iterative solution of linear systems in the 20th

century, J. Comput. Appl. Math. 123 (2000) 1e33. [29] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Pub. Co., Boston,

1996. [30] T. Chan, H. van der Vorst, Approximate and incomplete factorizations, in:

D.E. Keyes, A. Samed, V. Venkatakrishnan (Eds.), Parallel Numerical Algo- rithms, ICASE/LaRC Interdisciplinary Series in Science and Engineering, Vol. 4, Kluwer Academic, Dordecht, 1997, pp. 167e202.

[31] S. Jovic, D.M. Driver, Backward-facing step measurements at low Reynolds number, Re¼5000, NASA Technical Memorandum NO: 108807, 1994, pp. 1e24.

[32] P.R. Spalart, Direct simulation of a turbulent boundary layer up to Re¼1410, J. Fluid Mech. 187 (1988) 61e98.

[33] B. Panjwani, I.S. Ertesvag, A. Gruber, K.E. Rian, Large eddy simulation of backward facing step flow, in: 5th National Conference on Computational Mechanics, MekIT09, Trondheim, Norway, 2009.

[34] Z. Yang, Large-eddy simulation: past, present and the future, Chin. J. Aeronaut. 28 (2015) 11e24.

[35] H. Le, P. Moin, J. Kim, Direct numerical simulation of turbulent flow over a backward-facing step, J. Fluid Mech. 330 (1997) 349e374.

[36] B. Zhang, Research and application of filtering grid scale and meshing adaptive-control strategy for large eddy simulation, Shanghai Jiao Tong Uni- versity, 2011. Ph.D. thesis.

[37] B. Zhang, T. Wang, C. Gu, Y. Dai, An adaptive control strategy for proper mesh distribution in large eddy simulation, J. Hydrodyn. 22 (2010) 865e871.

[38] R.V. Westphal, J.P. Johnston, Effect of initial conditions on turbulent reat- tachment downstream of a backward-facing step, AIAA J. 22 (1984) 1727e1732.

[39] Y. Dubief, F. Delcayre, On coherent-vortex identification in turbulence, J. Turbul 1 (2000) 1e22.