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
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο