A large-eddy simulation method for low Mach number flows using preconditioning and multigrid

Computers & Fluids 35 (2006) 1126–1136 www.elsevier.com/locate/compfluid

A large-eddy simulation method for low Mach number flows using preconditioning and multigrid N. Alkishriwi *, M. Meinke, W. Schro¨der Aerodynamisches Institut, RWTH Aachen, Wu¨llnerstr. zw. 5 u. 7, 52062 Aachen, Germany Received 10 September 2004; received in revised form 20 January 2005; accepted 15 June 2005 Available online 19 September 2005

Abstract Large-eddy simulation (LES) has become one of the major tools to investigate the physics of turbulent compressible and incompressible flows. At low Mach numbers the performance of LES codes developed for the conservation equations of compressible fluids deteriorates due to the presence of two different time scales associated with acoustic and convective waves. In many subsonic turbulent flows low Mach number regions exist, which require large integration times until a fully developed flow is established. In such cases, the efficiency of algorithms for compressible flows can be improved considerably by low Mach number preconditioning methods. In this paper an efficient method of solution for low subsonic flows is developed based on an implicit dual time stepping scheme combined with low Mach number preconditioning and a multigrid acceleration technique. To validate the efficiency and the accuracy of the method, large-eddy simulations of turbulent channel flow at Res = 590 and cylinder flow at Re = 3900 are performed for several Mach numbers and the data are compared with numerical and experimental findings from the literature. The speedup compared to a purely explicit approach is in the range of 6–40.  2005 Elsevier Ltd. All rights reserved.

1. Introduction In many engineering problems compressible and nearly incompressible flow regimes occur simultaneously. In aerodynamics, for example, the flow over a multielement airfoil at high angle of attack, contains transonic flow regions even for small free stream Mach numbers. Other examples are low speed flows, which may be compressible due to surface heat transfer or volumetric heat addition. The numerical simulation of such flows requires to solve the conservation equations for compressible fluids to obtain the correct solution. When a compressible flow solver is applied to a nearly incompressible flow, its performance can deteriorate in terms of both speed and accuracy [2,3]. It is well *

Corresponding author. Tel.: +49 241 80 954 10; fax: +49 241 80 922 57. E-mail address: [email protected] (N. Alkishriwi). 0045-7930/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2005.06.002

known that most compressible codes do not converge to an acceptable solution when the Mach number of the flow field is smaller than Oð101 Þ. The main difficulty with such low speed flows arises from the large disparity between the wave speeds. The acoustic wave speed is ju ± cj, while entropy or vorticity waves travel at juj, which is quite small compared to ju ± cj. In explicit time-marching codes, the acoustic waves define the maximum time step, while the convective waves determine the total number of iterations such that the overall computational time becomes large for small Mach numbers. During the past years there has been a growing demand in the CFD community to solve such mixed flow problems by modifying the existing compressible flow solvers. Different methods have been proposed to address this problem [1–4]. One of most popular approaches to solve this problem, however, is to use low Mach number preconditioning methods for compressible codes [1–3,7–14]. The basic idea is to modify

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

the time marching behavior of the system of equations without altering the steady state solution. This is, however, only useful when the steady state is sought. Low Mach number preconditioners are evaluated using purely local information from the cell the residual of which will be preconditioned. Large-eddy simulations for low Mach number flows, i.e., flows that can be considered almost incompressible, have hardly been discussed in the context of solutions of the compressible Navier–Stokes. That is, generally large-eddy simulations of nearly incompressible fluids have been performed by incompressible LES solvers and not by especially formulated compressible LES solvers. The difficulty to design such a general LES solver lies in the strong susceptibility of the small scales to the dissipation level contained in the numerical scheme, which is changed, however, when an efficient and accurate flow solver for low Mach number flows with, e.g., preconditioning is developed. This means, the dilemma of a compressible LES solver is on the one hand, preconditioning is required to satisfy the physics of low Mach number flows while on the other hand, the artificial dissipation has to be well below the dissipation of the high wave number vortices, i.e., the small scales. Dailey and Pletcher [14] implemented multigrid for a preconditioned dual time step algorithm for the two-dimensional Navier–Stokes. They observed clear speedups for steady and unsteady simulations, however, the improvements for unsteady flows were problem dependent. Wang and Pletcher [35] performed a large-eddy simulation of a planar channel flow with significant heat transfer at a low Mach number using a method for the compressible governing equations to study effects of fluid property variations on the near-wall turbulence structure. They also used low Mach number preconditioning, but they did not elaborate on the impact on the rate of convergence. Chidambaram and Pletcher [36] applied a compressible LES solver with preconditioning and dual time stepping to isotropic decaying turbulence. Their research focused on the development of a colocated-grid algorithm with a remedy for the pressure-velocity decoupling problem. In other words, the bearing of preconditioning and dual time stepping on the numerical performance was not discussed. Lessani et al. [26] performed large-eddy simulations using preconditioning and multigrid for low Mach number turbulent channel and cylinder flows. They demonstrated for case studies at Mach numbers of 0.033 and 0.06 the impact of the dual time stepping approach and multigrid technique on the efficiency of the LES method by achieving a speedup factor of 4–7 compared to a basic explicit approach. In the present paper, a large-eddy simulation method is presented using an implicit dual time stepping scheme combined with preconditioning and multigrid, which is based on a different approximation of the inviscid terms and the subgrid scales compared to the Lessani et al.

1127

scheme [26]. This method leads to an efficiency increase up to a factor of 40 for problems at a Mach number 0.01, which is primarily due to the preconditioning approach. The paper is organized as follows. We start by a brief discussion of the LES formulation and filtering procedure of the governing equations. Then, the implementation of preconditioning in the LES context using the dual time stepping technique [18] is described. Subsequently, the numerical method is outlined. The time marching solution technique and the discretization within the dual time stepping approach are discussed in separate sections. Finally, the numerical results are analyzed followed by a thorough discussion of the impact of the dual time stepping approach and the multigrid technique on the overall efficiency of the LES method.

2. Governing equations The governing equations are the unsteady threedimensional compressible Navier–Stokes written in generalized coordinates ni, i = 1, 2, 3 oQ oðFci  Fvi Þ þ ¼ 0. ot oni

ð1Þ

The spatially filtered form of an arbitrary variable U(ni, t) is given by Z 1 Uðni ; tÞ ¼ Gðni ÞUðni ; tÞ dV ; ð2Þ V V where V is the entire domain and G(ni) is a filter function, which satisfies the normalization condition Z Gðni ÞdV ¼ 1. ð3Þ V

This function decomposes U(ni, t) into a large scale Uðni ; tÞ and a small scale or subgrid scale U0sgs ðni ; tÞ part Uðni ; tÞ ¼ Uðni ; tÞ þ U0sgs ðni ; tÞ.

ð4Þ

The problem of establishing the appropriate form of the averaged equations of a compressible fluid, which contains also temperature and density fluctuations, can be simplified by using the density-weighted averaging procedure suggested by Favre [16] e i ; tÞ ¼ qðni ; tÞUðni ; tÞ . Uðn ðni Þ q

ð5Þ

Thus, a variable is decomposed into its Favre-filtered e i ; tÞ and a fluctuating component U00 (ni, t). component Uðn The quantities Q, Fci , and Fvi are given by

1128

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

0 1 0 1 c 1 0 fi U q  q B oni r C B C B oxj ~j1 þ sj1 C Bq B C f c u þ oni p C B C ox B U i ~ C B q~u C B B C B C on 1 ~j2 þ sj2 C C; fi c~v þ oni p C; Fvi ¼ 1 B oxji r C; Fci ¼ 1 B q   ~ U Q¼ B q v oy C C C JB JB JB B oni C B C B C c on  ~ q w ~j3 þ sj3 C r B Bq C @ A i f  ox ~ U w þ p j i @ A @ A oz oni ~ ~e q b þj fi c q fi c p ~e þ U U oxj j 0

ð6Þ where the quantity Q represents the vector of conservative variables, J is the metric Jacobian, and Fci , F vi are inviscid and viscous flux vectors, respectively. The contravariant velocity in the i-direction is given by on on fi ¼ oni ~ ~. U u þ i ~v þ i w ox oy oz c

ð7Þ

~ij are expressed as The stress terms r   o~ ui onk o~ uj onk 2 ou~k onl e ~ij ¼ lð T Þ r þ  dij ; onk oxj onk oxi 3 onl oxk

In two-dimensions n = n1 and g = n2 Eq. (1) can be written as oQ oEc oFc oEv oFv þ þ þ þ ¼ 0. ð12Þ ot on og on og The eigenvalues of the q above system in the n-direction ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

2

are Un, Un, and U n  c ðnx Þ þ ðny Þ , where Un is the contravariant velocity and c is the speed of sound. When the Mach number approaches zero, these eigenvalues differ over several orders of magnitude. The objective of low Mach number preconditioning is to modify the wave speeds such that all eigenvalues possess the same order of magnitude. The rescaling of the characteristic wave speeds improves the convergence and leads to a more accurate solution. The preconditioned Navier–Stokes read

ð8Þ oQ ¼ CR; ot

ð13Þ

where lð Te Þ is the molecular viscosity based on the Favre-filtered static temperature Te and the subscripts define the coordinate directions. The quantities ~bj read ! e oni o T ~ ~ji þ kð Te Þ bj ¼ ~ ui r ; ð9Þ oni oxj

where R represents   oEc oFc oEv oFv R¼ þ þ þ on og on og

~ij Þ is neuj r in which the subgrid viscous work ðuj rij  ~ glected as in Knight et al. [17]. The pressure is calculated by    2 q ~e  ~ k sgs ; p ¼ ðc  1Þ q ð10Þ u q 2 i

and C is the preconditioning matrix, which is to be defined such that the new eigenvalues of the preconditioned system of equations are of similar magnitude. TurkelÕs preconditioning matrix [2] C transformed into conservative variables reads

k sgs is the subgrid scale kinetic energy per unit where q volume represented by the following relation 1 k sgs ¼ q ðg q ui ~ ui Þ. ui ui  ~ 2

ð11Þ c

c f uj Þ ð Ug The subgrid-scale stresses (SGS) sij ¼ q i uj  U i ~ c c g  ~e U fi Þ ð eU and the subgrid-scale energy term  i ¼ q

2 6 6 6 6 C¼6 6 6 4

3

M 2 ðc1ÞH þ1 2

uð1cÞH c2

vð1cÞH c2

ðc1ÞH c2

uM 2 ðc1ÞH 2

u2 ð1cÞH þ1 c2

uvð1cÞH c2

uðc1ÞH c2

vM 2 ðc1ÞH 2 c2 M 2 H 2

ð14Þ



M 2 ðc1Þ 2

þ1



2

v ð1cÞH þ1 c2   2  2 Hu M ð1cÞ  1 Hv M ð1cÞ 1 2 2 uvð1cÞH c2



vðc1ÞH c2 M 2 ðc1ÞH 2

þ

7 7 7 7 7. 7 7 5

ð15Þ

i

c

gc  ~ fi Þ can be modelled. Here, the large-eddy pU þð pU i simulations are carried out using the monotone integrated LES (MILES) approach [15] to represent the effect of non-resolved subgrid scales. The numerical algorithm defines the entire energy transfer between resolved and subgrid scales. This requires a discretization with a small dissipation. For that purpose a modified AUSM method with a central approximation for the pressure derivative is used, which has been validated for numerous test cases in [5]. 3. Preconditioned governing equations As mentioned before, preconditioning is required to remove the stiffness of the Navier–Stokes for compressible flow at low Mach numbers and to provide an efficient and accurate method of solution.

The new eigenvalues of the preconditioned conservation equations in the n-direction are U n; U n;

Un 1 ð1 þ Þ  2 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi H2 U 2n þ 4ðn2x þ n2y Þc2 ;

ð16Þ

where  is the steady state preconditioning parameter, which is defined by  ¼ MIN½MAXðM 2 ; M 2min Þ; 1

ð17Þ

and the quantity H is related to  via H =   1. The local Mach number is denoted by M and Mmin is a cutoff value to prevent the preconditioner from becoming singular. In our steady state analyses the best results are obtained by setting Mmin proportional to the free stream Mach number M1 M 2min ¼ M 21 =2.

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

The definition of the preconditioning parameter  ensures the wave speeds given in Eq. (16) to be of the same order of magnitude. Further below it will become clear that this steady state preconditioning formulation can be also applied when unsteady solutions are sought. 4. Numerical method The LES is carried out using a modified second-order accurate AUSM scheme using a centered 5-point low dissipation stencil to compute the pressure derivative in the convective fluxes. This scheme is described in detail and validated against several test problems such as channel, jet, and pipe flows in Meinke et al. [5]. The method is formulated for multiblock structured curvilinear grids and implemented on vector and parallel computers. The viscous stresses are discretized to secondorder accuracy using central differences, i.e., the overall spatial approximation is second-order accurate. 4.1. Basic time-stepping scheme A 5-step Runge–Kutta method is used to propagate the vector of solution Q in Eq. (13) from time level n to n + 1 Q0 ¼ Qn ; .. . Ql ¼ Q0  al DtR ðQl1 Þ; .. .

ð18Þ

Qnþ1 ¼ Q5 . The superscripts l, n denote, respectively, the step index l = 1, . . . , 5 and the time level, Dt = tn + 1  tn is the time step, and R* = CR represents the preconditioned inviscid and viscous residual. The

Runge–Kutta coefficients are 24 al ¼ 246 ; 244 ; 249 ; 12 ; . They are optimized for maxi24 24 mum stability of a centrally discretized scheme. To reduce the computational cost, local time stepping in conjunction with the implicit dual time stepping approach is used, where the temporal integration in pseudo-time for the preconditioned system is performed by exploiting the multigrid technique and local time stepping. 4.2. Dual-time stepping algorithm The low Mach number preconditioning introduced in Eq. (13) changes the time derivative of the Navier– Stokes without affecting the steady state solutions. For unsteady flow problems, a dual time stepping technique [18–20] for time accurate solutions is used. In this approach, the solution at the next physical time step is determined as a steady state problem to which preconditioning, local time stepping and multigrid is applied.

1129

Introducing of a pseudo-time s the unsteady governing equations with preconditioning are C1

oQ oQ þ þ R ¼ 0. os ot

ð19Þ

Note that only the pseudo-time terms are altered by the preconditioning, while the physical time and space derivatives retain their original form. This means, the acceleration techniques discussed in Section 3 for steady state computations can be immediately utilized to speed up the convergence within each physical time step to obtain an accurate solution for unsteady flow. Following the approach of Jameson [25], derivatives with respect to the real time t are discretized using a three-point backward difference scheme that results in an implicit scheme, which is second-order accurate in time  nþ1  oQ 3Q  4Qn þ Qn1 nþ1 ¼ RHS ¼ C þ RðQ Þ . os 2Dt ð20Þ Note that at s ! 1 the first term on left-hand side of Eq. (19) vanishes such that Eq. (1) is recovered. To advance the solution of the inner pseudo-time iteration, again a 5-step Runge–Kutta method in conjunction with local time stepping and multigrid is used. That is, in Eq. (18) RHS has to be substituted for R*. Arnone et al. [21] pointed out that the above multistage scheme becomes unstable when the physical time step Dt is of the order of the pseudo-time step Ds or smaller. This result was substantiated by Melson et al. [22] who demonstrated that the instability is caused by the nþ1 term 3Q2Dt , which becomes the dominant expression for small Dt. They suggested an implicit treatment of this term yielding the following formulation for the lth stage  1 3Ds l 0 al C RHS. Q ¼ Q  al Ds I þ 2Dt The additional term means that in smooth flows the development in pseudo-time is proportional to the evolution in t. In unsteady flows the preconditioning parameter  possesses the more general form [20,24]  ¼ MIN½MAXðM 2 ; M 2t ; M 2min Þ; 1.

ð21Þ

l The additional quantity M t ¼ pDtc , can be considered a non-dimensional signal propagation speed defined by the numerical approach and l is the maximum characteristic length scale of the flow problem.

4.3. Multigrid method The multigrid method [34] is based on a full approximation storage (FAS) scheme. Following Turkel [3], we transfer the residuals to the next grid based on the preconditioned system, since these residuals are more balanced than the non-preconditioned residuals. A 3-level sawtooth multigrid cycle is used in all calculations. To

1130

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

go from the fine to the coarse mesh, which contains every second grid point of the fine mesh, a residual based full weighting restriction operator is used. On the coarsest mesh the same iterative method of solution as on the finest mesh is applied. The corrections are transferred from the coarse to the fine grid by a trilinear prolongation operator. In general, 1 iteration is performed on the finest and each intermediate grid level, and 2 iterations on the coarsest mesh before the prolongation. More details of the multigrid approach are given in Meinke [6].

5. Results To show the improvement in the rate of convergence and the overall accuracy when an LES solver for the compressible Navier–Stokes is applied to flows at almost vanishing compressibility effects, the above described numerical method based on implicit dual-time stepping, multigrid, and preconditioning has been used to perform large-eddy simulations for two well-documented internal and external flow problems in the literature. First, the turbulent channel flow at a Reynolds number Res = 590 based on the friction velocity us is considered and then, the turbulent flow past a circular cylinder at a diameter based Reynolds number of ReD = 3900 is computed. For all flows the reference Mach number lies in the range of 0.01 6 Ma < 0.1. Finally, the impact of the dual-time stepping and the multigrid method on the computational efficiency is discussed.

direction the pressure and temperature fluctuations and the mass flow are assumed to be periodic. Using from Hussain and Reythe relation Res ¼ 0:1097Re0:911 cl nolds [33] the pressure difference dp ¼ qu2s L=H that drives the flow is determined by the maximum velocity. After a fully developed turbulent channel flow has been established the simulations are continued for another 200 dimensionless time units to compute the statistical results. Samples are stored each quarter of a dimensionless time unit. The profiles are generated by time averaging over a time interval of 200 dimensionless time units followed by combined streamwise and spanwise averaging over the entire channel. A direct numerical simulation for an incompressible flow at Res = 590 carried out by Moser et al. [31] on a mesh of 384 · 257 · 384 cells is used as reference solution for the LES results. The comparison of the DNS data and the LES findings on the fine grid at Ma = 0.01 and Dt = 0.01 in Figs. 1–3 show overall a convincing agreement for the mean velocity profile, turbulence intensity distributions, and resolved Reynolds stress component profile of the fine grid although the peak value in the streamwise fluctuations is slightly overpredicted. Furthermore, it is evident from Figs. 1–3 that besides the streamwise Reynolds stresses in Fig. 2 the results on the coarse and fine grids possess a convincing agreement such that only a slight grid dependence of the solution is observed. The over-

25

20

5.1. Channel flow

+

15 U

The channel flow computations at Res = usd/m = 590 correspond to Recl = ucld/m = 12448, where ucl is the maximum velocity on the centerline and d is the channel half-width. The Mach numbers based on the centerline velocity are set to Ma = 0.01, 0.05, 0.1. The size of the computational domain in the streamwise x, spanwise z, and normal y-direction is given by L = 2pd, B = pd, and H = 2d. The simulations are performed on two different grids, i.e., we use a fine grid with 2,063,205 points and a coarse grid with 887,841 points resulting in a grid refinement factor of approximately 2.3. The details of the meshes are summarized in Table 1. The boundary conditions consist of no-slip and isothermal conditions on the walls and periodic boundary conditions in the spanwise direction. In the streamwise

10

5

0

DNS MKM [31] LES-Dual Time Stepping-Fine Grid LES-Dual Time Stepping-Coarse Grid Log Law

5

30 y+

180

640

Fig. 1. LES of turbulent channel flow of the fine and coarse grids at Res = 590, Ma = 0.01, and physical time step Dt = 0.01; the data are given in inner wall units and compared with DNS results of [31]; mean velocity distribution.

Table 1 Mesh parameters of LES channel flow; the spatial steps are given in inner coordinates, e.g., Dy+ = Dyus/m LES channel flow

Dx+

Dy þ min

Dy þ max

Dz+

Nx · Ny · Nz

Total number of points

Fine grid Coarse grid

24.4 33.1

2.00 2.00

21.3 40.8

20.1 23.2

153 · 145 · 93 113 · 97 · 81

2,063,205 887,841

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136 3.5

20 2.5

15 +

2

U

σ(u’’)/uτ σ(w’’)/uτ σ(v’’)/uτ

25

DNS MKM [31] LES-Dual Time Stepping-Fine Grid LES-Dual Time Stepping-Coarse Grid

3

1131

1.5

10

1

LES-DT=0.01 Explicit LES-DT=0.001 LES-DT=0.025 LES-DT=0.1 Log Law

5 0.5 0

0

150

300 y+

450

0

600

Fig. 2. Distribution of the streamwise, spanwise, and normal Reynolds stresses. For details see caption of Fig. 1.

5

30 + y

180

640

Fig. 4. LES of turbulent channel flow at Res = 590 and Ma = 0.01; the data are given in inner wall units; mean velocity distribution at various time steps Dt.

0 0

-0.2

σ(u’’v’’)/uτ

σ(u’’v’’)/uτ

-0.2

-0.4

-0.4

-0.6

-0.6 LES-DT=0.01 LES-Explicit LES-DT=0.001 LES-DT=0.1

-0.8 -0.8 DNS MKM [31] LES-Dual Time Stepping-Fine Grid LES-Dual Time Stepping-Coarse Grid

-1 0

150

300 y+

450

-1

0

600

Fig. 3. Distribution of the Reynolds shear stress component (u00 v00 ). For details see caption of Fig. 1.

prediction of the maximum value in the streamwise fluctuations, which is emphasized in the coarse grid solution, is due to an insufficient turbulent energy redistribution from the streamwise to the spanwise direction, caused by slightly too weak a coherent near wall structure. We now turn to the discussion of the numerical results. It is well known that the backward differencing scheme has a dissipation error proportional to the time step [32]. To investigate the dependence of the solution on the size of the physical time step, several simulations at different values of Dt ranging from 0.001 to 0.1 are performed. Figs. 4–7 show the results at Ma = 0.01 for Dt = 0.001, Dt = 0.01, Dt = 0.025, and Dt = 0.1 to evidence the impact of an increasing time step on the accuracy of the solution.

150

300 y+

450

600

Fig. 5. Reynolds shear stress component (u00 v00 ) profile. For details see caption of Fig. 4.

Up to Dt = 0.01, which is 180 times larger than the maximum explicit time step, an accurate solution is obtained, i.e., hardly any deviations occur compared with the explicit or the Dt = 0.001 distributions of the Reynolds shear stress in Fig. 5 or of the normal and spanwise Reynolds stresses in Figs. 6 and 7. Furthermore, note that an acceptable agreement with the reference solution, i.e., the explicit data, can be obtained up to a physical time step of Dt = 0.025. This increase in Dt is achieved, since the modified AUSM based LES method possesses low overall dissipation [5]. 5.2. Cylinder flow The flow around a cylinder is performed at a Reynolds number ReD = 3900 based on the diameter D, a free stream Mach number M1 = 0.05, and a physical

1132

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

LES-DT=0.01 LES-Explicit LES-DT=0.001 LES-DT=0.025 LES-DT=0.1

1.4 1.2

σ(v’’)/uτ

1 0.8 0.6 0.4 0.2 0

0

150

300 y+

450

600

Fig. 6. Distribution of the normal Reynolds stress component. For details see caption of Fig. 4.

Fig. 8. C-type computational mesh for the flow over a circular cylinder.

2 LES-DT=0.01 LES-Explicit LES-DT=0.001 LES-DT=0.025 LES-DT=0.1

1.5

1.2 1

0.6 1

u/U∞

σ(w’’)/uτ

0.8

0.4 0.2

0.5

0

LES-Dual Time Stepping LES-Explicit LES-Incomp. solver [27] EXP. [29] EXP. [28]

-0.2 0

0

150

300 + y

450

600

Fig. 7. Distribution of the spanwise Reynolds stress component. For details see caption of Fig. 4.

time step Dt = 0.02. The computational domain extends 25D, 15D, and 1D in the streamwise, normal and spanwise direction. The C-type mesh, the near-wall resolution of which is shown in Fig. 8, possesses 105 · 273 · 49 in the circumferential, normal, and spanwise direction resulting in approximately 1.4 · 106 cells. The mesh is clustered near the wall such that the minimum normal spatial step is 0.013D. No-slip and isothermal conditions are specified on the walls and periodic boundary conditions are imposed in the spanwise direction. On all subsonic inflow and outflow boundaries non-reflecting boundary conditions with pressure relaxation are applied [23]. To show the quality of the results the LES findings are compared with experimental data from Lourenco and Shih [28] and Ong and Wallace [29] and LES results

-0.4 2

4

6

8

10

x/D Fig. 9. Streamwise velocity as a function of X/D on the centerline Y/D = 0, Z/D = 0 in the cylinder wake at ReD = 3900.

from Beaudan and Moin [27]. To further evidence the impact of the preconditioning on the solution several distributions computed with and without preconditioning are also presented. The statistics are based on a non-dimensional time interval of Te ¼ Tu1 =D ¼ 170 plus averaging in the spanwise direction. Fig. 9 shows the development of the streamwise velocity component on the centerline in the wake of the cylinder. Unlike the LES distribution computed using a pure explicit scheme without preconditioning at a Mach number M1 = 0.1 the new LES method captures the location and the minimum value in the u-distribution measured by Lourenco and Shih [28]. The agreement is even better than that of the LES of the incompressible Navier–Stokes from Beaudan and Moin [27].

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136 0.35

LES-Dual Time Stepping LES-Incomp. solver [27] Exp. [28]

0.3 0.25 2 uu/U∞

Such a convincing correspondence is also observed further downstream in the wake and in the normal distribution of the u-component at X/D = 1.06 presented in Fig. 10. Unlike the u-component profile from [27] and the LES simulation without preconditioning a satisfactory agreement with the experimental data from Lourenco and Shih [28] is still achieved at X/D = 2.02 in Fig. 11. The profiles of the velocity fluctuations of the streamwise and vertical components at X/D = 1.54 as a function of Y/D presented in Figs. 12 and 13 corroborate the correspondence with the measurements from [28]. The comparison with LES and experimental data from the literature with respect to the distribution of the pressure coefficient cp on the cylinder surface in Fig. 14 and the mean flow data comprised in Table 2 confirms the convincing quality of the results of the

1133

0.2 0.15 0.1 0.05 0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

2

y/D

Fig. 12. Streamwise velocity fluctuation as a function of Y/D in the cylinder wake at X/D = 1.54.

2 0.5

LES-Dual Time Stepping LES-Incomp. solver [27] Exp. [28]

1

0.4

0.5

0.3 2 vv/U∞

u/U∞

1.5

0

0.2

-0.5

LES-Dual Time Stepping LES-Explicit LES-Incomp. solver [27] Exp. [28]

-1 -1.5 -3

-2

-1

0

1

2

0.1

3

0 -2

-1.5

-1

y/D

0.5

1

1.5

2

Fig. 13. Vertical velocity fluctuation as a function of Y/D in the cylinder wake at X/D = 1.54.

LES method using preconditioning, dual time stepping, and multigrid.

2 1.5

5.3. Efficiency analysis

1 u/U∞

0 y/D

Fig. 10. Streamwise velocity as a function of Y/D in the cylinder wake at X/D = 1.06.

0.5 0 LES-Dual Time Stepping LES-Explicit LES-Incomp. solver [27] Exp. [28]

-0.5 -1 -1.5

-0.5

-1

-0.5

0

0.5

1

1.5

y/D Fig. 11. Streamwise velocity as a function of Y/D in the cylinder wake at X/D = 2.02.

The efficiency of the implicit dual time stepping scheme is defined as the ratio of the computational work Cwk required to advance the solution over a certain physical time interval of the explicit versus the implicit scheme Eim = Cwk,ex/Cwk,im. That is, the smaller Cwk,im compared to Cwk,ex, the larger the efficiency of the implicit scheme Eim. This ratio Eim depends mainly on the size of the physical time step, the maximum time step of the explicit scheme, and the number of inner iterations in the implicit scheme required to satisfy a convergence limit for the inner pseudo-time cycle. Let U1 = Dtphy/Dsmax,ex be the ratio of the physical time step to the maximum time step given by the stability

1134

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

1.5

0.01

LES-Dual Time Stepping LES-Explicit LES-Incomp. solver [27] EXP. [30]

1

0.001

0.5

0.0001

0

Residual

cp

Phy. Time Step DT =0.001 Phy. Time Step DT=0.005 Phy. Time Step DT=0.01 Phy. Time Step DT=0.025 Phy. Time Step DT=0.05 Phy. Time Step DT=0.075 Phy. Time Step DT=0.1

-0.5

1e-05

1e-06

-1 1e-07

-1.5 0

20

40

60

80

100 120 140 160 180

1e-08

θ Fig. 14. Pressure coefficient distribution on the cylinder surface.

1e-09

Table 2 Cylinder flow; experimental Strouhal number and bubble length from [37] and total drag from [30] ReD = 3900

Strouhal number

Total drag

Bubble length

Experimental data LES [38] Present LES

0.215 ± 0.005 0.210 0.217

0.98 ± 0.05 1.04 1.05

1.33 ± 0.2 1.35 1.31

5

10

15 20 25 30 35 Number of inner iterations

40

45

Fig. 15. Convergence of the inner iterations at Ma = 0.01.

M=0.01 M=0.05

70 60 50

Eim

analysis of the pure explicit scheme that defines the basis of the dual time stepping procedure. That is, in this analysis Dsmax,ex is the time step of the 5-step Runge–Kutta scheme. In other words, U1 describes the number of explicit time steps to proceed a time interval Dtphy. If the quantity U2 represents the total number of inner iterations, then the ratio U1/U2 defines the gain of the implicit dual time stepping method over the explicit scheme. This means, the efficiency Eim is equal to U1/U2. Furthermore, it is well known that the explicit time step has to be reduced at decreasing Mach numbers due to the discrepancy in the time scales of the acoustic and convective waves. For this reason, we will address in the following the impact of the size of the physical time step Dtphy and of the Mach number on the efficiency of the implicit scheme. Note that the size of the physical time step is determined only by the accuracy required, whereas the explicit time step is defined by the stability criteria and the required resolution, which is a direct function of the Reynolds number. We performed several computations at different time steps and various Mach numbers to evidence the impact on the overall efficiency Eim. Figs. 15 and 16 are based on channel flow solutions, whereas Fig. 17 represents the convergence behavior of the cylinder flow at M1 = 0.05. Fig. 15 shows the convergence of the inner iterations at Ma = 0.01. It is evident that as the size of the physical time step increases, the required number of iterations in the inner pseudo-time cycle grows. If a

40 30 20 10 0 0

0.02

0.04

0.06

0.08

0.1

Physical time step Δtphy Fig. 16. Efficiency Eim as a function of Dtphy at Ma = 0.01 and Ma = 0.05.

reduction of three-orders of magnitude suffices for a converged inner solution, and this is valid since the channel and cylinder flow results discussed before are based on such a convergence criterion, less than 25 inner iterations are necessary at Dtphy 6 0.025. Fig. 16 demonstrates the dramatic impact of the Mach number on the efficiency of the implicit scheme. The results, which are based on the fine grid solutions, show a speedup in the range of 9 to 60 compared to the explicit scheme for a reduction of two-orders of magnitude for a Ma = 0.01 flow. Note, that for a

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

MG-with Preconditioning SG-with Preconditioning MG-without Preconditioning SG-without Preconditioning

0.1 0.01

Residual

0.001 0.0001

1135

and thus of the overall solution method increases. However, the accuracy of the LES solution decreases due to the larger dissipation. To be more precise, especially the turbulent structures in the inertial subrange are affected by the multigrid method. From this, it can be concluded that preconditioning is the major component to improve the efficiency of a compressible LES method for very low Mach number flows and that the speedup due to multigrid is clearly limited due to physical constraints.

1e-05

6. Conclusion 1e-06 1e-07 20

40

60

80 100 120 140 160 180 200 Work units

Fig. 17. Convergence of single and multigrid methods with and without preconditioning.

reduction of three-orders of magnitude this speedup is lowered to a range of 6–26. Lessani et al. [26] propose for turbulent channel flow at Res = 180 a residual drop of 2.5 to be sufficient for accurate results. This leads in the cases considered to a maximum speedup of 40 compared to the explicit scheme. This efficiency gain drops dramatically when the compressibility is increased, i.e., when the Mach number is enlarged from Ma = 0.01 to Ma = 0.05. Note, however, that this result is also due to the relation between the Mach number and the explicit time step. It is quite clear that a large-eddy simulation based on the compressible conservation equations, which is used to seek solutions in the very low Mach number range, will deliver accurate results (Figs. 1–3 and Figs. 9–14) and can be considered efficient only if the preconditioning ansatz is implemented. Fig. 17 shows the impact of preconditioning and multigrid on the rate of convergence, i.e., the behavior of the maximum residual as a function of work units, where one work unit is defined as the computational work required to perform one time step on the finest mesh, is presented. The Mach number of the cylinder flow is M1 = 0.05 and the physical time step is set Dtphy = 0.02. It is evident that the main speedup is due to the preconditioning and that the best performance is achieved by the method based on preconditioning and multigrid. Note that without preconditioning the computation reveals a limit cycle behavior. At an inner residual drop of three-orders of magnitude the implicit dual time stepping method is approximately 7 times faster than the purely explicit approach. This result is comparable to the efficiency gain in the channel flow case. When artificial damping terms are added or the physical time step is increased, the efficiency of the multigrid

An efficient solution method to perform large-eddy simulations for nearly incompressible flows using the governing equations of compressible fluids has been developed. The method is based on an implicit time accurate dual time-stepping scheme in conjunction with low Mach number preconditioning and multigrid acceleration. To validate the scheme, large-eddy simulations of turbulent channel flow at Res = 590 and cylinder flow at ReD = 3900 for several Mach numbers and physical time steps have been performed. The results show the scheme to be effective in improving both the accuracy and the rate of convergence at low Mach number flows. Although the efficiency varies at different parameters such as the size of the physical time step, the required residual constraint for the inner pseudo-time iteration, the Mach number and the Reynolds number, the flow problems considered evidence the considerable benefit of the low Mach number preconditioning and the multigrid method. The new method is 6–40 times faster than the basic explicit 5-step Runge–Kutta scheme.

References [1] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2(12):12–26. [2] Turkel E. Preconditioning techniques in computational fluid dynamics. Annu Rev Fluid Mech 1999;31:385–416. [3] Turkel E. A review of preconditioning methods for fluid dynamics. Appl Numer Math 1993;12:257–85. [4] Wall C, Pierce C, Moin P. A semi-implicit method for resolution of acoustic waves in low Mach number flows. J Comput Phys 2002;181:545–63. [5] Meinke M, Schro¨der W, Krause E, Rister T. A comparison of second- and sixth-order methods for large-eddy simulations. Comput Fluids 2002;31:695–718. [6] Meinke M. Numerische Lo¨sung der Navier–Stokes Gleichungen fu¨r instationa¨re Stro¨mungen mit Hilfe der Mehrgittermethode. PhD thesis, Institute of Aerodynamics, RWTH Aachen University, Aachen, Germany, 1993. [7] Weiss JM, Smith WA. Preconditioning applied to variable and constant density flows. AIAA J 1995;33(11):2050–7. [8] Choi YH, Merkle CL. The application of preconditioning in viscous flows. J Comput Phys 1993;105:207–23. [9] Turkel E, Radespiel R, Kroll N. Assessment of preconditioning methods. DLR Report, 1995.

1136

N. Alkishriwi et al. / Computers & Fluids 35 (2006) 1126–1136

[10] Van Leer B, Lee WT, Roe P. Characteristics time stepping or local preconditioning of the Euler equations. AIAA Paper 91-1552, 1991. [11] Jespersen D, Pulliam T, Buning P. Recent enhancement to OVERFLOW. AIAA Paper 97-0644, 1997. [12] Pierce NA, Giles MB. Preconditioned multigrid methods for compressible flow calculations on stretched meshes. J Comput Phys 1997;136:425–45. [13] Unrau D, Zingg DW. Viscous airfoil computations using local preconditioning. AIAA Paper 97-2027, 1997. [14] Dailey LD, Pletcher RH. Evaluation of multigrid acceleration for preconditioned time-accurate Navier–Stokes algorithms. AIAA Paper 95-1668, 1995. [15] Boris JP, Grinstein FF, Oran ES, Kolbe RL. New insights into large eddy simulation. Fluid Dyn Res 1992;10:199–228. [16] Favre A. Equations des gaz turbulents compressibles. J Mec 1965;4:361–90. [17] Knight D, Zhou G, OkongÕo N, Shukla V. Compressible large eddy simulation using unstructured grids. AIAA Paper 98-0535, 1998. [18] Jameson A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. AIAA Paper 91-1596, 1991. [19] Schro¨der W, Keller HB. Wavy Taylor–Vortex flows via multigridcontinuation method. J Comput Phys 1990;91(1):197–227. [20] Venkateswaran S, Merkle CL. Dual time stepping and preconditioning for unsteady computations. AIAA Paper 95-0078, 1995. [21] Arnone A, Lion MS, Povinelli LA. Integration of Navier–Stokes equations using dual time stepping and multigrid method. AIAA J 1995;33(6):985–90. [22] Melson ND, Sanetrik MD, Atkins HL. Time-accurate Navier– Stokes calculations with multigrid acceleration. In: Proc 6th Copper Mountain Conf on Multigrid Methods, 1993. p. 423–39. [23] Poinsot TJ, Lele SK. Boundary conditions for direct simulations of compressible viscous flows. J Comput Phys 1992;101:104–29. [24] Pandya SA, Venkateswaran S, Pulliam TH. Implementation of preconditioning dual-time procedures in OVERFLOW. AIAA Paper 2003-0072, 2003.

[25] Hsu JM, Jameson A. An implicit–explicit hybrid scheme for calculating complex flows. AIAA Paper 2002-0714, 2002. [26] Lessani B, Ramboer J, Lacor Ch. Efficient large-eddy simulation of low Mach number flows using preconditioning and multigrid. Int J Comput Fluid Dyn 2004(3):221–33. [27] Beaudan P, Moin P. Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds numbers. Technical Report TF-62, Center Turb. Res., Stanford, 1994. [28] Lourenco LM, Shih C. Characteristics of the plane turbulent near wake of a circular cylinder. A particle image velocimetry study (Data taken from [27]), 1993. [29] Ong L, Wallace J. The velocity field of the turbulent very near wake of a circular cylinder. Exp Fluids 1996;20:441–53. [30] Norberg C. Effects of Reynolds number and low-intensity free stream turbulence on the flow around a circular cylinder. Publ. No. 87/2, Department of Applied Thermoscience and Fluid Mech., Chalmers University of Technology, 1987. [31] Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel flow up to Res = 590. Phys Fluids 1999;11(4): 445–943. [32] Hirsch C. Numerical computation of internal and external flows, vol. I. New York: Wiley; 1990. [33] Hussain AK, Reynolds WC. Measurement in fully developed turbulent channel flow. Trans ASME 1975:568–80. [34] Brandt A. Guide to multigrid development. Lectures notes in mathematics. Berlin: Springer; 1981. p. 220–312. [35] Wang W, Pletcher RH. On the large eddy simulation of a turbulent channel flow with significant heat transfer. Phys Fluids 1996;8(12):3354–66. [36] Chidambaram N, Pletcher RH. A colocated-grid, fully coupled algorithm for large eddy simulation of incompressible and compressible flows. Numer Heat Transfer Part B 2000;37: 1–23. [37] Cardell GS. Flow past a circular cylinder with a permeable splitter plate. PhD thesis, Graduate Aeronautical Laboratories, California Institute of Technology, 1993. [38] Kravchenko AG, Moin P. Numerical studies of flow over a circular cylinder at ReD 3900. Phys Fluids 2000;12(2):403–17.