Fourth order compact finite volume scheme on nonuniform grids with multi-blocking

Computers & Fluids 56 (2012) 1–16

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Fourth order compact finite volume scheme on nonuniform grids with multi-blocking M. Ghadimi, M. Farshchi ⇑ Aerospace Engineering Department, Sharif University of Technology, Azadi Av., Tehran, Iran

a r t i c l e

i n f o

Article history: Received 12 August 2010 Received in revised form 11 July 2011 Accepted 16 November 2011 Available online 7 December 2011 Keywords: Fourth order Compact scheme Finite volume Nonuniform grids Multi-block domain

a b s t r a c t We have developed a fourth order compact finite volume method for the solution of low Mach number compressible flow equations on arbitrary nonuniform grids. The formulation presented here uses collocated grid that preserves fourth order accuracy on nonuniform meshes. This was achieved by introduction of a new fourth order method for calculation of cell and face averaged metrics. A special treatment of nonlinear terms is used to guarantee the stability of the fourth order compact method. Moreover an approach for applying this method to multi-block domains is presented for complicated geometries and parallel processing applications. Several test cases including the flow in a lid-driven cavity, laminar flow over a flat plate, decay of the Taylor vortex, and transient flow behind a backward facing step show the accuracy, efficiency, and stability of the method in single and multi-block domains. To the best of our knowledge, this is the first symmetric compact method in the finite volume form that can achieve full fourth order accuracy in solving fluid flow equations on nonuniform collocated grids without any stability problems. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Direct numerical simulation (DNS) and large eddy simulation (LES) of turbulent flow fields require high accuracy numerical methods that can resolve small scale motion. High order accurate numerical methods with high frequency resolution have recently gained more popularity because of their application in turbulence simulations and computational aero-acoustics. The so-called compact schemes belong to this class of methods. Compact schemes, for spatial integration, in combination with a fourth order low storage Runge–Kutta time integration method can capture high frequency waves correctly with a few number of grid points. Moreover the computational cost of inverting matrices is low due to their three or five points’ stencil and resulting tri-diagonal or penta-diagonal matrices. Applicability of compact schemes to complex geometries makes these schemes more advantageous than spectral methods with the same or higher level of spectral resolution. Capturing flow fluctuations with high wave-number and small amplitude requires high order of truncation error and low numerical dissipation preserving these small scales. In particular large eddy simulations, in which sub-scale structures are modeled by a second order relation [1], require at least a third order numerical method to distinguish between the truncation errors and small scale fluctuations [2]. Also Since numerical methods with an even ⇑ Corresponding author. Tel.: +98 21 66164605; fax: +98 21 66022731. E-mail address: [email protected] (M. Farshchi). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.11.007

order of accuracy have less dissipation than odd order ones [3], then in order to have truncation errors that are an order of magnitude smaller than the sub-grid scale model errors with low wave amplitude dissipation a numerical method with at least fourth order of accuracy is required. In the present work a fourth order compact finite volume method for solving low Mach number compressible flows on nonuniform grids in multi-block domains is presented. Compact methods have been primarily presented in the finite difference form. However, for more complicated applications and simplicity in implementation of boundary conditions, finite volume formulation is more appropriate. Mattiussi [4] using a topological reasoning proves that integral methods, like finite volume methods, are more suitable for the numerical simulation of field problems. Moreover finite volume methods are conservative, which is an essential property of methods used for compressible fluid flow problems [5]. The paper by Gaitonde and Shang [6] appears to be the first work dealing with compact schemes in a finite volume formulation. They developed a number of fourth order compact finite volume schemes for one dimensional linear wave propagation. Kobayashi [7] formulated a wide range of Pade compact finite volume schemes and considered the related issues like stability, boundary conditions, and order of accuracy. There are a few researchers who have actually utilized compact finite volume methods for the solution of fluid flow problems on nonuniform grids [8–10]. For such applications finite volume approach can be formulated either in physical [8], or in transformed computational space [9,10]. The first approach does not require a specific

2

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

treatment of metrics and has a more insightful formulation; however its application results in an overall second order formulation due to an inherent second order approximation used in metrics determination. Smirnov et al. [11,12] used a fourth order compact finite volume method working in physical space with a collocated grid arrangement resulting in an overall second order method. Results of their simulations show that this method is both stable and accurate. To increase the order of metrics approximation in the physical space compact quadrature must be developed directly on the deformed mesh. Since the truncation error depends on all coordinate directions [12], in order to achieve a given formal order of accuracy, these derivations requires a larger number of neighboring nodes [9]. This drawback could become even more serious in three dimensions. Lacor et al. [8] attempted to overcome these obstacles and use this procedure for compressible flows; however they concluded that an overall fourth order of accuracy could not be achieved by this methodology. Results of their work include a third and a second order compact finite volume method, where the third order method had considerable stability problems. Formulation of a compact finite volume method in the computational space has been used by Pereira et al. [13] for incompressible flow simulations. They used a collocated grid system and reported fourth order accuracy on uniform grids. They also discussed fourth order accuracy on nonuniform grids, but did not present any results for this type of grid. Piller and Stalio [9] presented a successful fourth order accurate compact finite volume formulation on nonuniform staggered grids. To achieve this order of accuracy, a sixth order compact formulation was used. Later Piller and Stalio [10] presented a high order compact finite volume scheme on a curvilinear grid with collocated arrangement for solving a single scalar advection–diffusion equation. They used an asymmetric, fifth order upwind, compact scheme for advective fluxes in their work. Their results show fourth order formal accuracy. However, this level of accuracy was achieved for the solution of a single linear advection–diffusion equation and no results were provided for nonlinear flow equations. The nonlinear terms are a major source of instability in compact schemes and can also reduce the order of accuracy. Also, the application of upwind schemes suffers from addition of an inherent dissipation which is not suitable for LES and DNS applications and interferes with the subgrid scale model [14,15]. In this paper, we use a conservative form of the mass, momentum, and energy equations for low Mach number compressible flows on a nonuniform grid and enforce discrete conservation of these equations in the transformed computational coordinate with an overall fourth order accuracy. To do this we have drawn on the

work by Vasilyev [16], and his conclusion that grid-independent differentiation and interpolation schemes seem inevitable in order to extend the conservation properties to nonuniform grids which impose the use of a coordinate transformation. To achieve an overall fourth order accurate method we have developed a fourth order accurate algorithm for calculating cell and face averaged metrics in the computational space along with a symmetric fourth order compact method for integration of the nonlinear advective–diffusive conservation equations. The nonlinear treatment introduced by Pereira et al. [13] is extended in this work for evaluation of nonlinear advective fluxes. This treatment preserves the order of accuracy without introducing any stability problem. We have selected a collocated grid arrangement because of its simplicity, ease of book keeping, and low computational cost. Moreover, we have developed an algorithm for implementation of this scheme to multi-block domains in complicated geometry and applications with parallel processing. Several test cases; including flow in a lid-driven cavity, laminar flow over a flat plate, decay of a Taylor vortex, and transient flow behind a backward facing step are used to show the accuracy, efficiency, and stability of the method in single and multi-block domains.

2. Governing equations The governing equations for two dimensional compressible flows including continuity, momentum and energy equations can be represented in the following dimensionless form

  @U @F @G 1 @F v @Gv þ þ ¼ þ @t @x @y Re @x @y

ð1Þ

where U is the solution vector, F, G, Fv, and Gv, represent inviscid and viscous flux vectors in x and y direction respectively. These variables can be represented as

2

q

3

7 6 6 qu 7 7 6 U¼6 7 6 qv 7 5 4

ð2Þ

E 2

3

3

2

qv 6 qu2 þ p 7 6 quv 7 7 7 6 6 F¼6 7; G ¼ 6 7 4 quv 5 4 qv 2 þ p 5 qu

ðE þ pÞu

Fig. 1. Physical and computational space for nonuniform grids.

ðE þ pÞv

ð3Þ

3

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

2 6 6 Fv ¼ 6 6 4

0

3

2

0

3

sxx sxy

7 7 7; 7 5

6 6 Gv ¼ 6 6 4

sxy syy

7 7 7 7 5

1 usxx þ v sxy  Prðc1ÞM 2 qx 0

1 usxy þ v syy  Prðc1ÞM 2 qy 0

ð4Þ with P ¼ ðc  1Þ½E  12 qðu2 þ v 2 Þ Here q is density, u and v are velocity components, p is pressure, and E is total energy per unit volume. Ideal gas relations are used to relate pressure to density and total energy. Shear stress tensor components are represented by sxx, sxy, and syy. Heat flux vector components are represented by qx and qy in x and y direction respectively. Stokes stress law and Fourier heat transfer law are used for shear stress tensor and heat transfer vector determination, respectively. The parameters Re, Pr, M0, and c represent Reynolds number, Prandtl number, Mach number, and specific heat ratio respectively.

work a collocated grid arrangement is used due to its simplicity for achieving fourth order formulation and ease of required book keeping with minimal computational costs. We have introduced a special treatment of nonlinear terms to avoid stability problems and achieve robustness. This issue will be discussed in more details later in the paper. 3.2. Transformation to computational space As mentioned above, to achieve fourth order accuracy on a collocated grid, transformation to computational space is made before finite volume integration. The physical and computational space

3. Numerical method A time explicit finite volume approach is used to integrate Eq. (1). Time integration is performed using a fourth order low storage Runge–Kutta method. The low storage formulation is selected because of its low memory requirement. Traditionally in finite-volume formulations, flow equations (Eq. (1)) are integrated over a typical control volume (cell) of a grid with nonuniform meshes; however in order to maintain fourth order accuracy we shall first transform these equations from the physical space to a uniform grid computational space and then apply a finite volume integration method to the transformed equations. Before getting to the space transformation, let us consider physical space grid arrangements. 3.1. Collocated versus staggered grids Specification of solution variables over a computational cell can be done in two different ways. In the classical collocated arrangement all of the solution variables are stored at the same location on the cell. This arrangement requires a single grid and leads to a simple formulation and easy book keeping. An alternative arrangement is to store some variables (usually thermodynamic variables) at the cell center and others (usually kinematic variables) at cell edge centers. This arrangement, called staggered, requires multiple grids and was first introduced to avoid pressure–velocity decoupling in second order numerical methods solving incompressible flow fields. In the case of incompressible solvers with compact schemes, Pereira et al. [13] have shown that pressure and velocity fields are fully coupled in a fourth order compact finite volume scheme on a collocated grid provided that an odd number of cells is used in each coordinate direction. Although this was proven analytically only for flow fields with periodic or Dirichlet boundary conditions on uniform grids; the results obtained for other types of meshing and boundary conditions by Pereira et al. [13] and Smirnov et al. [12] show no spurious oscillations related to pressure–velocity decoupling. Staggered grid arrangement has also been used for compressible flow solvers with the aim of decreasing stability problems and increasing their robustness at high Reynolds numbers [17]. At high Reynolds numbers, aliasing errors associated with the nonlinear interaction of high frequency waves become important and can cause stability problems. This problem is more serious for high order and spectral-like resolution methods such as compact schemes. Use of a staggered grid arrangement is one way of reducing aliasing errors and related stability problem [17]. In the present

Fig. 2. Calculated (a) streamlines and (b) contour of pressure for Taylor vortex flow.

Table 1 Dependence of error norms on grid size in Taylor vortex flow for uniform grids. Grids

77 15  15 31  31

u

v

p

Error (L1)

Order

Error (L1)

Order

Error (L1)

Order

5.87E5 3.21E6 1.63E7

– 3.81 4.11

5.36E5 2.99E6 1.52E7

– 3.78 4.10

3.52E5 2.37E6 1.22E7

– 3.54 4.08

4

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Fig. 3. Effect of grid refinement on error norms in simulating Taylor vortex flow on a uniform mesh.

for a two dimensional domain, are shown in Fig. 1. Eq. (1) can be transformed to the computational space as:

" # ^ ^v ^ @ F^ @ G @U 1 @ F^v @ G þ þ ¼ þ @t @n @ g Re @n @g

ð5Þ

where n and g are computational space coordinates and ^ denotes the modified solution and flux vectors expressed as

^ ¼U U J

ð6Þ

   1 ^ ¼1 g Fþg G F^ ¼ nx F þ ny G ; G x y J J

ð8Þ

where nx , ny , gx and gy are metrics and J is the Jacobian of transformation. In the finite volume approach the cell averaged and face averaged values of variable u are defined as follows [9]:

Z

1 Dn Dg

g 1 ¼ u iþ ;j 2

1 Dg

n

n

Z

iþ1 2

Z

2

gj1

2

uðn; gÞdn dg

ð9Þ

gj1

i1 2

gjþ1

gjþ1

2

c

2

ð13Þ

2

 @ ug @ ug @ ug d  ng ng þ þ c ¼ u  u 1 1 3 iþ1;j i;j @n i2;j @n iþ2;j @n iþ2;j Dn

ð14Þ

Coefficients are computed using fourth (Eq. (13)) and fifth (Eq. (14)) order accurate truncated Taylor series expansion for parameter u. In Eq. (14) fifth order expansion is selected to achieve fourth order accuracy for u derivatives. The coefficients are given as:

a¼

1 ; 4

c¼

1 ; 10

ð7Þ

   1 ^ v ¼ 1 g F v þ g Gv F^v ¼ nx F v þ ny Gv ; G x y J J

 ni;jg ¼ u

  g g 1 þ u  g 1 þ au g 3 ¼ b u  niþ1;j  ni;jg au þ u i ;j iþ ;j iþ ;j

b¼

3 4

d¼

6 5

Convective fluxes at a boundary with a Neumann boundary condition and diffusive fluxes at a boundary with a Dirichlet boundary condition are obtained using following equations [9]:

@ ug  ng  ng 1 þ dui;j þ euiþ1;j @n i2;j

 g 1 þ bu  g 3 ¼ c Dn  g 1 þ au u i ;j iþ ;j iþ ;j 2

a¼

2

4 ; 3

b¼

2

1 1 23 7 ; c¼ ; d¼ ; e¼ 6 6 12 12

2





u niþ12 ; g dg

ð10Þ

 @ ug @ ug @ ug 1 g g g  ni;jg þ eu  niþ1;j  niþ2;j  1 þ du þfu cu 1 þa 1 þb 3 ¼ i2;j @n i2;j @n iþ2;j @n iþ2;j Dn ð16Þ

2

n 1 ¼ u i;jþ 2

1 Dn

Z n

n

iþ1 2

i1 2





u n; gjþ12 dn

ð11Þ

Integrating Eq. (5) over a sample cell and using Gauss’s theorem results in the following equation;

^ ng     @U i;j ^n 1  G ^ n 1 Dn DnDg þ F^giþ1;j  F^gi1;j Dg þ G i;jþ i;j 2 2 2 2 @t     ^v n 1  G ^ v n 1 Dn ¼ F^v g 1  F^v g 1 Dg þ G iþ2;j

ð15Þ

i2;j

i;jþ2

i;j2

a ¼ 52;

b¼

23 35 1033 767 14 ; c¼ ; d¼ ; e¼ ; f ¼ 2 2 12 12 3

The evaluation of coefficients of these equations is done by a fifth order Taylor series expansion chosen to guarantee fourth order accuracy for derivatives. 3.4. Deconvolution

ð12Þ

To advance the above equation in time convective and diffusive fluxes and metrics must be computed in a consistent fourth order formulation. Moreover nonlinear terms must be treated consistently. Appropriate Dirichlet, Neumann or mixed boundary conditions must be formulated at the boundaries. 3.3. Convective and diffusive fluxes The convective and diffusive fluxes are computed by a symmetric compact stencil suggested by Piller and Stalio [9] as follows:

The above control volume approach provides the cell averaged values at each time step. At the end of calculations, cell center val-

Table 2 Dependence of error norms on grid size in Taylor vortex flow for stretched grids. Grids

77 15  15 31  31

dc

0.518 0.237 0.114

u

v

p

Error (L1)

Order

Error (L1)

Order

Error (L1)

Order

4.51E5 2.65E6 1.43E7

– 3.62 3.99

4.72E5 2.76E6 1.45E7

– 3.63 4.02

3.31E5 2.16E6 1.12E7

– 3.49 4.04

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

5

Fig. 4. Effect of grid refinement on error norms in simulating Taylor vortex flow on a nonuniform mesh.

ues must be recovered from cell average ones. This will be done by a deconvolution procedure introduced by Pereira et al. [13] and Piller and Stalio [9]. The procedure used here for interior cells and the left boundary cells with Dirichlet or Neumann boundary conditions follows the formulation presented by Piller and Stalio [9]. In transient applications the initial conditions are given at cell centers. These values must be first converted to cell averaged ones. We have referred to this reverse deconvolution process as a convolution process. To do this a fifth order Taylor series expansion of the variable about the center of a cell is used in each direction. The coefficients are obtained by satisfying the following symmetric relationship: g g af ni1;j þ f ni;jg þ af niþ1;j ¼ bfi1;j þ cfi;j þ bfiþ1;j

the general idea consider the term fg ni;jg which is approximated by a fourth order relation with respect to the f ni;jg g ni;jg and cell or face averaged value of f and g in neighboring cells. This equation can be presented as:  ng ¼ f ng gng fg i;j i;j i;j h     i g g g g g g g g gniþ1;j þ f n;i;jþ1 gni;jþ1 þ a f niþ1;j  f ni1;j  gn;i1;j  f ni;j1  gni;j1

ð20Þ

ð17Þ

where f is the value of parameter at the cell center and indexes correspond to cell number. The coefficients are given as:

a¼

17 ; 206

b¼

27 93 ; c¼ 206 103

This sixth order accurate relation can be used for interior cells only. Similar formulations are required at the boundaries. We have developed the following relations for left Dirichlet and Neumann boundary conditions respectively:

f ng þ af ng þ bf ng þ cf ng ¼ df 1 þ ef i;j i ;j i;j iþ1;j iþ2;j iþ3;j 2

a¼

673 ; 5999

b¼

155 3 852 4608 ; c¼ ; d¼ ; e¼ 5999 857 5999 5999

f ng þ af ng þ bf ng ¼ cDn@f þ dfi;j þ efiþ1;j i;j iþ1;j iþ2;j @ni12;j a¼

155 ; 628

b¼

ð18Þ

ð19Þ

7 27 297 57 ; c¼ ; d¼ ; e¼ 628 628 314 314

where coefficients are found by a fifth order Taylor series expansion. We should point out that Eqs. (18) and (19) break the tri-diagonal structure of the linear system. We have used the Woodbury formulation [17], which breaks this linear system to three tri-diagonal matrices and applies the Thomas algorithm to each one separately. This does not affect the total computational cost, since the convolution procedure is done only once at the beginning of calculations. 3.5. Nonlinear terms To achieve a full fourth order accurate and stable scheme special treatment for nonlinear terms is necessary. Present treatment is an extension of the method used by Piller and Stalio [9]. To show

Fig. 5. Calculated (a) horizontal and (b) vertical velocity components over a flat plate.

6

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

where the coefficient is obtained by applying a fourth order Taylor series expansion of parameters f and g. In this case it results in 1 a ¼ 48 . To maintain fourth order accuracy at the boundaries with Dirichlet and Neumann boundary conditions we have developed the following equations respectively:

fg ni;jg ¼ f ni;jg gni;jg     i 1 hg f iþ1;j  f gi1;j ggiþ1;j  ggi1;j þ f ni;jþ1  f ni;j1 gni;jþ1  gni;j1 þ 12 2 2 2 2 2 2 2 2 ð21Þ

! @f g @f g @g g @g g þ þ Dn 2 1 1 1 1 @niþ2;j @ni2;j @niþ2;j @ni2;j ! # @f n @f n @g n @g n 2 þ þ D g þ 1 1 1 1 @ gi;jþ2 @ gi;j2 @ gi;jþ2 @ gi;j2

1 fg ni;jg ¼ f ni;jg gni;jg þ 48

"

aliasing errors appearing in high resolution schemes and provides improved stability as will be shown later in the paper.

3.6. Metrics Metrics and Jacobian determination must be consistently done in a manner that the total order of accuracy is preserved. We have developed a five-step method for calculation of face averaged metrics by using compact interpolations. Given the coordinates of cell vertices in the physical space, these steps are taken: 1. Calculate physical coordinates of midpoints of computational faces and cells by the following fourth order relations:

r i;j1 ¼

ð22Þ

2

We have also treated face averaged nonlinear terms which appear in Eq. (12) in a similar manner resulting in formulations given below:

fg giþ1;j ¼ f giþ1;j ggiþ1;j þ 2

f

f

2

2

 i 1 hng g g g f i;jþ1  f ni;j1  gni;j1 gni;jþ1 48

ð23Þ

@g g g @g g ¼ f iþ1;j iþ1;j @niþ12;j 2 @n 2  i 1 hng g g g g g f i;jþ1  f ni;j1 þ  gniþ1;j1  gni1;jþ1 þ gni1;j1 gniþ1;jþ1 96Dn ð24Þ @g g g @g g 1 ¼ f iþ1;j iþ1;j @ giþ2;j 2 @g 2  i 1 hng ng g g f i;jþ1  f ni;j1 þ  2gni;jg þ barg i;j1 gni;jþ1 24Dg

 1   9  r i1;j1 þ r iþ1;j1  ri3;j1 þ r iþ3;j1 2 2 2 2 2 2 2 2 16 16

 1   9  r i1;j þ r iþ1;j  r i3;j þ r iþ3;j 2 2 2 2 16 16

r i;j ¼

ð26Þ

ð27Þ

where r can be any of the physical coordinates, Fig 1. 2. Evaluate cell and face averaged values of physical coordinates by integrating these coordinates with a classical fourth order Runge–Kutta method and using Eqs. (9)–(11).

r n 1 ¼ i;j 2

r ni;jg ¼

 1 ri1;j1 þ 4r i;j1 þ riþ1;j1 2 2 2 2 2 6

ð28Þ

1  r 1 1 þ r iþ1;j1 þ 4ri;j1 þ 4r i1;j þ 16ri;j þ 4r iþ1;j 2 2 2 2 2 36 i2;j2 

þ4r i;jþ1 þ ri1;jþ1 þ r iþ1;jþ1 2

2

2

2

2

ð29Þ

ð25Þ 3. Calculate cell and face averaged values of Jacobian inverse

Appendix A contains additional formulations for triple products and products in three-dimensional space. This treatment reduces

(ð1=JÞni;jg ; ð1=JÞgi1;j ; . . .). The inverse of the transformation Jacobian 2

is defined as;

Table 3 Dependence of error norms on grid size in laminar flow over a flat plate for stretched grids. Grids

40  30 80  60 160  80

dc

1.69E2 8.47E3 4.99E3

u

v

p

Error (L1)

Order

Error (L1)

Order

Error (L1)

Order

2.34E4 1.68E5 2.04E6

– 3.81 3.98

3.41E5 2.48E6 3.07E7

– 3.79 3.95

2.14E3 1.60E4 2.00E5

– 3.75 3.93

Fig. 6. Effect of grid refinement on error norms in simulating a flow over a flat plate.

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

1 ¼ xn yg  xg yn J

ð30Þ

The integral of the Jacobian inverse over a cell is equal to the volume of that cell and can be evaluated in terms of cell vertices coordinates. Hence the cell averaged value of Jacobian inverse can be obtained from Eq. (9). The face averaged Jacobian inverse requires integration over a face. Two types of derivatives are involved; derivatives with respect to the face-parallel coordinate and those with respect to the face-perpendicular coordinate. The face-parallel integration can be performed explicitly by knowing the coordinates of face vertices; the face-perpendicular integration is calculated by a fourth order relation as follows;

@r 11 85 23 g 2 g Dn gi1;j ¼  rgi1;j þ r ni;jg  r niþ1;j þ rniþ2;j @n 2 3 18 18 9 2

7

Stokes equations [1,20]. Staggered grids can also lead to reduction of aliasing errors and improvement of robustness [18], however it adds formulation complexity and requires extensive book-keeping. We have used collocated grids to avoid additional computational cost and will show that the nonlinear treatments presented in previous sections, decrease aliasing error and make the scheme stable.

ð31Þ

where the cell and face averaged values of physical coordinates on the RHS of this equation are known from the previous step. g n;g n;g n;g  xi;j ; g  yi;j ) 4. Calculate cell averaged values of metrics ( nn; xi;j ; nyi;j ; g using compatibility equations: 8

< ny ¼ xg J :

gy J

¼ xn

(

;

nx J

gx J

¼ yg

ð32Þ

¼ yn

then use the nonlinear treatment procedure explained in the n;g previous section. For example nx can be computed by the nonlini;j ear relation:

h   ng 1 ng ng g g 1 nx i;j nx ð1J Þni;jg ¼ nx i;j ð J Þi;j þ 12  nx ni1;j ð1J Þni;jg  ð1J Þni1;j   i ng ng g  nx i;j1 ð1J Þni;jg  ð1J Þni;j1 þ nx i;j

ð33Þ

The LHS of this equation can be obtained by the cell averaged compatibility equation follows:

nx

  1 ng 1 ng DyðgÞi1;j þ 4DyðgÞi;j þ DyðgÞiþ1;j i;j ¼ yg i;j ¼ 2 2 J 6

ð34Þ

where a classical fourth order Runge–Kutta method is used in this equation and:

DyðgÞi;j ¼ yi;jþ1  yi;j1 2

ð35Þ

2

Midpoint values of coordinates in computational space are obtained from step 1. An Equation with face averaged values of Jacobian inverse is used instead of Eq. (33) for cells near boundaries.  gx 1 ; 5. Calculate the face averaged value of metrics,  ngx 1;j ;  ngy 1;j ; g i ;j i

2

i

2

2

g gyi1;j ; . . . ; using the compact scheme (Eq. (13)) and cell averaged 2

value of metrics computed in the previous step. 3.7. Stability The main drawback of compact schemes is their sensitivity to aliasing errors which affects their stability [18]. The aliasing error, which arises whenever nonlinear terms are approximated in the discrete physical space as a point-wise product [19], becomes important in schemes having high resolution at high wave numbers [20]. This type of error causes numerical instability when flow fields with high wave number waves, such as high Reynolds number flows, are simulated using high resolution methods such as compact schemes. Several approaches have been used to eliminate aliasing errors in compact schemes. A popular approach is the low-pass filtering method. In this method high frequency waves are removed by explicit low-pass filters [1]. Another approach, which reduces the aliasing error, is to utilize a modified form of nonlinear terms, i.e. the skewsymmetric form of the nonlinear convective term in the Navier–

Fig. 7. Calculated streamline of lid-driven cavity (a) Re = 100 (b) Re = 400 (c) Re = 1000.

8

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Fig. 8. Calculated and benchmark [23–25] velocity components in the lid-driven cavity; left figures: u in the vertical middle line, right figures v in the horizontal middle line. (a) Re = 100, (b) Re = 400 and (c) Re = 1000.

9

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Hence, the nonlinear treatment increases the robustness and stability of our scheme. 3.8. Multi-blocking One of the important contributions of this work is to extend our fourth order compact finite volume scheme to multi-block domains. This can help with the simulation of flow fields in complicated geometries or parallel processing applications. Our approach requires that the neighboring cells of different blocks have equally shared faces. This allows for fluxes on shared faces to be calculated using a recursive method. For a given block the variables on the shared faces are used as the boundary conditions allowing the governing equations to be solved explicitly for that block. Updating fluxes on shared faces using the compact relations (Eqs. (13) and (14)), this procedure is used for the next block. This recursive method is continued until an acceptable convergence is achieved. In the next section the solution of a transient flow over a backward facing step is used to show the stability and efficiency of the proposed approach. In this case convergence is achieved within 2 or 3 iterations. The presented approach is different from Chao et al. [27] method in two ways. First, we have used our compact scheme for computing the fluxes on shared faces whereas Chao et al. [27] have used an explicit formulation on shared points which is different from compact formulation of interior points. Hence the present approach leads to a consistent compact formulation of all fluxes at each time step. Second, the present approach is an iterative method at each time step making it suitable for time dependent problems.

Fig. 9. A stretched 60  06 grid used in simulation of lid-driven cavity.

To show this we consider a two dimensional discrete Fourier transform for a field variable f as: N

M

1 1 2 2 X X

fi;j ¼

n¼N2

^f nm e2pic ðinNþjm MÞ

ð36Þ

m¼M 2

If we integrate this equation over a cell according to Eq. (9), the cell averaged value is obtained as (see Appendix B): N

f ng ¼ i;j

M

1 1 2 2 X X n¼N2

Þ ^f e2pic ðinNþjm M nm

4. Results

ð37Þ

m¼M 2

In this section several two dimensional test cases have been simulated to ascertain the formal order of accuracy, stability and robustness of proposed method. Configurations have been selected in such a way that various boundary conditions and types of flow fields could be tested by the presented method. As the first test case, Taylor vortex flow is solved to check the fourth order formal accuracy of the method for both uniform and nonuniform grids. The exact solution of this transient flow with periodic boundary condition is available. For the second test case a laminar flow over a flat plate is simulated to evaluate the order of accuracy in a wallbounded flow with inflow and outflow boundary conditions. The laminar flow Reynolds number is set to 10,000 to assess the stability and robustness of the method in the presence of high frequency numerical noise. The third test case is the classical lid-driven cavity flow. In this test case the Richardson extrapolation method is used to obtain a nominal exact solution for the evaluation of the formal

where the cell averaged coefficients are related to the cell centered coefficients by the following relation:

pn  p m  ^f nm ¼ sin  N Dn sin M Dg ^f mn pn Dn pm Dg N M

ð38Þ

Here the term sinðkn Þ sinðkg Þ=ðkn kg Þ (kn and kg are wave numbers in the directions of computational coordinates) acts as a box-filter [21] and reduces the effect of high wave numbers. Recall that a nonlinear term such as fg ni;jg is replaced by the nonlinear treatment given by Eqs. (20)–(22). This treatment of the average of nonlinear terms, instead of using f ni;jg g ni;jg , reduces the effect of high wave numbers and subsequent aliasing errors. From a mathematical point of view, fg ni;jg in the Fourier space can be expressed as; 0 fg ni;jg

1

B X X X B B ¼ B B n¼N2 m¼M @ þ n n 2 1 2 ¼ n N 1 2

M 1 2

m1 þ m2 ¼ m

^f n m g^n m þ 1 1 2 2

X n1 þ n2 ¼ n þ N m1 þ m2 ¼ m

X

^f n m g^n m þ 1 1 2 2

^f n m g^n m þ 1 1 2 2

n1 þ n2 ¼ n m1 þ m2 ¼ m þ M

X n1 þ n2 ¼ n þ N m1 þ m2 ¼ m þ M

C C sin pn Dn sin pm Dg in jm ^f n m ^f n m C pNn pm M  e2pic ð N þ M Þ 1 1 2 2C C Dn M Dg N A ð39Þ

in which, all terms in the parentheses except the first one are the aliasing errors. Clearly, the contribution of aliasing errors is larger at high wave numbers. The box-filtering term in the above equation is the same as that appearing in Eq. (38) and plays the role of filtering and reducing high wave number effects. We must point out that this filtering effect is inherent in our cell averaging treatment of nonlinear terms and no additional explicit filtering has been used.

order of accuracy. The flow is then solved on 2, 4, 8 and 16 multiblock domains to evaluate the computational efficiency of the proposed multi-blocking method. The next test case considered is concerned with the initial flow development over a backward-facing step and dynamics of the recirculation region behind the step. This is to verify the capability of the proposed multi-blocking method in simulating transient characteristics of the flow.

10

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Table 4 Dependence of error norms on grid size in lid-driven cavity flow for uniform grids. Grids

25  25 50  50 100  100

u

v

p

Error (L1)

Order

Error (L1)

Order

Error (L1)

Order

7.87E3 9.10E4 5.81E5

– 3.11 3.97

6.97E3 8.63E4 5.46E5

– 3.01 3.98

5.14E03 6.16E04 3.75E05

– 3.06 4.04

4.1. Taylor vortex To evaluate the order of accuracy of the proposed scheme in an unsteady flow, a decaying Taylor vortex flow with an exact analytical solution is considered. The flow field variables are given as 2t

uðx; y; tÞ ¼  sinðxÞ cosðyÞ eRe

v ðx; y; tÞ ¼ cosðxÞ sinðyÞ Mach numbers of all above test cases are so low that they can be classified as incompressible flow fields. This choice was driven by the availability of benchmark test cases. The final test case considered is the subsonic compressible flow in a twosided lid-driven cavity. This final test case validates applicability of the present method to low Mach number compressible flow fields.

2t

eRe

ð40Þ

1 4t pðx; y; tÞ ¼ ðcosð2xÞ þ cosð2yÞÞ eRe 4 in the ½0; p2 domain. These equations are nondimensionalized and Re represents the Reynolds number of flow. The structure of the velocity field and the contours of pressure are shown in Fig. 2. The shape of streamlines is constant and only their values are changed as time evolves. We choose the Re = 100 for our study and results are considered at t ¼ 0:34657Re corresponding to the

Table 5 Dependence of error norms on grid size in lid-driven cavity flow for stretched grids. Grids

dc

u

30  30 60  60 120  120

5.64E2 1.88E2 9.41E3

v

p

Error (L1)

Order

Error (L1)

Order

Error (L1)

Order

8.70E3 6.38E4 3.87E5

– 3.76 4.05

7.67E3 5.27E4 3.32E5

– 3.86 3.99

2.14E2 1.45E3 9.65E5

– 3.87 3.92

Fig. 10. Effect of grid refinement on error norms in a lid-driven cavity flow on a uniform mesh (Re = 100).

Fig. 11. Effect of grid refinement on error norms in a lid-driven cavity flow on a nonuniform mesh (Re = 100).

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Fig. 12. A curvilinear 50  50 grid used in simulation of lid-driven cavity.

11

Fig. 15. A multi-block domain with 16 blocks for the simulation of lid-driven cavity flow.

Fig. 13. Calculated with curvilinear mesh and benchmark [23–25] results of velocity components in the lid-driven cavity; left figures: u in the vertical middle line, right figures v in the horizontal middle line (Re = 1000).

Fig. 14. Effect of grid refinement on error norms in a lid-driven cavity flow on a curvilinear mesh (Re = 1000).

12

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Fig. 16. Convergence History of a lid driven cavity flow (Re = 100) for 1, 2, 4, 8 and 16 blocks.

Fig. 18. Streamlines of the flow over a backward-facing step with Re = 304 at Ut ¼ 3:86 (a) present work (b) experimental result [28]. h

Table 6 Averaged iteration number in each step for 2, 4, 8 and 16 multi-block domains in lid-driven cavity flow. Number of blocks

Averaged iteration number

2 4 8 16

2.93 3.11 3.23 3.36

half-life of the velocity field. Table 1 lists the error norms of the velocity components and pressure obtained on three uniform grids with respect to the exact solution. The variations of velocity and pressure errors with respect to the grid size are shown on a logarithmic scale in Fig. 3. It is observed that the scheme has reached the fourth order of accuracy for both velocity and pressure. To create a nonuniform grid the square domain is divided into four equal squares. In each of the squares the following grid stretching is used at the squares outer corner,

8 1g 91 > > < bbx þ1 1= x 1 B C x ¼ L@1  bx  1g A > > : bbx þ1 ; þ1 1 x 8 1g 91 0 > > by þ1 < b 1 1= y B C y ¼ H@1  by  1g A > > : bby þ1 ; þ1 1 0

Fig. 19. Streamlines of the flow over a backward-facing step with Re = 304 at Ut ¼ 10:5 (a) present work (b) experimental result [28]. h

ð41Þ

y

1

y

0.75 0.5 0.25

Zone2

Zone1

0 0

0.5

1

x

1.5

2

Fig. 17. Multiblock domain and meshing used for the simulation of flow over a backward-facing step.

Fig. 20. Streamlines of the flow over a backward-facing step with Re = 304 at Ut ¼ 23:8 (a) present work (b) experimental result [28]. h

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

where symbols L and H represent the width and the height of each sub-domain. This formulation stretches grid near the boundaries by means of parameters bx and by : In this test case, these stretching parameters are set to 1.2, see Fig. 9. For such nonuniform grids, we define the characteristic grid size as:

dc ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 4 4 1 d n i i

ð42Þ

where di represents different grid sizes. The resulting error norms versus grid sizes are presented in Table 2. Fig. 4 shows this data on a logarithmic scale and compares it with a fourth order line. It verifies the fourth order accuracy of the solution for the nonuniform grid. 4.2. Flat plate The second test case considered is a laminar flow over a flat plate. This is an appropriate test case for evaluating the stability and accuracy of the method in a convective–diffusive flow field with various boundary conditions such as inflow-outflow and far-field boundary conditions. To simulate an incompressible flow field the Mach number is set to 0.01. At the inflow boundary the density and velocity components are set to their inlet values and pressure is extrapolated from the inside flow using a Neumann boundary condition. At the outflow and far-field boundaries the pressure is fixed and the density and velocity components are extrapolated from the inside flow using a Neumann boundary condition. Fig. 5 compares the profiles of vertical and horizontal components of the velocity at Re = 10,000 with the Blasius solution [22]. Here a 40  30 stretched grid with stretching away from the wall according to Eq. (41) is used. The stretching parameter, by , is set to 1.1. The results show the accuracy of the method in simulating high Reynolds number laminar flows and confirm the effectiveness of the proposed treatment of nonlinear terms in providing a stable solution. The formal order of accuracy is evaluated using the exact Blasius solution. Three stretched grids are used for this purpose. The stretching parameter, by for the 40  30 and 80  60 grids is set to 1.1, and in the 160  80 grid it is set to 1.003. To demonstrate the ability of the fourth order flow solver in handling grids with high aspect ratio cells, the grid stretching is performed in such a way to have cells with aspect ratio of 100 for the finest grid. The resulting error norms versus characteristic grid sizes are presented in Table 3 and are compared with the fourth order slope in Fig. 6. It is shown that the fourth order of accuracy is achieved for this case too. 4.3. Lid driven cavity The next test case considered is a lid driven cavity flow. This is a classical test case to validate the accuracy of the solution method for recirculating flow fields. This test case has been chosen for its simple wall boundary conditions, and having classical benchmark solutions [23,24]. Fig. 7 shows the simulation streamlines for three different Reynolds number. For these simulations, a uniform 50  50 grid is used and Mach number is set to 0.01. Small vortices are formed on the lower right and left corners of the cavity and evolve as the Reynolds number increases. These results are in accordance with the benchmark results [23,24]. Fig. 8 presents the resulting vertical and horizontal velocity components at the centerline of the cavity in x and y directions, respectively. An excellent agreement is observed with the benchmark results of Ghia et al. [23], Bottella and Peyret [24] and Erturk et al. results [25]. Ghia et al. applied 129  129 grid points and Bottella and Peyret required at least 96  96 grid points to obtain accurate results at

13

Re = 1000. We have been able to obtain as accurate a result using 50  50 grid points, indicating an important advantage of the full fourth order compact scheme used here. The formal order of accuracy of this test case can be evaluated using Richardson extrapolation [26]. This method uses numerical solutions at a given point in the domain obtain from two solutions at two levels of grid resolution. A nominal exact solution is then extrapolated from these solutions as follows:

E¼

mk N2  N1 mk  1

ð43Þ

where the symbols E, N1 and N2 represent the nominal exact solution, the numerical solution on the coarse grid and the numerical solution on the fine grid respectively. m is the ratio of the coarse grid size to the fine one and k is the order of accuracy in the numerical simulation. To use this method we have chosen the numerical solution on a 128  128 uniform grid as the reference fine grid for both uniform and nonuniform coarse grids. Since generally speaking a cell mid point in the coarse grid does not have a geometrically corresponding node in the fine grid; an interpolation scheme is required to find fine grid solutions at the computational nodes of the coarse grid. The interpolation scheme must have an order of accuracy higher than the order of simulation scheme. A sixth order interpolation scheme, using the nearest twenty points to evaluate the two dimensional Taylor series coefficients up to the fifth term, is used in the present work. The nonuniform grids used in this study are stretched according to Eq. (41) with stretching parameters ta-

Fig. 21. Time variation of the reattachment length for Re = 304.

Fig. 22. Calculated streamline of two-sided lid-driven cavity with Ma = 0.4 and Re = 700.

14

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

ken to be 1.2 in both directions. Fig. 9 shows a sample 60  60 nonuniform stretched grid used for the above case. The grid sizes and error norms of velocity components and pressure for the Reynolds number of 100 on uniform and nonuniform grids are shown in Tables 4 and 5 respectively. The formal order of accuracy for each grid is computed and presented in these tables. Results indicate that the fourth order of accuracy is reached on both uniform and nonuniform grids. Figs. 10 and 11 show the variation of these error norms with respect to the characteristic grid sizes in the logarithmic scale for uniform and nonuniform grids respectively. We have also demonstrated the solver competence on a curvilinear grid for this type of flow. Fig. 12 shows a fully curvilinear 50  50 grid. This grid has been generated by the following relations: Þðiimid Þ i;j Þ xi;j ¼ xi;j þ 0:02 ðiieimid sinð2py 2 mid Þ i;j þ 0:02 ðjje Þðjj sinð16pxi;j Þ yi;j ¼ y j2

ð44Þ

mid

i;j denote coordinates of the corresponding In which  xi;j and y uniform grid. The numbers imid and ie are the half number of grids in x direction and the edge number defined as follows:

ie ¼



1

ı 6 imid

nx

i > imid

ð45Þ

The same definition is used for jmid and je in the y direction. Fig. 13 shows vertical and horizontal velocity components at the centerline of the cavity in x and y directions compared with the benchmark results at Re = 1000. The variation of error norms with respect to the characteristic grid sizes in the logarithmic scale is also shown in Fig. 14; confirming fourth order of accuracy on this type of grid. 4.4. Multi-block lid-driven cavity The lid driven cavity flow is used to evaluate the additional computational cost associated with a multi-block simulation in compare to a single block simulation. The multi-blocking technique proposed here does not require any overlapping grids at the block interfaces, therefore there is no additional cost associated with the flux calculations at a block interfaces. However, there are a number of iterations in a given time step advancement to converge the value of the interface fluxes to a certain accuracy. The number of iterations increases if more convergence is imposed. It can be shown that the number of time steps to reach a stationary

Fig. 23. Calculated and benchmark [30] horizontal velocity component (left figures) and temperature (right figures) in two-sided lid-driven cavity with Ma = 0.4; (a) profiles along the horizontal centerline y = 0.5 (b) profiles along the vertical centerline x = 0.978.

15

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

solution can be increased to compensate for a low convergence requirement of the interface fluxes. The blocks are arranged in a serial manner here. Fig. 15 shows a 128  128 stretched grid and blocks arrangement for solving the flow in a cavity with 16 blocks. Fig. 16 compares the convergence history of the lid-driven cavity flow for a single block domain versus 2, 4, 8, and 16 multi-block domains for the case of Re = 100. Here the local interface flux convergence is reached when the percentage of error in interface fluxes reduces to about 5% of the local residual. Table 6 presents the average iteration numbers in each time step for the 2, 4, 8 and 16 blocks. These iteration numbers are calculated for a 128  128 grid. Adopting the Chao et al. [27] definition of computational efficiency of a multi-block approach for parallel processing, based on results presented in Table 6, we obtain an estimated computational efficiency of 79% for the 16 blocks, which is comparable to the their computational efficiency. 4.5. Laminar flow behind a backward-facing step The next test case is the flow behind a backward-facing step with the aim of testing applicability of the scheme in solving unsteady flows in multi-block domains. In this type of flow, fluid flow is separated from the wall due to suddenly expanded geometry and the resulted shear layer is reattached to the wall somewhere downstream of the step and develops a recirculating zone. The characteristic length of reattachment zone in the direction of flow and the number of vortices in this zone for unsteady flow are functions of the Reynolds number and time of evolution [28]. The multi-block grid used for the simulation of this flow is shown in Fig. 17. As observed in this figure, a relatively coarse grid (a 10  10 grid before and a 100  20 grid after step) is applied in this simulation. A symmetric boundary condition has been used for the upper boundary in which the vertical component of velocity is set to zero and the horizontal component of velocity, density and pressure are extrapolated from the inside flow using a Neumann boundary condition. The Reynolds number is set equal to 304 for comparison with the experimental results. Convergence of multi-block approach is achieved in each time step by only two or three iteration. The flow streamlines at three different nondimensionalized time intervals are compared with experimental results of Honji [28] Figs. 18–20. Time is nondimensionalized using the mean velocity of flow at the inlet, U, and the step height, h. These figures show an excellent accordance of vortices’ shape and evolution with the experimental results. The time variation of the reattachment location is compared with experimental data of Honji [28] and numerical results of Gerrard [29] in Fig. 21. The reattachment length xr is measured as the distance of reattachment location to the step. It can be seen that results of the present work have a good accordance with the numerical results of Gerrard [29]. However, there is a discrepancy between numerical and experimental results at later times which can be due to three-dimensionality, developing in the experiments at later times. 4.6. Compressible flow in a two-sided lid-driven cavity All of the test cases considered so far are very low Mach number flow cases. To show the presented method’s ability of simulating compressible flow fields we have selected a compressible flow in a two-sided lid-driven cavity [30]. The flow is generated by moving two opposing walls of a rectangular cavity in opposite directions. Its Prandtl number is 0.72 and its Mach number is 0.4. This was the highest Mach number used in the benchmark test case [30]. The ratio of the horizontal side of the cavity to its vertical side is 1.955. An 80  40 grid which is stretched near the walls according to Eq. (41) has been used in this test case. Stretching parameters in both directions are set to the 1.2. Fig. 22 shows streamlines for the

Reynolds number of 700 which is in accordance with the benchmark result [30]. Fig. 23 presents variations of horizontal velocity component and temperature along the horizontal and vertical centerline of the cavity. The present results have an excellent agreement with the benchmark results of Shah et al. [30]. 5. Conclusion We have presented a fully fourth order symmetric (central) compact finite volume method on nonuniform grids which can be used with multi-block domains. To achieve this we performed a coordinate transformation from the physical space to the computational space and used a new fourth order method for calculating cell and face averaged metrics. A new nonlinear treatment was used to preserve the accuracy and stability of the solution method. We have also extended our method to simulation of flows in multiblock domains making it appropriate for complicated geometries and parallel processing. Computational results for several test cases including; a decaying Taylor vortex, a laminar flow on a flat plate, and a lid-driven cavity flow, illustrate the fourth order formal accuracy and stability of this method for solving problems with different boundary condition. Moreover, the solution of liddriven cavity flow on multi-block domains and flow behind a backward facing step show the efficiency of presented approach for multi-block domains. The applicability of the present method to low Mach number compressible flows was demonstrated by considering a two-sided lid-driven cavity flow. Appendix A. Treatment of nonlinear terms The forth order formulation of nonlinear terms described in Section 3.5, can be extended to triple products as follows:

  1 h g  g g g g g g g g fghi1;j ¼ f i1;j g i1;j hi1;j þ hi1;j f i1;jþ1  f i1;j1 g i1;jþ1  g i1;j1 48 2 2 2 2 2 2 2 2 2    g gi1;j f gi1;jþ1  f gi1;j1 hgi1;jþ1  hgi1;j1 2 2 2  2  2 i þf gi1;j g gi1;jþ1  g gi1;j1 hgi1;jþ1  hgi1;j1 ðA:1Þ 2

2

2

2

2

This formulation can be derived explicitly or alternatively by twice using Eq. (23) for the term ðfgÞhgi1;j and then for the term 2 ðfgÞgi1;j . 2

Treatment of nonlinear terms in three-dimensions is similar to their two-dimensional counterpart, for example, Eq. (20) can be presented as follows in three-dimension.

  1 h ngf gf gf ngf gf gf gf fg ni;j;k ¼ f ni;j;k g i;j;k þ f iþ1;j;k  f ni1;j;k g niþ1;j;k  g ni1;j;k 48    gf gf gf gf þ f ni;jþ1;k  f ni;j1;k g ni;jþ1;k  g ni;j1;k   i gf gf gf gf  f ni;j;k1 g ni;j;kþ1  g ni;j;k1 þ f ni;j;kþ1

ðA:2Þ

Appendix B. DFT of a cell-averaged variable The Discrete Fourier transform of function f is defined by Eq. (36) of the paper: N

M

1 1 2 2 X X

fi;j ¼

^f e2pic ðinNþjm MÞ nm

ðB:1Þ

n¼N2 m¼M 2

On the other hand the cell averaged values of f is defined by Eq. (9) of the paper:

f ni;jg ¼

1 Dn Dg

Z n

n

iþ1 2

i1 2

Z

gjþ1 2

gj1 2

f ðn; gÞdn dg

ðB:2Þ

16

M. Ghadimi, M. Farshchi / Computers & Fluids 56 (2012) 1–16

Note that the computational space grid is uniform and the grid size is equal to unity. Hence at point (i, j) in the computational space grid (n, g)ij = (iDn, jDg), where Dn = 1 and Dg = 1, therefore (n, g)ij = (i, j). Transporting the origin of computational coordinates to the cell center, this integral can be stated as follows.

Z

1 DnDg

f ni;jg ¼

Z

Dn 2

Dg 2 Dg

D2n

f ði þ n; j þ gÞdn dg

ðB:3Þ

2

Substituting Eq. (B.1) in this relation, one can obtain:

f ni;jg

Z

1 ¼ DnDg

Z

Dn 2

Dg 2 Dg

D2n

2

0N 1 1 M 1 2 2 X X ðiþnÞn ðjþgÞm 2 p i þ ð Þ ^ c @ Adn dg N M f nm e

ðB:4Þ

n¼N2 m¼M 2

In this equation integrations can commute with summations and double integrations can be performed independently due to separability of the integrand. So one can obtain:

f ni;jg

! N 1 M1 Z Dn Z Dg 2 2 X 2 2 ðiþnÞn ðjþgÞm 1 X 2pic ð N þ M Þ ^ f nm ¼ e dndg Dg Dn Dg D2n  2 N M n¼ 2 m¼ 2

¼

! ! Z Dg N 1 M1 Z Dn 2 2 X 2 2 ðiþnÞn ðjþgÞm 1 X ^f nm e2pic ð N Þ dn e2pic ð M Þ dg Dg Dn Dg D2n 2 N M n¼ 2 m¼ 2

N1 2

M 1 2

X X

¼

1 Dn

^f nm e2pic ðinNþjm MÞ

n¼N2 m¼M 2

Z

Dn 2

! 2pic ðnn Þ N

e

1 Dg

dn

D2n

Z

Dg 2 Dg

gm

2pic ð M Þ

e

! dg

2

ðB:5Þ The integral terms can be computed as follows:

R D2n

nn

R Dg

gm

1 e2pic ð N Þ dn ¼ 2pi 1n Dn Dn Dn c ðN Þ 2

&

1 2 e2pic ð M Þ d Dg Dg 2



 Dnn Dnn epic ð N Þ  epic ð N Þ ¼



Dgm

Dgm

1

ðpNnÞDn



  sin pNn Dn

g ¼ 2pic ð1mÞDg epic ð M Þ  epic ð M Þ ¼ ðpm1ÞDg sin M M

p m M

 Dg

ðB:6Þ Substituting these relations in the Eq. (B.5), one can obtain: N

f ni;jg

¼

pn  Dn sinðpMm DgÞ 2pic ðinþjmÞ ^f nm sin N M pn N  pm  e Dn Dg N M m¼M M

1 1 2 2 X X n¼N2

ðB:7Þ

2

References [1] Meinke M, Schroder W, Krause E, Rister TH. A comparison of second- and sixthorder methods for large-eddy simulations. Comput Fluids 2002;31:695. [2] Ghosal S. An analysis of numerical errors in large eddy simulation of turbulence. J Comp Phys 1996;125:187.

[3] Hoffmann KA, Chiang ST. Computational fluid dynamics for engineers. Wichita (KS): Engineering Education System; 1993. [4] Mattiussi B. An analysis of finite volume, finite element, and finite difference methods using some concepts from algebraic topology. J Comp Phys 1997;133(2):289. [5] Hirsch C. Numerical computation of internal and external flows. West Sussex: Wiley; 1994. [6] Gaitonde D, Shang JS. Optimized compact-difference-based finite-volume schemes for linear wave phenomena. J Comp Phys 1997;138:617. [7] Kobayashi MH. On a class of Pade finite volume methods. J Comp Phys 1999;156:137. [8] Lacor C, Smirnov S, Baelmans M. A finite volume formulation of compact central schemes on arbitrary structured grids. J Comp Phys 2004;198:535. [9] Piller M, Stalio E. Finite-volume compact schemes on staggered grids. J Comp Phys 2004;197:299. [10] Piller M, Stalio E. Compact finite volume schemes on boundary fitted grids. J Comp Phys 2008;227:4736. [11] Smirnov S, Lacor C, Lessani B, Meyers J, Baelmans M. A finite volume formulation for compact schemes on arbitrary meshes with applications to RANS and LES. European congress on computational methods in applied sciences and engineering, September 11–14, Barcelona; 2000. [12] Smirnov S, Lacor C, Baelmans M. A finite volume formulation for compact scheme with applications to LES. In: AIAA Paper 2001-2546, 15th AIAA computational fluid dynamics conference, June 11–14, Anaheim, CA; 2001. [13] Pereira JMC, Kobayashi MH, Pereira JCF. A fourth-order-accurate finite volume compact method for the incompressible Navier–Stokes solutions. J Comp Phys 2001;167:217. [14] Beaudan P, Moin P. Numerical experiments on the flow past a circular cylinder at sub-critical Reynolds numbers. Technical report TF-62, Center of Turbulence Research; 1994. [15] Boersma BJ. A staggered compact finite difference formulation for the compressible Navier–Stokes equations. J Comp Phys 2005;208:675. [16] Vasilyev OV. High-order finite difference schemes on nonuniform meshes with good conservation properties. J Comp Phys 2000;157:746. [17] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the art of scientific computing. New York: Cambridge University Press; 2007. [18] Nagarajan S, Lele SK, Ferziger JH. A robust high-order compact method for large eddy simulation. J Comp Phys 2003;191:392. [19] Park N, Yoo JY, Choi H. Discretization errors in large eddy simulation: on the suitability of centered and upwind-biased compact difference schemes. J Comp Phys 2004;198:580. [20] Kravchenko AG, Moin P. On the effect of numerical errors in large eddy simulation of turbulent flows. J Comp Phys 1997;131:310. [21] Pope SB. Turbulent flows. Cambridge University Press; 2000. [22] Blasius H. Grenzschichten in Flussigkeiten mit kleiner Reibung. Z Math Phys 1908;56:1. [23] Ghia U, Ghia KN, Shin CT. High Re solutions for incompressible flow using the Navier Stokes equations and a multigrid method. J Comp Phys 1982;48:387. [24] Bottella O, Peyret R. Benchmark spectral results on the lid-driven cavity flow. Comput Fluids 1998;27:421. [25] Erturk E, Corke TC, Gokcol C. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Int J Numer Meth Fluids 2005;48:747. [26] Richardson LF. The approximate arithmetical solution by finite differences of physical problems including differential equations, with an application to the stresses in a masonry dam. Philos Trans Royal Soc Lond Ser A 1910;210:307. [27] Chao J, Haselbacher A, Balachandar S. A massively parallel multi-block hybrid compact–WENO scheme for compressible flows. J Comp Phys 2009;228:7473. [28] Honji H. The starting flow down a step. J Fluid Mech 1975;69(2):229. [29] Gerrard JH. The mechanics of the formation region of vortices behind bluff bodies. J Fluid Mech 1966;25:401. [30] Shah P, Rovagnati B, Mashayek F, Jacobs GB. Subsonic compressible flow in two-sided lid-driven cavity. Part I: Equal walls temperatures. Int J Heat Mass Transfer 2007;50:4206.