Direct numerical simulation of incompressible turbulent boundary layers and planar jets at high Reynolds numbers initialized with implicit large eddy simulation

Computers and Fluids 194 (2019) 104314

Contents lists available at ScienceDirect

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

Direct numerical simulation of incompressible turbulent boundary layers and planar jets at high Reynolds numbers initialized with implicit large eddy simulation Tomoaki Watanabe∗, Xinxian Zhang, Koji Nagata Department of Aerospace Engineering, Nagoya University, Nagoya 464-8603, Japan

a r t i c l e

i n f o

Article history: Received 22 June 2019 Revised 9 September 2019 Accepted 24 September 2019 Available online 25 September 2019 Keywords: Turbulent boundary layer Turbulent jet Direct numerical simulation Implicit large eddy simulation

a b s t r a c t A direct numerical simulation (DNS) initialized with an implicit large eddy simulation (ILES) is performed for temporally evolving planar jets and turbulent boundary layers. In the ILES, an initial laminar flow develops into a fully developed state of the planar jet or the boundary layer. Subsequently, the DNS is started from the flow field obtained by the ILES. This hybrid ILES/DNS methodology is tested for the planar jet and boundary layer by comparing the results with full DNS started from the initial laminar flow. The ILES results used as the initial conditions of the DNS do not possess small-scale fluctuations. However, the small-scale fluctuations in the DNS grow with time and develop well within an interval of the integral time scale, where the influences of initial conditions taken from the ILES disappear for an energy spectrum of velocity fluctuations. The DNS initialized with the ILES well reproduces small-scale characteristics of turbulence, such as Reynolds number dependence of skewness and flatness of velocity derivative and energy spectrum of velocity fluctuations in the inertial subrange and viscous range. The DNS initialized with the ILES predicts well statistics dominated by large scales, such as 1st- and 2ndorder statistics and longitudinal auto-correlation function, in agreement with previous experimental and numerical studies. Reynolds number dependence of the mean velocity, root-mean-squared velocity fluctuations, Reynolds stress, shape factor, and skin friction in the turbulent boundary layers in the present DNS are consistent with previous experimental studies. These investigations confirm advantages of applying the ILES at the transitional flow region in the DNS of turbulent shear flows at high Reynolds numbers. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Turbulence is an important phenomenon in fluid dynamics, and it has been studied with theories, experiments, and numerical simulations [1,2]. A direct numerical simulation (DNS) that numerically solves Navier–Stokes equations without turbulence models is frequently used for studying turbulence. The DNS can provide reliable three-dimensional flow data, which is difficult to obtain in experiments. The DNS needs to resolve all length scales in turbulent flows from the largest to the smallest. The characteristic length scale of the smallest scale of turbulent motions is the Kolmogorov length scale η, which decreases with Reynolds number. Therefore, the number of grid points used in the DNS, N, increases with the Reynolds number [2] as N ∼ Re9/4 . The DNS of turbulence at high Reynolds number is of great interest in turbulence researches. As

∗

Corresponding author. E-mail address: [email protected] (T. Watanabe).

https://doi.org/10.1016/j.compfluid.2019.104314 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

finer grid spacing is required at higher Reynolds number, the computational cost of DNS also increases with the Reynolds number. The DNS has been used to study transitional turbulent flows, such as mixing layers, boundary layers, and jets [3–5], where an initially laminar flow develops into turbulence. These simulations are often used to investigate a fully-developed turbulent region. For examples, statistics are often discussed in the self-similar region of free shear flows based on the DNS [3,6]. However, an important fraction of the computation time is spent for simulating the initial transitional region. For reducing the computational cost of DNS, some previous studies have combined a large eddy simulation (LES) with the DNS. Because the LES does not have to resolve small-scale motions, whose influences are modeled by subgrid scale (SGS) models, larger grid spacing can be used in the LES compared with the DNS. Cifuentes et al. [7] conducted hybrid LES/DNS of turbulent jet flame, where the LES was used in the regions close to the nozzle and the outflow boundary while a flow region of interests was simulated with the DNS. Hence, the LES was used to provide an inflow boundary condition of the

2

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

DNS region. Temporally evolving stably stratified wakes were also simulated with the DNS initialized with the LES [8,9]. When the DNS is initialized with the LES results, the simulated flow does not have small-scale fluid motions at the beginning of the DNS because the LES solves the low-pass filtered governing equations. However, small-scale motions grow with time in the DNS. This hybrid LES/DNS methodology might be justified by the behavior of very small scales: Lalescu et al. [10] found that motions of scales smaller than 20η slave to motions of larger scales, and small-scale fluctuations rapidly appear in the DNS initialized with artificial turbulence without small-scale motions, which is obtained by removing small-scale fluctuations from the flow field of the DNS. Their study implies that when the LES has spatial resolution smaller than 20η, smaller-scale fluctuations grow without changing the larger scales in the DNS initialized with the LES. If the LES can be used instead of the DNS in a part of the simulation, the total computational cost can be reduced. Several previous studies used LES results as the initial or inflow-boundary condition of DNS [7,9]. However, it is not clear yet how the transition from the LES to DNS solutions occurs even though this is related to the validity of the numerical method. The present paper investigates this issue in details in simulations of incompressible turbulent flows. Numerical simulations are performed for temporally evolving planar jets [11] and temporally evolving turbulent boundary layers [12,13]. These flows develop with time from a laminar state in a computational domain that is homogeneous in the streamwise direction. Implicit LES (ILES) that relies on an implicit SGS model is performed until the flow reaches the fully turbulent state. Then, the DNS is started by using the LES result as the initial condition. The DNS initialized with the ILES is compared with full DNS started from the initial laminar flow. Furthermore, the present method with the ILES and DNS is used to simulate the planar jet and boundary layer at high Reynolds numbers, and the statistics of these flows are compared with previous experimental and numerical studies. The paper is organized as follows. Section 2 describes the numerical methodology of the DNS and ILES. Section 3 presents the details of numerical simulations of temporal-developing planar jets and boundary layers. The results of the jets and boundary layers are presented in Sections 4.1 and 4.2, respectively. The statistics of small-scale quantities in both flows are compared with previous studies of turbulence in Section 4.3. Furthermore, computation time of the present ILES/DNS approach is discussed in Section 4.4. Finally, the conclusion is presented in Section 5. 2. Numerical simulations 2.1. Direct numerical simulations of incompressible turbulent flows This study considers an incompressible flow with a passive scalar transfer, whose evolutions are described by Navier–Stokes equations and scalar transport equation:

∂uj = 0, ∂xj

(1)

1 ∂p ∂ ui ∂ ui u j ∂ 2 ui + =− +ν , ∂t ∂xj ρ ∂ xi ∂x j∂x j

(2)

∂φ ∂φ u j ∂ 2φ + = Dφ , ∂t ∂xj ∂x j∂x j

(3)

denote x, y, and z directions, respectively, and the corresponding components of the velocity vector are denoted by u, v, and w. The DNS numerically solves Eqs. (1)–(3) to obtain temporal evolutions of ui , p, and φ with appropriate initial and boundary conditions. 2.2. Implicit large eddy simulations The LES is performed with a relatively small number of grid points compared with the DNS. Then, the computational grid does not resolve all length scales in turbulent flows, and the variables simulated in the LES are denoted with an overbar f , which can be defined as a low-pass filtered value of f. The governing equations of filtered variables are written as follows:

∂uj = 0, ∂xj

(4)

1 ∂p ∂ ui ∂ ui u j ∂ 2 ui + =− +ν + Rui , ∂t ∂xj ρ ∂ xi ∂x j∂x j

(5)

∂ φ ∂ φu j ∂ 2φ + = Dφ + Rφ , ∂t ∂xj ∂x j∂x j

(6)

where Rui and Rφ are related to unresolved scales. The ILES relies on the implicit model for Rui and Rφ that approximates the influences of the unresolved scales. For this purpose, the ILES often utilizes a low-pass filter because the most important feature in the unresolved scales is considered as the dissipation of turbulent kinetic energy and scalar variance [8,14]. The main difference in the numerical method between the ILES and DNS besides the spatial resolution is the low-pass filter used as Rui and Rφ in the ILES. The same numerical schemes can be used for spatial discretization and temporal advancement of the governing equations in the ILES and DNS. 2.3. Numerical methods The simulation code used in this paper is the same as the one used in our previous studies with DNS and ILES [15–17]. The code is based on the fractional step method. Variables are stored on staggered grids. Spatial derivatives are computed with the fullyconservative central difference schemes [18], where the fourthorder and second-order schemes are used respectively in homogeneous and inhomogeneous directions. The temporal advancement is based on the third-order Runge–Kutta method. Poisson equation for pressure is solved with the Bi-CGSTAB method [19]. The ILES results are interpolated onto the DNS grid points with a trilinear interpolation when the DNS is initialized with the ILES. For the ILES, the tenth-order low-pass filter [20] is applied to ui and φ as the implicit SGS model for Rui and Rφ . The computational time step t is determined by

 t = CFL × min



1 , | u | / x + | v | / y + | w | / z

(7)

with the constant CFL = 0.3. 3. Test problems

where ui is the ith component of the velocity vector, p is the pressure, ρ is the constant fluid density, φ is the passive scalar, such as dye concentration in a dilute solution, ν is the kinematic viscosity, and Dφ is the diffusivity coefficient for φ . Variables ui , p, and φ are given as a function of time t and position xi . Here, i = 1, 2, and 3

The DNS initialized with the ILES is tested for two flows: a temporally evolving planar jet [11,21,22] and a temporally evolving turbulent boundary layer [12,13,23]. These are often studied as an approximation of their spatially evolving counterparts because various statistics have similar profiles in spatial and temporal simulations. These flows are simulated in a computational domain that is periodic in the streamwise (x) and spanwise (z) directions. The statistics are computed with spatial averages  ·  taken on a x − z

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

Fig. 1. Schematic of numerical simulations of (a) temporally evolving planar jet and (b) temporally evolving boundary layer. Parameters of the flows are also shown in figures: H is the initial jet width, UJ is the initial jet velocity, D is the trip wire diameter, and UW is the speed of the moving wall.

plane as a function of the cross-streamwise position y and time. A fluctuation from the average is denoted by f  = f −  f . The size of the computational domain is Lx × Ly × Lz , which is represented by Nx × Ny × Nz grid points. Grid spacing is uniform in the x and z directions while non-uniform grids are used in the y direction. 3.1. Temporally evolving turbulent planar jet The temporally evolving turbulent planar jet [11,21] is characterized by the jet Reynolds number ReJ = UJ H/ν and Schmidt number Sc = ν /Dφ , where H is the initial jet width and UJ is the initial jet velocity. The initial condition is illustrated in Fig. 1(a). The initial velocity profile is given by

u = u  + u  , v = v , w = w  ,

1 2



1 2

u = UJ + UJ tanh



H − 2|y| , 4 θJ

(9)

where θJ = 0.01H is used in this study. The center of the planar jet is located at y = 0. For the turbulent transition from the initial state, u , v , and w are added to the mean velocity in Eq. (8), where the root-mean-squared (rms) value of u , v , and w is 0.025UJ . The velocity fluctuations have non-zero values only in the jet (|y| ≤ 0.5H). The method for generating the initial velocity fluctuations are presented in Appendix A. The initial scalar profile is given by

1 2

1 2



φ = φJ + φJ tanh



H − 2|y| , 4 θJ

(10)

without scalar fluctuations. Here, φ J is the initial scalar value in the jet. Free slip boundary conditions are applied at the boundaries in the y direction. The grid location in the y direction is given by



y( j ) = −

instantaneous velocity and scalar profiles in J01-LES at t = 17tr are used as the initial condition of J01b. Similarly, J04-LES and J10-LES are the ILES started from t = 0 (Eqs. (8)–(10)), and their results at t = 17tr are used as the initial condition of the DNS in J04 and J10. Table 2 summarizes the grid size in the simulations. Here, y is taken at y = 0 while the Kolmogorov scale at y = 0 and t = 20tr in the DNS (J01a, J04, and J10) is used for normalizing the grid size. The DNS has x = z ≈ 1.7η and y ≈ 1.3η. On the other hand, the ILES uses the grid with x = z ≈ 5η and y ≈ 3η. The Kolmogorov scale is slightly smaller at earlier time. However, the grid size is smaller than 10η throughout the simulation in all ILES. Therefore, most length scales of turbulence are resolved in the ILES. In the ILES of the planar jet, the low pass filter is applied as the implicit SGS model at every 20 computational time steps. Here, the filter is not applied at every time step because the time increment of the simulation is calculated with the CFL number using the instantaneous velocity field that contains the mean flow contribution while the appropriate filtering interval should be determined by the time scale of velocity fluctuations. Table 2 also shows the time interval of filtering τ f normalized by the Kolmogorov time scale τη = (ν /ε )1/2 . τ f /τ η is of the order of 100 , and the filter is applied at the time interval close to the time scale of the smallest turbulent motions. The filter in the ILES is expected to dissipate the velocity fluctuations at the scales close to the grid size, which is about 6η in the present ILES. The velocity fluctuations at these scales can grow at the time scale slightly larger than τ η . Therefore, the filter applied at every τ f ≈ τ η is expected to work well as the implicit SGS model. It is also noteworthy that the ILES results do not change even if the filter is applied at every time step.

(8)

where u , v , and w are velocity fluctuations, and the mean velocity u is given by

3





Ly j−1 atanh (tanhαy ) 1 − 2 2αy Ny − 1



,

(11)

with αy = 1.5 and an integer j = 1, . . . , Ny . Eq. (11) gives smaller grid spacing near the jet centerline. Table 1 summarizes parameters in numerical simulations of the planar jet. The simulations are performed until time tend in the table, where tr = H/UJ is used as the reference time scale of the jet. ReJ is 10,0 0 0, 40,0 0 0, or 10 0,0 0 0 while Sc = 1 is used in all simulations. The numerical simulation of the planar jet is called Jab with ab representing the jet Reynolds number ReJ = ab0 0 0 0 while the ILES is referred with “-LES.” The size of the computational domain is (Lx , Ly , Lz ) = (6H, 10H, 4H ). Three simulations are performed for ReJ = 10, 0 0 0. J01-LES and J01a are performed respectively as the ILES and DNS, both of which are initialized with Eqs. (8)–(10). The

3.2. Temporally evolving turbulent boundary layer The turbulent boundary layer developing on a plate moving at a constant speed UW is considered as shown in Fig. 1(b). Hereafter, the subscript W refers a value on the wall. The simulation of the temporally evolving boundary layer is initialized with the velocity field that approximates the flow induced by a tripping wire with diameter D [12]:

u = u  + u  , v = v , w = w  ,

(12)

with velocity fluctuations u , v , and w generated by the method shown in Appendix A. The mean velocity profile is given by

1 2

1 2

u = UW + UW tanh

 D   y 1− , 2 θS D

(13)

where θS = 0.03D is the initial shear layer thickness. The rms value of u , v , and w is 0.05UW . The initial scalar profile is given as a function of y without fluctuations:

1 2

1 2

φ = φW + φW tanh

 D   y 1− . 2 θS D

(14)

The flow is characterized by the Reynolds number based on the trip wire diameter ReD = UW D/ν and Schmidt number Sc = ν /Dφ . The reference time scale is defined as tr = D/UW in the case of the boundary layer. The grid location in the y direction is given by



y ( j ) = Ly 1 −

1 tanh tanhαy

   j−1 , αy 1 − Ny − 1

(15)

with αy = 3.5 and j = 1, . . . , Ny , where the wall is located at y = 0. The grid spacing is smaller near the wall. Free slip boundary conditions are applied at the boundary at y = Ly while non-slip boundary conditions with constant values of u = UW and φ = φW are applied at y = 0. Parameters in the numerical simulations are summarized in Table 3, where the Reynolds number based on the momentum

4

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314 Table 1 Computational parameters of simulations of temporally evolving turbulent planar jets. Case

Initial condition

ReJ

Lx /H

Ly /H

Lz /H

Nx

Ny

Nz

tend /tr

J01-LES J01a J01b

Eqs. (8)–(10) Eqs. (8)–(10) LES (t = 17tr )

10000 10000 10000

6 6 6

10 10 10

4 4 4

288 864 864

500 1200 1200

192 576 576

20 20 20

J04-LES J04

Eqs. (8)–(10) LES (t = 17tr )

40000 40000

6 6

10 10

4 4

768 2304

1600 3200

512 1536

17 20

J10-LES J10

Eqs. (8)–(10) LES (t = 17tr )

100000 100000

6 6

10 10

4 4

1536 4608

2500 5800

1024 3072

17 20

Table 2 Grid size in numerical simulations of temporally evolving turbulent planar jets. Time interval of filtering τ f in the ILES (t = 17tr ) is also compared with the Kolmogorov time scale τ η . η and τ η are the Kolmogorov length and time scales, respectively, on the centerline at t = 20tr in the DNS.

x / η

y / η

z / η

τ f /τ η

J01-LES J01a J01b

−2

2.1 × 10 6.9 × 10−3 6.9 × 10−3

−2

1.2 × 10 5.0 × 10−3 5.0 × 10−3

−2

2.1 × 10 6.9 × 10−3 6.9 × 10−3

5.0 1.7 1.7

2.9 1.2 1.2

5.0 1.7 1.7

3.8 – –

J04-LES J04

7.8 × 10−3 2.6 × 10−3

3.8 × 10−3 1.9 × 10−3

7.8 × 10−3 2.6 × 10−3

5.2 1.7

2.5 1.3

5.2 1.7

2.4 –

J10-LES J10

3.9 × 10−3 1.3 × 10−3

2.4 × 10−3 1.0 × 10−3

3.9 × 10−3 1.3 × 10−3

5.3 1.8

3.2 1.4

5.3 1.8

1.8 –

Case

x / H

y / H

z / H

Table 3 Computational parameters of simulations of temporally evolving turbulent boundary layers. Case

Initial condition

L x /D

Ly /D

Lz /D

Nx

Ny

Nz

Reθ end

tend /tr

B02-LES B02a B02b B02c B02d

Eqs. (12)–(14) Eqs. (12)–(14) LES (Reθ = 1540) LES (Reθ = 1750) LES (Reθ = 1900)

56 56 56 56 56

76 76 76 76 76

28 28 28 28 28

288 864 864 864 864

300 400 400 400 400

288 512 512 512 512

2000 2000 2000 2000 2000

235 239 235 235 234

B04-LES B04

Eqs. (12)–(14) LES (Reθ = 3300)

120 120

162 162

60 60

648 1536

648 864

648 1024

3300 4000

746 1052

B06-LES B06

Eqs. (12)–(14) LES (Reθ = 5200)

170 170

230 230

85 85

864 2048

864 1152

864 1536

5200 6000

1587 1967

B08-LES B08

Eqs. (12)–(14) LES (Reθ = 6400 )

220 220

298 298

110 110

1152 2916

1152 1728

1152 1458

6400 8000

2340 2831

B10-LES B10

Eqs. (12)–(14) LES (Reθ = 8800 )

276 276

360 360

138 138

1458 3456

1458 2048

1458 2592

8800 10000

3247 3846

B13-LES B13

Eqs. (12)–(14) LES (Reθ = 11800 )

360 360

480 480

180 180

1728 4608

1944 2700

1728 3456

11800 13000

4690 5382

thickness is defined as Reθ = UW θ /ν with the momentum thick ∞ 2 dy. Turbulent transition of the temness θ = 0 u(UW − u )/UW porally evolving boundary layer successfully occurs for ReD ≥ 500 in previous DNS study [12]. As ReD increases, the transition occurs at earlier time. However, we have found that ReD = 40 0 0 results in an over-tripped boundary layer, in which the mean streamwise velocity profile does not agree with experiments of zero-pressure gradient turbulent boundary layers. Therefore, all simulations use ReD = 20 0 0 as in previous studies of temporally evolving turbulent boundary layers [12,13,23]. All simulations are conducted with Sc = 1. The size of the computational domain (Lx , Ly , Lz ) is determined based on the target value of Reθ . The 99% boundary layer thickness δ is defined as the vertical location of (UW − u )/UW = 0.99. In the fully developed state, δ ≈ (8Reθ /ReD )D can be used to estimate δ /D from Reθ and ReD [12]. According to previous studies [12,24,25], the computational domain (Lx , Ly , Lz ) should satisfy Lx > 2π δ , Ly > 3δ , and Lz > π δ . Substituting δ = (8Reθ /ReD )D yields Lx /D > 16π Reθ /ReD , Ly /D > 24Reθ /ReD , and Lz /D > 8π Reθ /ReD . The target value of Reθ is 20 0 0, 40 0 0, 60 0 0, 80 0 0, 10,0 0 0, or 13,0 0 0. The numerical simulation of the boundary layer is called

Bab with the target value of Reθ = ab0 0 0. Reθ monotonically increases with time, and the simulations are performed until Reθ reaches Reθ end at t = tend shown in Table 3. Five simulations for Reθ = 20 0 0 are performed for testing the DNS initialized with the ILES. B02-LES and B02a are respectively the ILES and DNS, both of which are initialized with Eqs. (12)–(14). B02b, B02c, and B02d are the DNS initialized with the B02-LES at different Reθ , and time is advanced until Reθ = 20 0 0 in the DNS. Similarly, the DNS of B04, B06, B08, B10, and B13 are performed by using the ILES results as the initial condition as summarized in Table 3. The low pass filter in the ILES is applied at every 500 time steps. The filtering interval is determined such that the filter does not cause overdamping of velocity fluctuations, which affects the initial turbulent transition of the boundary layer. In this case, the mean velocity profile normalized by the viscous scales does not agree with previous experiments and DNS of spatially developing turbulent boundary layers. It turned out that the overdamping occurs when the filtering time interval τ f divided by the viscous time scale ν /u2τ is O(10−2 ). The present ILES is performed with τ f /(ν /u2τ ) = 5–10. The ILES with the filtering interval of 200 steps or 700 steps was also performed with the same condition as B02-LES. However, influences of the

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

5

Table 4 Spatial resolution in the ILES of turbulent boundary layers. Grid size divided by the viscous = i /(ν /uτ ) is evaluated at the end of the simulation. Table also shows the grid length + i size divided by the Kolmogorov length scale η and the time interval of filtering τ f divided by τ η at y = 0.35δ, where δ , η, and τ η are obtained at t = 41tr (Reθ = 1200) in the DNS (B02a). Case

+x

+yW

+z

[x /η]0.35δ

[y /η]0.35δ

[z /η]0.35δ

[τ f /τ η ]0.35δ

B02-LES B04-LES B06-LES B08-LES B10-LES B13-LES

16.8 14.9 14.9 14.2 13.9 14.8

0.28 0.26 0.26 0.25 0.23 0.22

8.4 7.5 7.5 7.1 7.0 7.4

10.7 10.2 10.8 10.5 10.4 11.4

2.7 1.4 1.1 0.8 0.7 0.6

5.3 5.1 5.4 5.2 5.2 5.7

3.5 2.9 2.6 2.5 2.1 1.9

Table 5 Computational domain size divided by the 99% boundary layer thickness, grid size divided by the viscous length +i = i /(ν /uτ ), and grid size divided by the Kolmogorov length scale at y = 0.35δ or y = 0.5δ . These results are taken at the end of the DNS. Case

Lx /δ

Ly /δ

Lz /δ

+x

+yW

+z

[x /η]0.35δ

[y /η]0.35δ

[z /η]0.35δ

[y /η]0.5δ

B02a B02b B02c B02d B04 B06 B08 B10 B13

6.6 6.7 6.7 6.7 7.1 6.9 6.2 6.5 6.4

9.0 9.1 9.1 9.1 9.6 9.3 8.4 8.4 8.5

3.3 3.4 3.4 3.4 3.6 3.4 3.1 3.2 3.2

5.5 5.5 5.5 5.5 6.1 6.2 5.5 5.7 5.5

0.21 0.21 0.21 0.21 0.19 0.19 0.16 0.16 0.16

4.7 4.7 4.7 4.7 4.6 4.1 5.5 3.8 3.7

1.8 1.8 1.9 1.9 1.7 1.6 1.3 1.3 1.2

1.5 1.5 1.5 1.5 1.1 1.0 0.9 0.8 0.8

1.5 1.5 1.6 1.6 1.3 1.1 1.3 0.9 0.8

1.9 1.9 1.9 1.9 1.3 1.3 1.1 1.1 1.0

filtering interval were not observed for statistics presented in this paper. Table 4 shows the spatial resolution of the ILES, where the grid size divided by the viscous length ν /uτ is denoted by √ + = i /(ν /uτ ), uτ = τW /ρ is the friction velocity, and τW = i −ρν (∂ u/∂ y )W is the wall shear stress. DNS of wall turbulence + is often conducted under the condition of + x < 9.7, yW < 0.2, and + z < 4.8 in previous study [26]. The grid size in the present ILES is slightly greater than these values. The table also shows the grid size divided by the Kolmogorov length scale η at y = 0.35δ, where η and δ are obtained at t = 41tr (Reθ ≈ 1200) in the DNS (B02a). Here, η and δ in B02a are used because all simulations are performed with the same ReD and the main difference among the simulations is the size of the computational domain, which hardly affects η and δ at t = 41tr . All simulations have x /η ≈ 10 and z /η ≈ 5 in the ILES. Smaller y /η compared with x /η and z /η are used in the ILES such that the mean velocity profile near the wall is well resolved throughout the simulations. Table 4 also shows the time interval of filtering τ f normalized by the Kolmogorov time scale τη = (ν /ε )1/2 . As in the ILES of the planar jet, τ f /τ η is of the order of 100 . Table 5 shows the computational domain size and the grid size in relation to the length scales of the flow at the end of the DNS. The grid size is small enough for performing the DNS in terms of the viscous length (ν /uτ ) and Kolmogorov scale η. 4. Results and discussions 4.1. Temporally evolving turbulent planar jet Fig. 2 shows passive scalar φ on a x − y plane of the planer jet in J10-LES. The turbulent planar jet develops with time and spreads in the y direction. The turbulence appears near the edge of the jet in Fig. 2(a) and (b), where the flow near the centerline has φ = φJ . In later time in Fig. 2(d), the entire jet region becomes turbulent as attested by the patterns of small-scale scalar fluctuations in the jet. The ILES of the planar jet is validated by comparing J01a and J01-LES. Fig. 3(a) and (b) show the cross-streamwise profiles of

mean streamwise velocity u and rms streamwise velocity fluctuation urms at different time instances. The ILES results of u and urms agree well with the DNS results, and the present numerical methods and conditions are suitable for the ILES of the planar jet. The flow does not have small-scale fluctuations at the beginning of the DNS initialized with the ILES. To understand the evolution of small-scale velocity fluctuations in the DNS, an energy spectrum of the streamwise velocity fluctuation Eu is computed with the Fourier transform in the x direction (kx : streamwise wavenumber). Fig. 4 shows the temporal evolution of Eu in J01b, which is initialized with the ILES at t = 17tr . Fig. 4 also shows Eu in the ILES, which rapidly decreases with kx around kx H = 80 although the inertial subrange with Eu ∼ k−5/3 is well resolved in the ILES. In the DNS, Eu for kx > 80 rapidly grows with time in J01b while Eu in lower wavenumbers hardly changes with time. Fig. 5 compares Eu between the DNS performed from t = 0 (J01a) and the DNS initialized with the ILES (J01b) at the same time instance for t ≥ 17.5tr , where the ILES result at t = 17tr is also shown for comparison. In Fig. 5(a), Eu in J01b is slightly smaller than that in J01a for kx H ≥ 80, and Eu at these large wavenumbers is still growing with time in J01b. At t = 18tr and 19tr , Eu hardly differs between J01a and J01b even at the large wavenumbers. The growth of Eu indicates that small-scale motions that are missing at the initial state of J01b have developed within an interval of less than tr = H/UJ . Convergence of the transition from the ILES to DNS solutions is further studied here. Previous numerical simulations of channel flows often used a linear profile of total shear stress as a criterion of convergence from the initial condition to the fully developed turbulent state [27]. The transition from the ILES to the DNS occurs mainly at small scales even though the total shear stress in the planar jet is dominated by large scales. Therefore, the convergence criterion is established by comparing the energy spectrum with a model spectrum. Here, we use a model of a three-dimensional energy spectrum E(k) in isotropic turbulence [28] given by





Em (k ) = C ε 2/3 k−5/3 exp −β (kη )4 + cη4

1/4



+ β cη ,

(16)

6

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

Fig. 2. Passive scalar φ on a x − y plane of the planer jet in J10-LES: (a) t/tr = 2; (b) t/tr = 5; (c) t/tr = 10; (a) t/tr = 15.

Fig. 3. Cross-streamwise (y) profiles of (a) mean streamwise velocity u and (b) rms streamwise velocity fluctuation urms at t = 5tr , 10tr , and 15tr in the DNS (J01a) and LES (J01-LES).

10-2 10-3 10-4 10-5 10-6 10-7 J01-LES t = 17.0tr 10-8 J01b -9 t = 17.2tr 10 t = 17.3tr 10-10 t = 17.4tr -11 10 1 10

tained by the following relation [2]:

Eu(kx) ~ kx-5/3

Eu,m (k ) =

k

∞



1 − (k/p)2

 Em ( p ) p

dp.

(17)

The contribution of the wavenumbers unresolved in the ILES,

∞ k ≥ kLES , to u2rms can be written as k Eu (k )dk. From the sharp LES drop in Eu in the ILES caused by the filtering, kLES is estimated as 2π /(4x ). We define the ratio of the kinetic energy for k ≥ kLES between the DNS and the model as

102

Fig. 4. Energy spectra of streamwise velocity fluctuations, Eu , on the jet centerline in the LES (J01-LES) at t = 17tr and the DNS initialized with the ILES (J01b) for 17.2tr ≤ t ≤ 17.4tr .

where C = 1.5, β = 5.2, and cη = 0.4, and the subscript m denotes a quantity given by the model. Eq. (16) can be applied at large wavenumbers at large Reynolds numbers. The one-dimensional energy spectrum corresponding to the model spectrum can be ob-

∞ Eu (k )dk k R = ∞LES . kLES Eu,m (k )dk

(18)

Note that kLES is assumed to be large enough for the spectral shape for k ≥ kLES to be almost independent of flows. For such large wavenumbers, the model is expected to be valid for various turbulent flows including turbulent boundary layers and turbulent jets [28]. Even though the model spectrum for isotropic turbulence is used, the well-converged DNS is expected to have R = O(100 ). Because Eu in the ILES is much smaller than the model at large wavenumbers, R is much smaller than 1 when the DNS is started from the ILES result. Fig. 6 shows the temporal evolution of R in J01a and J01b, where kLES in J01-LES is used for calculating R. Here,

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

J01-LES (t = 17tr) 10-2

10-3

10-3

10-3

10

-4

10

-4

10-4

10

-5

10

-5

10-5

10

-6

10

-6

10-6

10

-7

10

-7

10-7

10-8

10-8

10-8

10-9

10-9

10-9

10

-10

101

10

102

J01b

-5/3

-2

10

J01a

-10

101

102

(b)

(a)

7

Eu(kx) ~ kx 10-2

10-10

101

102

(c)

Fig. 5. Comparison of Eu on the jet centerline between the DNS performed from t = 0 (J01a) and the DNS initialized with the ILES (J01b): (a) t = 17.5tr ; (b) t = 18tr ; (c) t = 19tr . Eu in the LES (J01-LES) at t = 17tr is also shown for comparison.

Fig. 8. Longitudinal auto-correlation function Ru of the streamwise velocity plotted against streamwise separation distance rx divided by the jet half width δ u . The present DNS results (J01b, J04, and J10) at t = 20tr are compared with DNS (Klein et al. [29]) and experiment (Namer & Ötügen [30]) of spatially evolving planar jets.

Fig. 6. Temporal evolutions of R on the jet centerline in J01a and J01b.

R is not equal to 1 even in the DNS because of difference between the model for isotropic turbulence and genuine turbulent flows. When the DNS is started in J01b at t/tr = 17, small-scale fluctuations are missing and R takes a small value R ≈ 0.2. As the smallscale fluctuations grow with time, R increases until R ≈ 0.7, which is close to the value in J01a. R tends to be independent of time, indicating that small-scale fluctuations have grown well. Time independence of R can be a good measure of the convergence of the transition from the ILES to the DNS. The integral time scale τI = (u2rms + v2rms + w2rms )/3ε and Kolmogorov time scale τ η are compared with the time scale at which small-scale motions grow in the DNS initialized with the ILES. Fig. 7(a) shows the temporal evolution of τ I and τ η on the centerline in the DNS (J01a). Once the turbulent planar jet has developed, τ I and τ η gradually increase with time, where τ η is of the order of O(10−1 tr ) at ReJ = 10, 0 0 0. In Figs. 4 and 5, Eu rapidly grows within

an interval of tr , which is about 10 times of τ η . The time scale of the growth of Eu , O(10τ η ), is consistent with the time scale of the chaotic synchronization of small scales in turbulence, where development of small scales also occurs within a time interval of O(10τ η ) [10]. Fig. 7(b) plots τ I , τ η , and the integration time of the DNS, τDNS = 3tr , against ReJ for J01b, J04, and J10. τ DNS is close to τ I in the present DNS, and is much larger than τ η . τ DNS ≈ τ I is long enough for Eu in small scales to develop. The spectral shape in small scales in the DNS is further investigated in Section 4.3. The DNS initialized with the ILES is validated by comparing the statistics of the temporally evolving planar jets with previous DNS and experiments of spatially evolving planar jets. Hereafter, the statistics at the end of the simulation are presented. Fig. 8 shows a longitudinal auto-correlation function of streamwise velocity

τI τη

102

101

101

100

100

10-1

10-1 0 (a)

10

10-2 0

20

50000

100000

(b)

Fig. 7. (a) Temporal evolution of integral time scale τ I and Kolmogorov time scale τ η on the jet centerline in J01a. (b) Comparison among τ I , τ η , and integration time of DNS (τ DNS ) in J01b, J04, and J10, where τ I and τ η are taken at the end of the simulations.

8

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

Present DNS: J01a; J04; J10 Stanley et al. (x/H = 12); Klein et al. (x/H = 15); Watanabe et al. (x/H = 20) DNS: Experiments: Deo et al. (x/H = 50); Gutmark & Wygnanski (x/H = 143); Watanabe et al. (x/H = 20)

(a)

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 -3

-2

-1

0

1

2

3

(b)

0 -3

0.4

0.04

0.3

0.02

0.2

0

0.1

-0.02

0 -3

-2

-1

0

1

2

3

(c)

-0.04 -3

-2

-1

0

1

2

3

-2

-1

0

1

2

3

(d)

Fig. 9. Cross-streamwise profiles of rms velocity fluctuations, (a) urms , (b) vrms , and (c) wrms , and (d) Reynolds stress u v , normalized by the mean centerline velocity uC . The DNS results of J01a, J04, and J10 are compared with previous DNS by Klein et al. (ReJ = 40 0 0) [29], Stanley et al. (ReJ = 30 0 0) [31], and Watanabe et al. (ReJ = 2200) [32] and experiments by Deo et al. (ReJ = 16, 500) [33], Gutmark and Wygnanski (ReJ = 30, 0 0 0) [34], and Watanabe et al. (ReJ = 2200) [35].

Fig. 10. Cross-streamwise profiles of mean scalar φ and rms scalar fluctuation φ rms normalized by the mean scalar on the centerline φC . The DNS results of J01a, J04, and J10 are compared with previous DNS by Watanabe et al. (ReJ = 2200 and Sc = 1) [36] and da Silva et al. (ReJ = 60 0 0 and Sc = 0.7) [67] and experiments (ReJ = 2200 and Sc ≈ 600) [35].

fluctuations, Ru , calculated with streamwise separation distance rx on the jet centerline. When rx is normalized by the jet half width δ u based on the mean velocity profile, Ru at different streamwise locations tends to collapse onto a single curve [29]. Ru in the present DNS agrees well with those in spatially evolving planar jets. Ru reaches 0 around rx /δu = 2, which is about 1/4 of the computational domain size in the x direction. Fig. 9 shows rms velocity fluctuations urms , vrms , and wrms and Reynolds stress u v  normalized by the mean centerline velocity uC , where the coordinate y is normalized by the jet half width δ u . Fig. 10 shows the mean scalar φ and rms scalar fluctuation φ rms normalized by the mean scalar on the centerline φC , where the coordinate is normalized by the jet half width of the mean scalar profile, δ φ . The profiles of these statistics in the present DNS agree

with previous studies of spatially evolving planar jets in a selfsimilar region. The present DNS exhibits weak Reynolds number dependence while the spatially evolving planar jets also show a similar level of variations of these statistics in Figs. 9 and 10. These variations can also be related to the inflow condition at the nozzle used in the experiments [37]. The 1st and 2nd order statistics shown in Figs. 9 and 10 are dominated by large-scale motions of the jet. Further comparison of statistics related to small-scale characteristics is presented in Section 4.3 along with the results from the DNS of turbulent boundary layers. Here, how the second-order statistics depend on the Reynolds number is different between experiments and DNS. The flow inside the jet nozzle is also influenced by the Reynolds number, while the numerical simulations are often performed without including the

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

9

Fig. 11. Passive scalar φ on a x − y plane of the turbulent boundary layer in B13-LES: (a) t/tr = 12 (Reθ = 501 ); (b) t/tr = 25 (Reθ = 1007); (c) t/tr = 230 (Reθ = 2001); (d) t/tr = 1024 (Reθ = 4016); (e) t/tr = 1920 (Reθ = 60 0 0); (f) t/tr = 3378 (Reθ = 90 0 0).

nozzle by assuming the velocity profile at the nozzle exit. This difference can be the reason why the Reynolds number dependence differs between simulations and experiments. In temporal simulations, the size of the computational domain can affect the statistical convergence, which is investigated by repeating the simulation of J01a with different sets of random numbers used for generating the initial velocity fluctuations. Ensemble average and standard deviation of urms /UJ at y = 0 and t = 20tr among four simulations are 0.11 and 0.0064, respectively, where the standard deviation is about 5% of the ensemble average. The rms velocity computed from a single simulation can contain statistical uncertainty of about 5%. The statistics obtained by simulations might be influenced by various factors: numerical schemes, the ILES field used in the DNS, and degree of statistical convergence. It seems that switching from the ILES to the DNS hardly affects the statistics as shown in Fig. 5(c). The impact of the initialization method on the statistics might be quantified and further studied by uncertainty quantification techniques [38].

4.2. Temporally evolving turbulent boundary layer Fig. 11 shows passive scalar φ on a x − y plane in the temporally evolving boundary layer in B13-LES. The turbulent boundary layer develops with time, and larger-scale structures appear in later time. Even though Fig. 11 is taken from the ILES, one can see that the most length scales of the turbulent boundary layer are resolved by the ILES, and only very small scales are in the subgrid scale. Fig. 12 (a) compares temporal evolutions of Reθ and the viscous length ν /uτ between the DNS (B02a) and the LES (B02LES). The LES results are almost identical to the DNS results. Once Reθ exceeds 10 0 0, ν /uτ tends to be constant. Fig. 12(b) plots the 2 /2 ) as a function of Re skin friction coefficient C f = τW /(ρUW θ in B02a and B02-LES. For Reθ ≥ 1200, Cf follows the friction law

/4 C f = 0.024Re−1 , indicating that the turbulent boundary layer is θ fully developed in both DNS and ILES. Fig. 13 (a) shows the mean streamwise velocity profile U + = u/uτ in B02a (DNS) and B02-LES, where y+ = y/(ν /uτ ) is the distance from the wall normalized by the viscous length. The DNS results of spatially evolving turbulent boundary layer [39] at Reθ = 1986 are also shown in the figure. U + follows the velocity profiles of the viscous sublayer and the log-low region. The present DNS and ILES provide the mean velocity profile similar to that in the spatially evolving turbulent boundary layer. Fig. 13(b) shows the energy spectrum of streamwise velocity fluctuation, Eu , as a function of streamwise wavenumber kx in B02a and B02-LES at two vertical locations y/δ = 0.04 and 0.5, which are equivalent to y+ = 28 and 357, respectively. Both DNS and ILES yield a similar profile of Eu at low wavenumbers while Eu in the LES rapidly decreases for kx D ≥ 10, which is not resolved well in the ILES. The effects of the low-pass filter applied as the implicit SGS model are examined by conducting the LES without filtering, where the numerical parameters and conditions except the filter are the same as in B02-LES. Fig. 13(b) also shows Eu at y/δ = 0.04 in the LES without filtering. Because the energy dissipation rate is underestimated in the LES, the energy accumulates at large wavenumbers if the low-pass filter is not applied. Therefore, the spectrum in the LES without filtering has larger values for 5 ≤ kx D ≤ 10 than those in B02a and B02-LES. The filter introduces the numerical energy dissipation that can suppress unphysical growth of Eu at large wavenumbers. Fig. 14 (a) compares the integration time of the DNS, τ DNS , integral time scale τ I , and Kolmogorov time scale τ η as a function of Reθ in all DNS initialized with the LES. The results are taken at y/δ = 0.35 at the end of the simulations. τ DNS is much larger than τ η , and the integration time is comparable to the integral time scale similarly in the DNS of turbulent planar jets. Fig. 14(b) shows Eu in B02a-d at the end of the simulations, where B02a is performed from t = 0. The DNS initialized with the LES (B02b, B02c,

10

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

0.04

0.008

0.03

0.006

0.02

0.004

0.01

0.002

B02a (DNS) B02-LES

2000

1000 B02a (DNS) B02-LES (a

0 0

100

200

0

0 b)

500

1000

1500

2000

Fig. 12. (a) Temporal evolution of the Reynolds number based on momentum thickness Reθ and the viscous length ν /uτ in B02a and B02-LES. (b) Skin friction coefficient /4 . plotted against Reθ in comparison with a friction law of turbulent boundary layers C f = 0.024Re−1 θ

Fig. 13. (a) Mean streamwise velocity at Reθ = 20 0 0 in B02a and B02-LES. The present numerical simulations are compared with the DNS results of spatially developing turbulent boundary layer (Reθ = 1986) [39]. The figure also shows the velocity profiles for the viscous sublayer, U + = y+ , and for the log-law region, U + = (1/κ )lny+ + A with κ = 0.41 and A = 5.1. (b) Energy spectra of streamwise velocity fluctuations, Eu , at y/δ = 0.04 and 0.5 in B02a and B02-LES. Eu at y/δ = 0.04 obtained from the LES without filtering is also shown for comparison.

10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10

103

102

101

100

0

5000

10000

15000

(a

b)

B02a B02b B02c B02d

100

101

Fig. 14. (a) Integration time of the DNS, τ DNS , integral time scale τ I , and Kolmogorov time scale τ η at the end of the DNS. The results are shown against Reθ at the end of the simulation, and τ I and τ η are taken at y/δ = 0.35. (b) Energy spectra of streamwise velocity fluctuations, Eu , at y/δ = 0.04 and 0.5 in B02a, B02b, B02c, and B02d.

and B02d) yields Eu consistent with B02a even at large wavenumbers. Although B02b, B02c, and B02d are started at different time instances of the ILES, Eu at large wavenumbers is identical at the end of the simulations. It is noteworthy that τ DNS is smaller than the integral time scale in B02d, and Eu at large wavenumbers is recovered in the DNS within the integral time scale. Thus, the integration time in the DNS is long enough for small-scale motions to develop from the initial state. R defined as Eq. (18) is also computed at y/δ = 0.5. At Reθ = 20 0 0, R = 0.63, 0.64, 0.63, and 0.62 are obtained for B02a, B02b, B02c, and B02d, respectively. Thus, R in B02a (DNS) is also close to those in other cases, and R is not affected by the initial condition provided by the ILES. R hardly changes with the integration time in B02b, B02c, and B02d. These results also confirm that the transition from the LES to the DNS has converged in these simulations. The smaller value of R in the boundary layer than in the jet can be related to accuracy of the

model spectrum because the turbulent boundary layer (B02) is more anisotropic and has a lower turbulent Reynolds number. The DNS initialized with the ILES is further validated by comparing the statistics with experiments and numerical simulations of spatially evolving turbulent boundary layers conducted at similar Reθ . Hereafter, the statistics at the end of the simulations in B02c, B04, B06, B08, B10, and B13 are presented. Fig. 15 shows the mean streamwise velocity U + in comparison with the experiments at Reθ = 2260 (Erm and Joubert [40]) and Reθ = 130 0 0 (De Graaff and Eaton [41]). U + in these experiments is similar to those in B02c and B13. There is a range of y+ for which the DNS results follow the log law. The log-law region is extended to larger y+ for higher Reθ as expected from previous studies [41,47]. Figs. 16 (a) and (b) plot Cf and the shape factor H = δ ∗ /θ against Reθ , respectively, while Fig. 17 shows the relation between Reτ

∞ and Reθ , where δ ∗ = 0 u/UW dy is the displacement thickness.

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

11

Fig. 15. Mean streamwise velocity U + = u/uτ in the DNS (B02c, B04, B06, B08, B10, and B13) at the end of the simulations. The DNS results are compared with experiments by Erm and Joubert [40] (Reθ = 2260) and De Graaff and Eaton [41] (Reθ = 13, 0 0 0). The figure also shows the velocity profiles for the viscous sublayer, U + = y+ , and for the log-law region, U + = (1/κ )lny+ + A with κ = 0.41 and A = 5.1.

Fig. 16. (a) Skin friction coefficient Cf and (b) shape factor H plotted against Reθ in the DNS (B02c, B04, B06, B08, B10, and B13) at the end of the simulations. In (a), the

/4 , and the present results are compared with experiments by De Graaff and Eaton [41] and Österlund et al. [42], a friction law of turbulent boundary layers C f = 0.024Re−1 θ Coles–Fernholz best fit from Nagib et al. [43], C f = 2[(1/0.384 )ln(Reθ ) + 4.127]−2 . (b) also shows results from experiments by De Graaff and Eaton [41], zonal detached eddy simulations by Deck et al. [44], and the relation by Coles [45].

Fig. 17. Relation between Reτ and Reθ in the DNS (B02c, B04, B06, B08, B10, and B13) and experiments by De Graaff and Eaton [41] and Smith [46].

The present DNS results are compared with previous experimental and numerical studies as shown in the captions of the figures [41–45]. At high Reθ , Cf tends to be larger than the friction law

/4 C f = 0.024Re−1 . This tendency can be found in both present DNS θ and previous experiments. The present DNS initialized with the ILES yields the Reθ dependence of these quantities consistent with previous studies. Reτ reaches 4 × 103 in the present DNS with the highest Reθ (B13). + Figs. 18, 19, and 20 compare the 2nd-order statistics, u+ rms , vrms ,  v + , between the present DNS and previous experiw+ , and  u rms 2 1/2 /u , v+ = ments. Here, the superscript + represents u+ τ rms = u  rms +  2 1 / 2  2 1 / 2   +   v  /uτ , wrms = w  /uτ , and u v  = u v /u2τ . Profiles of

the 2nd-order statistics at different Reθ are well reproduced in the present DNS. In order to examine the statistical convergence, the simulation of B02c is repeated four times. Ensemble average and standard deviation of the peak value of u+ rms are 2.80 and 6.99 × 10−3 among four simulations, respectively. On the other hand, en+ semble average and standard deviation of u+ rms at y = 400 (the −2 outer region) are 1.53 and 6.02 × 10 , respectively. The standard −3 deviation divided by the ensemble average for u+ rms is 2.50 × 10 + −2 + for the peak value of urms and is 3.93 × 10 at y = 400. + A peak value of u+ rms , denoted by [urms ]max , appears around y+ = 16. It is known that [u+ ] increases with the Reynolds max rms number [49]. Fig. 21 shows [u+ rms ]max as a function of Reθ in comparison with experimental results by De Graaff and Eaton [41] and the logarithmic curve fit [u+ rms ]max = 1.84 + 0.29log10 Reθ obtained by Metzger et al. [49] [u+ rms ]max in the present DNS increases with Reθ as also found in previous studies. These comparisons with spatially evolving turbulent boundary layers confirm that the hybrid ILES/DNS methodology captures well the Reynolds number dependence of turbulent boundary layers. The equation for the turbulent kinetic energy kT = ui ui /2 can be written as

∂ kT ∂ kT + u j  = Pk + Tk + k + Dk + ε ∂t ∂xj

(19)

with

Pk = −ui uj 

1 ∂ ∂ui  ∂    ,T = − u u u /2, k = −  p uj , ∂xj k ∂xj j i i ρ ∂xj (20)

12

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

Fig. 18. Vertical profiles of rms velocity fluctuations and Reynolds stress at (a) Reθ = 20 0 0 in B02c and (b) Reθ = 40 0 0 in B04. Experimental results at Reθ = 2100 and 4400 (Osaka et al. [47]) are also shown in (a) and (b), respectively.

Fig. 19. Vertical profiles of rms velocity fluctuations and Reynolds stress at (a) Reθ = 60 0 0 in B06 and (b) Reθ = 80 0 0 in B08. Experimental results at Reθ = 6040 (Osaka et al. [47]) and Reθ = 8100 (Carlier and Stanislas [48]) are also shown in (a) and (b), respectively.

Fig. 20. Vertical profiles of rms velocity fluctuations and Reynolds stress at (a) Reθ = 10, 0 0 0 in B10 and (b) Reθ = 13, 0 0 0 in B13. Experimental results at Reθ = 11, 500 (Carlier and Stanislas [48]) and Reθ = 13, 0 0 0 (De Graaff and Eaton [41]) are also shown in (a) and (b), respectively.



Dk = ν

∂ u

∂ u



∂ kT i i , ε = −ν , ∂ x j∂ x j ∂xj ∂xj

(21)

which are called the production, turbulent diffusion, pressure diffusion, viscous diffusion, and dissipation terms, respectively. Fig. 22 shows the turbulent kinetic energy budget in B02d in comparison with previous DNS study [50]. The DNS initialized with the ILES well predicts the turbulent kinetic energy budget near the wall. The dissipation term in B02a is also shown in the figure, and the effects of the initial field provided by the ILES do not appear for the turbulent kinetic energy dissipation rate. The dissipation rate is related to the energy spectrum of the velocity derivative ∂ u/∂ x,

which is written as k2x Eu . Small-scale velocity fluctuations make large contribution to the variance of velocity derivative. Even for large wavenumbers, Eu hardly differs between B02a and B02d in Fig. 14(b), indicating that the DNS initialized with the ILES predicts well the small-scale contribution to the dissipation rate. This is also the case in the simulations of the planar jet, where the effects of the initial conditions generated by the ILES do not appear in Eu at small scales in Fig. 5(c). 4.3. Statistics related to small-scale motions of turbulence Statistical features of small scales are examined for both planar jets and boundary layers. The results in this subsection are

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

4

3

Present DNS De Graaff & Eaton 2 103

104

+ Fig. 21. Peak value of u+ rms , [urms ]max , plotted as a function of Reθ . Experimental results by De Graaff and Eaton [41] and the logarithmic curve fit [u+ rms ]max = 1.84 + 0.29log10 Reθ by Metzger et al. [49] are shown for comparison.

13

Reynolds number Reλ = urms λx /ν, where λx = urms /(∂ u /∂ x )2 1/2 is the Taylor microscale. The present DNS results agree well with previous studies. Fig. 24(a) also includes the relation of (Reλ , −S∂ u/∂ x ) obtained in the ILES. Mi et al. computed −S∂ u/∂ x from time series data of streamwise velocity with Taylor’s hypothesis in a turbulent round jet, and found that −S∂ u/∂ x computed with low-pass filtered velocity is smaller than that of unfiltered velocity [53]. The ILES also does not contain small-scale fluctuations. Therefore, the ILES yields smaller skewness than the DNS in Fig. 24(a). However, (Reλ ,-S∂ u/∂ x ) in the DNS initialized with the ILES exhibits good agreement with other studies. Fig. 24(b) shows the relation of (Reλ , F∂ u/∂ x ). F∂ u/∂ x increases with Reλ . The present DNS well predicts this tendency. Furthermore, flatness of the passive scalar derivative, F∂ φ /∂ x = (∂ φ  /∂ x )4 /(∂ φ  /∂ x )2 2 , is plotted as a function of Reλ in Fig. 25. F∂ φ /∂ x has stronger dependence on Reλ than the flatness of the velocity derivative. The present DNS data also follow the trend obtained from the previous studies. The DNS initialized by the ILES well recovers the features of these small-scales statistics. These results confirm that the smallscale fluctuations, which are missing at the beginning of the DNS, fully develop within τ DNS comparable to large-scale eddy turnover time τ I . 4.4. Computation time of DNS initialized by ILES

Fig. 22. Turbulent kinetic energy budget near the wall in B02d compared with DNS results by Spalart [50] and the dissipation term in B02a. The terms in Eq. (19) are normalized by the viscous scales (u4τ /ν ).

105 104 10

Present DNS J01b J04 J10 B02c B04 B06 B08 B10 B13 Experiments Sadeghi et al.

3

102 101 100 10-1 10-2 10-3 -3 10

10-2

10-1

100

Fig. 23. Normalized energy spectra of the streamwise velocity fluctuation Eu /(εν 5 )1/4 plotted as a function of kx η. The DNS results are taken on the centerline of the jet and at y/δ = 0.5 of the boundary layer at the end of the simulations. The experimental result of a turbulent round jet by Sadeghi et al. [51] is also shown /3 with C1 = 0.49. for comparison. The dashed-dotted line represents Eu = C1 ε 2/3 k−5 x

shown for the centerline of the jet and for y/δ = 0.5 of the boundary layer at the end of the simulations. Fig. 23 shows normalized energy spectra of the streamwise velocity fluctuation, Eu /(εν 5 )1/4 , as a function of kx η. The normalized spectra were shown to be universal except for large scales [1]. The present DNS results are consistent with the experimental result of a turbulent round jet [51] and the theoretical prediction for the inertial subrange [1], /3 Eu = C1 ε 2/3 k−5 with C1 = 0.49. The spectral shape in small scales x is well predicted by the present DNS. Derivative skewness and flatness are computed as S∂ u/∂ x = (∂ u /∂ x )3 /(∂ u /∂ x )2 3/2 and F∂ u/∂ x = (∂ u /∂ x )4 /(∂ u /∂ x )2 2 , respectively. Fig. 24(a) shows −S∂ u/∂ x as a function of the turbulent

CPU time of the DNS initialized by the ILES is computed as hL,D = hL + hD from the CPU time spent for the ILES and DNS denoted by hL and hD , respectively. hD is also used for estimating CPU time required for the full DNS as hF D = hD (tend /τDNS ), where the CPU time is assumed to be proportional to the integration time of simulations. Table 6 summarizes hL,D , hFD , and hL,D /hFD . The initialization with the ILES can save more than 80% of the CPU time in the case of the planar jet. For the turbulent boundary layer, hL,D /hFD tends to decrease with the Reynolds number. This is because longer computation is required until the turbulent boundary layer reaches a higher target value of Reθ . 2–3 times more grid points in each direction are used in the DNS compared with the ILES. Therefore, the total number of grid points is different by a factor of 8 − 27 between the DNS and ILES. Hence, even though the ILES used in this study resolves most length scales of turbulence, usage of ILES in initial parts of the simulations can significantly reduce the CPU time. 4.5. Application hybrid ILES/DNS to spatial simulations and comparison with synthetic inflow methods. The present approach can also be used for spatial simulations. A computational domain can be divided into upstream and downstream regions, which are simulated by LES and DNS, respectively. The flow simulated by the LES is used as the inflow boundary condition at the boundary between the DNS and LES regions. The present study found that small-scale fluctuations grow in the DNS within an interval of the integral time scale τ I . Therefore, the transition from the LES to the DNS solutions is expected to occur within streamwise distance of Uc τ I , where Uc is the convective velocity of the flow. An inflow-boundary or initial condition of the DNS can be given by a synthetic velocity, which may be generated by superimposing artificial velocity fluctuations on a mean velocity profile [63]. The synthetic velocity generators have a lower computational cost than LES because they often convert random noise into artificial turbulent velocity fluctuations without temporal advancement of Navier–Stokes equations [64–66]. However, the well-resolved LES yields more realistic velocity fluctuations for a wide range of length scales than most synthetic velocity generators. Therefore, it

14

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314 Table 6 Computation (CPU) time (hours) of full DNS (hFD ) and DNS initialized with ILES (hL,D ). Case

J01

J04

J10

B02

B04

B06

B08

B10

B13

hL,D [h] hFD [h] hL,D /hFD

109 681 0.16

5138 29161 0.18

65678 391711 0.17

174 459 0.33

8145 20912 0.39

27809 88346 0.31

112766 396043 0.25

275886 1104058 0.25

804716 3627271 0.22

Fig. 24. (a) Skewness and (b) flatness of velocity derivative ∂ u /∂ x plotted as a function of Reλ . The DNS results are taken on the centerline of the jet and at y/δ = 0.5 of the boundary layer at the end of the simulations. The results from the LES of J01-LES, J04-LES, J10-LES, B02-LES, B10-LES, and B13-LES are also shown in (a) for comparison. Figure includes the results by Kerr [52], Mi et al. [53], Antonia and Chambers [54], Mydlarski and Warhaft [55], Sreenivasan and Antonia [56], and Kuo and Corrsin [57]. Some of these results were compiled by Sreenivasan and Antonia [56].

Fig. 25. Flatness of scalar derivative ∂ φ  /∂ x plotted as a function of Reλ . Figure includes the results by Antonia and Chambers [54], Antonia and Danh [58], Antonia and Van Atta [59], Kerr [52], Gibson et al. [60], Gibson and Masiello [61], and McConnell [62]. These results were also compiled by Sreenivasan and Antonia [56].

is expected that the influence of initial/inflow conditions disappears within a shorter time in the DNS combined with the LES than with the synthetic velocity generator. 5. Conclusion Direct numerical simulations (DNS) initialized with implicit large eddy simulations (ILES) are tested for the temporally evolving planar jets and turbulent boundary layers. The ILES is performed with a very fine resolution, where most scales of turbulent motions are resolved except for very small scales close to the Kolmogorov length scale. The ILES is performed until the flow develops into a fully turbulent state. Subsequently, the ILES results are used as the initial condition of the DNS. As time is advanced in the DNS, the energy spectrum in small scales, which are not resolved in the ILES, rapidly increases with time without changing the spectrum at large scales. Within an interval of the integral time scale, the spectrum becomes almost identical for the DNS initialized with the ILES and the DNS started from t = 0. The DNS initialized with the ILES also well reproduces well-known small-scale characteristics, such as Reynolds number dependence of skewness and flatness of velocity and scalar derivatives and energy spectrum of ve-

locity fluctuations in the inertial subrange and viscous range. The 1st and 2nd order statistics and longitudinal correlation function of the planar jets are also well predicted by the present DNS. The DNS initialized with the ILES yields Reθ dependence of the mean velocity, rms velocity fluctuation, Reynolds stress, shape factor, and skin friction in the turbulent boundary layers, which agrees well with previous studies. For example, a peak value of the rms streamwise velocity fluctuation increases with Reθ as observed in experiments. Even though most length scales of turbulence are resolved in the ILES, the hybrid ILES/DNS significantly reduces the computational cost. These investigations confirm that DNS combined with ILES is useful to study turbulent flows at high Reynolds numbers at a limited computational resource. Acknowledgment The numerical simulations presented in this manuscript were carried out on the high-performance computing system (NEC SXACE) in the Japan Agency for Marine-Earth Science and Technology. This work was partially supported by “Collaborative Research Project on Computer Science with High-Performance Computing in Nagoya University” and by JSPS KAKENHI grant number 18K13682 and 18H01367. Appendix A. Generation of initial velocity fluctuations In LES and DNS of turbulent boundary layers and planar jets, the successful transition from an initial laminar state to turbulence requires the initial or inflow velocity to contain artificial fluctuations or turbulent fluctuations. The simplest method for generating velocity fluctuations is based on random numbers that give each component of velocity fluctuating vector. However, the velocity fluctuations generated in this way have a small spatialcorrelation length scale, and a large amount of turbulent kinetic energy is contained in very small scales [64]. This results in a rapid decay of velocity fluctuations. Therefore, some methods for generating velocity fluctuations for inflow boundary or initial conditions rely on spatially correlated random fluctuations [64,65]. In the present study, a relatively simple approach is used for

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

Fig. 26. Premultiplied energy spectrum kx Eu of initial velocity fluctuations in the simulations of the planar jet, where kx is the wavenumber in the streamwise direction. Eu normalized by the rms streamwise velocity fluctuation is plotted against the wavelength x = 2π /kx divided by H.

generating spatially correlated fluctuations. The velocity fluctuations generated by random numbers are spatially correlated over the length of grid size. Therefore, a coarse grid is defined for generating the initial fluctuation with Nxi × Nyi × Nzi points. The number of the grid points Nxi × Nyi × Nzi are determined such that the grid spacing is in ≈ 0.02H near the centerline of the planar jets and in ≈ D near the wall in the boundary layers. The random numbers between −1 and 1 are assigned onto three components of the velocity fluctuation vector at all grid points. Then, the fluctuations are interpolated onto the grid used for simulations. Finally, the fluctuations are normalized to fulfill the requirement for rms velocity fluctuations in simulations. This method yields the velocity fluctuations with a spatial correlation over the length of in . Fig. 26 shows the premultiplied energy spectrum of initial u in one of the simulations of the planar jets. A peak appears around the wavelength of 0.1H, which is much larger than the grid size used in the simulation in Table 3. References [1] Pope SB. Turbulent flows. Cambridge Univ. Pr.; 20 0 0. [2] Davidson PA. Turbulence: an introduction for scientists and engineers. Oxford Univ. Pr.; 2004. [3] Rogers MM, Moser RD. Direct simulation of a self-similar turbulent mixing layer. Phys Fluids 1994;6(2):903–23. [4] da Silva CB, Métais O. On the influence of coherent structures upon interscale interactions in turbulent plane jets. J Fluid Mech 2002;473:103–45. [5] Wu X, Moin P. Direct numerical simulation of turbulence in a nominally zero– pressure-gradient flat-plate boundary layer. J Fluid Mech 2009;630:5–41. [6] Moser RD, Rogers MM, Ewing DW. Self-similarity of time-evolving plane wakes. J Fluid Mech 1998;367:255–89. [7] Cifuentes L, Dopazo C, Martin J, Domingo P, Vervisch L. Effects of the local flow topologies upon the structure of a premixed methane-air turbulent jet flame. Flow Turbul Combust 2016;96(2):535–46. [8] Diamessis PJ, Spedding GR, Domaradzki JA. Similarity scaling and vorticity structure in high-Reynolds-number stably stratified turbulent wakes. J Fluid Mech 2011;671:52–95. [9] Watanabe T, Riley JJ, de Bruyn Kops SM, Diamessis PJ, Zhou Q. Turbulent/non-turbulent interfaces in wakes in stably stratified fluids. J Fluid Mech 2016;797:R1. [10] Lalescu CC, Meneveau C, Eyink GL. Synchronization of chaos in fully developed turbulence. Phys Rev Lett 2013;110(8):084102. [11] da Silva CB, Pereira JCF. Invariants of the velocity-gradient, rate-of-strain, and rate-of-rotation tensors across the turbulent/nonturbulent interface in jets. Phys Fluids 2008;20(5):055101. [12] Kozul M, Chung D, Monty JP. Direct numerical simulation of the incompressible temporally developing turbulent boundary layer. J Fluid Mech 2016;796:437–72. [13] Watanabe T, Zhang X, Nagata K. Turbulent/non-turbulent interfaces detected in DNS of incompressible turbulent boundary layers. Phys Fluids 2018;30(3):035102. [14] Bogey C, Bailly C. Turbulence and energy budget in a self-preserving round jet: direct evaluation using large eddy simulation. J Fluid Mech 2009;627:129–60. [15] Watanabe T, Nagata K. Gradients estimation from random points with volumetric tensor in turbulence. J Comput Phys 2017;350:518–29. [16] Watanabe T, Riley JJ, Nagata K, Onishi R, Matsuda K. A localized turbulent mixing layer in a uniformly stratified environment. J Fluid Mech 2018;849:245–76.

15

[17] Tanaka S, Watanabe T, Nagata K. Multi-particle model of coarse-grained scalar dissipation rate with volumetric tensor in turbulence. J Comput Phys 2019;389(15):128–46. [18] Morinishi Y, Lund TS, Vasilyev OV, Moin P. Fully conservative higher order finite difference schemes for incompressible flow. J Comput Phys 1998;143(1):90–124. [19] Van der Vorst HA. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J Sci Stat Comput 1992;13(2):631–44. [20] Kennedy CA, Carpenter MH. Several new numerical methods for compressible shear-layer simulations. Appl Numer Math 1994;14(4):397–433. [21] van Reeuwijk M, Holzner M. The turbulence boundary of a temporal jet. J Fluid Mech 2014;739:254–75. [22] Watanabe T, da Silva CB, Sakai Y, Nagata K, Hayase T. Lagrangian properties of the entrainment across turbulent/non-turbulent interface layers. Phys Fluids 2016;28(3):031701. [23] Zhang X, Watanabe T, Nagata K. Turbulent/non-turbulent interfaces in high resolution direct numerical simulation of temporally-evolving compressible turbulent boundary layers. Phys Rev Fluids 2018. Accepted for publication [24] Schlatter P, Örlü R. Assessment of direct numerical simulation data of turbulent boundary layers. J Fluid Mech 2010;659:116–26. [25] Lozano-Durán A, Jiménez J. Effect of the computational domain on direct simulations of turbulent channels up to Reτ = 4200. Phys Fluids 2014;26(1):011702. [26] Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel flow up to Reτ = 590. Phys Fluids 1999;11(4):943–5. [27] Vinuesa R, Prus C, Schlatter P, Nagib HM. Convergence of numerical simulations of turbulent wall-bounded flows and mean cross-flow structure of rectangular ducts. Meccanica 2016;51(12):3025–42. [28] Pope SB. PDF methods for turbulent reactive flows. Prog Energy Combust Sci 1985;11(2):119–92. [29] Klein M, Sadiki A, Janicka J. Investigation of the influence of the Reynolds number on a plane jet using direct numerical simulation. Int J Heat Fluid Flow 2003;24(6):785–94. [30] Namer I, Ötügen MV. Velocity measurements in a plane turbulent air jet at moderate Reynolds numbers. Exp Fluids 1988;6(6):387–99. [31] Stanley SA, Sarkar S, Mellado JP. A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation. J Fluid Mech 2002;450:377–407. [32] Watanabe T, Sakai Y, Nagata K, Ito Y, Hayase T. Vortex stretching and compression near the turbulent/nonturbulent interface in a planar jet. J Fluid Mech 2014;758:754–85. [33] Deo RC, Nathan GJ, Mi J. Similarity analysis of the momentum field of a subsonic, plane air jet with varying jet-exit and local Reynolds numbers. Phys Fluids 2013;25(1):015115. [34] Gutmark E, Wygnanski I. The planar turbulent jet. J Fluid Mech 1976;73:465–95. [35] Watanabe T, Sakai Y, Nagata K, Terashima O. Experimental study on the reaction rate of a second-order chemical reaction in a planar liquid jet. AIChE J 2014;60(11):3969–88. [36] Watanabe T, Sakai Y, Nagata K, Ito Y, Hayase T. Enstrophy and passive scalar transport near the turbulent/non-turbulent interface in a turbulent planar jet flow. Phys Fluids 2014;26(10):105103. [37] Wu N, Sakai Y, Nagata K, Suzuki H, Terashima O, Hayase T. Analysis of flow characteristics of turbulent plane jets based on velocity and scalar fields using DNS. J Fluid Sci Tech 2013;8(3):247–61. [38] Rezaeiravesh S, Vinuesa R, Liefvendahl M, Schlatter P. Assessment of uncertainties in hot-wire anemometry and oil-film interferometry measurements for wall-bounded turbulent flows. Eur J Mech-B/Fluids 2018;72:57–73. [39] Jiménez J, Hoyas S, Simens MP, Mizuno Y. Turbulent boundary layers and channels at moderate Reynolds numbers. J Fluid Mech 2010;657:335. [40] Erm LP, Joubert PN. Low-Reynolds-number turbulent boundary layers. J Fluid Mech 1991;230:1–44. [41] De Graaff DB, Eaton JK. Reynolds-number scaling of the flat-plate turbulent boundary layer. J Fluid Mech 20 0 0;422:319–46. [42] Österlund JM, Johansson AV, Nagib HM, Hites MH. A note on the overlap region in turbulent boundary layers. Phys Fluids 20 0 0;12(1):1–4. [43] Nagib HM, Chauhan KA, Monkewitz PA. Approach to an asymptotic state for zero pressure gradient turbulent boundary layers. Phil Trans R Soc A 2007;365(1852):755–70. [44] Deck S, Renard N, Laraufie R, Weiss P-É. Large-scale contribution to mean wall shear stress in high-Reynolds-number flat-plate boundary layers up to Reθ = 13650. J Fluid Mech 2014;743:202–48. [45] Coles D. The turbulent boundary layer in a compressible fluid. Tech Rep; 1962. R–403–PR. United States Air Force Project RAND. [46] Smith RW. Effect of Reynolds number on the structure of turbulent boundary layers. PhD thesis, Princeton University 1994. [47] Osaka H, Kameda T, Mochizuki S. Re-examination of the Reynolds-number-effect on the mean flow quantities in a smooth wall turbulent boundary layer. JSME Int J, Ser-B 1998;41(1):123–9. [48] Carlier J, Stanislas M. Experimental study of eddy structures in a turbulent boundary layer using particle image velocimetry. J Fluid Mech 2005;535:143–88. [49] Metzger MM, Klewicki JC, Bradshaw KL, Sadr R. Scaling the near-wall axial turbulent stress in the zero pressure gradient boundary layer. Phys Fluids 2001;13(6):1819–21.

16

T. Watanabe, X. Zhang and K. Nagata / Computers and Fluids 194 (2019) 104314

[50] Spalart PR. Direct simulation of a turbulent boundary layer up to reθ = 1410. J Fluid Mech 1988;187:61–98. [51] Sadeghi H, Lavoie P, Pollard A. The effect of Reynolds number on the scaling range along the centreline of a round turbulent jet. J Turbulence 2014;15(6):335–49. [52] Kerr RM. Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence. J Fluid Mech 1985;153:31–58. [53] Mi J, Xu M, Zhou T. Reynolds number influence on statistical behaviors of turbulence in a circular free jet. Phys Fluids 2013;25(7):075101. [54] Antonia RA, Chambers AJ. On the correlation between turbulent velocity and temperature derivatives in the atmospheric surface layer. Boundary-Layer Meteorol 1980;18(4):399–410. [55] Mydlarski L, Warhaft Z. On the onset of high-Reynolds-number grid-generated wind tunnel turbulence. J Fluid Mech 1996;320(1):331–68. [56] Sreenivasan KR, Antonia RA. The phenomenology of small-scale turbulence. Annu Rev Fluid Mech 1997;29(1):435–72. [57] Kuo AY-S, Corrsin S. Experiments on internal intermittency and fine-structure distribution functions in fully turbulent fluid. J Fluid Mech 1971;50(2):285–319. [58] Antonia RA, Danh HQ. Structure of temperature fluctuations in a turbulent boundary layer. Phys Fluids 1977;20(7):1050–7. [59] Antonia RA, Van Atta CW. On the correlation between temperature and velocity dissipation fields in a heated turbulent jet. J Fluid Mech 1975;67(2):273–88.

[60] Gibson CH, Stegen GR, Williams RB. Statistics of the fine structure of turbulent velocity and temperature fields measured at high Reynolds number. J Fluid Mech 1970;41(1):153–67. [61] Gibson CH, Masiello PJ. Observations of the validity of dissipation rates of turbulent velocity and temperature fields. Lecture Notes Phys Stat Models Turbul 1972;12(427):53. [62] McConnell SO. The fine structure of velocity and temperature measured in the laboratory and the atmospheric marine boundary layer. University of California; 1976. PhD thesis. [63] Martín MP. Direct numerical simulation of hypersonic turbulent boundary layers. part 1. initialization and comparison with experiments. J Fluid Mech 2007;570:347–64. [64] Klein M, Sadiki A, Janicka J. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations. J Comput Phys 2003;186(2):652–65. [65] Kempf A, Klein M, Janicka J. Efficient generation of initial-and inflow-conditions for transient turbulent flows in arbitrary geometries. Flow Turbul Combust 2005;74(1):67–84. [66] Wu X. Inflow turbulence generation methods. Annu Rev Fluid Mech 2017;49:23–49. [67] da Silva CB, Lopes DC, Raman V. The effect of subgrid-scale models on the entrainment of a passive scalar in a turbulent planar jet. J. Turbulence 2015;16(4):342–66.