High Accurate Finite Volume Method for Large Eddy

fingersfieldMécanique

22 févr. 2014 (il y a 3 années et 1 mois)

65 vue(s)




1

High Accurate

F
inite
V
olume
Method
for L
arge
E
ddy
S
imulation of Complex Turbulent Flows

Lan Xu

a,
d

, Guixiang Cui

b
*
, Chunxiao Xu
b
,
Zhishi Wang

c
,

Zhaoshun Zhang

b

and
Naixiang Chen

a

a

Department of thermal energy engineering, Tsinghua university, Beijin
g, 100084, China

b

Department of Engineering Mechanics, Tsinghua university, Beijing, 100084, China

c

Department of Civil and Environmental Engineering, Faculty of Science and Technology,

University of Macau , Macao ,

China

d

College of Science, Donghua U
niversity, Shanghai
,

200051, China


______________________________________________________________________

Abstract

This paper proposes a
finite volume
method with
compact

fourth order accuracy scheme
for large
eddy simulation (LES)
. T
wo
-
dimensional lid
-
dr
iven

cavity

flow

and
a
flow

over
a
n oscillating
plate

are used as examples to
verify

both t
he accuracy

and
convenience

of the proposed scheme. A

turbulent channel
flow

and a turbulent

flow
over a

back
ward facing
step

are numerically tested for
its effectiv
eness
by LES with dynamic Smagorinsky subgrid model.
I
mmersed boundary method
(IBM) is

applied in this paper t
o deal with flows with
complex
configuration

so that
the boundary
condition on the rigid wall

can be
satisfied

well
.

A
curved

channel

flow

and
a f
low
around a
n airfoil

of NACA0012

are computed with
immersed boundary method
, and
comparison

with experimental
data is also made, showing that t
he higher
accurate

finite volume method
for LES
is proved to be

a
promising numerical method.


Keywords
:
l
arge
eddy simulation,
f
inite volume method, compact scheme,
immersed

boundary

method

_______________________________________
_______________________________
_________


1.
Introduction


Turbulence is the most difficult problem in the fluid mechanics, and it has n
ot yet
fully solved, it appears everywhere even in nano flows[
1
-
3
].
Large eddy simulation is
believed to be a potential prediction method for complex turbulent flows in foreseeable
future
[
4
]
, h
owever
,

a
number of
problems are still needed to be further
res
olved before
it can be used in practice. On the physical aspect the subgrid model is the
most
significant issue in LES and great efforts have been made in the development of LES
since it was proposed in 1970’s. A number of subgrid models are available, amo
ng
which

the dynamic Smagroinsky model is considered to be
the
best

one
,

which will be
used in this paper
. Wall model for subgrid stress is another important issue in practical



*
Corresponding author. Address:

Department of
Engineering Mechanic
s
, Tsinghua university, Beiji
ng
,

100084, China
.

Tel.: +86
-
10
-
62785551; Fax: +86
-
10
-
627
1824.

E
-
mail addresses
:
cgx@mail.tsinghua.edu.cn

(
Guixiang Cui
).




2

use of LES for complex wall
-
bounded turbulent
flows

[
5
]
. A number of wall model
s
have been proposed, e.g.
,

equilibrium layer model
,

two
-
layer model and
others
.
Although the wall model is not

a

fully resolved issue in LES
,
the equilibrium layer
model based on the logarithm law in near
-
wall turbulence is easy to use in numerical
calcul
ation and
will be used
in the paper.

On

contrast
to the physical issues the numerical issues are more serious in the
practical development of LES. Since LES is a
n

unsteady solver it needs powerful
computer for long time calculation to obtain reliable stati
stics. Both accuracy and
efficiency

of the

numerical method are the
major
key
s

to
practical use of LES.
Comparing with other
various numerical solution
s,
the finite volume method (FVM)

[6,7]

is
more
robust and has
better

conservative property of mass, mome
ntum and
ene
rgy transfer and also it is eas
ier

to accommodate the wall subgrid model. However
the conventional FVM, e.g.
,
SIMPLE with up
-
wind scheme, is in low accuracy
and
is
not suitable for LES. As we know that the subgrid stress is proportional to the
square of
filter length


, which is in the same order of magnitude as the grid mesh length
[
8
]
.
When

a

low
er

accuracy numerical scheme
is used
,

the numerical errors may exceed the
subgrid
stress and one can not obtain a real LES. To keep the advantage of FV
M we try
to use higher order accuracy scheme in the interpolation between volume and surface
average
s

of flow variables. The
Padé
finite volume
compact
method

[
9
], which is
successfully
applied
in numerical computation of two dimensional Navier
-
Stokes
equa
tion

[
10
], is
extended to

LES of three dimensional unsteady flow.
To effectivel
y
resolve

wall
-
bounded turbulence
,

the non
-
uniform grids is needed in the wall region and
we
deduce an
higher order accuracy FVM compact scheme
on

non
-
uniform grids.

The
propose
d method is tested in
a three
-
dimensional lid
-
driven

cavity

flow
, a
turbulent
channel

flow

and a turbulent flow over a backward facing step

with satisfaction.

To deal
with flows with
complex
configuration the immersed boundary method (IBM) is used to
satis
fy the boundary condition on the rigid wall. A
curved

channel

flow

and
a flow
around a
n airfoil

of NACA0012

are computed
successfully

with
immersed boundary
method.

The paper
is

arranged as follows. In section 2
a high accurate scheme is
formulate
d,

and
numerical verification is made; I
n section 3

the accuracy and
effectives

of the
proposed scheme are numerically illustrated;

A
pplications of the numerical method to
the simulation of
complex turbulent flows are
presen
ted

in section 4.

Conclusions

are

summ
a
rized

in section 5.


2. Formulation of the numerical method


2
.
1
.

Governing equation
s

of LES



The governing equation
s

of large eddy simulation
can be

obtained by filtering
Navier
-
Stokes equation
which can be
written as:




3


( )
0
i
i
u
y





(
1
)


2
( ) ( )
( )
i j ij
i i
j i j j j
u u
u u
p
t y y y y y
 


 
 

    
     

(2)

in which variables with upper bar stand for the filtered quantities or so
-
called resolved
scale turbulent quantities and
ij i i j
j
u u u u

 
is
the
subgrid stress. We apply eddy
vi
scosity type model in the paper
:




2
1
2,
3
ij t ij kk ij t s
S C S
   
   

(3)

in which


1
2
1
2,
2
j
i
ij ij ij
j i
u
u
S S S S
y y
 


  
 
 
 
 

and the eddy viscosity will be
determined by dynamical procedure
[
11
]
.


2.
2
.

The
spatial
di
scretization of

the

equation
s


The conservative law
s

of mass and momentum

for

an element of incompressible fluid
can be written
as follows


( ) ( ) ( ) 0,
e w n s t b
C C y z C C x z C C x y
           

(4)


( ) ( ) ( )
( )[ ]
e e w w n n s s t t b b
P
t e w n s t b
x y z C C y z C C x z C C x y
t
y z x z x y S
x x y y z z

     
     
 

 
              
 

 
 
     
   
             
 
   
     
   
 

(5)

where the Cartesian coordinates are used with
x
-
axis in west
-
east direction,
y
-
axis
in

south
-
north direction and
z
-
axis
in

bottom
-
top direction.
( ),( ),
e e w w
C u C u
 
 

( ),( ),( ),( )
n n s s t t b b
C v C v C w C w
   
   

are the mass fluxes on the element surfaces

,


stands for the velocity components and


is the fluid viscosity. The subscripts
e, w, n, s,
t, b

in equation (
5
) denote the surface of the element and
P

is the center of the element
as illustrated in Figure 1 and
S
is the external force terms in momentum

equation, such
as the pressure gradient.





4



(a) the control element

(b) the grids

Figure 1 Illustration for the non
-
staggered grids used in computation


We use non
-
staggered grids in computation since it is easy in programming. To
avoid the dec
oupling between velocity and pressure in non
-
staggered grids we use
momentum interpolation with fourth order compact scheme. The accuracy of FVM
depends on the discretization formula used in connection between primary variables at
grid points W, E, P and
f
luxes on element surfaces e, w

as shown in Figure 1 (b). We use
the fourth order Padé type compact scheme proposed by

Kobayashi

(1999) as the
discretization formula. The one dimensional formulae for interior points in
x

direction
can be written as follows



4
/2/2
( )
i i i i
yz yz yz xyz xyz
i h i h i h i h
a b c d O h
       
 
    

(
6
)

with


2 2
1 1
2 2
1
2 3
1 1 1
2 3
1 1
/( )
/( )
2( 2 )/( )
2( 2 )/( )
i i i i
i i i i
i i i i i i
i i i i i i
a h h h
b h h h
c h h h h h
d h h h h h
 

  
 

 

 


  


  



in which
/2/2/2
/2/2/2
(,,)
x y z
xyz
x y z
x y z d d d x y z
     
  
  
      
  

is the average of


on
the volume, i.e. the variable at center of the element,
/2/2
/2/2
(,,)
y z
yz
y z
x y z d d y z
    
 
 
    
 

is the average of


on the element surfaces, i.e.
the variables at center of surfaces e and w;
1
( 1,2,,1)
i i i
h x x i N

   
are
mesh
lengths

of the grids;



x

is a shif
t operator defined as follows






s
f x f x s

 

(
7
)

Similar formula can be derived for the boundary points as follows


1 2 3/2 5/2 7/2
4
1 1 1 1
( )
yz yz xyz xyz xyz
a b c d O h
    
    

(
8
)


1 1/2 3/2 5/2
4
( )
N N N N N
yz yz xyz xyz xyz
N N N N
a b c d O h
    
   
    

(
9
)

with




5


1 2 1 2 3
1
2 3 2
2 2
2 2 3 1 2 1 3 1
1
1 2 1 2 3
2 2 2
1 1 2 3 2 2 3 1 2 3 1 3
1
2
2 1 2 2 3 1 2 3
2
1 1 2
1
2
2 3 1 2 3
( )( )
( )
2 2 6 3 4
( )( )
( 2 )(2 2 2 )
( )( ) ( )
( )
( ) ( )
h h h h h
a
h h h
h h h h h h h h
b
h h h h h
h h h h h h h h h h h h
c
h h h h h h h h
h h h
d
h h h h h
  






   


  


     



   

 



  




1 2 1 2 3
2 3 2
2 2
2 2 3 1 2 1 3 1
1 2 1 2 3
2 2 2
1 1 2 3 2 2 3 1 2 3 1 3
2 1 2
( )( )
( )
2 2 6 3 4
( )( )
( 2 )(2 2 2 )
(
N N N N N
N
N N N
N N N N N N N N
N
N N N N N
N N N N N N N N N N N N
N
N N N
h h h h h
a
h h h
h h h h h h h h
b
h h h h h
h h h h h h h h h h h h
c
h h h
    
  
       
    
           
  
  


   

  
     


2
2 3 1 2 3
2
1 1 2
2
2 3 1 2 3
)( ) ( )
( )
( ) ( )
N N N N N
N N N
N
N N N N N
h h h h h
h h h
d
h h h h h
    
  
    









  

 



  



When
( 1,2,,1)
i
h x i N
   

t
he formulas
on

uniform grids are as
follows


4
/2/2
1 1 3 3
( )
4 4 4 4
yz yz yz xyz xyz
x x x x
O h
       
   
    

(
10
)


4
1 2 3/2 5/2 7/2
3 (17 8 )/6 ( ),
yz yz xyz xyz xyz
O h
    
    

(
11
)


4
1 1/2 3/2 5/2
3 (17 8 )/6 ( )
yz yz xyz xyz xyz
N N N N N
O h
    
   
    

(
12
)

The derivatives of variable


can also be discretized with fourth order accuracy, for
instance in interior points


4
/2/2
( )
i i i i
yz yz yz
xyz xyz
i h i h i h i h
a b c d O h
x x x
  
     
 
  
    
  

(
13
)

with


2 2
1 1 1
2 2
1 1 1
2 2
1 1
2 2
1 1 1
1
2 2
1 1 1
1
2 2
1 1 1
( )
( )( 3 )
( )
( )( 3 )
12
( )( 3 )
12
( )( 3 )
i i i i i
i
i i i i i i
i i i i i
i
i i i i i i
i i
i
i i i i i i
i i
i
i i i i i i
h h hh h
a
h h h hh h
h h hh h
b
h h h hh h
hh
c
h h h hh h
hh
d
h h h hh h
  
  
 
  

  

  

 


  


 


  





  





  



and on boundary




6


1 3/2 5/2 7/2
1
4
1 1 1 1 1
2
( )
yz yz
yz xyz xyz xyz
a b c d e O h
x x
 
   
 
     
 

(
14
)


1/2 3/2 5/2
1
4
( )
N N N N
N N
yz yz
yz xyz xyz xyz
N N N N N
a b c d e O h
x x
 
   
  

 
     
 

(
15
)

with


1 2 1 2 3
1
2 3 2
2 2 2
1 1 2 1 2 1 3 1 3 2 2 2 3 2 3
1
1 1 2 1 2 3
4 3 3 3 3 2 2 2 2
1 1 2 1 2 1 3 1 3 1 2 1 2
2 2 2 2 2
1 2 3 1 2 3 1 3 1 3
1
2( )( )
( )
2( 3 4 2 2 )
( )( )
6 8 16 4 8 9 15
9 15 3 3
2
h h h h h
a
h h h
h h h ah h ah h h h ah h ah h h h
b
h h h h h h
h ah h h h ah h h h ah h h h
ah h h h h h ah h h h
c
  


        

  
      
   
 
2 3 2
1 2 1 2 3
2 4 4 3 3 2 2 2 2
1 2 3 2 2 2 3 2 3 2 3 2 3
2 2
1 1 2 1 2 3
2 2 2 2
1 2 3 2 3 1 3 1 2 2 3 1 2 3
3 2 3 3 3 3 3 2
3 1 3 2 2 1 3 1 1 2
1
1 1
6 9
3 2 2
( ) ( )
9 8 3 9 12 9
3 4 8 2 6
9
2
h h h h h
h h h ah h ah h h h ah h h h
h h h h h h
h h h ah h ah h ah h ah h ah h h
h h h h ah ah ah h h h
h
d h
 
 
 
 
 
      
 
  
    
       

 
2 2 2 2
2 1 3 3 2 2 3
2 2
2 3 1 2 1 2 3
1 1 1 2 2
1
2
2 3 1 2 3
3 4 6
( )( ) ( )
2 ( 2 )
( )( )
h h h h h h h
h h h h h h h
h ah h ah h
e
h h h h h














 

 

 

 
  

 

   


  


  




1 2 1 2 3
2 3 2
2
1 1 2 3 2 2 3
1 1 2 1 2 3
4 3
1 2 3 1 2 3
2
1 2
2( )( )
( )
3 ( 2) (2 ) ( 1) ( )
2
( )( )
6 4( 2)(2 ) (9 15)( )
2
N N N N N
N
N N N
N N N N N N N N N
N
N N N N N N
N N N N N N N N
N N
N
h h h h h
a
h h h
h a h h h a h h h
b
h h h h h h
h a h h h a h h
h h
c
    
  
      
     
     
 
  


      
 
  
      


2 2 2
1 3 1 2 2 2 3
2 2 2
3 2 2 3
2 2
1 1 2 1 2 3
2 2 3
1 2 2 3 3 2
3 2 2
3 2 3 2
1
3( 1) 3 (2 3
) ( 1) ( )
( ) ( )
3(1 ) (3 3 ) (1 2 )(4
4 6
2
N N N N N N N N
N N N N N
N N N N N N
N N N N N N N N
N N N N
N N
a h h h h h h h
h a h h h
h h h h h h
a h h h h h a h
h h h h
d h
      
   
     
     
   

 
 
  
 
 
   
 
  
    
  

3
3 1
2
2 3 1
2 2
2 3 1 2 1 2 3
1 1 1 2 2
2
2 3 1 2 3
) ( 1)
3(2 )
( )( ) ( )
2 ( 2 )
( )( )
N N N
N N N
N N N N N N N
N N N N N N N
N
N N N N N
h a h
h h h
h h h h h h h
h a h h a h h
e
h h h h h
 
  
      
    
    














 

 

 
 

 
 

 

   


  
 

  






7

T
he formulas
on

uniform grids are as
follows


4
/2/2
1 1 6 6
( )
10 10 5 5
yz yz yz
xyz xyz
x x x x
O h
x x x x x
  
     
   
  
     
    

(
1
6
)


4
1 3/2 5/2 7/2
1 2
1 5 89 127 4
6 ( )
3 18 18 9
yz yz
yz xyz xyz xyz
O h
x x x
 
   
 
 
      
 
  
 

(
1
7
)


4
1/2 3/2 5/2
1
1 5 89 127 4
6 ( )
3 18 18 9
yz yz
yz xyz xyz xyz
N N N N
N N
O h
x x x
 
   
  

 
 
     
 
  
 

(
1
8
)

For brevity all interpolation formulas can be writ
ten in matrix form as



C yz C xyz
x x
A B
 


(
1
9
)


yz
D D xyz
x x
A B
x






(
20
)

Inserting the compact interpolation into equation (4) and (
5
) we have the
semi
-
discrete form of equation (
2
).


P P E E W W N N S S T T B B
P
x y z a a a a a a a b
t

      

 
           
 

 

(
21
)

in

which
(,,,,,....), =,,,,,,
D D C C D D
i i x x x x y y
a a A B A B A B i P E W N S T B


an
d

b

is a source term
as in usual FVM
.



2.
3
.

Time advancement


Fourth order Runge
-
Kutta integration is used in time marching, the governing
equation is written as


( )( )
C D Gp
t



   


(
22
)

in which
C
and
D

are the convective and diffusive operators r
espectively and discretized
by the compact scheme as stated above.
Gp

stands for the
pressure gradient
which is
solved by coupling the momentum equation with continuity equation as follows.


1 1
1 2 3 4
1
( 2 2 )
6
0
n n n
n
t
k k k k tGp
M
 

 



     






(
2
3
)

in which
M

stands for the operator

of continuity equation and
( )
( )( )
i
i
k C D

  
,
1,2,3,4
i


with


n



)
1
(

(
2
4
)




8













0
2
)
2
(
)
2
(
1
)
2
(



M
tGp
k
t
n

(
2
5
)













0
2
)
3
(
)
3
(
2
)
3
(



M
tGp
k
t
n

(
2
6
)











0
)
4
(
)
4
(
3
)
4
(



M
tGp
tk
n

(
2
7
)


2.
4
.

The solution of pressure


The pressure i
s calculated by continuity equation. To avoid the decoupling between
velocity and pressure
on

non
-
staggered grids we use momentum interpolation with
fourth order accuracy

to determine the velocity at the surface of the FVM cell. For
instance the momentum e
quation in
x

direction
in the second step of
Runge
-
Kutta
integration

is


(2) (2)
u
xyz
u H tGp
 

(
2
8
)

in which
( )
2
xyz n n
t
H u C D u

   
,
(2)
( )
yz yz
u u x
Gp d p p


 
,

d
u

is geometry
parameter for FVM cell.

From equation
(
16
)

the momentum interpolation can be ob
tained as
:


1
[ ]
yz C C xyz
x x
H A B H



(
2
9
)

and velocity on the cell surface can be calculated by the following equation



(2)
u
( )
yz xyz xyz
face x
u H d p p t


   

30

After the solution of cell surface velocity the pressure will be computed by strong
implicit procedure (SIP
)

[
12
]
.


2.5.
Immersed boundary method


The immersed boundary method (IBM) was originally introduced in the pioneering
work of Peskin

(
1972)

[
13
]
. The basic idea of IBM lies in the definition of solid (either
moving or
stationary
) boundaries. Traditionall
y most computational methods use
complicated boundary fitted grids to define the geometrical configuration and the flow
regions. The IBM actually
simulate
s the presence of solid bodies by means of suitably
defined body forces applied on the momentum equati
ons. The N
-
S equations allow the
specification of such a body force
f
i
, inserted as a source term i.e.,




9


2
( ) ( )
( )
i j ij
i i
i
j i j j j
u u
u u
p
f
t y y y y y
 


 
 

     
     

(
3
1
)

The body force is computed at every time step, so that the velocity field on an
arbitrary surface is driven to a spec
ified
value
b
V
. In general, that surface can move and
does not necessarily have to coincide with a grid line. As shown by Mohd
-
Yusof (1997)

[
14
]
, considering the time
-
discretized version of Eq. (
3
1
), one can write




1
n n
i i i
u u t RHS f

  

(
3
2
)

where the term RHS contains all the pressure gradient, advection and diffusion terms. In
order to drive the velocity field at the next time step
1
n
i
u

, to the desired level
bi
V
, it is
sufficien
t to formulate the source term f of the Eq. (
3
2
) as


n
bi i
i
V u
f RHS
t

  


(
3
3
)

which is imposed appropriately in the discretized form of the conservation equations
through the boundary conditions of the the solid walls. In general the
boundary

of t
he
region where we want
1
n
i bi
u V



is

not coincide with a coordinate line. So the value of
f
i

on the cell node closest to the surface but outside the body is linearly interpolated.


3.

The check of numerical accuracy



3
.
1
.
T
wo
-
dimensional

lid
-
driven

cavity

flow with gravity


The real accuracy of proposed method is checked by a cavity flow
with gravity at

Re
=
1
. The flow is driven by upper cover with particular motion and
an analytical
solution can be obtained for
the flow field inside cavit
y. The details of the flow case can
be
found in

Pereira (2001)
[
10
]
. The accuracy of proposed method is presented in Table
1 and 2 with
non
-
uniform and uniform grids respectively
.

The errors are the maximum
deviations between the numerical solution and anal
ytical results and the accuracy is
estimated by the exponent of power law.
It is clear that the proposed method is of fourth
order accuracy approximately

and error is smaller on non
-
uniform grids than
that
on
uniform grids

with same size
.


Table 1

Errors a
nd accuracy of proposed method with non
-
uniform grids

Grid
points

Longitudinal velocity

Vertical velocity

Pressure

errors

accuracy

errors

accuracy

errors

accuracy




10

8×8

16×16

32×32

64×64

3.24E
-
3

2.59E
-
4

2.18E
-
5

1.60E
-
6



3.54

3.45

3.70

7.32E
-
3

6.13E
-
4

4.3
7E
-
5

3.52E
-
6



3.45

3.75

3.54

2.95E
-
2

2.51E
-
3

2.12E
-
4

1.63E
-
5



3.44

3.44

3.60


Table 2

Errors and accuracy of proposed method with uniform grids

Grid
points

Longitudinal velocity

Vertical velocity

Pressure

errors

accuracy

errors

accuracy

errors

accurac
y

8×8

16×16

32×32

64×64

3.43E
-
3

2.48E
-
4

1.86E
-
5

1.43E
-
6



3.79

3.74

3.70

6.71E
-
3

4.72E
-
4

3.46E
-
5

2.61E
-
6



3.83

3.77

3.73

4.55E
-
2

3.43E
-
3

2.74E
-
4

2.21E
-
5



3.73

3.65

3.63


3
.
2
.
Oscillating

plate


A

flow
above
an oscillating plate

is calculated in order t
o
check

t
he accuracy of
the

method

for unsteady flows
.
This
coordinates

are

set up
with

the
x
-
axis along the plate,
and the
y
-
axis normal to it. The velocity of the plate is
prescribed as

cos
u t



and
the

exact

solution can be obtained
f
or this problem
.
T
he

errors between our computational
results and
exact

solution
s

are listed in Table
3
.
A
nd the accuracy is estimated by the
exponent of power law
. It is clear that the method has approximately fourth order
accuracy.


Table
3

Errors and a
ccuracy of proposed method with
different
numerical resolutions

Grid
points

Longitudinal velocity

Vertical velocity

Pressure

errors

accuracy

errors

accuracy

errors

accuracy

25×25

50
×
50

100×100

8.90E
-
3

6
.
21
E
-
4

4.43E
-
5



3.84

3.81

7.41E
-
3

5.43

E
-
4

3.82E
-
5



3.77

3.83

1.11E
-
2

1.
03
E
-
3

8
.
09
E
-
5



3.43

3.67


4
. Applications


4
.1
.

Case 1 Turbulent channel flow


The first testing case is
a
turbulent flow through
channel;

the geometry is illustrated
in Figure 2.
The Reynolds number is
Re 395
u H
 

 

and the half width of channel
H
=1.




11


Figure 2 The geometry of turbulent flow in channel


The streamwise, normal and spanwise directions are
x
,
y
,
z

res
pectively

and
u
,
v
,
w

are the velocity components in correspondent directions. Computational
domain is composed of a rectangular box with the normal height equaling 2, streamwise
l
ength
4


and spanwise width

2

. The periodic conditions are posed in both streamwise
and spanwise directions; and a wall model of power law (
Werner, 1991)

[
1
5
]
is given at
the channel wall such that

2
1
2
1
1
1
1
2 | |/| | 0.5/
| |
1 1
| |,otherwise
2
B
P P P P
B B
w
B
B
B
P
P P
u u A
B B
A u
A
   

 
 











 
   
 

 

   

 
   
 

,



(

3
4
)

with

8.3,1/7
A B
 
.

In equation (
3
4
)

u
P

is the
tangential
velocity component to a wall
at the
first
grid point next to the wall and
P


is the vertical distance of the
first
grid
point to the wall.


4
.1.1
.

Comparison of re
sults between
fourth
-
order
accuracy

and second
-
order

accuracy

The results are shown in Figure 3

in wall units.
Uniform grids are used in t
he
fourth order accuracy scheme with
32
×
64×32

while

finer
non
-
uniform grids with
64×64×64

are adopted in
second
-
order
-
accura
cy

scheme.

The results are compared with
the direct numerical simulation performed by
Moser et al. (1991)

[
1
6
]
with fine grids of
256
3
.

Although the resolution is higher in second order accuracy scheme

near the wall

the mean velocity deviates from th
e DNS results remarkably
; the
turbulence intensities

and Reynolds stress are much lower than the DNS results. This indicates that second
order accuracy scheme has much numerical dissipation.

Note that the first grid point to
the wall is located at
12
y



in the fourth order accuracy scheme while the first grid
point is around
5
y



in the second order accuracy scheme. The results indicate that
the higher order accuracy is more effective than use of higher resolutio
n. Of course
combination of higher accuracy scheme with higher resolution of non
-
uniform grids is
much more effective
. T
he results will be shown below.


y
z
x



12

y
+
u
+
1
0
0
1
0
1
1
0
2
1
0
3
0
5
1
0
1
5
2
0
D
N
S
L
E
S
(
f
o
u
r
t
h
o
r
d
e
r
)
L
E
S
(
s
e
c
o
n
d
o
r
d
e
r
)
y
+
U
r
m
s
,
V
r
m
s
,
W
r
m
s
0
1
0
0
2
0
0
3
0
0
4
0
0
0
0
.
5
1
1
.
5
2
2
.
5
3
U
r
m
s
W
r
m
s
V
r
m
s
D
N
S
L
E
S
(
f
o
u
r
t
h
o
r
d
e
r
)
L
E
S
(
s
e
c
o
n
d
o
r
d
e
r
)
y
+
-
u
'
v
'
0
1
0
0
2
0
0
3
0
0
4
0
0
0
0
.
5
1
D
N
S
L
E
S
(
f
o
u
r
t
h
o
r
d
e
r
)
L
E
S
(
s
e
c
o
n
d
o
r
d
e
r
)


(a) Mean velocity profile


(b)
Turbulence intensities



(c) Reynolds stress

Figure 3

The comparison between fourth order and second order accuracy s
cheme in turbulent channel
flow


(─): DNS Mansour et al [
14
], (□): fourth order accuracy results, (
-

-

or

): second order accuracy
results.

4
.1.2
.

Comparison of results between uniform grids
and
non
-
uniform grids

The results are shown in Figure
4

in wall
unit
s
. The uniform grids are
32×64×32
while the non
-
uniform grids are
32×
48
×32
. The results are compared with the direct
numerical simulation performed by
Moser

et al. (1991)

with fine grids o
f
256
3
.

The first grid above the wall is set up at

12
y



in both uniform and non
-
uniform
grid systems. It can be seen that almost same accuracy is obtained in coarse
non
-
uniform girds and finer uniform grids.

Note that the computing c
ost is reduced by
1/3
when the

non
-
uniform coarse grids
are

used.


y
+
u
+
1
0
0
1
0
1
1
0
2
1
0
3
0
5
1
0
1
5
2
0
D
N
S
L
E
S
(
n
o
n
u
n
i
f
o
r
m
)
L
E
S
(
u
n
i
f
o
r
m
)
y
+
U
r
m
s
,
V
r
m
s
,
W
r
m
s
0
1
0
0
2
0
0
3
0
0
4
0
0
0
0
.
5
1
1
.
5
2
2
.
5
3
U
r
m
s
W
r
m
s
V
r
m
s
D
N
S
L
E
S
(
n
o
n
u
n
i
f
o
r
m
)
L
E
S
(
u
n
i
f
o
r
m
)
y
+
-
u
'
v
'
0
1
0
0
2
0
0
3
0
0
4
0
0
0
0
.
5
1
D
N
S
L
E
S
(
n
o
n
u
n
i
f
o
r
m
)
L
E
S
(
u
n
i
f
o
r
m
)


(a) Mean velocity profile (b)
Turbulence intensities

(c) Reynolds stress

Figure 4 The comparison between uniform grids and non
-
uniform grids in turbulent channel flow

(─): DNS Mansour et al [
14
]
,
(

): fourth order accuracy results
,

(
-

-

or

): second order accuracy
results
.


4
.2
.

Case 2 Turbulent flow over back
-
step


To examine proposed method in complex turbulent flows, e.g. flow with separation,
we take the flow over

back
-
step as
a
testing case.

The

flow geometry is shown in Figure
5

and the Reynolds number equals 5100
(
0
U H

), in which
U
0

is free stream speed and
H

is the height of the step. The streamwise,
normal and spanwise directions are
x
,
y
,
z

respectively and
u
,
v
,
w

are



13

correspondent velocity components
.
The computational domain is com
posed of a
rectangular box with streamwise length of
30
H
, spanwise width of
4
H

and normal
height of
6
H
. The inlet velocity condition is composed of a mean velocity profile plus
random fluctuations with prescribed intensit
ies

profiles. The length of inflow
to
back
-
step is 10
h

in order to have sufficient length
in

the development of turbulence.
Periodic condition is posed in spanwise direction and non
-
reflection condition is posed
at outlet and upper boundary. Wall function with power law
[
1
5
]

is given at the
rigid
wall as in the case of turbulent channel flow.

Uniform grids are used in spanwise
direction while non
-
uniform grids are prescribed in normal and streamwise direction
and the grids near rigid wall are

shown in Figure
5
b, the total grids are 172×83×32
in
streamwise, normal and spanwise directions respectively.


After the flow reaches statistically stationary the computation is continued to take
200 samples. The time interval of the samples is 0.1
H/U
0
. The statistics are taken
average over samples and s
panwise direction. The results are compared with the DNS
by Hung et al (1997)

[1
7
] with grid points 768×192×64.



X
/
H
Y
/
H
-
1
0
1
2
3
0
0
.
5
1
1
.
5
2
2
.
5
3
3
.
5


(a) computational domain


(b) local grid design

Figure 5 The illustration of the geometry for flow over b
ack
-
step and local grid design

The mean streamlines are shown in Figure
6

and the
mean
reattach
ment

length is
estimated as
6.
4
3
H

by the location at which the mean velocity
U
=0 at the first grid point
away from the wall
which is in good agreement with DNS r
esult of
6.28
.


Figure 7 presents the comparison between computational LES results and DNS
results for t
he mean
streamwise
velocity profiles

at three representative locations in the
recirculation, reattachment and recovery regions
with
X/
H

=
4
,
6

and
15
. G
ood
agreement

between computational LES and DNS results is obtained at all locations.



X
/
H
Y
/
H
0
5
1
0
1
5
0
5
1
0
1
5

Figure
6

The mean streamlines with reattached length
Xr
=6.
4
3
H





14


U
A
/
U
0
Y
/
H
-
0
.
2
5
0
0
.
2
5
0
.
5
0
.
7
5
1
1
.
2
5
0
0
.
5
1
1
.
5
2
2
.
5
3
3
.
5
4
D
N
S
L
E
S
X
/
H
=
4
U
A
/
U
0
Y
/
H
-
0
.
2
5
0
0
.
2
5
0
.
5
0
.
7
5
1
1
.
2
5
0
0
.
5
1
1
.
5
2
2
.
5
3
3
.
5
4
D
N
S
L
E
S
X
/
H
=
6
U
A
/
U
0
Y
/
H
-
0
.
2
5
0
0
.
2
5
0
.
5
0
.
7
5
1
1
.
2
5
0
0
.
5
1
1
.
5
2
2
.
5
3
3
.
5
4
D
N
S
L
E
S
X
/
H
=
1
5


(a)
X/
H
=4


(b)
X/
H
=6


(c)
X/
H
=15

Figure 7 Velocity profiles after
back
-
step
,
(─): DNS Hung et al [
1
7
]
,
(

):
LES

results
.


The
turbulence intensities

profiles normalized by
0
U
and Reynolds stress
es

profile
s

normalized by
2
0
U

are given in Figure
8
, 9
, 10
and 1
1

at
x
-
stations
with
X/H
= 4, 6 and
15
respectively
. All statistics are in fairly good agreement with DNS results.

The

near
wall peak of streamwise turbulence intensit
y

appears to develop earlier in DNS than in
LES in the recover
y

zone

(
X/H

15).

T
he differences

of t
he nor
mal and spanwise

turbulence intensit
ies between DNS and LES
are

noticed mainly at the vicinity of the
wall
.

The difference
s

may be caused by the wall model used in LES and it indicates that
the wall model should be further investigated.

U
r
m
s
/
U
0
Y
/
H
0
0
.
0
5
0
.
1
0
.
1
5
0
.
2
0
.
2
5
0
.
3
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
4
U
r
m
s
/
U
0
Y
/
H
0
0
.
0
5
0
.
1
0
.
1
5
0
.
2
0
.
2
5
0
.
3
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
6
U
r
m
s
/
U
0
Y
/
H
0
0
.
0
5
0
.
1
0
.
1
5
0
.
2
0
.
2
5
0
.
3
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
1
5


(a)
X/H
=4

(b)
X/H
=6 (c)
X/H
=15

Figure 8 Root mean square of streamwise fluctuations
,
(─): DNS Hung et al [
1
7
]
,
(

):
LES

results
.


V
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
4
V
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
6
V
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
1
5


(a)
X/H
=4


(b)
X/H
=6


(c)
X/H
=15

Figure 9 Root mean square of normal fluctuations
,
(─): DNS Hung et al [
1
7
]
,
(

):
LES

results
.





15

W
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
.
2
2
5
0
.
2
5
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
4
W
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
6

W
r
m
s
/
U
0
Y
/
H
0
0
.
0
2
5
0
.
0
5
0
.
0
7
5
0
.
1
0
.
1
2
5
0
.
1
5
0
.
1
7
5
0
.
2
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
1
5


(a)
X/H
=4


(b)
X/H
=6



(c)
X/H
=15

Figure 10 Root mean square of spanwise fluctuations
,
(─): DNS Hung et al [
1
7
]
,
(

):
LES

results
.



u
v
/
(
U
0
*
U
0
)
Y
/
H
-
0
.
0
2
-
0
.
0
1
5
-
0
.
0
1
-
0
.
0
0
5
0
0
.
0
0
5
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
4
u
v
/
(
U
0
*
U
0
)
Y
/
H
-
0
.
0
2
-
0
.
0
1
5
-
0
.
0
1
-
0
.
0
0
5
0
0
.
0
0
5
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
6
u
v
/
(
U
0
*
U
0
)
Y
/
H
-
0
.
0
2
-
0
.
0
1
5
-
0
.
0
1
-
0
.
0
0
5
0
0
.
0
0
5
0
1
2
3
4
D
N
S
L
E
S
X
/
H
=
1
5


(a)
X/H
=4


(b)
X/H
=6


(c)
X/H
=15

Figure 11 Reynolds stress distributions
,
(─): DNS Hung et al [
1
7
]
,
(

):
LES

results
.



4.
3
. Case
3

Curved t
urbulent

channel

flow


A mildly curved turbulent channel flow has been simulated by the above numerical
method

with orthogonal grids
. T
he total grids are
64
×
48
×32
. The curvature parameter
R
h

is

chosen to be 1/79=0.0127, and the Reynolds number is around 2800 based on
mean velocity of the inlet and half width of the channel.

The

flow geometry is shown in
Figure
12. The computational domain is 12.6
h
, 2
h

and 4.2
h

in streamwise, normal and
spanwise
direction
s

respectively

and the dynamical Smagoringsky model is used for the
closure of subgrid stress
.





16


x
/
h
y
/
h
0
1
2
3
0
1
2
3


(a) computational domain (b) local grid design

Figure
12

The illustration of the geom
etry for
curved turbulent channel flow

and local grid design


The computed flow fields are used to study the effects of streamline curvature by
comparing the concave and convex sides of the channel. Observed effects are shown in
Figure 13 and they are cons
istent with the DNS case completed by Moser and Moin
(
1987)

[
1
8
]
.

As expected, turbulence intensities near the concave wall are higher than
those near the convex wall.

In addition, it is found that stationary Taylor
-
Gortler
vortices are present and that th
ey have a significant effect on the flow which is shown in
Figure 14.


y
+
u
+
1
0
0
1
0
1
1
0
2
5
1
0
1
5
2
0
t
o
p
w
a
l
l
(
D
N
S
)
b
o
t
t
o
m
w
a
l
l
(
D
N
S
)
t
o
p
w
a
l
l
(
L
E
S
)
b
o
t
t
o
m
w
a
l
l
(
L
E
S
)
2
.
5
l
n
(
y
+
)
+
5
.
5
y
+
U
r
m
s
0
5
0
1
0
0
1
5
0
0
1
2
3
t
o
p
w
a
l
l
(
D
N
S
)
b
o
t
t
o
m
w
a
l
l
(
D
N
S
)
t
o
p
w
a
l
l
(
L
E
S
)
b
o
t
t
o
m
w
a
l
l
(
L
E
S
)
y
+
V
r
m
s
0
5
0
1
0
0
1
5
0
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
0
.
7
0
.
8
0
.
9
1
t
o
p
w
a
l
l
(
D
N
S
)
b
o
t
t
o
m
w
a
l
l
(
D
N
S
)
t
o
p
w
a
l
l
(
L
E
S
)
b
o
t
t
o
m
w
a
l
l
(
L
E
S
)


(a) Mean velocity profile

(b) Root mean square of streamwise fluctuations

(c) Root mean square of
normal

fluctuations






y
+
W
r
m
s
0
5
0
1
0
0
1
5
0
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
0
.
7
0
.
8
0
.
9
1
1
.
1
1
.
2
t
o
p
w
a
l
l
(
D
N
S
)
b
o
t
t
o
m
w
a
l
l
(
D
N
S
)
t
o
p
w
a
l
l
(
L
E
S
)
b
o
t
t
o
m
w
a
l
l
(
L
E
S
)

y
+
-
u
v
0
5
0
1
0
0
1
5
0
-
0
.
2
-
0
.
1
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
0
.
7
t
o
p
w
a
l
l
(
D
N
S
)
b
o
t
t
o
m
w
a
l
l
(
D
N
S
)
t
o
p
w
a
l
l
(
L
E
S
)
b
o
t
t
o
m
w
a
l
l
(
L
E
S
)


(
d
) Root mean square of
spanwise

fluctuations



(
e
) Reynolds stress

Figure
13

The
results

of curved turbulent channel flow with LES





17

z
y
0
1
2
3
4
0
.
5
1
1
.
5
2
2
.
5
3
3
.
5
4


a

Velocity vector of the Taylor
-
Gortler vortices

Z
+
T
a
o
0
1
0
0
2
0
0
3
0
0
4
0
0
5
0
0
0
.
0
0
3
0
.
0
0
4
0
.
0
0
5
0
.
0
0
6
0
.
0
0
7
0
.
0
0
8
0
.
0
0
9
0
.
0
1
0
.
0
1
1
0
.
0
1
2


b

Variation of the wall shear stress in the spanwise direction (

: concave wall

---
: convex wall)

Figure 1
4

Taylor
-
Gortler

vortices


In Figure 14(b) the spanwise variation of the wall shear stress is sho
wn for both walls.
On the concave wall there is a very sharp minimum in shear stress located between the
vortices.


4.
4
. Case
4

F
low

around an airfoil of NACA0012


A flow
around a
n airfoil

of NACA0012

has been simulated by the proposed
numerical method

wit
h orthogonal grids
.
The NACA 0012 airfoil was simulated
by
proposed method respectively at two different Reynolds numbers

based on
free stream

velocity and chord length
C

of the airfoil
, Re =
8
00 and Re = 660,000.
T
he geometry
and the computational domain

is illustrated in Figure
15
.

The computational domain is
12
c
, 4
c

and 4
c

in streamwise, normal and spanwise direction
s

respectively

and the
dynamical Smagoringsky model is used for the closure of subgrid stress
. T
he total grids
are
132
×
72
×
18.


x
/
c
y
/
c
0
5
-
2
-
1
0
1
2
3
4
5
6
7
8




18

Figure
15

The illustration of the geometry

and

grid design for
flow

around an airfoil



4.
4
.1.
Re=800

The wing is placed in a uniform flow
with attack angle

20

and t
he Reynolds
number is 800. The time evolution of the drag coefficient

(
C
D
)

a
nd the lift coefficient

(
C
L
)

shows the periodicity of the phenomenon shown in Figure 1
6

and

t
he streamlines

at
four time steps with in the
period

are shown in
Figure 1
7
. They are consistent with the
DNS case com
uted

by Hoarau et al
(
2003)

[1
9
]
.

t
C
L
,
C
D
5
1
0
1
5
2
0
2
5
3
0
0
0
.
5
1
C
L
(
L
E
S
)
C
D
(
L
E
S
)
C
L
(
D
N
S
)
C
D
(
D
N
S
)

Figure 1
6

The
results

of curved turbulent channel flow with LES

X
Y
-
0
.
5
0
0
.
5
1
1
.
5
-
1
-
0
.
5
0
0
.
5
1

X
Y
-
0
.
5
0
0
.
5
1
1
.
5
-
1
-
0
.
5
0
0
.
5
1

(a)T/4 (b) T/2




19

X
Y
-
0
.
5
0
0
.
5
1
1
.
5
-
1
-
0
.
5
0
0
.
5
1

X
Y
-
0
.
5
0
0
.
5
1
1
.
5
-
1
-
0
.
5
0
0
.
5
1

(c)3T/4 (d) T

Figure 1
7

T
he streamlines
at four time steps within one

peri
od


4.
4
.
2
.
Re=660,000

The wing is placed in a uniform flow upstream at
various attack angles

and the
Reynolds number is
660,0
00.
Simulations for a range of attack angles from
-
7
o

up to
7

o

were performed. Fig.
18 and Table 4

show the lift

coefficient

and dr
ag

coefficient

at
various attack angles. These are compared with experimental data and with results from
PowerFLOW
[
20
]
. It is known that both

C
L

and

C
D

are Reynolds number dependent as
shown in Fig.
19

(based on available experimental data of closest possi
ble Reynolds
numbers).
I
t can be seen that
the method

gives accurate prediction for all attack angles.



A
n
g
l
e
C
L
-
1
0
-
5
0
5
1
0
-
1
.
2
-
1
-
0
.
8
-
0
.
6
-
0
.
4
-
0
.
2
0
0
.
2
0
.
4
0
.
6
0
.
8
1
1
.
2
C
o
m
p
u
t
e
E
x
p
e
r
i
m
e
n
t
C
L
C
D
-
1
-
0
.
5
0
0
.
5
1
0
.
0
1
0
.
0
1
5
0
.
0
2
0
.
0
2
5
0
.
0
3
C
o
m
p
u
t
e
E
x
p
e
r
i
m
e
n
t

Figure 18 The l
ift
coefficient

at various attack angles

Figure 19 The l
ift
coefficient

and drag
coefficient


Table
4

Summary of results for simula
tions and experiment for flow past a

NACA0012 airfoil



P
roposed method


Experiment

PowerFLOW

d
C

l
C

d
C

l
C

d
C

l
C




20

7



0.01381

-
0.74829

0
.
0
1480

-
0
.7
12
0
1





3


0.01251

-
0.33765

0.0
12
31

-
0.
32499





0

0.01167

-
3.5742E
-
4

0.01
128

-
1.2258
E
-
3

0.0
2010

-
3.518
E
-
3

3

0.01247

0.33584

0.0
12
3
8

0.33
247

0.01
32
7

0.
18
8
28

7

0.01402

0.75472

0.0
150
3

0.7
1
46
49

0.01574

0.7449


In summary, the numerical examples performed in the thesis give strong evidence
that the proposed meth
od can be used to compute complex
turbulent
flows
.


5. Concluding remarks


A fourth
-
order
-
accurate finite volume method is
propos
ed for Large Eddy Simulation.
The Padé type compact scheme is used in discretization for both interior and boundary
grids and t
he accuracy of
the
proposed method is checked by a two dimensional cavity
flow

and a
flow

over
a
n oscillating
plate
. It is evaluated that the proposed method is of
fourth
-
order accuracy approximately
.

The numerical simulations of a
turbulent channel

flow

and a turbulent flow over a
backward facing step are
performed successfully
by the proposed method. R
esults are in
good agreement with DNS data
in coarse girds
.
It is concluded that
the proposed
method
has the advantages of higher accuracy, less grids and

lower
comput
ing

cost
.

To deal with flows with
complex
configuration the immersed boundary method (IBM)
is used to satisfy the boundary condition on the rigid wall. A
curved

channel

flow

and a
flow
around a
n airfoil

of NACA0012
are computed with
immersed b
oundary method.
The

results are consistent with
experiment

and DNS data

and
valid
ate the feasibility

of
immersed boundary method.

It is concluded that the proposed method
is a powerful
approach to precisely simulate the flow fields with complex configurati
ons.


Acknowledgements


The authors would like to thank National Science Foundation of China
(Grant
10
5
720
73

and 50
4
790
06
)

for
finance

support
.


Reference

[1]

H
e
J
h,
W
u
Y
,
Z
uo
W
w.
critical length of straight jet in electrospinning

,
P
olymer 46 (26):
12637
-
12640, 12 2005




21

[2]

H
e
J
h,
W
an
YQ
,
Y
u
J
y.
allometric scaling and instabi
lity in electrospinning

,
I
nternational
journal of nonlinear sciences and numerical simulation, 5 (3): 243
-
252, 2004

[3]

H
e
J
h,
W
u
Y
,
P
an
N
.
a mathematical model f
or ac
-
electrospinning

,
I
nternational journal of
nonlinear sciences and numerical simulation, 6 (3): 243
-
248, 2005

[4]

Piomelli U. Large
-
eddy simulation: achievements and challenges. Prog. Aerosp. Sci. 1999; 35:
335
-
362.

[5]

Piomelli U, and Balaras E. Wall
-
layer
models for large
-
eddy simulations. Annual Review of
Fluid Mechanics 2002; 34: 349
-
374.

[6]

Ma HM, Fan SQ, Chen HP
.
Numerical study of unsteady flow in thrust vector
ing nozzle

,
I
nternational

J
ournal of
T
urbo &
J
et
-
E
ngines

22 (1): 31
-
40
,
2005

[7]

Bigerelle M, Iost A
.
Physical interpretations of the numerical instabilities in
diffusion equations via
statistical thermodynamics
,

I
nternational
J
ournal of
N
onlinear
N
ciences and
N
umerical
S
imulation,
5 (2): 121
-
134, 2004

[8]

Pope S. B. Turbulent flows. Cambridge: the United Kingdom at the University Press; 2000.

[9]

Kobayashi M H. On a cla
ss of Pad´e finite volume methods
.

Journal of Computational Physics
1999
; 156:
137

180.

[10]

Pereira J M C, Kobayashi M H, Pereira I C F. A fourth
-
order
-
accurate finite volume compact
method for the incompressible Navier
-
Stokes solutions . Journal of Computatio
n

Physics 2001;
167: 217

243.

[11]

Germano M, Piomelli U, Moin P. A dynamic subgrid scale eddy viscosity model. Physics of
Fluid A 1991; 3(7): 1760
-
1765.

[12]

Stone H. L. Iterative solution of implicit approximation of multidimensional partial differential
equations.

SIAM J. on Num. Analysis 1968; 5: 530
-
558.

[13]

Peskin

C

S. Flow patterns around heart valves:

A numerical method. Journal of Computational
Physics 1972
;
10
:
252
-
271
.

[14]

Mohd
-
Yusof

J
.

Combined Immersed Boundaries/B
-
Splines Methods for Simulations of Flows
in Comp
lex Geometries
.

CTR Annual Research Briefs
.

NASA Ames/Stanford University, 1997
;

317
-
328
.

[15]

Werner H, Wengle H.
Large
-
eddy simulation of turbulent flow over and around a cube in a plate
channel
.

The Eighth International Symposium on Turbulent Shear Flows
. Be
rlin:
Springer
-
Verlag; 1991,
1941
-
1946
.

[16]

Moser R D, Kim J, Mansour N N. Direct numerical simulation of turbulent channel flow up to
Re
=590. Physics of Fluid 1999
;

11:

943
-
945
.

[17]

Hung Le, Parviz Moin, John Kim.
Direct numerical simulation of turbulent flow ov
er a
backward
-
facing step. Journal of Fluid Mechanics 1997
;

330:

349

374
.

[18]

Moser R D, Moin P. The effects of curvature in wall
-
bounded turbulent flows.

Journal of Fluid
Mechanics

1987
;

175:

479
-
510
.

[19]

Hoarau Y, Faghani D
,
et

al
. Direct
n
umerical
s
imulation of

the
t
hree
-
d
imensional
t
ransition to
t
urbulence in the
i
ncompressible
f
low around a
w
ing.
Flow, Turbulence and Combustion

2003
;

71:

119

132
.

[20]

David P Lockard, Li
-
Shi Luo, Bart A Singer. Evaluation of the Lattice
-
Boltzmann equation
solver PowerFLOW for aerodynamic applications. ICASE Report, 2000
-
40.




22