Gas kinetic scheme for turbulence simulation

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.1 (1-14)

Aerospace Science and Technology ••• (••••) •••–•••

1

Contents lists available at ScienceDirect

2

67 68

3

Aerospace Science and Technology

4 5

69 70 71

6

72

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13 14 15 16

Gas kinetic scheme for turbulence simulation

77

Shuang Tan, Qibing Li ∗ , Zhixiang Xiao, Song Fu

79 80 81

AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China

82

17 18

83

a r t i c l e

i n f o

84

a b s t r a c t

19 20 21 22 23 24 25 26 27 28 29 30

78

85

Article history: Received 3 November 2017 Received in revised form 19 March 2018 Accepted 13 April 2018 Available online xxxx Keywords: GKS RANS model LES model Hybrid RANS/LES model Multiscale method Hypersonic turbulence

31 32 33

The extended gas-kinetic scheme (GKS) for turbulence simulations is developed based on the generalized BGK equation with the effective relaxation time. This relaxation time can be computed from the turbulent viscosity, through which turbulence models can be directly combined. For engineering lowcost simulations of high Reynolds number flows, common-used RANS models are applied, while the LES and hybrid RANS/LES (DES and IDDES) models as well as the minimized dispersion and controllable dissipation (MDCD) reconstruction are adopted in high fidelity turbulence simulations. In addition, the turbulent transport equations are solved in a strongly coupled way by using GKS with scalar transport. The extended GKS is applied in typical turbulent flow predictions including the RANS simulation of hypersonic compression ramp flow and the detailed simulation of multiscale turbulent structures in lowspeed cylinder flow with hybrid models. The predicted results agree well with existing experimental measurements and numerical studies, which shows the good accuracy, resolution and robustness of the extended GKS and reveals the wide prospects in turbulence simulations on different model scales, including the multiscale models such as hybrid RANS/LES methods. Furthermore, the multiscale turbulence simulation methods are worthy of further studies based on the generalized BGK model. © 2018 Elsevier Masson SAS. All rights reserved.

86 87 88 89 90 91 92 93 94 95 96 97 98 99

34

100

35

101

36

102

37 38

103

1. Introduction

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

The gas-kinetic scheme (GKS or GKS-NS) is a new CFD method based on the approximate solution of mesoscopic BGK equation on Navier–Stokes (NS) level [1]. In this scheme, the viscous and inviscid transports are coupled automatically and with consistent inherent dissipation mechanism, which guarantees its good performances in various flow problems, such as hypersonic viscous flows [2,3], magnetohydrodynamics [4], as well as the high accuracy and high resolution schemes for flow fine simulations [5–7]. The most existing applications focus on laminar flows so far. Turbulence is a typical flow problem in engineering applications. Due to the multiscale features, it is a challenge to balance the accuracy requirements and computational costs [8,9] in numerical simulations, especially for high-Reynolds-number flows. Among turbulence models or simulation methods, the direct numerical simulation (DNS) solves the NS equations on grid cell size near the smallest (dissipation) scale. DNS method can capture fluctuations on all the scales, but the computational cost is too high for engineering high-Reynolds-number flow simulations. On the contrary, with the help of turbulence models, the Reynolds averaged Navier–Stokes (RANS) equations can be cheaply solved with

60 61 62 63 64 65 66

*

Corresponding author. E-mail address: [email protected] (Q. Li).

https://doi.org/10.1016/j.ast.2018.04.022 1270-9638/© 2018 Elsevier Masson SAS. All rights reserved.

coarse grid only to capture the very large (integral) scale structures. Different from DNS and RANS, large eddy simulation (LES) [10] handles the filtered NS equations, and is closed by additional sub-grid models based on the filtering scale or grid cell size. Thus LES can be regarded as a multiscale method. To further reduce the grid cost of LES in near wall regions where turbulent fluctuations are strong anisotropic and typically small scale, the hybrid RANS/LES method is proposed and become the hot topic in turbulence high fidelity simulations [11]. The hybrid method is also the multiscale or multilevel method [12], which can keep good balance between flow resolution accuracy and computational cost. GKS can be directly applied to DNS of low-Reynolds-number turbulent flows, such as the DNS of mixing layer [13,14], the homogeneous turbulence [15,16] and the channel turbulent flow [17]. For high-Reynolds-number turbulence, the extended BGK equation with the effective relaxation time τe can be used to describe the relaxation of turbulent fluctuations [18]. τe can be determined by traditional turbulence models. The corresponding scheme is the extended GKS, and will be introduced in details in the following. In our previous studies, typical turbulence models such as Baldwin– Lomax (BL), k–ω SST and RNG k–ε models were implemented in GKS and showed good performances [19,20]. Similar conclusions are also got in other studies [21,22]. In recent research [23], k–ω model is considered in GKS via τe in renormalized form [18,24], and the turbulent ‘rarefaction’ effect is proposed and discussed. In

104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.2 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

turbulence fine simulations, LES simulation with GKS has also been conducted in complex flows around a car side mirror [25]. Despite many existing studies on the extended GKS, there are still some problems worthy of further investigations. The first one is the influences of numerical schemes and turbulence models. The self-adapting numerical dissipation [1] in the extended GKS is of more interest which comes from the cross-scale evolving solution of the extended BGK equation. The performance in hypersonic turbulent flows also requires evaluation. The second one is the capacity of the extended GKS in high fidelity simulations with multiscale models such as the hybrid RANS/LES methods. This work is expected to be a basis for further researches of turbulence simulating and modeling based on gas kinetic theory. The rest of this paper is arranged as follows. The construction of GKS is briefly introduced in section 2. The construction and analysis of the extended GKS for turbulence simulations are in section 3, while followed by the numerical simulations and discussions in section 4. And the final section is the conclusions.

19 20

2. Gas-kinetic scheme

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

GKS which describes macroscopic flows from the mesoscopic BGK-Boltzmann equation is established by Xu et al. [1,3]. The scheme is briefly introduced in follows. The BGK equation is,

38 39 40 41

and g is the equilibrium state. f and g are both functions in the high-dimensional phase space (x, t , u, ξ ) where ξ is the internal degree of freedom. And g is the Maxwell distribution,

 g=ρ

λ

( K +3)/2 e

π

44 45 46

−λ(|u−U|2 +ξ 2 )

(2)

.

Here λ = ρ /(2p ), ξ 2 = ξ12 + ξ22 + ... + ξ K2 . K is the total number of internal degrees of freedom, with value K = 2 for the diatomic gas molecule in three-dimensional (3D) flows. If f is known, the macroscopic conserved quantities and fluxes at the cell interface in the finite volume scheme can be got,



42 43

T

Q = (ρ , ρ U, ρ E ) =

 f d ,



Fm =

(3) um  f d ,

m = 1, 2, 3.

47 48 49 50 51 52 53 54

During the particle collision processes, f and g satisfy the conservation constraint,



( g − f )d = 0,

(4)



T where d = du 1 du 2 du 3 dξ1 dξ2 ...dξ K ,  = 1, u, 12 (|u − U|2 + ξ 2 ) . BGK equation Eq. (1) has a general solution,

55 56

59 60 61 62 63 64 65 66

68



=

1

t





g (x , t , u, ξ )e

τ

−(t −t  )/τ



dt + e

−t /τ

f 0 (x − ut , u, ξ ).

(5)

0

Here f 0 is the initial distribution function at t = 0, in which x = x − u(t − t  ) is the trajectory of particle motions. The characteristic scale of f 0 is the kinetic scale, while the characteristic scale of g is the hydrodynamic scale. Thus, the general solution Eq. (5) naturally contains the multiscale transition between these two extreme scales, and the transition process is controlled by



f 0 (x, u, ξ ) = 1 + alm xm − τ (alm um + A l ) (1 − H[x1 ]) g l

  r r + 1 + am xm − τ (am um + A r ) H[x1 ] g r ,

g (x, t , u, ξ ) =

(1 + (1 − H[x1 ])¯alm xm l,r l,r am , a¯ m ,

l,r

r + H[x1 ]¯am xm

¯ ) g0 . + At

71 72 73 74 75

78 79 80

¯ come from the derivatives of and A

81 82 83 84 85 86 87 88 89 90 91

+ (−τ + C 1 + C 2 )(¯alm um H[u 1 ]

92

r um (1 − H[u 1 ])) g 0 + a¯ m (7)   l l l + C 0 − (C 1 + C 2 )am um − C 1 A H[u 1 ] g   r + C 0 − (C 1 + C 2 )am um − C 1 A r (1 − H[u 1 ]) g r ,

93

where C 0 = e−t /τ , C 1 = τ C 0 and C 2 = tC 0 . With f , the macroscopic numerical fluxes can be calculated according to Eq. (3), and the macroscopic conserved quantities are updated through the finite volume method. Besides, the second-order accuracy in both time and space can be achieved with a single step. Meanwhile, the distribution function f in Eq. (7) is a combination of Maxwellian function through which the integration to compute macroscopic quantities is very simple. So the computational cost is competitive with traditional schemes directly based on macroscopic governing equations. The usual reconstruction techniques can be directly applied to obtain the macroscopic quantities and their slopes at the cell interface. More details of GKS can be found in the literature [1,3,5,26] and will not be repeated here. It should be noted that in the above distribution function, the particle movements in both the normal and tangential directions of the cell interface are taken into account, thus the resulted GKS is truly multidimensional which is difficult for traditional CFD schemes based on macroscopic Euler equations. If neglect the tangential coefficients, the directional splitting scheme can be recovered. Besides, the viscous effect is controlled by the collision time τ which is computed by

p

70

77

(6)

¯ g0 f (0, t , u, ξ ) = (1 − C 0 ) g 0 + (t − τ + C 1 ) A

μ

69

76

The coefficients A Maxwellian distribution, and can be computed by the gradients of macroscopic conserved quantities. The superscripts l, r represent the quantities on the left or the right side of the cell interface, respectively. H[x] is the Heviside function to account for the discontinuity at a cell interface. So the time evolving distribution function f at a cell interface (x = 0) can be explicitly expressed in terms of local Maxwellian distributions,

τ=

f (x, t , u, ξ )

57 58

67

by the first-order Chapman–Enskog expansion instead of being directly discretized in the high dimensional phase space, while the Taylor expansion is adopted to achieve the target accuracy order. Therefore, huge amounts of computational costs can be avoided. In the second-order accurate GKS, the first order Taylor expansion is adopted, and then f 0 and g in Eq. (5) are constructed as, (taking the normal direction x1 of the cell interface as an example),

∂f ∂f g− f +u· = , (1) ∂t ∂x τ where τ is the relaxation time. f is the gas distribution function

36 37

t /τ . When simulating flows on NS level, f 0 can be constructed

+ C t

|p − p | . | pl + p r | l

r

94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

(8)

Here μ is the molecular viscosity, while pl,r is the reconstructed pressure at the corresponding side of a cell interface. The constant C = 1 is simply set and the time step t is constrained by the Courant–Friedrichs–Lewy condition. The last term in the above equation takes into account the effect of the discontinuity at a cell interface. The free transport of particle in f 0 leads to the upwind characteristics of the corresponding scheme, such as the kinetic flux vector splitting method (KFVS) [27]. The equilibrium state g corresponds to a central scheme. It should be noted that if f 0 and

122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.3 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

3

 1 2 3 4 5 6 7 8 9 10 11

g are constructed on the same scale, such as the NS scale or a RANS one, the transition from f 0 to g in Eqs. (5) and (7) is more like a switching from upwind scheme to a central one. This means a dynamic adjustment of numerical dissipation which is controlled by the ratio t /τ . That is, the numerical dissipation in GKS flux is included not only in an explicit way through the modification of τ , but also an implicit method through a dynamic adjustment of the weight of upwind scheme. This mechanism leads to a good balance between strong robustness and high accuracy, which is especially important for numerical schemes in high-speed viscous flow simulations [1–3].

12 13

3. Extended GKS for turbulence simulation

14 15

3.1. The extended GKS

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

To simulate turbulence effectively the RANS or the filtered NS equations can be adopted along with the corresponding turbulence models. These governing equations are similar to the NS equations if the eddy viscosity model is implemented, except that the molecular viscosity is replaced by the effective turbulent viscosity. Thus GKS can be extended for turbulence simulation which is based on the generalized relaxation model to describe the relaxation of turbulent fluctuations [18],

∂f ∂f g− f +u· = . ∂t ∂x τe

Here the effective relaxation time τe uniformly describes the effects of particle thermal motions and turbulent fluctuations,

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

(9)

τe =

μ + μt p + 2k/3

+ C t

|p − p | , | pl + p r | l

r

(10)

where k is the turbulent kinetic energy and μt is the turbulent viscosity, while the second part on the right side is for the artificial viscosity. Through the first-order Chapman–Enskog expansion, which is similar to the derivation of NS equations, and after the integration of the distribution function in the phase space, the corresponding macroscopic RANS or filtered NS equations can be recovered. In the extended GKS to approach the solution of Eq. (9), the construction of the distribution function f is similar to Eqs. (5) and (6), while the relaxation time is modified as τe . The corresponding coefficients are computed by the quantities and their gradients of the simulated turbulent flow fields. To determined the turbulent eddy viscosity μt , the turbulence model should be adopted. For LES models such as the Smagorinsky model, the turbulent viscosity can be directly calculated by the shear rate S and grid scale  as μt = C S ρ S 2 , where C s is the model constant. For RANS models, the transport equations of turbulence quantities, such as those for k and ω in the k–ω model should be solved. The transport equations of present considered turbulence models such as S-A and k–ω are given in the Subsection 3.2. When solving these equations, one common-used method is the decoupled treatment [21,28] in which other numerical schemes, such as FVS or Roe scheme, can be directly used to solve the turbulence transport equations. In the present study, a coupled method is applied [19,20]. With a operator splitting to handle the source terms separately, the turbulence transport equations can be solved via GKS with scalar transport [29]. In the scheme, the turbulence quantities are handled as the new internal degrees of freedom. Taking the one-equation S-A model as an example in the extended GKS, the macroscopic variables are calculated with f as,

Q = (ρ , ρ U, ρ E , ρ ν˜ ) T =



 f d

(11)

68

T

69

Here  = 1, u, 12 (|u − U|2 + ξ 2 ), θν , where θν is the corresponding microscopic degree of freedom of the turbulence quantity ρ ν˜ , and d = du 1 du 2 du 3 dξ1 dξ2 ...dξ K dθν . It should be noted that the Favre-averaged or spatially filtered macroscopic quantities hereafter are represented by the same symbols for the original ones for simplicity. The Maxwell distribution g in the extended GKS is,

 g=ρ

λ

π

( K +3)/2

67

70 71 72 73 74 75 76 77

2 2 2 e−λ(|u−U| +ξ +(θν −ν˜ ) ) .

(12)

The numerical fluxes at the cell interface are calculated with f which has the same form as Fm in Eq. (3). Furthermore, in the scheme with scalar transport [29], the computations of the scalar fluxes are split from the fluxes of the conservative variables, so that the scheme can simulated flows with arbitrary Schmidt number and thus can be directly used in the computations of the turbulence transport equations. In such treatments, the numerical accuracy and robustness of turbulence quantities can keep consistent with the conservative variables. More details of the scheme can be found in the literature [29] and will not be repeated here. In the construction of the extended GKS, besides the basic conceptions of the effective relaxation time and the solution of the turbulence model equations, two important points require further discussions. Firstly, when the Chapman–Enskog expansion is used to solve the extended BGK equation, the scale of f 0 and g are both determined by the specific scale of utilized turbulence models, which means that the multiscale transition in the extended BGK equation from f 0 to g still performs like the automatic adjustment of numerical dissipation from upwind to central type in the current frame. Thus, the extended scheme cannot achieve multiscale simulations rooted from the generalized relaxation model. The analysis is validated to both RANS and LES models, as well as the hybrid RANS/LES methods. On the other hand, due to that the scale of turbulence model itself changes with the grid scale (filtering scale) in LES and hybrid RANS/LES models, the extended GKS coupled with these models can still performance as a multiscale method. Secondly, Chen et al. [18] has concluded some advantages in turbulence simulations with the extended BGK equation. And it is considered that the projection of high-dimensional space to lowdimensional will enable the scheme to represent rich and complex physics. However, in the present frame, the extended BGK equation is not discretized in the high-dimensional phase space, but approximated with the first-order Chapman–Enskog expansion. These two points are important in the analysis and evaluations of the extended GKS. In fact, the physical contents of the extended BGK equation are far more than that in current extended GKS. It is worthwhile to realize the multiscale descriptions of turbulence from the cross-scale solution of the extended BGK equation in the future studies. It should be mentioned that in turbulence fine simulations with LES or hybrid RANS/LES models, the high accuracy and high resolution of the numerical scheme are very important to capture multiscale turbulent structures with limited grid resolution. In this study, GKS with the minimized dispersion and controllable dissipation (MDCD) reconstruction [30,31] is used in the LES and hybrid RANS/LES simulations. MDCD-GKS has good performances in complex flow field simulations which is discussed in detail in our previous work [7]. In subsonic flow simulations, the linear MDCD reconstruction can be directly used without the hybrid with WENO, thus the maximum dissipation is strictly limited by the dissipation parameter γdiss of MDCD.

78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.4 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

4

1

3.2. Turbulence models

2 3 4 5 6 7

ρk

μt =

The S-A, k–ω and SST k–ω RANS models are used in the simulation of high-Reynolds-number flows. The Vreman–Smagorinsky LES model and hybrid DES-SA/IDDES-SA are applied in turbulence fine simulations. The models are briefly introduced in the following.

max(ω, C lim

10 11 12 13 14 15 16 17

3.2.1. S-A model The S-A model is a typical one-equation model with the transport equation [32],

∂ ρ νˆ U j ∂ ρ νˆ + = C b1 (1 − f t2 ) ρ νˆ ∂t ∂xj   C b1 ((1 − f t2 ) f v2 + f t2 ) +ρ − C f w1 w 2  2 νˆ

18

∂ ν ρ νˆ 2 − × d σ ∂xj    ∂ νˆ  1 ∂ . + μ + (1 + C b2 ) ρ νˆ σ ∂xj ∂xj

19 20 21 22 23 24

2ˆ

C b2

(13)

27 28

μt = ρ νˆ f v1 .

(14)

The related parameters and constants in Eqs. (13)–(14) are,

29 30 31 32

χ˜ 3

νˆ 2 f ν1 = , χ˜ = , f t2 = C t3 e−C t4 χ˜ , ν χ˜ 3 + C 3v1

33 34 35 36 37 38 39 40 41 42 43 44 45 46

fw = g

1 + C 6w3 g6

+

1

6

, g = r + C w2 (r 6 − r ),

C 6w3

νˆ νˆ f v2 χ˜ r= , Sˆ =  + 2 2 , f v2 = 1 − ; ˆS κ 2 d2 1 + χ˜ f v1 κ d C b1 = 0.1355,

(15)

2

σ = , C b2 = 0.622, κ = 0.41, C w3 = 2.0, C b1

κ2

+

1 + C b2

49 50 51

σ

54 55 56 57 58 59 60 61 62 63 64 65 66

72

2

S¯ i j =

α=

1 2

13 25



3

∂U j ∂Ui + ∂xj ∂ xi

, β∗ =

9



73 74

−

1 ∂ Uk 3 ∂ xk

7

75

8

76

δi j , C lim = ,

1

3

2

5

77

, σ = , σ∗ = ,

100 ∂k ∂ ω

78 79

≤0 ∂xj ∂xj 1 σd = , σd0 = , ∂ k ∂ ω ⎪ 8 ⎪ ⎩ σd0 , >0 ∂xj ∂xj 1 + 85χω β = β0 f b eta, β0 = 0.0708, f β = , 1 + 100χω        Sˆ  ∂Ui 1 ∂ Uk 1 ∂ Um  i j jk ki  ˆ − χ ≡ S = + δki . ,  ki ∗ 3  (β ω)  2 ∂ xi ∂ xk 2 ∂ xm 0,

80 81

(18)

82 83 84 85 86 87 88 89

To improve the computational efficiency the wall function [33] is adopted in the above S-A model. 3.2.2. k–ω model The k–ω (Wilcox06) model is a typical two-equation model in which the transport equations for the turbulence quantities k and ω are [9],

∂ ρ k ∂ ρ kU j + ∂t ∂xj    ∂ ∂Ui ρk ∂ k , = ρ Tij − β ∗ ρ kω + μ + σ∗ ∂xj ∂xj ω ∂xj ∂ ρω ∂ ρω U j + ∂t ∂xj ω ∂Ui ρ ∂k ∂ ω = α ρ Tij − β ρω2 + σd k ∂xj ω ∂xj ∂xj    ∂ ρk ∂ ω . + μ+σ ∂xj ω ∂xj The turbulence viscosity computed by,

The SST k–ω model is a combination of k–ω model and k– model. The details of the model can be found in the literature [34], and will not be repeated here. In hypersonic flows, the contribution of turbulent kinetic energy k usually cannot be ignored. k should be explicitly taken into account in the Reynolds stress term ρτi j of the momentum equation as shown in Eq. (18). Besides, k also need to be explicitly included in the total energy ρ E = ρ (e + 12 U i U i + k), while the flux of ρ E need to be modified with the viscous flux of k correspondingly [9,35],

∂ F˜ ρ E = F ρ E + ∂xj

 

 ∂k μ + σ ∗ μt ∂xj



μt = min(μt , ρκ d a1k).

52 53

71

ρ T i j = 2μt S¯ i j − ρ kδi j ,

∗

.

47 48

68 69



(16)

91 92 93 94 95 96 97 98 99 100 101 102

.

(19)

103 104

The compressibility correction of k–ω model is also considered, in which the length scale correction is employed to reduce the magnitude of heat transfer near the reattachment point [36,35],

3 C v1 = 7.1, C t3 = 1.2, C t4 = 0.5, C w2 = 0.3, C w1 =

)

67

90

The turbulence viscosity is,

25 26

β∗

(17)

.

70

⎧ ⎪ ⎪ ⎨

κ

2 S¯ i j S¯ i j

The related parameters and constants in Eqs. (16)–(17) are,

8 9



105 106 107 108

(20)

Here d is the wall distance, and a1 = 0.31, κ = 0.41. The corresponding model is labeled as k–ω –LS. No other compressibility corrections are applied in the turbulence models. According to Eq. (13) and Eq. (16), it can be seen that except for the convection and diffusion terms, the sources terms such as the production and dissipation terms should be computed when solving the turbulence transport equations. The source terms are usually related with the gradients of the flow field quantities. In super/hypersonic flow fields, the gradients near a shock are usually nonphysical, not only because of the capability of NS equations to resolve the real shock structure, but also because the grid resolution is always limited. It often leads to artificial sudden increase of the turbulent kinetic energy k. In computations of the convection and diffusion terms, such numerical effects will be restrained via the slope limiter in the reconstruction. Therefore, the gradients terms in the source terms should also be limited reasonably. In recent study, a damping function based smooth measurement factor of WENO scheme is used to eliminate the fake amplification of k [37]. In the present study, a simple limiter is developed based on the pressure shock detector [38] (taking the i-direction as an example),

109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.5 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10

s=

| p i −1 − 2p i + p i +1 | . | p i −1 + 2p i + p i +1 |

(21)

The maximal value smax is chosen among the seven neighboring computational cells (three cells along each direction). Then the gradient term of macroscopic quantities ∇ Q is limited as,

 ∇ Q = (1 − smax )∇ Q .

(22)

With the modification, the stability of the turbulence equations can be greatly improved.

11 12 13 14 15 16 17 18 19 20 21 22



μt = ρ c

25 26 27 28 29 30 31 32

Bβ ai j ai j

(23)

.

Here c = is computed by the model constant C s = 0.1 in the original Smagorinsky model. The related parameters are, 2.5C s2

∂U j ai j = , βi j = 2aim amj , ∂ xi

(24)

2 2 2 B β = β11 β22 − β12 + β11 β33 − β13 + β22 β33 − β23 .

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

3.2.4. Hybrid RANS/LES model The RANS and filtered NS equations can be uniformly described with the integration of Eq. (9). Thus, the hybrid model can be directly implemented in GKS via the hybrid turbulent viscosity in Eq. (10). Among the hybrid RANS/LES methods, DES (detached eddy simulation) class models [40] only need minor changes from the corresponding RANS model, and have been widely used. In the present study, the DES/IDDES (improved delayed DES)-SA models are both considered in the extended GKS to conduct efficient high fidelity turbulence simulations. The influence of turbulence model is discussed in the extended GKS. The DES-SA model [41,42] is developed with the modification of length scale in the destruction term of S-A model,

d˜ = min(d, ϕ S A C D E S ).

56

(25)

Here d is the distance to the closest wall and the grid scale  = max is chosen as the maximum grid spacing in all three directions. C D E S = 0.65 is the model constant. ϕ S A is the correction function which prevent μt from declining to 0 in LES regions,

⎛

54 55

ϕ

2 SA

= min ⎝10 , 2

1−

C b1

∗ C w1 κ 2 f w

( f t2 + (1 − f t2 ) f v2 )

f v1 max(10−20 , 1 − f t2 )

⎞

⎠.

(26)

57 58 59 60 61 62 63 64 65 66

In LES and hybrid RANS/LES methods, the turbulence model changes with the computational grid resolution (grid scale). These multiscale methods can achieve good balance between resolution accuracy and computational cost. However, it should be emphasized again that, the multiscale transitions in the simulations come from the corresponding basic principle of LES, or partition modeling concept of hybrid models. In current frame of the extended GKS, the inherent multiscale mechanism of the gas-kinetic method has no contributions on the multiscale simulations of turbulence.

The related parameters are same as S-A model in Eq. (15). DES model has inherent defects such as the modeled stress depletion (MSD), the logarithmic-layer mismatch (LLM) and the grid-induced separation (GIS) problems. In order to settle these problems, IDDES [43] is developed, and is recognized as one of the best hybrid RANS/LES model. IDDES combines the so-called wall-modeled LES and the DDES (delayed DES) method. The performance of the model is determined by the upstream flow features. The details of the IDDES model can be found in the literature [43].

67 68 69 70 71 72 73 74 75 76 77

4. Numerical tests and discussions 3.2.3. LES model Other than the RANS simulations, the unsteady fine simulations are also conducted by using LES method. A typical eddy viscosity LES model is the Smagorinsky model. However, its dissipation is too strong in the near wall regions, and need to be modified with the wall function. A modified Smagorinsky model developed by Vreman [39] is implemented in the present study. In the nearwall regions, the dissipation is naturally small, which is helpful to capture turbulent boundary. The specific expression of the model is,

23 24

5

78

Several typical turbulent flows are simulated with the present extended GKS. Of great interest is the performance of the present scheme in the RANS simulation of hypersonic flow, in which the accuracy and robustness are both very important. In addition, the grid convergence is emphasized for hypersonic turbulence simulations which is also fundamental but limited in existing studies. Special attention is also paid to the accuracy and resolution on multiscale structures in the fine simulation of turbulent flows with LES and hybrid RANS/LES models. Without special remarks, the van Leer limiter is used in the reconstruction based on conservative variables. For steady flows, the time-implicit LU-SGS scheme and local time steps are applied to accelerate the convergence. The converged result is achieved after the 2-norm of the total energy residual drops by four orders of magnitude in the transonic M6 wing flow, and three orders of magnitude in the hypersonic ramp flow when the computation starts from a uniform initial flow field. The viscosity follows the Sutherland’s law in RANS simulations and the exponential law in LES and DES/IDDES simulations. The turbulent Prandtl number is set as 8/9 for k–ω model and 0.9 elsewhere. In unsteady simulations, the time averaged quantities are symbolized by the angle brackets ‘< >’, and the corresponding fluctuating quantities are denoted with ‘  ’.

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

4.1. M6 wing flow

103

The transonic flow around the ONERA M6 wing is a standard 3D benchmark for engineering simulations. Besides the 3D geometry, the flow structures are complex which include the interaction of shock and turbulent boundary and thus it is a good candidate to test the performance of the present extended GKS. The experimental and numerical data can be found in the NPARC Alliance Validation Archive project [44]. The computational conditions are, ◦

104 105 106 107 108 109 110 111 112

(27)

113

where the Reynolds number Rec is based on the mean aerodynamic chord (C ) of the wing. The standard grid marked as G2 is also provided in the NPARC project, which contains 4 structured blocks with 48 × 48 × 32 cells for each block. With the same turbulence model and grid, the advantages of the extended GKS can be measured clearly. An additional grid (G1) with 32 × 33 × 21 cells for each block is also tried to grid convergence test. In the flux computations, the directional splitting scheme is applied. The subsonic inflow, outflow and far field boundaries are all set according to the locally 1D Riemann invariants. The solid wall is adiabatic. Considering the excellent performance in aerodynamics simulations, S-A model is used in this case. Besides, according to the first grid size normal to the wing surface the wall function [33] is adopted which has simple algebraic forms, and can cover large y + ranges. The global CFL number is 50, while the local time step is also adopted to speedup the computation. The simulation with grid G2 takes less than an hour using 4 cores of Intel Xeon X5670 CPU.

115

6

Ma∞ = 0.84, AOA = 3.06 , Rec = 11.72 × 10 ,

114 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

6

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.6 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Fig. 1. Normalized pressure on M6 wing surface and the symmetry plane. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

21 22 23 24 25 26 27

Fig. 1 shows the normalized pressure on the wing surface and at the symmetry plane predicted with the G2 grid. The pressure shows an evident increase at the leading edge, which is followed by a decrease on the leeward side along with velocity acceleration until a shock wave formation. Due to the 3D geometry of the wing, the λ shock can be clearly identified at the leeward side. The pres-

sure coefficient at several spanwise sections are detailedly shown in Fig. 2, and compared with the experimental data and numerical results also computed with S-A model from the NPARC project. In the figures, b is the wing span length and C z is the chord at the local section. The present results of extended GKS-SA agree well with the experimental data, which show a bit better than the referenced numerical ones. Especially, at the section of z/b = 0.8, the λ shock (double-step) structure is well captured by the present GKSSA, while it is not resolved in the referenced numerical study with the same computational grid and turbulence model. The coarser grid G1 gives similar results except that the predicted shock wave seems a little more dissipative when compared with the standard grid. The wall friction coefficients at three spanwise locations are shown in Fig. 3. The friction profiles are similar on the windward side at different spanwise locations. However, there are large differences among these profiles on the leeward side, which is induced by the λ structure. And two additional dramatic decreases in the friction coefficient are triggered by the shock as well. In fact, the prediction of grid G1 also agrees well with G2 on the windward side, while distinct difference can be observed on the leeward side, mainly due to the different resolution of the shock structure. When compared with the referenced data computed with the same grid but a different numerical scheme in the WIND code, the present results show distinct improvements using the extended

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65 66

131

Fig. 2. Pressure coefficients at different spanwise locations of the M6 wing. The referenced experimental and S-A data are from the NPARC project [44].

132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.7 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

Fig. 3. Wall friction coefficient at different spanwise locations of M6 wing.

20 21 22 23 24 25

86 87

GKS. It should be pointed out that details of the used wall function for S-A model are not given in the referenced simulations which will also contributes to the differences between the simulated results. In brief, the feasibility of the extended GKS in high efficiency turbulence simulations is verified in this standard case.

88 89 90 91

26 27

92 93

4.2. Hypersonic compression ramp flow

28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

94

Due to the strong interaction among the shock wave, boundary layer and compressible turbulence, the hypersonic flow around a 2D compression ramp is very challenging for CFD methods, especially to predict the wall heat flux as well as the flow separation accurately, which requires both the strong robustness and high accuracy. Benefited from the automatic adjustments dissipation, GKS shows good performances in hypersonic laminar flows [2,3]. It is of great importance to validate the present extended GKS in this typical hypersonic turbulent flow. Besides, the grid convergence is paid more attention to which is also fundamental but limited in existing studies. The flow conditions are [45],

Ma∞ = 9.22, Re∞ = 4.7 × 107 m−1 , T ∞ = 64.5 K, T wall = 295.0 K.

(28)

The working gas is nitrogen and the ramp angle is 38◦ . As no stagnation or free-stream pressure is reported in the experiment, p ∞ is set as 2,300 Pa [46]. The experiment is carried out in the hypersonic gun tunnel and the length of wedge flat plate is L = L 1 = 0.43 m for pressure measurement [47], while L = L 2 = 0.56 m for heat transfer [48]. As the difference in experimental condition is distinct, the computations of pressure and heat flux are conducted with different geometry conditions respectively. The computational domain is shown in Fig. 4, with L 2 condition for an example. The supersonic inflow condition is used on the left and top boundaries with turbulence intensity Tu = 0.1% and turbulent viscosity μt /μ = 10. The linear extrapolation condition is used on the right boundary. On the bottom boundary, the symmetry condition is taken when x < − L, and the isothermal wall is taken when x ≥ − L. The grid parameters are shown in Table 1, where three grids are tried for the grid independent study. The primary variables are used in the reconstruction for better positivity of the reconstructed temperature or pressure. The computed pressure field is displayed in Fig. 4, in which the complex shock structures, including the leading edge shock wave, the separation shock wave and the reattachment shock wave, along with their interactions with the turbulent boundary layer bring

95 96 97 98

Fig. 4. The pressure contours in the hypersonic compression ramp flow computed by SST model with L = 0.56 m.

99 100 101

Table 1 The grid parameters in the hypersonic compression ramp flow. xmin is the grid cell size at both the leading edge and the corner. Grid

Nx L1

G1 G2 G3

– 312 –

Ny

 y min

xmin

102 103 104 105 106

L2

flat plate

ramp

168 324 648

3.17 × 10−6

2.64 × 10−7

2.0 × 10−4

107

7.61 × 10−7

6.24 × 10−8

5.0 × 10−5

109

75 150 300

1.54 × 10−6

1.27 × 10−7

1.0 × 10−4

108 110 111

great challenges for numerical simulations. The separation bubble can be clearly observed near the corner. The strong interaction between the separation shock and the wall after the bubble results in a strong reattachment shock and a dramatic increase of pressure and heat flux. Fig. 5 shows the computed wall heat flux with k–ω , k–ω –LS and SST model on different grids. It is interesting that with the coarse grid G1, SST model predicts heat flux which agrees well with the experimental data, but the results show obvious deviations when using the fine grid. The differences between G2 and G3 grid are small for k–ω –LS and SST. Thus, the results on the G2 grid of these two models can be considered to be grid converged in which the maximum y + is about 0.5, while the results of k–ω on the G2 grid are also included to qualitatively analyze the effect of the length scale correction in k–ω –LS. The computation on G2 with L 2 takes about 11 hours using the Intel Xeon X5670 CPU with 12 cores. Fig. 6 shows the wall pressure and heat flux predicted with different turbulence models, as well as different plate length L 1 and L 2 . Firstly it can be obviously observed that the length of flat plate L indeed has obvious influences on the flow fields. In the for-

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

8

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.8 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19 20

Fig. 5. The normalized wall heat flux in the hypersonic compression ramp flow computed with different grids and compared with the experimental data [45], q∞ = 6.56 W/cm2 .

85 86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40 41

Fig. 6. The normalized wall pressure and heat flux in the hypersonic compression ramp flow predicted with k–ω , k–ω –LS and SST model and compared with the experimental data [45]. The left figure is for L = 0.43 m and the right figure is for L = 0.56 m. The referenced results of k–ω –LS and SST are computed with L = 0.56 m [46].

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

106 107 108

mer studies, the same condition L = 0.56 m is used both in pressure and heat flux predictions [46,49]. Thus, the conclusions about the numerical results and the assessments about the compressibility corrections of turbulence model may be debatable. In fact, the study based on coarse-grid results may also be questionable. The referenced numerical results [46] are computed with the same turbulence models and plate length but with a much coarser grid than the current used G2 grid. Thus, these two factors may both contribute to the differences between the extended GKS and the referenced numerical data [46]. Secondly, differences of the considered turbulence models are evident in the figure. The peak positions of p w and q w predicted by the original k–ω model get the best agreement with the experimental data, while the peak value of q w is much higher. The modification of the length scale correction can reduce the heat transfer near the reattachment point, but it also enlarges the separation bubble. Meanwhile, the compressibility corrections usually consider the additional dissipation effects of compressibility on turbulent kinetic energy [49–51], which will further decrease the turbulent kinetic energy and enlarge the separation. Therefore, the compressibility corrections such as the dilatation–dissipation and pressure–dilatation terms for k–ω model are not further used in the extended GKS. Fig. 7 shows the temperature and kinetic energy near the corner, computed by k–ω and SST model with L = 0.56 m. The rapid

strain feature in SST model enhances the turbulent dissipation in regions of rapid compression, so that the turbulent kinetic energy and eddy viscosity decrease in the separation region, which leads to much larger separation bubble. Based on this case, two conclusions can be got. Firstly, the extended GKS shows good robustness and accuracy in the hypersonic flows. Secondly, the influences of geometry conditions, grid convergence and turbulence models are comprehensively considered, and the reliable simulation results are obtained for this challenging hypersonic case which agree well with the experimental data both in the pressure and heat flux. The simulation results are meaningful in the verification, assessment and improvement of turbulence models. It should be mentioned that the constant turbulent Prandtl number is adopted which is a strong simplification in hypersonic flows where the heat transfer is important. For further analysis of the differences between experimental data and the numerical simulations, a variable turbulent Prandtl number can be considered [52,53]. Besides, according to the studies of Babinsky [54], the instability and small structures in the flow field may also be the reason to cause the differences between the RANS simulations and the experiments, especially in the heat flux. Thus, the developments and improvements of turbulence model still require further studies, in which the robust GKS for high speed simulation is with more important application prospects.

109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.9 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14 15

80

Fig. 7. The normalized temperature (T / T ∞ , left) and kinetic energy (k/(γ R T ∞ ), right) near the corner computed by k–ω (solid lines) and SST (dashed lines) models.

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34 35

100

Fig. 8. The Reynolds stress profiles in channel flow, which are normalized by the friction velocity U τ2 . DNS data are from [55].

36 37

102

4.3. Turbulent channel flow

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

101

For better prediction of the turbulent fluctuation, it is necessary to do fine simulations. The channel flow between two infinite parallel flat plates has simple geometry but abundant multiscale turbulent structures, and thus can be adopted to validate the extended GKS in turbulence high fidelity simulation. The accuracy and resolution is of great interest which play a dominant role in the capturing of multiscale fluctuations. Besides, there exist a lot of studies on this benchmark wall-bounded turbulence. In multiscale turbulence simulations with LES and hybrid RANS/LES models, the MDCD reconstruction is adopted to enhance the resolution of GKS, which has been validated in our previous study [7]. The dissipation parameter in MDCD is set as 3 × 10−4 in this case. Two turbulence models are considered, including Vreman–Smagorinsky [39] LES model (c = 0.025) and SA-IDDES [43] model. SA-DES [41] model is not taken into account due to its MSD and LLM problems in channel flows. The Reynolds number is 3300 based on the volume-averaged velocity and channel half-width. It is about 203 for IDDES and 205 for LES, respectively, when based the computed friction velocity (Reτ ). The mean Mach number is 0.3. The computational domain is 2π × 2 × π , which is divided into 33 × 64 × 32 cells. The dimensionless grid size is about x+ ≈ 34,  z+ ≈ 17 and  y + ≈ 0.5, min which stratifies the grid resolution demands in LES simulation. The solid walls are isothermal and periodic conditions are adopted in the streamwise and spanwise directions. A external force is imposed on the flow to hold a constant mass flow rate. The CFL number is set as 0.2.

More than 150 time units based on the volume-averaged velocity and the channel length are taken up to compute the statistical quantities. Fig. 8 shows the predicted Reynolds stress. The DNS data for incompressible flow with Reτ = 178, using 128 × 129 × 128 computational cells [55] are provided for comparison. The IDDES prediction agrees well with the DNS simulation which shows the good performance of present numerical methods for turbulence fine simulation. In LES simulations, the modeled shear stress of Vreman–Smagorinsky model seems insufficient especially in the buffer layer, which results in stronger turbulent fluctuations when compared with DNS and IDDES results. The predicted transverse intensity < V  V  > and the total shear stress are also little larger which leads to the deviation of normalized mean velocity from U + = 2.5 ln y + + 5.0 in log-law region as shown in Fig. 9. With the help of RANS model in near wall region, GKS-IDDES shows effective improvement. In this typical wall-bounded turbulence, the good performance of extended GKS is validated in the LES and hybrid RANS/LES simulation. The simulation results validate the τe model in BGK-type equation for the multi-scale simulation of turbulence. With the verification of turbulent boundary, the extended GKS with multiscale turbulence models can be further used in more complex and practical turbulence fine simulations.

103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

4.4. Cylinder flow

127 128

The flow past a long, circular cylinder flow is a typical benchmark with massive separations, which is difficult for not only RANS simulations but also hybrid RANS/LES methods to predict the flow separation and multiscale turbulent fluctuations.

129 130 131 132

JID:AESCTE

10

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.10 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19 20 21

85

Fig. 9. Mean streamwise velocity profiles of turbulent channel flow compared with DNS data [55].

86 87

22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

88

In the present study, the flow with subcritical Reynolds number Re D = 3900 (where D is the cylinder diameter) is considered. There are a lot of studies on this case, both in experiments [56–58] and numerical simulations which cover DNS, LES and hybrid RANS/LES simulations [59,60,58,61]. In these studies, two quantities are mainly focused on, which are the recirculation length (L r / D) and the mean streamwise velocity profile at the station x/ D = 1.06 whether it is the ‘U’ shape or the ‘V’ shape. In the discussions, the consensus cause is the shear layer transition. An earlier shear layer transition will lead to smaller L r / D and ‘V’ shape velocity distribution. Kravchenko and Moin [60] pointed out that the inflow disturbance in the experiment of Lourenco and Shih [57] was the reason to cause the shear layer earlier transition. Parnaudeau et al. [58] further supplemented that the lack of statistical convergence in the PIV measurement may be another factor, and also indicated that L r / D is sensitive to conditions such as the aspect ratio (L / D), where L is the cylinder length in the experiment. In numerical studies, the scheme dissipation and dispersion, the grid resolution, the sub-grid stress model used in LES and the computational domain in the spanwise direction all have influences on the shear layer transition. Kravchenko and Moin also concluded that, ‘the agreement of the calculations with experiments that inhibits early shear layer transition due to external disturbances is fortuitous’. Among the numerical results, in the DNS simulation of Ma et al. [59], the grid resolution and the spanwise domain are all carefully discussed, while the computational conditions are also easily to be followed. Thus, this DNS data is considered to be more suitable to assess the present simulation with GKS-DES/IDDES, rather than the experimental and other numerical data. In the present simulation, the Mach number of the free stream is set as 0.2 to approach the incompressible experimental conditions and the turbulent viscosity is chosen as ten percent of the molecular one. The computational grid is O-type with 15D in the radial direction of the cross-section and π D in the spanwise direction. The grid size is 160 × 160 × 32, and is uniform in spanwise direction and non-uniform in the cross-section. The first cell size off the solid wall is 7 × 10−4 D, while the mean y + for the first grid is about 0.12. The grid radial stretch ratio is 1.068 in the near-wall region within a distance 3D from the wall, and 1.0 else [62]. The grid in the shear layer and separation region is relatively coarse. The flow quantities at far field boundary are set according to the locally 1D Riemann invariants, and the wall is isothermal with T wall = T ∞ . MDCD-GKS is used and the dissipation parameter is

89 90 91 92 93

Fig. 10. The instantaneous flow structures identified by Q-criterion. The iso-surface is colored by streamwise velocity.

94 95 96

set as 0.005. The CFL number is set as 0.37. Fig. 10 shows the instantaneous flow structures of the cylinder flow predicted by GKS-DES and GKS-IDDES, respectively. The large vortex shedding can be clearly observed in which small vortices twist each other and form complicated flow structures. The mean flow field quantities are obtained by averaging within 200D /U ∞ time units which are about 40 shedding cycles. The parameters such as the base pressure coefficient (< C pb >), Strouhal number (St), the separation angle (θsep ), recirculation length L r / D and the minimum < U >min on the centerline are shown in Table 2. The predicted separation bubble size by GKS-DES is significantly smaller, which is due to the GIS problem. The GIS problem leads to the earlier transition of the hybrid model to LES in the shear layer. Along with the relatively inadequate grid resolution for LES simulations in this region, the earlier shear layer instability occurs. At the same time, the good agreements of GKS-IDDES with the experimental data well illustrate the good performance of IDDES in solving the GIS problem. The IDDES results in the literature [61] are also presented for comparison, with different numerical scheme and grid topology, where the grid size is about 5 times of the present grid. With the same model, the accuracy of GKSIDDES is better for the sensitive quantity L r / D when compared with both the experimental and DNS data. The precision of St got by GKS-IDDES and the referenced IDDES is similar when compared with the DNS data, while the < C pb > in the referenced study is better. Overall, the < C pb > and L r / D computed by the IDDES model still show distinct differences from the DNS/experimental data. This is probably from the different flow features predicted in the upstream attached boundary layer which is laminar in the DNS/experiments, and turbulent in the hybrid model simulation. Fig. 11 shows the mean pressure coefficient on the cylinder surface and the mean streamwise velocity at the centerline in the wake region, which are compared with experimental [56,58,57] and DNS [59] data. Fig. 12 shows the mean velocity profiles in the wake region. In these figures, the difference in the separation length between GKS-DES and GKS-IDDES can be obviously

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.11 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 5

< C pb >

6

St

7

θsep Lr / D < U >min /U ∞

9

67

Table 2 The mean flow parameters for Re D = 3900 cylinder flow predicted by different numerical methods and compared with experimental and DNS data.

4

8

11

68 69

Exp. [56,58]

DNS [59]

IDDES [61]

GKS-DES

GKS-IDDES

70

−0.90 ± 0.05 0.208 ± 0.002 85.0 ± 2◦ 1.51 −0.34∼−0.252

−0.84

−0.8780

−1.0271

−0.8904

71

0.219 – 1.59 –

0.222 87.00◦ 1.4270 –

0.207 88.39◦ 1.1142 −0.2945

0.214 86.92◦ 1.5014 −0.3259

72 73 74 75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29 30 31

Fig. 11. The mean pressure coefficient on the cylinder surface (left) and streamwise velocity at the centerline of wake region (right). The stagnation point is located at θ and the cylinder center is at x/ D = 0. The experimental data of < C p > are from [56] and < U > /U ∞ from [57,58], while the DES/IDDES data are from [61].

= 0◦ ,

95 96 97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50 51 52

116

Fig. 12. The mean streamwise (left) and transverse (right) velocities at different locations in the wake region of cylinder flow. The solid line is for GKS-DES, and the dashed line for GKS-IDDES. ◦ is for experiment of Parnaudeau et al. [58].  is for experiment of Lourenco and Shih [57]. × is for DNS of Ma et al. [59].

53 54 55 56 57 58 59 60 61 62 63 64 65 66

117 118 119

observed. In addition, influences of the GIS problems in DES are intuitively shown by the mean vorticity magnitude (< z >) in Fig. 13 where the shear layer is obviously shorter in GKS-DES simulation. Meanwhile, in Fig. 12, the ‘V’ shape velocity profile at x/ D = 1.06 is obtained. In IDDES, the grid scale is modified as [43],

 = min[max(C w d w , C w max ,  wn ), max ].

(29)

Here max is the maximum grid spacing in all three directions which also is the grid scale of DES and DDES. C w is 0.15. d w is the wall distance, and  wn is the grid spacing in the wall-normal direction. Based on Eq. (29), the grid scale becomes smaller in the near wall region where d w and  wn are both smaller than max .

The grid scale modification is capable of the faster generation of ‘new’ turbulence near the interface of RANS and LES, which is helpful to solve the MSD problem. The effects on the flow field can be seen in Fig. 14 in which more smaller vortex structures in the near wall wake region can be captured by IDDES. As a branch of IDDES, DDES is not considered independently in the present study. According to the assessments of Spalart [40], the overmuch empirical functions have made IDDES ‘less readable than that of DES97’. Nevertheless, IDDES is still recognized as one of the best hybrid method so far. In Figs. 11 and 12, the results of GKS-IDDES show good agreements with the experimental data of Parnaudeau et al. [58], as well as the DNS data. Besides the first-order statistics, the second-order statistics shown in Fig. 15 also agree well

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

12

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.12 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21 22 23 24

87

Fig. 13. The mean spanwise vorticity magnitude of DES and IDDES simulations.

88

Fig. 14. Instantaneous spanwise vorticity field. The red solid line is for positive value and blue dashed line is for negative value.

89 90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64 65 66

130

Fig. 15. The mean resolved streamwise (top, left), transverse (top, right) and shear (bottom) Reynolds stress at different locations in the wake region. For more details, see caption of Fig. 12.

131 132

JID:AESCTE AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.13 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Fig. 16. The power spectra of transverse velocity fluctuations at the center line of wake region. The solid line is for GKS-DES, dashed line for GKS-IDDES, and line with ‘×’ symbol for the PIV experimental data [58]. The spectra at x/ D = 1.55, 3.0, 5.0 are rescaled by 10−3 , 10−6 and 10−9 , respectively.

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

with experimental and DNS data, which is much better than GKSDES. To further validate the resolution of multiscale structures, the power spectra of the transverse velocity fluctuations at four locations in the wake region are calculated and shown in Fig. 16, where f vs is the Kármán vortex shedding frequency. At locations of x/ D = 0.60 and 1.55, the predicted spectra by GKS-DES and GKS-IDDES are similar, while the low frequency energy is higher in GKS-DES. The strong large scale fluctuations increase the transverse Reynolds stress component in GKS-DES as shown in Fig. 15. The difference of grid scale between turbulence models may be responsible for the results. In the near wall wake region, the grid scale is max for DES and C w max for IDDES as in Eq. (29). Furthermore, the base frequency with a peak value in the spectra f / f vs = 1 is well captured by both GKS-DES and GKS-IDDES. And the second harmonic frequency as f / f vs = 3 can also be accurately predicted at locations x/ D = 3.0 and 5.0, which agrees well with the experimental data [58]. In the wake region far from the wall, GKS-DES and GKS-IDDES describe the flow field on LES scale, thus the resolved small scales by these two models are similar. However, due to the limited grid resolution without special refinement in the wake region, the inertial region where the power spectra follows the −5/3 law is much shorter than the experimental data. In this high fidelity turbulence simulation, the developed DES/IDDES-GKS obtains satisfactory results compared with the experiments, which demonstrate that turbulence models on different scale can be nicely combined via τe in the extended BGK equation. Furthermore, from the numerical results it can also be revealed that there are still some risks or overmuch empirical corrections in current hybrid models. So the demands are still strong for developing a universal hybrid method with more solid physical basis.

57 58

5. Conclusions

59 60 61 62 63 64 65 66

In order to satisfy the engineering demands of efficient simulation of high-Reynolds-number turbulent flows, the extended GKS is developed based on the extended BGK equation with the effective relaxation time model. In the extended GKS, different turbulence models are employed, including several common-used RANS models for engineering simulation and LES, RANS/LES hybrid models for high fidelity simulation. The turbulent transport equations are

13

solved by the GKS with scalar transport in a coupled way, so that the numerical accuracy and robustness keep the same as those of the conservative variables. For super/hypersonic flows, a simple method is proposed to limit the gradient in the source terms of turbulence transport equations to suppress the artificial production near a shock wave. In all the cases, the extended GKS can get good results which agree well with the experimental data. When compared with existing numerical studies with same turbulence models, the present results also show improvements, especially in the shock capturing. This is mainly benefited from the self-adapting dissipation mechanism of GKS. The unified numerical scheme for the conservative variables and turbulence quantities maybe another factor. The numerical performance of extended GKS combined with different turbulence models adequately proves that the effective relaxation time model can uniformly describe turbulence on different scales. From this work, it is observed that the performances of numerical schemes will directly affect the assessments of turbulence models, which will further influence the development or improvement of the models. Thus it is important to develop a numerical scheme with high accuracy, high resolution and strong robustness for turbulence simulation. Besides the RANS simulations, the current extended GKS can simulate turbulence on multiscale with the help of the corresponding turbulence models. However, the multiscale mechanism is from the turbulence model, not the cross scale evolution in the general solution of the generalized BGK model. In the future work, it is very worthwhile to develop the multiscale turbulence model or simulation method by taking advantages of the multiscale evolution mechanism of the BGK equation.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

Conflict of interest statement

97 98

None declared.

99 100

Acknowledgements

101 102

This work is supported by the National Natural Science Foundation of China (11672158, 11172154, 10932005), National Key Basic Research and Development Program (2014CB744100) and Special Program for Applied Research on Super Computation of the NSFCGuangdong Joint Fund (the second phase) (U1501501). We also would like to acknowledge the technical support of PARATERA and the “Explorer 100” cluster system of Tsinghua National Laboratory for Information Science and Technology. References [1] K. Xu, A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method, J. Comput. Phys. 171 (1) (2001) 289–335, https://doi.org/10.1006/jcph.2001.6790. [2] Q.B. Li, S. Fu, K. Xu, Application of gas-kinetic scheme with kinetic boundary conditions in hypersonic flow, AIAA J. 43 (10) (2005) 2170–2176, https://doi. org/10.2514/1.14130. [3] K. Xu, M.L. Mao, L. Tang, A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow, J. Comput. Phys. 203 (2) (2005) 405–421, https:// doi.org/10.1016/j.jcp.2004.09.001. [4] H.Z. Tang, K. Xu, C.P. Cai, Gas-kinetic BGK scheme for three dimensional magnetohydrodynamics, Numer. Math. Theory Methods Appl. 3 (4) (2010) 387–404, https://doi.org/10.4208/nmtma.2010.m9007. [5] Q.B. Li, K. Xu, S. Fu, A high-order gas-kinetic Navier–Stokes flow solver, J. Comput. Phys. 229 (19) (2010) 6715–6731, https://doi.org/10.1016/j.jcp.2010. 05.019. [6] L. Pan, K. Xu, Q.B. Li, J.Q. Li, An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Euler and Navier–Stokes equations, J. Comput. Phys. 326 (2016) 197–221, https://doi.org/10.1016/j.jcp.2016.08.054. [7] S. Tan, Q. Li, A high-resolution gas-kinetic scheme with minimized dispersion and controllable dissipation reconstruction, Sci. China, Phys. Mech. Astron. 60 (11) (2017) 114713. [8] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [9] D.C. Wilcox, Turbulence Modeling for CFD, 3rd edition, DCW Industries, Inc., La Canada, CA, 2006.

103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

AID:4527 /FLA

[m5G; v1.235; Prn:18/04/2018; 10:34] P.14 (1-14)

S. Tan et al. / Aerospace Science and Technology ••• (••••) •••–•••

[10] C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Annu. Rev. Fluid Mech. 32 (1) (2000) 1–32, https://doi.org/10.1146/ annurev.fluid.32.1.1. [11] J. Fröhlich, D. von Terzi, Hybrid LES/RANS methods for the simulation of turbulent flows, Prog. Aerosp. Sci. 44 (5) (2008) 349–377, https://doi.org/10.1016/ j.paerosci.2008.05.001. [12] P. Sagaut, S. Deck, M. Terracol, Multiscale and Multiresolution Approaches in Turbulence: LES, DES and Hybrid RANS/LES Methods: Applications and Guidelines, World Scientific, 2013. [13] Q.B. Li, Numerical Study of Compressible Mixing Layer with BGK Scheme, Ph.D. thesis, Tsinghua University, Beijing, 2002. [14] S. Fu, Q.B. Li, Numerical simulation of compressible mixing layers, Int. J. Heat Fluid Flow 27 (5) (2006) 895–901, https://doi.org/10.1016/j.ijheatfluidflow. 2006.03.028. [15] G. Kumar, S.S. Girimaji, J. Kerimo, Gas-kinetic schemes for direct numerical simulations of compressible homogeneous turbulence, Phys. Rev. E 80 (2009) 046702, https://doi.org/10.1103/PhysRevE.80.046702. [16] G. Kumar, S.S. Girimaji, J. Kerimo, WENO-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence, J. Comput. Phys. 234 (2013) 499–523, https://doi.org/10.1016/j.jcp.2012.10.005. [17] Q.B. Li, S. Fu, High-order accurate gas-kinetic scheme and turbulence simulation (in Chinese), Sci. China, Ser. G, Phys. Mech. Astron. 44 (2014) 278–284, https://doi.org/10.1360/132013-62. [18] H.D. Chen, S. Kandasamy, S. Orszag, R. Shock, S. Succi, V. Yakhot, Extended Boltzmann kinetic equation for turbulent flows, Science 301 (5633) (2003) 633–636, https://doi.org/10.1126/science.1085048. [19] Q.B. Li, S. Tan, S. Fu, K. Xu, Numerical simulation of compressible turbulence with gas-kinetic BGK scheme, in: The Proceedings of ACFM13, 2010, pp. 1109–1112. [20] S. Tan, Q.B. Li, S. Fu, S. Zeng, Engineering simulation of turbulence with gaskinetic BGK scheme, in: Proceeding of ICFM6, vol. 1376, AIP Publishing, 2011, pp. 78–80. [21] S.W. Xiong, C.W. Zhong, C.S. Zhuo, K. Li, X.P. Chen, J. Cao, Numerical simulation of compressible turbulent flow via improved gas-kinetic BGK scheme, Int. J. Numer. Methods Fluids 67 (12) (2011) 1833–1847, https://doi.org/10.1002/fld. 2449. [22] J. Jiang, Y. Qian, Implicit gas-kinetic BGK scheme with multigrid for 3D stationary transonic high-Reynolds number flows, Comput. Fluids 66 (2012) 21–28, https://doi.org/10.1016/j.compfluid.2012.04.029. [23] M. Righi, A gas-kinetic scheme for turbulent flow, Flow Turbul. Combust. 97 (1) (2016) 121–139, https://doi.org/10.1007/s10494-015-9677-2. [24] S. Succi, O. Filippova, H.D. Chen, S. Orszag, Towards a renormalized lattice Boltzmann equation for fluid turbulence, J. Stat. Phys. 107 (1–2) (2002) 261–278, https://doi.org/10.1023/A:1014570923357. [25] M.D. Su, J.D. Yu, A parallel large eddy simulation with unstructured meshes applied to turbulent flow around car side mirror, Comput. Fluids 55 (2012) 24–28, https://doi.org/10.1016/j.compfluid.2011.10.017. [26] Q.B. Li, K. Xu, S. Fu, A new high-order multidimensional scheme, in: Computational Fluid Dynamics 2010, Springer, 2011, pp. 629–635. [27] K. Xu, Gas-Kinetic Schemes for Unsteady Compressible Flow Simulations, Von Karman Institute for Fluid Dynamics Lecture Series 1998-03, 1998. [28] J. Li, C. Zhong, D. Pan, C. Zhuo, A gas-kinetic scheme coupled with SST model for turbulent flows, Comput. Math. Appl. (2018), https://doi.org/10.1016/ j.camwa.2016.09.012, submitted for publication. [29] Q.B. Li, S. Fu, K. Xu, A compressible Navier–Stokes flow solver with scalar transport, J. Comput. Phys. 204 (2) (2005) 692–714, https://doi.org/10.1016/j.jcp. 2004.10.026. [30] Z.S. Sun, Y.X. Ren, C. Larricq, S.y. Zhang, Y.c. Yang, A class of finite difference schemes with low dispersion and controllable dissipation for DNS of compressible turbulence, J. Comput. Phys. 230 (12) (2011) 4616–4635, https:// doi.org/10.1016/j.jcp.2011.02.038. [31] Q.J. Wang, Y.X. Ren, Z.S. Sun, Y.T. Sun, Low dispersion finite volume scheme based on reconstruction with minimized dispersion and controllable dissipation, Sci. China, Phys. Mech. Astron. 56 (2) (2013) 423–431, https://doi.org/10. 1007/s11433-012-4987-z. [32] P.R. Spalart, S.R. Allmaras, A One-Equation Turbulence Model for Aerodynamic Flows, AIAA Paper 92-0439, 1992, https://doi.org/10.2514/6.1992-439. [33] S.L. Krist, R.T. Biedron, C.L. Rumsey, CFL3D user’s manual (version 5.0), 1998. [34] F. Menter, C. Rumsey, Assessment of Two-Equation Turbulence Models for Transonic Flows, AIAA Paper 94-2343, 1994, https://doi.org/10.2514/6.1994-2343. [35] C.L. Rumsey, Compressibility considerations for k–ω turbulence models in hypersonic boundary-layer applications, J. Spacecr. Rockets 47 (1) (2010) 11–20, https://doi.org/10.2514/1.45350. [36] T.J. Coakley, C.C. Horstman, J.G. Marvin, J.R. Viegas, J.E. Bardina, P.G. Huang, M.I. Kussoy, Turbulence Compressibility Corrections, NASA Technical Memorandum 108827, 1994. [37] Z. Zhang, Z. Gao, C. Jiang, C.-H. Lee, A RANS model correction on unphysical over-prediction of turbulent quantities across shock wave, Int. J. Heat Mass

[38]

[39]

[40] [41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50] [51]

[52]

[53]

[54] [55] [56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

Transf. 106 (2017) 1107–1119, https://doi.org/10.1016/j.ijheatmasstransfer. 2016.10.087. A. Jameson, W. Schmidt, E. Turkel, et al., Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge–Kutta Time-Stepping Schemes, AIAA Paper 81-1259, 1981, https://doi.org/10.2514/6.1981-1259. A. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications, Phys. Fluids 16 (10) (2004) 3670–3681, https:// doi.org/10.1063/1.1785131. P.R. Spalart, Detached-eddy simulation, Annu. Rev. Fluid Mech. 41 (1) (2009) 181–202, https://doi.org/10.1146/annurev.fluid.010908.165130. P. Spalart, W.H. Jou, M. Strelets, S.R. Allmaras, Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach, in: Advances in DNS/LES, vol. 1, 1997, pp. 4–8. P.R. Spalart, S. Deck, M.L. Shur, K.D. Squires, M.K. Strelets, A. Travin, A new version of detached-eddy simulation, resistant to ambiguous grid densities, Theor. Comput. Fluid Dyn. 20 (3) (2006) 181–195, https://doi.org/10.1007/s00162006-0015-0. M.L. Shur, P.R. Spalart, M.K. Strelets, A.K. Travin, A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities, Int. J. Heat Fluid Flow 29 (6) (2008) 1638–1649, https://doi.org/10.1016/j.ijheatfluidflow.2008.07.001. J.W. Slater, J.C. Dudek, K.E. Tatum, The NPARC Alliance Verification and Validation Archive, NASA Technical Report NASA/TM-2000-209946, ASME-2000-FED11233, 2000. J.G. Marvin, J.L. Brown, P.A. Gnoffo, Experimental Database with Baseline CFD Solutions: 2-D and Axisymmetric Hypersonic Shock-Wave/Turbulent-BoundaryLayer Interactions, NASA Technical Report NASA/TM-2013-216604/ARC-E-DAATN11460, 2013. T. Coratekin, J.V. Keuk, J. Ballman, Performance of upwind schemes and turbulence models in hypersonic flows, AIAA J. 42 (5) (2004) 945–957, https:// doi.org/10.2514/1.9588. G.M. Elfstrom, Turbulent hypersonic flow at a wedge-compression corner, J. Fluid Mech. 53 (Part 1) (1972) 113–129, https://doi.org/10.1017/ s0022112072000060. G.T. Coleman, J.L. Stollery, Heat transfer from hypersonic turbulent flow at a wedge compression corner, J. Fluid Mech. 56 (04) (1972) 741–752, https://doi. org/10.1017/S0022112072002630. F. Grasso, D. Falconit, High-speed turbulence modeling of shockwave/boundary-layer interaction, AIAA J. 31 (7) (1993) 1199–1206, https:// doi.org/10.2514/3.49062. O. Zeman, A New Model for Super/Hypersonic Turbulent Boundary Layers, AIAA Paper 93-0897, 1993, https://doi.org/10.2514/6.1993-897. Y.B. Suzen, K.A. Hoffmann, Investigation of Supersonic Jet Exhaust Flow by Oneand Two-Equation Turbulence Models, AIAA Paper 98-0322, 1998, https://doi. org/10.2514/6.1998-322. T.P. Sommer, R.M.C. So, H.S. Zhang, Near-wall variable-Prandtl-number turbulence model for compressible flows, AIAA J. 31 (1) (1993) 27–35, https:// doi.org/10.2514/3.11314. X.D. Xiao, J.R. Edwards, H.A. Hassan, R.L. Gaffney, Role of turbulent Prandtl numbers on heat flux at hypersonic Mach numbers, AIAA J. 45 (4) (2007) 806–813, https://doi.org/10.2514/1.21447. H. Babinsky, J.K. Harvey, Shock Wave-Boundary-Layer Interactions, vol. 32, Cambridge University Press, 2011. R.D. Moser, J. Kim, N.N. Mansour, Direct numerical simulation of turbulent channel flow up to Reτ = 590, Phys. Fluids 11 (4) (1999) 943–945. C. Norberg, Effects of Reynolds Number and a Low-Intensity Freestream Turbulence on the Flow Around a Circular Cylinder, Technological Publication 87/2, Chalmers University of Technology, 1987. L.M. Lourenco, C. Shih, Characteristics of the plane turbulent near wake of a circular cylinder, a particle image velocimetry study, 1993 (data published in Refs. [58,63]). P. Parnaudeau, J. Carlier, D. Heitz, E. Lamballais, Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900, Phys. Fluids 20 (8) (2008) 085101, https://doi.org/10.1063/1.2957018. X. Ma, G.S. Karamanos, G.E. Karniadakis, Dynamics and low-dimensionality of a turbulent near wake, J. Fluid Mech. 410 (2000) 29–65, https://doi.org/10.1017/ S0022112099007934. A.G. Kravchenko, P. Moin, Numerical studies of flow over a circular cylinder at ReD = 3900, Phys. Fluids 12 (2) (2000) 403–417, https://doi.org/10.1063/1. 870318. V. D’Alessandro, S. Montelpare, R. Ricci, Detached-eddy simulations of the flow over a cylinder at Re = 3900 using OpenFOAM, Comput. Fluids 136 (2016) 152–169, https://doi.org/10.1016/j.compfluid.2016.05.031. S. Chen, Z. Xia, S. Pei, J. Wang, Y. Yang, Z. Xiao, Y. Shi, Reynolds-stressconstrained large-eddy simulation of wall-bounded turbulent flows, J. Fluid Mech. 703 (2012) 1–28, https://doi.org/10.1017/jfm.2012.150. P. Beaudan, P. Moin, Numerical Experiments on the Flow Past a Circular Cylinder at Sub-Critical Reynolds Number, Report TF-62, Thermosciences Division, Stanford University, 1994.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132