Below are the first 10 and last 10 pages of uncorrected machineread text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapterrepresentative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 285
Large Eddy Simulation by Using FiniteDifference Method
Y. Doi
Hiroshima University
Hiroshim, Japan
Abstract
Turbulent channel f low is studied
numerically by using Large eddy simulation
(LES). Finite difference method is employed in
the LES. The simulation is stably executed by
using the 3rd order upwind difference scheme
which dissipate numerical errors. First,
a computing test is performed for a laminar
f low past a circular cylinder. After that,
several pilot tests are performed with respect
to a turbulent channel f low, in order to
investigate the ef feet of the numerical
dissipation terms on the turbulent f low
structures. Turbulent f low structures in
turbulent channel flows can be well simulated
by using the 3rd order upwind difference
scheme .
1. Introduction
Large eddy simulation is a relatively new
approach to the calculation of turbulent
f lows. Since the first application of LES was made
by Deardrof f [ 1 ], many attempts on LES have
been made. However the typical LES
calculations are limited to the flows in the
simple conf figurations in order to resolve the
precise structure of the turbulent f lows and
to calculate the turbulent f lows stably. For
the practical engineering applications, the
application for the f lows in or around the
more complicated configurations is necessity.
For the convenience to simulate the f low
around an arbitrary conf iguration for general
purposes, boundary fitted coordinate systems
and regular grid systems are usually
introduced, although they cause numerical
errors which diverge the calculations. So that
there needs some numerical damping.
Some attempts to calculate LES have been
done for the f lows around the complicated
configurations by use of the f inite dif ference
method( [ 2 ], [ 3 ] ). In these methods, the 4th
derivative numerical dissipation terms are
285
introduced artif icially in order to damp
numerical disturbances with short wave length.
The numerical dissipation, which damps
numerical disturbances caused by numerical
truncations, must be distinguished from SGS
model which expresses the subgrid scale
ef feats.
In this paper, several numerical
experiments are performed in order to
investigate the ef feet of numerical
dissipation on turbulent f lows. First,
computing test is performed for a laminar f low
past a circular cylinder. After that, several
pilot tests are performed for a turbulent
paral lel channel f low.
2. Formulation
2.1 Governing Equations for the Large Scale
Field
We consider an incompressible f low whose
time evolution i s given by the NavierStokes
and continuity equations for the velocity
component s Hi; ( i =1, 2, 3 ) and pre s s ure p,
auk a tutus' ap
= _ _
at axe axe
a u k/ a Xk=0
a2u'
+ (l/Re*) (1)
ark axk
(2)
where i( j)=1,2,3 correspond to x1,x2,x3
respectively, where Oxlx2x3 is the Cartesian
coordinate(x1 is in the downstream direction,
X2 is the direction normal to boundary,x3 is
the lateral direction, see Fig.1). Subscripts
denote the partial dif ferentiation with
respect to the referred variables. All
variables and coordinates have been made
dimensionless by means of the length scale D,
distance of parallel channel, and the friction
OCR for page 285
velocity of the parallel channel(u
*
1 ).
2.2 Filtered Momentum and Continuity Equations
In LES, each flow variable "f" is
decomposed as follows,
f=f~f '
(3)
where f is the resolvable grid scale(GS)
component and f' is the residual subgrid
scale(SGS) component. We denote the grid scale
field as,
where ~ ij is Kronecker's delta. K is SGS
eddy coefficient, which is modeled by
Smagorinsky et al.[4] as follows,
K=(c ^)2{aul/axd(au,/ax~
+3uJ/axI)~/2 (10)
c is Smagorinsky's constant, which is
statistically analysed by Yoshizawa[5]. Here c
is taken as 0.1, which is a standard value
used for a parallel channel flow.
Leonard term and Cross term are
approximated as follows by use of Taylor
series expansion,
f (x, x2, x3)=~ G(xlx1 ~x2x2', X3X3')
·f(x1',x2',X3')dx1'dx2'dx3' (4) ~2 a 2(U~Uj) (11)
24 axkaxk
In the present study, TopHat filter is used
as G.
G(xl, X2, X3)= 1/^ 3 ; ~ x,x, ' 1 ~ Axi/2
(A 3=^Xt ~X2 ~X3)
= 0 ; other (5)
After applying the filtering operator to the
momentum equations and the continuity
equation, we get the following equations for
the filtered field.
auk a (gum) ap
= _ +
at axe ax,
a2u;
v ",
+ (lyre*) (6)
~ Xk a Xk
In the present study, following assumption is
applied.
u,u ,=u,u~+ R l a+ L I ~+C I ~ (1)
R,~=u, UJ'
L ~ ~=u~u~—U~uJ
Cl3=UIU] +UI UJ
(8)
The terms Rij, Lij and Cij are SGS Reynolds
Stress term, Leonard term and Cross term which
must be modeled in terms of grid scale
components to proceed the calculation.
2.3 Representation of SGS Reynolds Stresses,
Leonard Term and Cross Term
SGS Reynolds stresses are modeled as
follows introducing SGS eddy viscosity K.
Ut Us ~ ~ Auk Uk /3
=—K(au,/ax,+ au,/ax~) (9)
^2 ague abut
Ctu=   (ut ~ ~ (12)
24 aXkBxk EXkAXk
~2 auk auk
L ~ ~+C ~ a= (13)
12 a Xk a Xk
Substituting eqs.(7),(8) to eq.(6), we get
the following momentum equations.
au, a [Gaul—K(au~/ax,+ au,/ax~)]
_ = _
at axe
a EA2/l2. au~/axk. au~/axk)]
axe
a(~2x~) a2u,
— + (1/Re*) (14)
axe axk axk
F=p+uk us /3 +2x1 (15)
Here the quantity a (2xl)/axl in eq.(14) is
the dimensionless gross downstream pressure
gradient in case of parallel channel flows.
The term uk'uk'/3 is subtracted from the
Reynolds stress terms and added to the
pressure term.
2.4 Wall Damping Function
Near the walls/ Smagorinsky model is not
appropriate because of the inhomogeneous
effect of the wall turbulence. Therefore the
filter length ~ is multiplied by the Van
Driest exponential damping function fD after
MoinKim[6] as follows,
fD=( 1exp(—y+/A+) ) (16)
y+=x2Re* O X0.5
286
OCR for page 285
y+=(lX2)Re* 0.5
Cyclic boundary conditions are imposed for
velocity and pressure, in downstream and
spanwise directions. Noslip boundary
conditions are imposed at walls. Zero gradient
of the pressure normal to the wall is imposed.
4.Computed Results
4.1 Computing Test for Laminar Flow
Before the LES calculation, computing test
is performed for a laminar flow past a
circular cylinder. Although the flow for this
case is 2 dimensional one, the computation is
performed by 3 dimensional code in order to
check the code. In this case, all variables
have been made dimensionless by means of the
cylinder diameter(D) and the uniform flow
velocity(U). The number of grid points are 80
in circumferential direction, 40 in radial and
5 in axial of the cylinder. The diameter of
the computational domain is about 28 times of
that of the circular cylinder. The minimum
spacings in the radial direction, in the
circumferential direction and in the axial
direction are 0.01D, 0.008D and 2.0D,
respectively. In this calculation, Euler
explicit method is employed for the time
derivative. The flow around the circular
cylinder is accelerated from U=0 to the steady
speed(U). The acceleration is set to be
1.0U2/d. The time increment is set to be
0.001U/D when the flow is accelerated, and
after the acceleration the time increment is
set to be 0.0015U/D.
Noslip boundary condition is imposed on
walls. Zero gradient of the pressure normal to
the wall is imposed. On the outer boundary,
the pressure is fixed to zero. With respect to
the velocity on the outer boundary, a linear
extrapolation is given except a upstream
point at which the velocity is given by the
inviscid flow solution. Cyclic boundary
conditions are imposed for the velocity and
pressure in the axial direction.
In this case, the flow is laminar so that
the Smagorinsky's constant 'c' is set to be
zero. The artificial dissipation factor a is
set to be 1.0. The convergence condition of
the Poisson eq.(25) is set as the divergence
of the velocity is less than 10 3.
Thecalculated pressure distribution on the
circular cylinder is shown in Fig.2 together
with the experimental data[7] and the
calculated result by Kodama[8]. Reynolds
number(Re) based on the cylinder diameter and
the uniform flow velocity is 40. The computed
result shows good agreement with the
experimental data. In Fig.3 drag coefficient
is compared with the Tritton's experimental
curved and Kodama's calculated results. The
calculated result agrees well with them.
In Kodama's calculation, 4th order
artificial dissipation terms were also added
to the momentum equations in order to
stabilize the calculation. Kodama analyzed
the truncation errors and indicated that the
artificial dissipation terms were so small
that the accuracy of the solution was not
degraded by them.
Without the 4th order artificial
dissipation terms, the calculation diverges
always. Although the 2nd order artificial
dissipation terms distort the calculated
solutions, the distortion of the 4th order
artificial dissipation terms is not so large.
4.2 Computed Results for a Channel Flow
Simulations are performed for a parallel
channel flow. Reynolds number(Re ) based on
the channel breadth(D) and the friction
velocity(u ) is 500. The calculated region has
a downstream length of 3D and a lateral width
of 0.75D. The downstream length and the
lateral width are subdivided into 30 equal
grid intervals. Therefore, the intervals in a
downstream direction(x1) and lateral
direction(x3) are as follows,
~Xt=O. 1 ~ X3=0. 025.
Nonuniform grid systems are used in the
normal direction. The minimum and maximum
intervals in the normal direction are as
follows,
~xamt n=O. 0033 Axamax=O. 0728*
In this calculation, AdamsBashforth and Euler
explicit methods are employed for the time
derivative.Smagorinsky's constant 'c' is set
to be 0.1. The convergence condition of the
Poisson eq.(25) is set as the divergence of
the velocity is less than 10 2.
Dependency of the artificial dissipation
factor a on the calculated results is shown
in Figs.4,5. Figs.4a,b show the turbulence
intensities normalized by the friction
velocity. From top to bottom, they are
1/2 1/2, / ,
respectively. The angular brackets denotes
the x1x3 plane average of a quantity and the
double prime " denotes the deviation of a
quantity from the xlx3 plane average. They
are compared with the experimental data by
Clark[10] and HussainReynolds[11] and with
the calculated result by MoinKim[6]. The
calculation diverges when ~ is less than
0.03. Therefore artificial dissipation is
required for the stable calculation. However,
the serious effect of the artificial
dissipation on the turbulence intensities is
2~
OCR for page 285
shown in case of FEZ =0.06.
The prof i le of the mean velocity non
dimensionalized by the friction velocity is
shown in Figs.5a,b. The solid line
represents the law of the wall and the Coles'
log law. The experimental data by Hussain
Reynolds [ 11 ] (Re =1280) is also plotted in
Fig.5. The calculated velocity in the log law
region i s higher than the experimental one.
This discrepancy partially depends on the
dif ference of Reynolds number. A noticeable
discrepancy does not exist between such a
change of artificial dissipation factor
a=0.03 and 0.06. However, if we take a large
artificial dissipation factor, for instance
a=l.o( so called Kawamura's scheme[l2]), the
velocity profile is seriously deformed[ 13 ].
The order of magnitude of each term in the
streamwise momentum equation is compared in
Figs.6, 7. The contour lines shown in Figs.6, 7
are the convection terms, artificial
dissipation terms, the SGS model
terms(Leonard, Cross and SGS Reynolds terms)
and the viscous terms on xlx3 plane. The
convection terms are dominant in the middle of
channel, whereas the convection and the
viscous terms are dominant near the wall. The
SGS model terms are too smal 1 compared with
the convection terms or the vi scous terms. In
order to implement the LES calculation, the
artificial dissipation terms must be smaller
than the SGS model terms. In case of ar=0.06,
the artificial dissipation terms are greater
than or the same order as the SGS model terms
in places on y =22.2 plane. Because of thi s
artificial dissipation terms, the turbulence
intensities are excessively damped in the
middle of the channel. As a result each term
is damped as shown in Fig.7b.
Even for the case of a =0.03, the
artificial dissipation terms are the same
order as the SGS model terms in the middle of
the channe 1. So that, there sti 11 remains room
to improve the present computing scheme.
When we use Euler explicit method instead
of the AdamsBashforth method, the scheme of
the time derivative is only first order
accuracy. To proceed the stable calculation,
the artificial dissipation factor must be
greater than 0.15 for Euler explicit method,
which is 5 times as large as the factor
required for AdamsBashforth method.
The contour lines of each term calculated
by Euler explicit method are shown in Fig.8.
The artificial dissipation terms are greater
than the SGS model terms and also greater than
the results calculated by AdamsBashforth
method. This is caused by the large
di s s ipation factor.
Fig.9 shows the turbulence intensities
calculated by Euler explicit method. Although
the artificial dissipation terms are greater
than the SGS model terms, the calculated
turbulence intensities are not so much damped
but agree with the experiments. In Fig.10,
the profile of the mean velocity is shown. The
velocity profile calculated by Euler method is
almost same with that of AdamsBashforth
method except the middle of channel. When we
take a larger dissipation factor, the
turbulence intensities are damped in the
middle of channel as shown in Fig.11. This
tendency is the same with the case of Adams
Bashforth method(Fig.4b). Figs.12a,b,c show
the instantaneous velocity vectors on x2x3
plane. Although the velocity f fluctuation is
attenuated in the middle of the channel for
a =0.06(AdamsBashforth), the velocity
fluctuation is not so much attenuated for
a=0.15(Euler). When we chose a dissipation
factor suitable for the calculating method,
even for Euler method, the global feature of
the turbulence f low can be obtained.
Conclusions
Turbulent channel f low is studied
numerically by using finite difference method.
The simulation is stably executed by using the
3 rd orde r upw i nd di f f e ren ce scheme. The
numerical dissipation terms attenuate the
turbulence intensities. When we chose a
dissipation factor suitable for the
calculating method, turbulent f low structures
in a turbulent channel f low can be well
simulated. However a problem still remains
how to determine the dissipation factor.
The present study is partially undertaken
as a student pro ject at Hiroshima University.
Author i s indebted to Mr. K. Murakami who
discussed and helped the author in the course
of this study.
The computations were executed by apollo
DN10000 and HITAC M680H at Hiroshima
University.
Ref erences
[1] Deardroff,J.W.; A Numerical Study of
Three Dimensional Turbulent Channel Flow
at Large Reynolds Numbers, J. Fluid
Mech. vol.41, (1970), p.453.
[2] Miyata, H., Kajitani, H., Zhu, M.,
Kawano, T., Takaki, M.: Numerical Study
of Some WaveBreaking Problems by a
FiniteDifference Method, J. of the
Kansai Soc. of Naval Architects, (1987),
pp.1124
289
OCR for page 285
[3] Sato, T., Miyata, H., Baba, N., Kajitani,
H.: FiniteDifference Simulation Method
for Waves and Viscous Flows about a Ship,
J. of the Soc. of Naval Architects of
Japan, (1986), pp.l420
[4] Smagorinsky,J., Manabe,S. & Holloway,J.L.
; Numerical Results from a NineLevel
General Circulation Model of the Atmo
sphere, Mon. Weath. Rev. vol.93,
(1965), p.727.
[5] Yoshizawa,A.; A StatisticallyDerived
Subgrid Model for the LargeEddy Simula
tion of Turbulence, Phys. of Fluids.
25(9), (1982), p.1532.
[6] Moin,P., Kim,J.; Numerical Investigation
of Turbulent Channel Flow, J. Fluid Mech.
vol.118, (1982), p.341.
[7] Grove, A.S., Shair, E.S., Petersen, E.E.
and Acrivos, A.: An experimental
investigation of the steady separated
flow past a circular cylinder, JFM
vol.19, (1964), pp.6080.
[8] Kodama Y. : Computation of the Two
Dimensional Incompressible NavierStokes
Equations for Flow Past a Circular
Cylinder Using an Implicit Faetored
Method, Papers of Ship Research
Institute, Vol.22, No.4, July (1985),
pp.335377.
[9] Tritton, D.J.: Experiments on the flow
past a circular cylinder at low Reynolds
numbers, JFM vol.6, (1959), pp.547567.
[10] Clark, J.A.; A Study of Incompress
ible Turbulent Boundary Layers in Chan
nel Flow, Trans. ASME, J. Basic Eng.
vol.90, (1968), p.455.
[11] Hussain,A.K.M.F., Reynolds,W.C.; Measure
ment in Fully Developed Turbulent Chan
nel Flow, Trans. ASME, J. Fluid Eng.
vol.97, (1975), p.568.
[12] Kawamura, T. et al.:Computation of High
Reynolds Number Flow around a Circular
Cylinder with Surface Roughness, AIAA
paper, 840340 (1984)
[13] Doi, Y., Kimura, T.: LargeEddy
Simulation of Turbulent Channel Flow,
J. of the Soc. of Naval Architects of
Japan, (1988), pp.2330
, ,
~ // U2
/ ~ U3
2/ ~ X3
X' ~

o
. _
o
~3
o_
o
o
— _
1.0
290
/'
/
/ 1
/ ' ~
o
Fig.1 Coordinate system
Fore :~\\ Aft
Exp.(Grove) ~ ~,k,' Pe ~nt Cal.
~ 1 1
0.00 30.00 60.00 90.00 1 20.00 1 50.00 1 80.00
Fore
(deg.) Aft
Fig.2 Pressure distribution on a eireular
cylinder at Re=40
\
Tritton
· Kodama
Fig.3 Drag eoeffieient
OCR for page 285
a) $1
=1 ~ ~
g
c~
b) O
o
o
8
~^p(^
O ~ ~~
_
8 oo ~o oo 80 00 t20.oo 160.oo 200.oo 240.oo
Y+
+ Exp.(Clark Re =1300)
* Exp.(HussainReynolds Re =1280)
ca l ( MoinKim R~ =1280 )
_ t.~
_ /
00 40. 00 80. oo 1~. 00 1 60. oo =0. 00 240.
Y+
b)
Fig.4a Turbulence intensities
t=10.5, Re =500
a 1/2, b 1/2
c: (U3) >
a) O
J —
~_


8
. ~ ~ ~ ~ ~ , , , , . ,
o.oo ~o.oo 80.oo t20.oo I eo oo 200.oo 240.00
Yt
+ Exp.(Clark Re =1300)
* Exp.(HussainReynolds Re =1280)
Cal.(MoinKim Re =1280)
8 r~
O
O. 00 40. 00 80. 00 120. 00 1 60. ~ 200. 00 240.
Y~
A
~.

o
o
o
OC
+

\
_ /
,_
~O 00 40. 00 80. 00 120. 00 1 80. ~ 200. 00 240.
Y~
a =0.03 Fig.4b Turbulence intensities
t=10, Re =500
a l/2' b 1/2
C: 1/2
o
o
u~—
* Exp.(HussainReynolds Re =1280)
Present Cal.
/~
L
_ ~
~ ,'
_^
A~
1 1 1 1
 10 1o2 y+ 1
Fig.5a Mean velocity profile Cl =0.03
t=10.5, Re =500
291
~ =0.06
* Exp.(HussainReynolds Re =1280)
Present Cal.
t. —
/\ Yd_
_
_ ^
a ,^
,'
i.'
I _
/~
~_
/~

1 1 1 1 1 +
10 1o2 y
Fig.5b Mean velocity profile ~ =0.06
t=10, Re =500
OCR for page 285
b)
a)
b)
d)
b)
c'
d)
~ ~ d)
Fig.6a Contour lines of each term in the stream Fig.6b Same as Fig.6a
wise momentum equation on x1x3 plane, a =0.06, y+=22.2, t=10, Re =500,
a =0.03, y =22.2, t=9.5, Re =500, e: convection terms interval 50,
e: convection terms interval 50, b:artificial dissipation terms interval 2
b:artificial dissipation terms interval 2 c:SGS model terms interval 2, d:viscous
c:SGS model terms interval 2, d:viscous terms interval 10
terms interval 10
Fig.7a Same as Fig.6a
a=0.03, y =250, t=9.5, Re =500,
e: convection terms interval 50,
b:artificial dissipation terms interval
c:SGS model terms interval 2, d:~7~scous
terms interval 10
292
b
r:~~^O~:
C) ~
it== =
Fig.7b Same as Fig.6a
a =0.06, y+=250, t=10, Re =500,
a:convection terms interval 50,
b:artificial dissipation terms interval 1
c:SGS model terms interval 1, d:viscous
terms interval 1
OCR for page 285
c)
d)
h)
d
.~
~=f~=
f~ r
Fig.8a Contour lines of each term in the stream
wise momentum equation calculated by
Euler method,
~ =0.15, y =22.2, t=9.5, Re =500, ~
a:convection terms interval 50, o
b:artificial dissipation terms interval 5
c:SGS model terms interval 5, d:viscous
terms interval 10
o 1
1. >'~ : 
b)
. _
~,

o
o
o
_ ~
°O 00 40. 00 80. 00 ~ 20. 00 1 60. 00 200. 00 24 0. 00
Yt
/~+~
/
o
o~
Fig.8b Same as Fig.8a
a =0.15, y =250, t=9.5, Re =500,
e: convection terms interval 50,
b:artificial dissipation terms interval 2
c:SGS model terms interval 2, d:viscous
terms interval 10
80. 00 1 20. 00 1 80. 00 200. 00 240.
Y~
U
W
~
40. 00 80. 00 120. 00 I 00. 00 200. 00 240. 00
Yt
Fig.9 Turbulence intensities calculated by
Euler method, ~=0.15
*
t=9.5, Re =500
a l/2' b 1/2
c:1/2 2
o
o
u'_

o
o
* Exp.(HussainReynolds Re =1280)
Present Cal.
/#
3
~a
3~
~
Fig.10 Mean velocity profile calculated by
Euler method, ~ =0.15
293 t=9.5, Re =500
OCR for page 285
U2 ~ US
, ~~ ~ · 2 ~ ~ ~ · ~ ~ `~ '"in
 44 ~ .` ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ .
^~ ~ ~ ~ ~ ~ , , ~ ~ ~ ~ ~ SO ;
,~ ~ ~ ~ . , , ~ ~ . ~ ~ ~ ~ t
a) ~ , , . ~ . .
Fig.12a Velocity vectors on x2x3 plane
oo 00 40.00 90.00 t20.00 lS0.00' 200.00' 240.00 calculated by AdamsBashforth method
YE a =o. 03, t=9.5, x1=1.5, Re =500
+ Exp.(Clark Re =1300)
b)~ * Exp.(HussainReynolds Re =1280)
_ ~ Cal.(MoinKim Re =1280)
'I _ '2_ ·P ·, ·, 2 ~  :~]
t~ ~ ~ :1
Fig.12b Same as Fig.12a
O / ~=0.06, t=10.0, xl=1.5, Re =500
o
._ , . . . . . . . . . .
=0 0 40. 00 80. 00 120. 00 1 60. 00 200. 00 240. 00
Y~
Fig.ll Turbulence intensities calculated by ;~;ttit;~,' ~ ~ ~ ~ ~ ~ {< .~`t. .;
Euler method ~ =0.6 · ~ ~ ~ ~ _ ~ ~ ~ ~ <~' ~`
a ~/2, b:<~O '2~/2 ~ ~ ,
Fig.12c Velocity vectors on x2x3 plane
calculated by Euler method
~=0.15, t=9.5, x1=1.5, Re =500
294