A parallel, finite-volume algorithm for large-eddy simulation of turbulent flows

Computers & Fluids 29 (2000) 877±915

www.elsevier.com/locate/compfluid

A parallel, ®nite-volume algorithm for large-eddy simulation of turbulent ¯ows Trong T. Bui NASA Dryden Flight Research Center Edwards CA, USA Received 4 February 1999; received in revised form 1 July 1999; accepted 9 September 1999

Abstract A parallel, ®nite-volume algorithm has been developed for large-eddy simulation (LES) of compressible turbulent ¯ows. This algorithm includes piecewise linear least-square reconstruction, trilinear ®nite-element interpolation, Roe ¯ux-di€erence splitting (FDS), and second-order MacCormack time marching. A systematic and consistent means of evaluating the surface and volume integrals of the control volume is described. Parallel implementation is done using the message-passing programming model. To validate the numerical method for turbulence simulation, LES of fully developed turbulent ¯ow in a square duct is performed for a Reynolds number of 320 based on the average friction velocity and the hydraulic diameter of the duct. Direct numerical simulation (DNS) results are available for this test case, and the accuracy of this algorithm for turbulence simulations can be ascertained by comparing the LES solutions with the DNS results. For the ®rst time, a ®nite volume method with Roe FDS was used for LES of turbulent ¯ow in a square duct, and the e€ects of grid resolution, upwind numerical dissipation, and subgrid-scale dissipation on the accuracy of the LES are examined. Comparison with DNS results shows that the standard Roe FDS adversely a€ects the accuracy of the turbulence simulation. For accurate turbulence simulations, only 3±5% of the standard Roe FDS dissipation is needed. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Unstructured; Parallel; Finite-volume CFD algorithm; Compressible

E-mail address: [email protected] (T.T. Bui). 0045-7930/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 9 9 ) 0 0 0 4 0 - 7

878

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Nomenclature A ^ jAj As c C CFL CS cv d d‡ D DH d~s dv Et f ~ F ~i F ~v F G I j J k K n~ N p Pg q2 ql R S skl t T u u‡ uave ÿu 0 v 0 ut2

Area Roe ¯ux-di€erence splitting matrix Area of duct side walls Local speed of sound SGS model constant Courant number Smagorinsky constant speci®c heat at constant volume Normal distance from a solid wall Normal distance from a solid wall in wall units, d ‡ ˆ rumt d Entire ¯ow domain Hydraulic diameter Elemental surface area on the boundary of a control volume Elemental volume of a control volume Total energy/unit volume Normal component of the inviscid ¯ux vector Flux vector of the Navier±Stokes equations Inviscid ¯ux vector Viscous ¯ux vector Spatial ®lter used in the LES equations Total number of grid points in the streamwise direction Jacobian determinant Total number of grid points in the wall-normal direction Conduction heat-transfer coecient Total number of grid points in the spanwise direction Normal unit vector Total number of parallel processing nodes static pressure Mean pressure gradient Trace of the SGS Reynolds stress tensor SGS term in the LES energy equation Speci®c gas constant Total area on the boundary of a control volume Velocity gradient tensor Time Temperature x-component velocity Velocity in wall units, u‡ ˆ uut Mean streamwise velocity Mean Reynolds stress

T.T. Bui / Computers & Fluids 29 (2000) 877±915

ut U v V w x, y, z y‡ dkl D Dts e1 m x, Z, z r skl tkl tw

879

q Friction velocity, ut ˆ trw State vector of the Navier±Stokes equations y-component velocity Total volume of a control volume z-component velocity Coordinates of the physical space y-coordinate in wall units,y‡ ˆ rumt y Kronecker delta Width of ®lter used in LES equations Sampling time Scaling factor for Roe ¯ux-di€erence splitting Molecular viscosity coecient Coordinates of the computational space Static density SGS term in the LES momentum equations (the SGS Reynolds stress tensor) Viscous stress tensor Wall shear stress

Subscripts a Node index L ¯ow conditions to the left of a cell face R ¯ow conditions to the right of a cell face Superscript n time level ± cell-averaged quantities in the Navier±Stokes equations, or ®ltered or large-scale quantities in the LES equations ~ Favre-®ltered (density-weighted) variables 4 vector quantity 1. Introduction Although direct numerical simulation (DNS) and large-eddy simulation (LES) of turbulence have became important tools to study ¯ow physics, their application to turbulent ¯ows in complex engineering applications is presently not possible because of the tremendous computational resources required. However, technologies that could make such a feat possible in the future are available today. These technologies include LES of turbulent ¯ows, unstructured computational ¯uid dynamics (CFD) algorithms, and parallel computer systems. The above three technologies have been the subjects of ongoing intensive research, and a large body of knowledge has been accumulated on each of these subjects. Yet there have been relatively few studies into the use of unstructured CFD algorithms for DNS and LES of turbulent ¯ows. This is because unstructured CFD algorithms, DNS, and LES are very expensive in terms of computer memory and processing time.

880

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Many popular unstructured CFD codes use the ®nite volume method. Two principal ingredients in a ®nite volume method are reconstruction and ¯ux evaluation. Reconstruction is the process of computing the spatial gradients of the ¯ow variables. In structured methods, this is done using ®nite di€erence formulas and the (i, j, k ) grid indices. The (i, j, k ) block structure is not used in unstructured methods, and other techniques have been devised to compute the gradient of the ¯ow variables. One successful technique is the least-square reconstruction algorithm proposed by Barth [1], where the ¯ow gradients in a cell are computed through a least-square ®tting procedure. Although this technique has been used quite successfully for the Euler and RANS simulations, it has not yet been applied for LES or DNS of turbulent ¯ows. Flux evaluation is another essential ingredient of structured as well as unstructured ®nite volume methods. In ®nite volume methods, the numerical approximation for the spatial variation of the ¯ow variables is discontinuous across a cell boundary. Therefore, a ¯ux formula is needed to compute a single ¯ux at a cell boundary given the two di€erent ¯ow states on the left and right sides of the boundary. If the physical propagation of information of the inviscid ¯ow is taken into account in computing the inviscid ¯ux, then this results in a family of numerical methods known as upwinding. Two main approaches in upwinding methods are ¯ux vector splitting (FVS) [2] and ¯ux di€erence splitting (FDS) [3]. The FVS approach has been used in an upwind-biased ®nite di€erence method developed by Rai and Moin [4] for the DNS of spatially evolving compressible turbulent boundary layer over a ¯at plate. On the other hand, FDS has not been used in DNS and LES of turbulent ¯ows. Yet there is a potential for improvement in accuracy with the FDS approach. Van Leer et al. [5] have found that using the FVS method resulted in excessive dissipation in the boundary layer region for viscous computations, whereas the FDS method provided good results, even for coarse grids. The CFD algorithm developed in this paper incorporates both the least-square reconstruction technique and the FDS method mentioned previously. This will allow, for the ®rst time, an examination of these techniques for turbulence simulation. Since both of these approaches are very expensive in terms of computational resources, parallel computer systems are used in the present investigation. In this paper, a new ®nite volume method is described. This new method applies to both structured and unstructured grids, and it includes a systematic and consistent means of evaluating the surface and volume integrals of the control volume. The parallel implementation and parallel eciency of the method are discussed. The accuracy of this method for turbulence simulations are examined, in detail, using the test case of LES of fully developed turbulent ¯ow in a straight square duct.

2. Numerical algorithm Development of the numerical algorithm has previously been described in detail in Ref. [7]. This algorithm has been validated for time-accurate inviscid Euler simulations [8] and threedimensional viscous Navier±Stokes simulations [9] with good results. To describe the numerical algorithm, the Navier±Stokes equations are used in this section. These equations can be written

T.T. Bui / Computers & Fluids 29 (2000) 877±915

881

in vector form as @U ~ ~ ‡r
Fˆ0 @t

…1†

~ is the ¯ux vector of the Navier±Stokes equations. where U is the state vector and F While the ®nite di€erence and ®nite element methods start from the di€erential form of the governing equations, the ®nite volume method discretizes the Navier±Stokes equations directly in the integral form, ensuring the conservation of mass, momentum, and energy both locally at the discrete cell level and globally over the entire ¯ow domain. Conservation is important for capturing shocks and other ¯ow discontinuities accurately in high speed compressible ¯ow simulations, and it is the strength of ®nite volume methods. Since many aerospace engineering applications involve high speed compressible ¯ows, the ®nite volume method is used in this investigation. Using the ®nite volume method, Eq. (1) is integrated over a ®nite volume. Assuming the grid does not change with time and using the Gauss divergence theorem, the resulting equation is … d Å 1 ~
d~s F …2† Uˆÿ dt V S where Å ˆ 1 U V

… U dv

…3†

~
d~s ˆ F ~
n~ ds F

…4†

V

and

To numerically solve Eq. (2), the major steps of the solution procedure are reconstruction, ¯ux computation, and evolution. This standard, ®nite-volume solution procedure has been used in previous works and has been described in detail by Barth [1]. The steps are given below. 2.1. Step one: reconstruction For the ®rst step, reconstruction, a cell-centered scheme is used. The piecewise linear, leastsquare reconstruction procedure used here is similar to those used by Barth [1] and Coirier [10]. Each of the ®ve primitive variables r, u, v, w, and p is assumed to linearly vary within a ®nite volume as: Å ‡ Ux …x ÿ x † ‡ Uy …y ÿ y†  ‡ Uz …z ÿ z † U…x, y, z † ˆ U

…5†

where U can be any of the above variables. The bars in Eq. (5) denote cell-averaged values as de®ned in Eq. (3). When used for high-speed compressible ¯ow simulations, a gradient limiter is normally used in Eq. (5) to ensure that the reconstruction polynomial does not produce new extrema near a ¯ow discontinuity such as a shock wave. In this paper, the gradient limiter is not used because the test case is low-Mach number turbulent ¯ow in a square duct.

882

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Following Coirier [10], the gradients Ux , Uy , and Uz in the target cell are computed using a least-square procedure that minimizes the sum of the squares of the di€erences between the values computed using the reconstruction polynomial from the target cell and the cell averages of the support set. For a three-dimensional hexahedron cell, the support set is the six neighboring cells that share their faces with the target cell. Algebraically, the minimization statement above can be expressed as: 32 3 2 3 2 Ux c1 a1 b1 b2 4 b1 a2 b3 54 Uy 5 ˆ 4 c2 5 …6† b2 b3 a3 Uz c3

a1 ˆ

6 X

…x i ÿ x 0 † 2

iˆ1

a2 ˆ

6 X ÿ
2 yi ÿ y0 iˆ1

a3 ˆ

6 X

…zi ÿ z0 † 2

iˆ1

b1 ˆ

6 X iˆ1

b2 ˆ

6 X

ÿ
…x i ÿ x 0 † yi ÿ y0

…x i ÿ x 0 †…zi ÿ z0 †

iˆ1

6 X ÿ
yi ÿ y0 …zi ÿ z0 † b3 ˆ iˆ1

c1 ˆ

6 ÿ X


U i ÿ U 0 …x i ÿ x 0 †

iˆ1

c2 ˆ

6 ÿ X iˆ1


ÿ
U i ÿ U 0 yi ÿ y0

T.T. Bui / Computers & Fluids 29 (2000) 877±915

c3 ˆ

6 ÿ X


U i ÿ U 0 …zi ÿ z0 †

883

…7†

iˆ1

where i ˆ 1±6 denotes the neighboring cells, and i ˆ 0 denotes the target cell. 2.2. Step two: ¯ux computation With the piecewise linear reconstruction, the unknown variables are continuous and assumed to linearly vary within a ®nite volume. However, no guarantee exists that the variables will be continuous across adjacent ®nite volumes because a di€erent polynomial is used in each ®nite volume. As a result, a ¯ux formula is needed to compute a single ¯ux at a ®nite-volume boundary using ¯uxes from the adjacent volumes. In the numerical solution of the Navier± Stokes equations, splitting the total ¯ux vector into the inviscid ¯ux vector and viscous ¯ux vector is convenient: ~ v
n~ ~
n~ ˆ F ~ i
n~ ‡ F F

…8†

For the viscous ¯ux, a simple arithmetic average is used. The normal component of the ~ i
n~, is approximated using the Roe FDS method [3] without the entropy inviscid ¯ux vector, F correction. The entropy correction is normally used to remove the nonphysical expansion shock at the sonic transition point and is not needed here because of the low-Mach number test case. De®ne the normal component of the inviscid ¯ux vector as ~ i
n~ fˆF

…9†

Then f can be computed using the Roe FDS method as 1 1 à f ˆ …f L ‡ f R † ÿ jAj…U R ÿ UL † 2 2

…10†

where the subscripts L and R denote the ¯ow conditions to the left and right sides of the cell face, respectively. Fig. 1 shows the de®nitions used for the left and right states. Consider a cell, A, and its neighbor, B, sharing a common face, 1±2. When the ¯ux across face 1±2 is computed for cell A, the left state (L) of the face 1±2 is on the side of cell A, and the right state (R) is on the

Fig. 1. De®nition of the left and right states of a face.

884

T.T. Bui / Computers & Fluids 29 (2000) 877±915

side of cell B. This de®nition of the left and right states of a face is used because the face normal unit vector n~, which also serves as the locally one-dimensional coordinate system for the wave propagation across face 1±2, points from cell A to cell B. The L and R states are reversed when the ¯ux is computed for cell B. 2.3. Step three: evolution The two-stage, second-order, MacCormack time-marching algorithm [6] is used to advance the solution in time. This explicit predictor-corrector time-marching method is accurate, ecient, and simple to implement on parallel computer systems. Eq. (2) can be rewritten as d Å UˆR dt

…11†

When applied to Eq. (11), the MacCormack time-marching method gives Å n‡1 ˆ U Å n ‡ DtRn U

…12†

  1 Å n‡1 Å n n‡1 n‡1 Å U ˆ ‡ U ‡ DtR U 2

…13†

In addition to the major solution steps outlined above, it is important to describe how the volume and surface integrals in Eqs. (2) and (3) are evaluated. Present unstructured ®nite volume CFD algorithms lack a consistent and systematic means to approximate the spatial coordinates as well as evaluate the volume and surface integrals. To complete the development of this ®nite volume algorithm, a means of approximating the spatial coordinates is borrowed from the ®nite element methods. The spatial coordinates x, y, and z of a ®nite volume are approximated by a continuous trilinear hexahedral element [11]. All of the integrations are then numerically evaluated using the one-point Gauss quadrature formula. ~ y, z† space is mapped to a trilinear hexahedral element in the Each cell in the physical x…x,

Fig. 2. Local mapping between the physical ®nite-volume and the trilinear hexahedral element.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

885

~ Z, z† space as shown in Fig. 2. The nodes, indexed 1 to 8, have the nodal coordinates in x…x, the x~ space shown in Table 1. This type of mapping is di€erent from the mapping used in generalized curvilinear ®nitedi€erence methods. In the ®nite-di€erence methods, the mapping applies to the entire computational block. In this algorithm, the local mapping applies to the local cell only. Each cell has its own mapping function, which is similar to a ®nite-di€erence multiblock method in which each cell is its own block. This ability gives considerable ¯exibility in the grid topology that can be used and is one of the strengths of this unstructured ®nite-volume algorithm. The mapping from the x~ space to the x~ space is given by:     X 8 ~ Na x~ x a x x ˆ aˆ1

    X 8 y x~ ˆ Na x~ ya aˆ1

    X 8 ~ z x ˆ Na x~ za

…14†

aˆ1

  ÿ
1 Na x~ ˆ …1 ‡ xa x† 1 ‡ Za Z …1 ‡ za z† 8

…15†

where the subscript a denotes the node index, ranging from 1 to 8. In the physical (x, y, z ) coordinate system, node a has the coordinate …x a , ya , za †: In the computational …x, Z, z† space, node a has the coordinate …xa , Za , za †: The coordinates …x a , ya , za † vary from cell to cell, depending on the physical grid. The coordinates …xa , Za , za † are the same for every cell and are shown in Table 1. To evaluate the volume integral, the following relation [11] is used:

Table 1 Nodal coordinates in the ~z space Node index (a )

za

Za

za

1 2 3 4 5 6 7 8

ÿ1 1 1 ÿ1 ÿ1 1 1 ÿ1

ÿ1 ÿ1 1 1 ÿ1 ÿ1 1 1

ÿ1 ÿ1 ÿ1 ÿ1 1 1 1 1

886

T.T. Bui / Computers & Fluids 29 (2000) 877±915

… V

f…x, y, z † dv ˆ

…1 …1 …1 ÿ1 ÿ1 ÿ1

f…x…x, Z, z†, y…x, Z, z†, z…x, Z, z††j…x, Z, z† dx dZ dz

where j is the Jacobian determinant, de®ned as 2 3 ! xx xZ xz @ x~ ˆ det4 yx yZ yz 5 j ˆ det @ ~z zx zZ zz

…16†

…17†

The partial derivatives x x , yZ , zz , and so forth can be obtained by di€erentiating Eq. (14). Evaluating the determinant in Eq. (17) results in the following: j ˆ x z yx zZ ÿ x x yz zZ ÿ x z yZ zx ‡ x Z yz zx ‡ x x yZ zz ÿ x Z yx zz

…18†

To evaluate the integral in Eq. (16), the one-point Gauss quadrature formula is used. In one dimension, the Gauss quadrature formulas are optimal, which means that accuracy of order (2n ) is achieved using (n ) integration points. Gaussian rules for integrals in several dimensions are constructed by employing the one-dimensional Gaussian rules on each coordinate separately. In three dimensions, the one-point Gaussian rule is given as …1 …1 …1 ÿ1 ÿ1 ÿ1

f…x, Z, z† dx dZ dz ˆ 8f…x ˆ 0, Z ˆ 0, z ˆ 0†

Using the tools described above, the volume of the cell is computed as … V ˆ 1 dv V

Vˆ

…1 …1 …1 ÿ1 ÿ1 ÿ1

j…x, Z, z† dx dZ dz

…19†

…20†

…21†

Using Eq. (19), V ˆ 8j…x ˆ 0, Z ˆ 0, z ˆ 0†

…22†

so that the cell volume is approximately eight times the Jacobian determinant evaluated at the center of the cell …x ˆ 0, Z ˆ 0, z ˆ 0†: Note that Eq. (22) contains two approximations: the physical coordinates …x, y, z† in the cell are approximated by Eq. (14), and the one-point Gaussian rule given by Eq. (19) is used for the numerical integration. Better approximation of the cell volume can be obtained using a higher-order approximation for the physical coordinates and a Gaussian rule with more points. Computing the centroid of the cell also requires the evaluation of the volume integral. The  y,  z†:  The numerical approximation for the x coordinates of a cell centroid are given by …x, coordinate of the centroid is developed below. Approximations for the y and z coordinates are made exactly the same way.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

… x ˆ

v

x dv

x ˆ

…23†

V

…1 …1 …1 ÿ1 ÿ1 ÿ1 …1 …1

887

x…x, Z, z†j…x, Z, z† dx dZ dz …1 j…x, Z, z† dx dZ dz

…24†

ÿ1 ÿ1 ÿ1

x ˆ

8x…x ˆ 0, Z ˆ 0, z ˆ 0†j…x ˆ 0, Z ˆ 0, z ˆ 0† 8j…x ˆ 0, Z ˆ 0, z ˆ 0†

…25†

or x ˆ x…x ˆ 0, Z ˆ 0, z ˆ 0† y ˆ y…x ˆ 0, Z ˆ 0, z ˆ 0† z ˆ z…x ˆ 0, Z ˆ 0, z ˆ 0†

…26†

From Eq. (14), x…x ˆ 0, Z ˆ 0, z ˆ 0† ˆ

8 1X xa 8 aˆ1

y…x ˆ 0, Z ˆ 0, z ˆ 0† ˆ

8 1X ya 8 aˆ1

z…x ˆ 0, Z ˆ 0, z ˆ 0† ˆ

8 1X za 8 aˆ1

…27†

With the approximations in this algorithm, the coordinates of the cell centroid are simply the averages of the coordinates of the eight nodes de®ning the three-dimensional hexahedron cell. The task of evaluating the surface integrals is described next. The surface integrals on the „ ~
d~s: To develop the procedure, one face of the right side of Eq. (2) are of the form F hexahedron shown in Fig. 2 is considered. The results for the other faces can be obtained using the same procedure. Consider face 6-2-3-7 of the hexahedron in Fig. 2. Fig. 3 shows this face redrawn for convenience. Because x ˆ 1:0 for the nodes 6, 2, 3, 7 as well as any other point on this face, Eq. (14) reduces to the following:

888

T.T. Bui / Computers & Fluids 29 (2000) 877±915

xˆ


1 X ÿ 1 ‡ Za Z …1 ‡ za z†x a 4 aˆ6, 2, 3, 7

yˆ


1 X ÿ 1 ‡ Za Z …1 ‡ za z†ya 4 aˆ6, 2, 3, 7

zˆ


1 X ÿ 1 ‡ Za Z …1 ‡ za z†za 4 aˆ6, 2, 3, 7

…28†

so that x, y, and z are functions of Z and z only. Following the development outlined by Greenberg [12], the tangent vectors d~s1 and d~s2 are tangents on the plane 6-2-3-7. d~s1 is de®ned to be along the Z ˆ constant curve on the face, and d~s2 is along the z ˆ constant curve. The vector d~s1 may be expressed as follows: d~s1 ˆ dxi~‡ dyj~‡ dzk~ dx ˆ x x dx ‡ x Z dZ ‡ x z dz dy ˆ yx dx ‡ yZ dZ ‡ yz dz dz ˆ zx dx ‡ zZ dZ ‡ zz dz

…29†

where dx ˆ 0 because the entire plane 6-2-3-7 is an x-constant plane, and dZ ˆ 0 because the vector is de®ned to be along the Z ˆ constant curve. So   ~ ~ ~ d~s1 ˆ x z i ‡ yz j ‡ zz k dz …30† and similarly

Fig. 3. Coordinate systems for surface integral evaluation.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

 d~s2 ˆ x Z i~‡ yZ j~‡ zZ k~ dZ 

889

…31†

The elemental area vector d~s, denoted by the shaded parallelogram, can be computed as

or

d~s ˆ d~s2  d~s1

…32†

i h l~s ˆ …yZ zz ÿ zZ yz †i~‡ …zZ x z ÿ x Z zz †j~‡ …x Z yz ÿ yZ x z †k~ dZ dz

…33†

Note that the order of the cross product in Eq. (32) is chosen so that the elemental area vector ~ d~s is positive pointing out of the cell and negative pointing into the cell. With the vector F de®ned as ~ ˆ Fx i~‡ Fy j~‡ Fz k~ F

…34†

the dot product is ~
d~s ˆ …yZ zz ÿ zZ yz †Fx ‡ …zZ x z ÿ x Z zz †Fy ‡ …x Z yz ÿ yZ x z †Fz  dZ dz F

…35†

Finally, using the one-point Gaussian rule, the surface integral can be evaluated as … ~
d~s ˆ 4…yZ zz ÿ zZ yz †Fx ‡ …zZ zz ÿ x Z zz †Fy ‡ …x Z yz ÿ yZ x z †Fz jxˆ1, Zˆ0, zˆ0 F

…36†

…6237 †

The surface integrals for the other faces can be approximated in an analogous fashion. The results are given below. … ~
d~s ˆ 4…zZ yz ÿ yZ zz †Fx ‡ …x Z zz ÿ zZ x z †Fy ‡ …yZ x z ÿ x Z yz †Fz jxˆÿ1, Zˆ0, zˆ0 …37† F …1584 †

… …8734 †

… …1265 †

… …5678 †

… …2143 †

~
d~s ˆ 4…zx yz ÿ yx zz †Fx ‡ …x x zz ÿ zx x z †Fy ‡ …yx x z ÿ x x yz †Fz jxˆ0, Zˆ1, zˆ0 F

…38†

~
d~s ˆ 4…yx zz ÿ zx yz †Fx ‡ …zx x z ÿ x x zz †Fy ‡ …x x yz ÿ yx x z †Fz jxˆ0, Zˆÿ1, zˆ0 F

…39†

~
d~s ˆ 4…yx zZ ÿ zx yZ †Fx ‡ …zx x Z ÿ x x zZ †Fy ‡ …x x yZ ÿ yx x Z †Fz jxˆ0, Zˆ0, zˆ1 F

…40†

~
d~s ˆ 4…zx yZ ÿ yx zZ †Fx ‡ …x x zZ ÿ zx x Z †Fy ‡ …yx x Z ÿ x x yZ †Fz jxˆ0, Zˆ0, zˆÿ1 F

…41†

Symbolically, the surface integrals of the opposite faces are the negative of each other (for example, face x ˆ 1 as given by Eq. (36) and face x ˆ ÿ1 as given by Eq. (37)). However, for

890

T.T. Bui / Computers & Fluids 29 (2000) 877±915

actual numerical values, each of the integrals will need to be separately evaluated because the integrands depend on the coordinates of the nodal points …x a , ya , za † and the Gaussian point coordinates. To calculate the Jacobian determinant of the coordinate transformation, Eqs. (14) and (15) are used. For example, the derivative x x can be evaluated as follows:     X 8 ~ Na x~ x a x x ˆ aˆ1

xx ˆ

8 X Na, x x a

…42†

aˆ1

Although the above algorithm can be used for both structured and unstructured grids, a structured grid was used in this study. This allows the evaluation of the numerical accuracy for turbulence simulation without the complication of possible grid e€ects. Also, the structured grid facilitates direct comparison with the DNS results. Finally, the structured grid provides considerable savings in computer memory and CPU time for this simple duct geometry.

3. Parallel implementation Details of parallel implementation and the results of the parallel performance studies can be found in Refs. [7] and [9]. The code was implemented on parallel computer systems using the message-passing programming model and message-passing libraries Message-Passing Interface (MPI) and Parallel Virtual Machine (PVM). A simple but e€ective domain decomposition strategy was used. The CFD mesh of the duct was partitioned in the streamwise direction. For a mesh of I  J  K cells, each parallel processing node is responsible for advancing (I/N )  J  K in time, where N is the total number of parallel processing nodes available. For distributed memory parallel computers, it is important that the communication cost is small compared to the computation e€ort. In this parallel algorithm, the communication cost is minimized by restricting the data exchange to the left and right neighboring nodes only. All necessary data are sent and received in one single message passing call as one large contiguous binary data string. This reduces the amount of work in packing, unpacking, sending, and receiving data. Furthermore, the communication cost for this algorithm scales well, i.e., it does not increase with increasing number of nodes. Regardless of the number of parallel processing nodes, each node only needs to exchange data with its two adjacent neighbors, and this can be arranged to occur simultaneously for all nodes. Finally, dedicated communication channels between the adjacent processors were used whenever possible to accelerate the communication speed. To determine the parallel speedup of this algorithm, Navier±Stokes simulations of laminar ¯ow in an S-duct were performed using three di€erent mesh sizes [9]. The parallel speedup factor is plotted in Fig. 4. As can be seen in Fig. 4, the actual parallel speedup using both the

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 4. Parallel speedup results for the IBM SP2 and Cray T3D.

891

892

T.T. Bui / Computers & Fluids 29 (2000) 877±915

T3D and the SP2 systems is very good for most cases, and the speedup is nearly ideal for large (I/N ) ratios. For a ®xed grid size, increasing the number of processing nodes eventually causes the actual parallel speedup curve to drop below the ideal speedup line. This is caused by an increase of the parallel overhead inherent in the algorithm. The parallel overhead is the additional work that is required for parallel processing. In this algorithm it includes the communication cost to exchange data and the extra ¯ux evaluations at the boundaries of each grid partition assigned to a parallel node. For a ®xed grid size, the communication cost is constant, but the overhead in the extra ¯ux evaluations increases as the number of processors increases. But even for the worse case, where each node has only two streamwise cells, the parallel eciency is still quite good. For example, when 40 nodes are used to solve the problem with an 81  32  32 grid, the Cray T3D has a parallel eciency of 90%. It is informative to compare the computing performance of a parallel computer system and a serial supercomputer. Such a comparison is shown in Fig. 5. As discussed previously, better parallel performance is obtained for larger grids. For this CFD code, it can be seen that the T3D is approximately 10±30 times faster than the YMP. Note that parallel speedup is highly dependent on the CFD algorithm and code implementation. Signi®cant speedup is possible in this case, because an explicit algorithm is used. Also, the code was carefully implemented and optimized on the T3D. For other CFD algorithms and code implementations, the parallel speedup may vary. 4. Governing equations for large-eddy simulation In LES, the large scale of turbulence is computed directly in the numerical simulation, and the e€ects of the small scale stresses are modeled using a subgrid-scale (SGS) model. The governing equations for LES of turbulent ¯ows can be obtained from ®ltering (local volumeaveraging) the compressible Navier±Stokes equations. From Moin et al. [13] the LES equations

Fig. 5. Performance comparison between the Cray T3D and YMP.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

893

for compressible ¯ows (using tensor notation) are given by the following: @ @ r u~k ˆ 0 r ‡ @t @x k

…43†

@ @ @ @ t~ kl @skl r u~k u~l ‡ p ÿ r u~k ‡ ‡ ˆ0 @t @x l @x k @x l @ xl

…44†

  @ E~ t @ ~ @ @ @ @ ~ @skl @q1 ~ pu~ l ÿ ‡ cv ˆ0 ‡ Et u~l ‡ T ‡ u~k k u~k t~ kl ÿ @x l @x l @x l @x l @x l @t @x l @x l

…45†

The bar in the LES equations (43)±(45) denotes a ®ltered or large-scale ¯ow quantity, de®ned as … 0 0 0 f ˆ G…x~ ÿ x~ †f…x~ † dx~ D

where G is a spatial ®lter and the integral is over the ¯ow domain, D. The tilde in the LES equations denotes a Favre-®ltered (density-weighted) variable, de®ned as rf f~ ˆ r The ®ltered ideal gas equation of state is given by  T~ p ˆ rR Moin et al. [13] used an internal energy equation in their derivation of the LES equations. The LES total energy equation, Eq. (45), is obtained from adding the dot product of the LES momentum equation, Eq. (44), and the ®ltered velocity ®eld u~k to the LES internal energy equation in Moin et al. [13] The LES equations given by Eqs. (43)±(45) are essentially the Navier±Stokes equations written for the ®ltered variables plus the additional subgrid terms in the momentum and total energy equations. Thus, the numerical algorithm developed in the last section can be used to solve the LES equations. The treatment of the subgrid terms are discussed in the next section.

5. Subgrid-scale models Detailed studies have previously been performed to assess the relative importance of the subgrid terms in the ®ltered total energy equation for compressible turbulent shear ¯ows at di€erent Mach numbers. [14,15]. These studies led to the conclusion that the energy subgrid terms may be neglected if the Mach number of the simulation is low. Because of the low-Mach number of the turbulent square duct test case, this assumption is used. The subgrid term in the momentum equations, Eq. (44), is

894

ÿ
~ k u~l skl ˆ r ug k ul ÿ u

T.T. Bui / Computers & Fluids 29 (2000) 877±915

To close the system of LES equations, this term needs to be modeled. In Moin et al. [13] this term is approximated as follows:   1~ 1 2 ~ ~  jSj Skl ÿ Smm dkl ‡ q 2 dkl skl ˆ ÿ2CrD …46† 3 3 where q 2 ˆ sii is the trace of the SGS Reynolds stress tensor. The ®ltered velocity gradient tensor is   1 @ u~k @ u~ l ~ Skl ˆ ‡ 2 @x l @x k and
ÿ ~ ˆ 2S~kl S~kl 1=2 jSj In Eq. (46), C is a constant to be determined according to the particular SGS model used. For LES of turbulent channel and duct ¯ows using the Smagorinsky SGS model [16] a value of C ˆ 0:01 is commonly used with good results. Note that the constant C in Eq. (39) is the square of the Smagorinsky constant CS ˆ 0:1: Unlike LES of isotropic turbulence, C is not constant in wall-bounded ¯ows and varies according to distance from the wall. The dynamic SGS model developed by Germano et al. [17] would correctly determine the value for C using a dynamic procedure; however, this model is computationally expensive because of the extra ®ltering operations that must be done. Also, a question currently exists on the mathematical well-posedness of this model [18]. Finally, the dynamic model has been known to compensate for the e€ects of numerical dissipation by automatically varying the magnitude of the constant C. One of the main objectives of this research is to quantify the e€ects of the upwind numerical dissipation on the accuracy of the turbulence simulations. Also, the simpler Smagorinsky model has been found to work as well as the dynamic SGS model for this simple test case [19]. As a result, the Smagorinsky SGS model is used in this study with the constant C given by the following:  ‡ 3 !! d …47† C ˆ 0:01 1:0 ÿ exp ÿ 25 where d ‡ is the normal distance from the wall in wall units, de®ned as d‡ ˆ

rut d m

and the friction velocity is de®ned as

r tw ut ˆ r

T.T. Bui / Computers & Fluids 29 (2000) 877±915

895

Because the turbulent ¯ow in the corner of the square duct encounters walls in two di€erent directions, d is taken to be dˆ

2yz ÿ
1=2 y ‡ z ‡ y2 ‡ z2

…48†

Eq. (48) is frequently used in the turbulence modeling of ¯ows in the vicinity of a wall corner. The variables y and z are the normal distances to the nearest walls in the y and z directions. Note that d tends to y as ( y/z ) tends to 0, and d tends to z as (z/y ) tends to 0, which are the intended results. In LES, the width of the ®lter used in the process of volume-averaging the Navier±Stokes equations, D, pis typically chosen to be the grid spacing size. This study de®ned the grid spacing size to be 3 V, where V is the cell volume. The term q 2 in Eq. (46) is the isotropic part of the SGS Reynolds stress tensor. Like the rest of this tensor, the term cannot be calculated directly in an LES and has to be modeled. A number of di€erent models for q 2 has been proposed [20±22]. However, results from recent studies indicated that this term is not important for accurate LES of low-Mach number, lowReynolds number compressible turbulent ¯ows. An evidence in support of the above conclusion was presented by Squires [23] who compared two di€erent models of q 2 in addition to setting q 2 ˆ 0: Squires found essentially no di€erence in the results of LES of compressible isotropic turbulence at a low-Mach number and, in fact, observed that neglecting q 2 slightly improved the agreement between the LES and DNS results. Vreman et al. [14] con®rmed the above results with their simulations of a three-dimensional temporal compressible mixing layer at a mean convective Mach number of 0.2. In a priori tests, the SGS model that neglects q 2 was found to give a better correlation with the DNS results. Furthermore, LESs conducted with the dynamic SGS model for q 2 were unstable for the cases that were studied. For low-Mach number turbulence LES, neglecting the term q 2 will not introduce large errors in the results and is actually desirable in some cases, as the above ®ndings showed. As a result, the term q 2 is neglected in this study. This assumption is analogous to the Stokes assumption for the viscous stress tensor in the Navier±Stokes equations. With the term q 2 omitted, the SGS stresses can be included in the Navier±Stokes equations by simply replacing the laminar ~  2 jSj: viscosity coecient m with meff , where meff ˆ m ‡ mSGS and mSGS ˆCrD 6. Large-eddy simulation of turbulent ¯ow in a square duct To validate the numerical method for turbulence simulations of duct ¯ows, LES of fully developed turbulent ¯ow in a square duct is performed. Although this test case has a simple geometry, the turbulent ¯ow physics is quite complex, and the accuracy of the numerical algorithm can be evaluated without the complication of complex grid topologies. For the purpose of comparison, a low-Reynolds number, square duct DNS solution is

896

T.T. Bui / Computers & Fluids 29 (2000) 877±915

available. This DNS database was used by Mompean et al. [24] to evaluate nonlinear k±e turbulence models. Another DNS solution of the fully developed turbulent square duct ¯ow at a slightly lower Reynolds number is also available.[25] Fig. 6 shows the coordinate system and geometry for the square duct ¯ow. In this test case, the Reynolds number based on the mean streamwise velocity and hydraulic diameter is 4800. Based on the friction velocity and hydraulic diameter, the Reynolds number is 320. Table 2 summarizes the ¯ow properties for the test case, assuming an average Mach number of 0.3 and standard sea level properties for air. The computational domain size used in the LES is 12DH  DH  DH : In choosing the size of the computational domain, care must be taken to ensure that the length of the computation domain is large enough to adequately contain the largest turbulence structure. Two-point velocity correlations for three di€erent cross-stream positions in the duct were computed from the DNS solution by Gavrilakis [25]. The correlations for all three velocity components become essentially 0 at a duct length of approximately 6DH , so that a length of 12DH should be adequate to capture the streamwise turbulence structures. Two di€erent grids are used in the present LES, and Table 3 shows the simulation parameters. The sampling time for the turbulent statistics is large compared to the time step size, but is small compared to the eddy turnover time in order to capture the unsteadiness of the turbulent ¯ow. The periodic boundary condition used for the in¯ow and exit boundary of the square duct is similar to the one used by Coleman et al. [26]. With this boundary condition, all of the ¯ow conditions are periodic at the duct inlet and exit planes. The driving pressure gradient in the duct is speci®ed in the ¯ow equations as an extra body force term. An isothermal no-slip wall boundary condition is used at the solid walls. Conservation of energy in the duct requires some heat transfer through the wall to balance the work of the body force term, and the isothermal wall boundary condition allows the simulation to set the necessary heat transfer rate automatically. To reduce the number of iterations required for convergence, the initial conditions for the

Fig. 6. Coordinate system and geometry for the square duct ¯ow.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

897

Table 2 Flow properties for the turbulent square duct test case Flow properties

Values

Average Mach number, uave c Average streamwise velocity, uave Average friction velocity, ut Re ˆ ruavem DH Ret ˆ rutmDH Hydraulic diameter, DH Mean pressure gradient, Pg H Eddy turnover time, 0:5D ut

0.3 102.4 m/s 6.83 m/s 4800 320 7.37  10ÿ4 m ÿ289,646 Pa/m 5.4  10ÿ5 s

large eddy simulations were obtained from interpolating a DNS solution provided by Gavrilakis. The DNS was done using a 768  127  127 grid. The current simulations were done using two di€erent grid sizes, 129  90  90 (grid A) and 257  100  100 (grid B). The ®ner LES grid B, which gave good results, is approximately 20% of the total size of the DNS grid. The convergence of the LES is determined by monitoring the time history of the total wall shear stress. For fully developed turbulent ¯ow in a straight square duct, conservation of the mean streamwise momentum shows that the mean driving pressure gradient and the total wall shear stress are related by the following: … tw dA ˆ ÿVPg …49† As

The surface integral is over the four side walls of the square duct, so that As ˆ 4  12  DH2 : Pg 3 is the mean driving „ pressure gradient, and V is the total volume of the duct, given by 12  DH : 1 De®ning t w ˆ As As tw dA and Pg ˆ dp=dx, the familiar relation between the mean pressure gradient and wall shear stress in fully developed ¯ow in a square duct can be recovered.   DH dp …50† t w ˆ ÿ 4 dx Table 3 Parameters for the LES Parameters

Grid A

Grid B

Grid size (I  J  K ) Number of cells Minimum resolution (in wall units) Maximum resolution (in wall units) Sampling time, Dts (s) Time step size, Dt (s) CFL number

129  90  90 1,013,888 30  1.88  1.88 30  4.86  4.86 6.0  10ÿ7 3.0  10ÿ9 0.98

257  100  100 2,509,056 15  1.69  1.69 15  4.37  4.37 5.0  10ÿ7 2.5  10ÿ9 0.93

898

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 7 shows the time history of the total side wall shear force for the LES using grid B. The instantaneous side wall shear force level from the LES, shown as a solid line, is seen to ¯uctuate about a mean value, indicating that ¯ow equilibrium has been reached in the current simulation. The time average of the computed side wall shear force is 0.001387 N. This value is in excellent agreement with the exact value of 0.001391 N, computed from Eq. (50) and shown as the horizontal dashed line in Fig. 7. Because the time step is constant, the number of iterations is directly proportional to the time elapsed. The simulation was conducted for 218,600 time iterations, which is approximately 10 eddy turnover times (as de®ned in Table 2). Using 128 processors on the T3D computer (Cray Research; Eagan, Minnesota), the simulation of these 10 eddy turnover times (5.5  10ÿ4 in physical ¯ow time) took 772 h or approximately one month of central processing unit (CPU) time. The same simulation would have taken approximately 150 h on an SP2 (IBM Corporation, Austin, TX) or T3E (Cray Research) computer with the same number of processors. The processing rate in ¯oating point operations per second (¯ops) of the Cray machines was not obtained. On a 64-node SP2, the code achieved 2.2 GFlops. The LES results shown in Figs. 7±17 are obtained using grid B and a modi®ed Roe FDS with an e1 value of 0.03. The modi®cation to the Roe FDS and the de®nition of the e1 parameter will be discussed below. Fig. 8 shows the mean streamwise velocity pro®le along a wall bisector. The LES solution (solid line) is compared with the DNS solution (diamond dots) supplied by Gavrilakis. The mean velocity pro®le in the LES was averaged both in time and space. The agreement can be seen to be very good. The mean velocity pro®le in wall units are plotted in Fig. 9. In this ®gure, the discrete DNS and LES solutions at the cell centroid are plotted. The LES grid is not as dense as the DNS grid, but the agreement is very good. Unlike the 2D turbulent boundary layer and channel ¯ows, the mean wall shear stress in a square duct actually varies around the perimeter of the side wall at a given streamwise location. In Fig. 9, the mean wall shear stress value from Eq. (49) was used to calculate the friction velocity scale. A log-law region can clearly be seen for y‡ values greater than 20. However, the log-law constants are di€erent from the standard

Fig. 7. Convergence history for the total wall shear stress.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

899

Fig. 8. Mean streamwise velocity pro®le for fully developed turbulent ¯ow in a square duct along the wall bisector.

values for 2D turbulent boundary layer and channel ¯ows. Gavrilakis [25] attributed this di€erence to the in¯uence of the secondary ¯ow ®eld in the square duct. Fig. 10 shows the mean secondary velocity vectors from the LES. In straight ducts of noncircular cross-sections, turbulence-driven secondary ¯ows are known to exist. These ¯ows are di€erent from the pressure-driven secondary ¯ows found in curved ducts. In straight square ducts, the turbulence-driven secondary ¯ows are directed from the center of the duct toward the corners along the corner bisectors, and have been found to be produced by the anisotropy of the Reynolds stresses in the cross-sectional plane of the square duct [27]. Although the magnitudes of these secondary velocities are extremely small compared to the mean average streamwise velocity (on the order of 2% in this simulation), these velocities have been found to be important features of this ¯ow.

Fig. 9. Mean streamwise velocity pro®le in wall units.

900

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Some asymmetry is still evident in Fig. 10, but the overall features of the secondary ¯ows are well captured by the current simulation. To determine the accuracy of the simulation in capturing turbulence-driven secondary ¯ows, the mean secondary velocity pro®les along the lines z=…0:5DH † ˆ 0:15, 0.30, 0.50, 0.70, and 0.80 are compared with the DNS results in Figs. 11±15. Good agreement is obtained between the LES and DNS mean secondary velocity pro®les. In Fig. 16, the mean Reynolds stress pro®le along a wall bisector is compared with the DNS solution. In Fig. 17, the mean turbulence intensities urms , vrms , and wrms are plotted. These results have been quadrant-averaged as well as averaged in space and time. Good agreement is obtained for these turbulence statistics. Although the agreements between the LES and DNS solutions are very good for this simulation, the accuracy of the LES in capturing the turbulence velocity ¯uctuations was found to be highly dependent on the numerical dissipation and the grid size used. The e€ects of the Roe FDS upwinding term and grid size on the computed turbulence velocity ¯uctuations are examined next. 7. E€ect of the Roe ¯ux-di€erence splitting term When this new ®nite volume algorithm was ®rst used to perform LES of turbulent ¯ow in the square duct, incorrect levels of turbulent velocity ¯uctuations were obtained. As a result, the mean Reynolds stress pro®le was signi®cantly lower than expected, and other turbulent mean ¯ow properties were in poor agreement with the DNS solution. A systematic search was conducted for the cause of this problem. Finer mesh, smaller time steps, higher-order time

Fig. 10. Mean secondary velocity vectors from LES with a 257  100  100 grid.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

901

marching methods, better wall boundary conditions, and other SGS models have all been tried without success. Finally, it was discovered that removing the Roe FDS term from the LES was the only way to raise the mean Reynolds stress pro®le. However, without the numerical provided by the FDS term, the LES was numerically unstable. Although the Roe FDS implemented in this code gave good results for Euler [8] and laminar Navier±Stokes [9] test cases, the full Roe FDS term turned out to be too dissipative for LES. Yet without the FDS term, the LES was numerically unstable. This dilemma was resolved with a simple modi®cation to the Roe FDS algorithm. In Eq. (10), the inviscid ¯uxes normal to a cell boundary are approximated as 1 1 ^ f ˆ …f L ‡ f R † ÿ jAj…U R ÿ UL † 2 2 This approximation can be interpreted to state that the normal component of the inviscid ¯ux at a cell boundary is the sum of the central di€erence of the ¯uxes on the left and right states, ^ 1=2…f L ‡ f R †, plus the Roe upwinding dissipation term, ÿf1=2jAj…U R ÿ UL †g: If this interpretation is used, then the amount of Roe upwinding dissipation can be controlled using a

Fig. 11. Mean secondary velocity pro®les along z=…0:5DH † ˆ 0:15:

902

T.T. Bui / Computers & Fluids 29 (2000) 877±915

multiplying factor in front of the Roe FDS term, such that   1 1 ^ f ˆ …f L ‡ f R † ÿ e1 jAj…UR ÿ UL † 2 2

…51†

where e1 can range between 0 and 1. e1 ˆ 0 corresponds to central di€erencing only, and e1 ˆ 1 corresponds to the full Roe FDS. Note that Lin et al. [28], using the same interpretation of the Roe upwinding term as Eq. (51), also concluded that the normal Roe upwinding term produces too much numerical dissipation for computational aeroacoustics applications. They found that using e1 values of approximately 0.1 gave good results for acoustics computations. In the LES conducted here, even smaller values of e1 are required. It should be emphasized that the selection of e1 is not arbitrary. Omitting the Roe FDS term altogether …e1 ˆ 0† causes all calculations to be unstable, and the best turbulence solutions are obtained using the smallest possible values of e1 that can still provide stable calculations. In general, it was found that smaller minimum values of e1 can be used for ®ner grids. For the 129  90  90 grid, the

Fig. 12. Mean secondary velocity pro®les along z=…0:5DH † ˆ 0:3:

T.T. Bui / Computers & Fluids 29 (2000) 877±915

903

minimum value of e1 for stable calculation is 0.05, and for the 257  100  100 grid, the minimum value drops to 0.03. The LES results presented in the preceding section were obtained using e1 ˆ 0:03 and the ®ner grid. To study the e€ect of the Roe upwinding term on the turbulence simulation, LES are made for the same grid size of 129  90  90 but with di€erent values of e1 : Fig. 18 shows the e€ect of Roe FDS on the mean streamwise velocity pro®le. Near the wall, using the full Roe FDS term produces a mean velocity gradient that is much less than both the DNS solution and the LES solution with the reduced Roe FDS. A similar e€ect is observed in Fig. 19, where the mean Reynolds stress pro®le obtained with e1 ˆ 1:0 is much lower than expected. Fig. 20 also shows the excessive numerical dissipation of Roe FDS in the turbulence solution. In this ®gure, the solution with e1 ˆ 1:0 gives a signi®cantly higher level of urms and lower levels of vrms and wrms : Although they did not use the Roe FDS, Wang and Pletcher [29] reported the same problem in their LES of fully developed turbulent channel ¯ow using an upwind CFD algorithm. From studying the results shown in Figs. 18±20, the full Roe FDS upwinding dissipation can be seen to be detrimental to the turbulence solution. Although an improvement in the solution

Fig. 13. Mean secondary velocity pro®les along z=…0:5DH † ˆ 0:5:

904

T.T. Bui / Computers & Fluids 29 (2000) 877±915

quality can be obtained by reducing the contribution of the Roe upwinding term, getting a good solution is not possible with coarse grids. For the 129  90  90 grid, the minimum e1 value for numerical stability is 0.05, Yet the agreement with the DNS solution is still not satisfactory even with this minimum e1 value. To improve the accuracy of the turbulence simulation, it is necessary to use a ®ner grid which in turn allows a smaller minimum e1 value. In the next section, the e€ect of a grid size on the quality of the turbulence solution will be studied. 8. E€ect of the grid size In this section, the LES solutions from two di€erent grid sizes, 129  90  90 and 257  100  100, are compared. Both calculations used the minimum values o e1 that are required for stability. Fig. 21 shows a comparison of the mean streamwise velocity pro®les. Using the ®ner grid in the LES gives an almost perfect agreement with the DNS solution. Figs. 22 and 23 show the

Fig. 14. Mean secondary velocity pro®les along z=…0:5DH † ˆ 0:7:

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 15. Mean secondary velocity pro®les along z=…0:5DH † ˆ 0:8:

Fig. 16. Mean Reynolds stress pro®le along the wall normal bisector.

905

906

T.T. Bui / Computers & Fluids 29 (2000) 877±915

same improvement in the pro®les of the turbulence intensities and the mean Reynolds stress, respectively. Another e€ect of grid size can be observed by examining the mean secondary velocity vectors in a cross section of the square duct. Fig. 24 shows the solution obtained using grid A and e1 ˆ 0:05: The secondary turbulence-induced velocity ®eld is captured by the coarser grid. However, comparing this result with the ®ner grid B result in Fig. 10 shows that the corner vortices in Fig. 10 are somewhat smaller than those shown in Fig. 24. Fig. 25 shows the streamwise-averaged instantaneous secondary velocity vectors from the DNS solution for

Fig. 17. Turbulent intensities along the wall normal bisector.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

907

Fig. 18. E€ect of Roe FDS on the mean streamwise velocity pro®le.

reference. For the DNS grid, the near wall vortical structures are even smaller than either of the LES grids. This comparison shows that the near wall turbulent structures are better resolved with ®ner grids. 9. E€ect of the subgrid-scale model The previous section showed that the LES solution with the ®ne grid gives the best agreement with the DNS solution. Although the LES ®ne grid is only 20% of the size of the DNS grid, the LES grid density in the cross¯ow plane is 60% of the DNS cross¯ow plane

Fig. 19. E€ect of Roe FDS on the mean Reynolds stress pro®le.

908

T.T. Bui / Computers & Fluids 29 (2000) 877±915

density. As the LES grid resolution in the cross¯ow plane approaches the DNS resolution, it is interesting to see the e€ect of the SGS model on the turbulence solution. An additional simulation was performed using the ®ner LES grid, with no SGS model. This simulation is e€ectively a coarse grid DNS. Fig. 26 shows a comparison of the mean streamwise velocity pro®les. The use of the SGS model makes essentially no di€erence in the mean streamwise velocity solution. Small di€erences are also observed in the turbulent velocity ¯uctuations

Fig. 20. E€ect of Roe FDS on the turbulent intensities.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

909

Fig. 21. E€ect of grid size on the mean streamwise velocity pro®le.

shown in Fig. 27. Fig. 28 shows a comparison of the mean Reynolds stress pro®les. The biggest di€erence is at the peak of the Reynolds stress pro®le, where the LES solution with no SGS model predicts a slightly higher peak than the DNS, and the LES solution with the SGS model predicts a slightly lower peak than the DNS. Note that the LES turbulence statistics were computed for the resolved velocity ®eld only. As a results, the LES Reynolds stress pro®le with the SGS model should be lower than the DNS solution because the SGS contribution was not included. These results indicate that an SGS model is not needed for an accurate simulation of this test case. As discussed earlier, the grid resolution in the near wall region has to be very ®ne to resolve the small energy-producing structures there. So that to resolve the near wall region, the ®ne grid LES has e€ectively approached the DNS in the near wall limit.

10. Conclusion A new, parallel, ®nite-volume computational ¯uid dynamics algorithm was developed for large-eddy simulation (LES) of turbulent ¯ows using parallel computer systems. Major components of the algorithm included piecewise linear least-square reconstruction of the unknown variables, trilinear ®nite-element interpolation for the spatial coordinates, Roe ¯uxdi€erence splitting (FDS), and second-order MacCormack explicit time marching. A systematic and consistent means of evaluating the surface and volume integrals of the control volume was described. Although the test cases were computed using structured meshes, this algorithm could be used for unstructured meshes as well. The parallel implementation was accomplished using the message-passing programming model.

910

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 22. E€ect of grid size on the turbulent intensities.

For the ®rst time, a parallel, ®nite-volume numerical algorithm with Roe FDS was used for LES of turbulent ¯ow in a square duct, and several conclusions have been drawn regarding the accuracy and eciency of this numerical algorithm. Comparison with the direct numerical simulation (DNS) solution showed that the piecewise linear least-square reconstruction algorithm and the MacCormack time marching method work well. However, the standard Roe FDS upwind dissipation adversely a€ects the accuracy of the turbulence simulations. A

T.T. Bui / Computers & Fluids 29 (2000) 877±915

911

Fig. 23. E€ect of grid size on the mean Reynolds stress pro®le.

modi®cation to the standard Roe FDS method was proposed in which the inviscid ¯ux is computed as the arithmetic average of the right and left ¯uxes plus the product of the Roe FDS dissipation term and a reduction factor. For accurate turbulence simulations, only 3±5% of the normal Roe FDS dissipation was found to be needed. The ®ner, 257  100  100 LES grid required less Roe FDS upwind dissipation for stability

Fig. 24. Mean secondary velocity vectors from LES with a 129  90  90 grid.

912

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 25. Streamwise-averaged secondary velocity vectors from the DNS solution.

and produced a more accurate solution than the 129  90  90 LES grid. The near wall vortical structures were better simulated by the ®ner grid LES, and the e€ect of the subgridscale model on the accuracy of the results was found to be small for the ®ne grid LES, which is nearly as ®ne as the DNS grid in the near wall region.

Fig. 26. E€ect of the SGS model on the mean streamwise velocity pro®le.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 27. E€ect of the SGS model on the turbulent intensities.

913

914

T.T. Bui / Computers & Fluids 29 (2000) 877±915

Fig. 28. E€ect of the SGS model on the mean Reynolds stress pro®le.

Acknowledgements The computational resources were provided by the High Performance Computing and Communications Program (HPCCP) /Computational Aerosciences (CAS) at the NASA Ames Research Center, the HPCCP/Earth and Space Sciences (ESS) at the NASA Goddard Space Flight Center and the NASA Jet Propulsion Laboratory, and the NASA Glenn Research Center. The author would like to thank Dr. Spyros Gavrilakis of the Ecole Polytechnique Federale de Lausanne, Switzerland, for graciously providing the DNS solution used in this study and giving helpful insights into the ¯ow physics of the turbulent ¯ow in a square duct. The author also thanks Dr. Tom Lund of the University of Texas at Arlington for his helpful advice on the topic of LES. Finally, the author would like to thank Prof. Robert MacCormack of Stanford University for his advice on the numerical algorithm development and his support as well as encouragement throughout this research project.

References [1] Barth Timothy J. Recent developments in high order k-exact reconstruction on unstructured meshes, AIAA-930668, January 1993. [2] Steger JL, Warming RF. Flux vector splitting of the inviscid gas dynamic equations with application to ®nitedi€erence method. Journal of Computational Physics 1981;40:263±93. [3] Roe PL. Approximate Riemann solvers, parameter vectors, and di€erence schemes. Journal of Computational Physics 1981; 43:357±72. [4] Rai MM, Moin P. Direct numerical simulation of transition and turbulence in a spatially evolving boundary layer. Journal Computational Physics 1993;109:169±92. [5] Van Leer B, Thomas JL, Roe PL, Newsome RW. A comparison of numerical ¯ux formulas for the Euler and Navier±Stokes equations, AIAA Paper 87±1104, 1987.

T.T. Bui / Computers & Fluids 29 (2000) 877±915

915

[6] MacCormack RW. The e€ect of viscosity in hypervelocity impact cratering, AIAA Paper 69±354. [7] Bui Trong Tri. A parallel ®nite volume algorithm for large-eddy simulation of turbulent ¯ows, Ph.D. Thesis, Department of Aeronautics and Astronautics, Stanford University, CA, March 1998. [8] Bui TT, Mankbadi RR. Numerical simulation of acoustic waves interacting with a shock wave in a quasi-1D convergent-divergent nozzle using an unstructured ®nite volume algorithm. International Journal of Computational Fluid Dynamics 1998;10(4):281±90. [9] Bui Trong T. Numerical simulation of 3D low mach number viscous duct ¯ows using an explicit method on massively parallel computer systems, AIAA-97-0335, January 1997. [10] Coirier William John, An adaptively-re®ned, cartesian, cell-based scheme for the Euler and Navier±Stokes equations, NASA TM-106754, 1994. [11] Hughes Thomas JR. The ®nite element method: linear static and dynamic ®nite element analysis. Englewood Cli€s, NJ: Prentice-Hall, 1987. [12] Greenberg Michael D. Foundations of applied mathematics. Englewood Cli€s, NJ: Prentice-Hall, 1978. [13] Moin P, Squires K, Cabot W, Lee S. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A 1991;3:2746±57. [14] Vreman AW, Geurts BJ, Kuerten JGM. Subgrid modelling in LES of compressible ¯ow. In: Voke PR, et al., editors. Direct and large-eddy simulation I. Dordrecht, The Netherlands: Kluwer Academic Publishers, 1994. p. 133±44. [15] Vreman Albertus Willem. Direct and large-eddy simulation of the compressible turbulent mixing layer, Ph.D. Thesis, Universiteit Twente, The Netherlands, December 1995. [16] Smagorinsky J. General circulation experiments with the primitive equations. Monthly Weather Review 1963;91:99±164. [17] Germano M, Piomelli U, Moin P, Cabot WH. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A 1991;3:1760±5. [18] Ghosal S, Lund TS, Moin P, Akselvoll K. A dynamic localization model for large-eddy simulation of turbulent ¯ows. Journal of Fluid Mechanics 1995;286:229±55. [19] Breuer M, Rodi W. Large-eddy simulation of turbulent ¯ow through a straight square duct and a 180 degree bend. In: Proc. of the First ERCOFTAC Workshop on DNS and LES, Guildford, Surrey, UK. Dordrecht, The Netherlands: Kluwer Academic Publishers, 1994. [20] Erlebacher G, Hussaini MY, Speziale CG, Zang TA. Toward the large-eddy simulation of compressible turbulent ¯ows. Journal of Fluid Mechanics 1992;238:155±85. [21] Squires Kyle, Zeman Otto. On the subgrid-scale modeling of compressible turbulence, Proceedings of the 1990 Summer Program. NASA/Stanford Center for Turbulence Research, December 1990, pp. 47±59. [22] Yoshizawa A. Statistical theory for compressible turbulent shear ¯ows, th the application to subgrid modeling. Physics of Fluids A 1986;29:2152±64. [23] Squires Kyle D. Dynamic subgrid-scale modeling of compressible turbulence, Annual Research Briefs, Stanford University, Center for Turbulence Research, December 1991, pp. 207±223. [24] Mompean G, Gavrilakis S, Machiels L, Deville MO. On predicting the turbulence- induced secondary ¯ows using nonlinear k±e models. Physics of Fluids 1996;8(7):1856±68. [25] Gavrilakis S. Numerical simulation of low-Reynolds-number turbulent ¯ow through a straight square duct. Journal of Fluid Mechanics 1992;244:101±29. [26] Coleman GN, Kim J, Moser RD. A numerical study of turbulent supersonic isothermal-wall channel ¯ow. Journal of Fluid Mechanics 1995;305:159±83. [27] Madabhushi RK, Vanka SP. Large eddy simulation of turbulence-driven secondary ¯ow in a square duct. Physics of Fluids 1991;A 3:2734±45. [28] Lin San-Yih, Yu-Fene Chen, Sheng-Chang Shih. Numerical Study of MUSCL Schemes for Computational Aeroacoustics, AIAA-97-0023, January 1997. [29] Wang WP, Pletcher RH. Evaluation of some coupled algorithms for large eddy simulation of turbulent ¯ow using a dynamic SGS model, AIAA-95-2244, June 1995.