Below are the first 10 and last 10 pages of uncorrected machine-read 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, chapter-representative 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 Finite-Difference 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 sub-grid 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 Navier-Stokes 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 O-xlx2x3 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 sub-grid 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(xl-x1 ~x2-x2', X3-X3') ·f(x1',x2',X3')dx1'dx2'dx3' (4) ~2 a 2(U~Uj) (11) 24 axkaxk In the present study, Top-Hat 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 Moin-Kim[6] as follows, fD=( 1-exp(—y+/A+) ) (16) y+=x2Re* O X0.5 286

OCR for page 285
y+=(l-X2)Re* 0.5 OCR for page 285
Cyclic boundary conditions are imposed for velocity and pressure, in downstream and spanwise directions. No-slip 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. No-slip 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. The-calculated 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, 4-th 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 4-th order artificial dissipation terms, the calculation diverges always. Although the 2-nd order artificial dissipation terms distort the calculated solutions, the distortion of the 4-th 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. Non-uniform 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, Adams-Bashforth 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.4-a,-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 x1-x3 plane average of a quantity and the double prime " denotes the deviation of a quantity from the xl-x3 plane average. They are compared with the experimental data by Clark[10] and Hussain-Reynolds[11] and with the calculated result by Moin-Kim[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.5-a,-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 stream-wise 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 xl-x3 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.7-b. 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 Adams-Bashforth 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 Adams-Bashforth 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 Adams-Bashforth 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 Adams-Bashforth 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.4-b). Figs.12-a,b,c show the instantaneous velocity vectors on x2-x3 plane. Although the velocity f fluctuation is attenuated in the middle of the channel for a =0.06(Adams-Bashforth), 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 M-680H 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 Wave-Breaking Problems by a Finite-Difference Method, J. of the Kansai Soc. of Naval Architects, (1987), pp.11-24 289

OCR for page 285
[3] Sato, T., Miyata, H., Baba, N., Kajitani, H.: Finite-Difference Simulation Method for Waves and Viscous Flows about a Ship, J. of the Soc. of Naval Architects of Japan, (1986), pp.l4-20 [4] Smagorinsky,J., Manabe,S. & Holloway,J.L. ; Numerical Results from a Nine-Level General Circulation Model of the Atmo- sphere, Mon. Weath. Rev. vol.93, (1965), p.727. [5] Yoshizawa,A.; A Statistically-Derived Subgrid Model for the Large-Eddy 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.60-80. [8] Kodama Y. : Computation of the Two- Dimensional Incompressible Navier-Stokes Equations for Flow Past a Circular Cylinder Using an Implicit Faetored Method, Papers of Ship Research Institute, Vol.22, No.4, July (1985), pp.335-377. [9] Tritton, D.J.: Experiments on the flow past a circular cylinder at low Reynolds numbers, JFM vol.6, (1959), pp.547-567. [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, 84-0340 (1984) [13] Doi, Y., Kimura, T.: Large-Eddy Simulation of Turbulent Channel Flow, J. of the Soc. of Naval Architects of Japan, (1988), pp.23-30 , , ~ // 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,' P-e ~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.(Hussain-Reynolds Re =1280) ca l ( Moin-Kim R~ =1280 ) -_ t.~ _ / 00 40. 00 80. oo 1~. 00 1 60. oo =0. 00 240. Y+ b) Fig.4-a 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.(Hussain-Reynolds Re =1280) Cal.(Moin-Kim 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.4-b Turbulence intensities t=10, Re =500 a l/2' b 1/2 C: 1/2 o o u~— * Exp.(Hussain-Reynolds Re =1280) Present Cal. /~ L _ ~ ~ ,' _^ A~ 1 1 1 1 - 10 1o2 y+ 1 Fig.5-a Mean velocity profile Cl =0.03 t=10.5, Re =500 291 ~ =0.06 * Exp.(Hussain-Reynolds Re =1280) Present Cal. t. -— /\ Yd_ _ _ ^ a ,^ ,' i.' I _ /~ ~_ /~ - 1 1 1 1 1 + 10 1o2 y Fig.5-b Mean velocity profile ~ =0.06 t=10, Re =500

OCR for page 285
b) a) b) d) b) c' d) ~ ~ d) Fig.6-a Contour lines of each term in the stream Fig.6-b Same as Fig.6-a wise momentum equation on x1-x3 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.7-a Same as Fig.6-a 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.7-b Same as Fig.6-a 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.8-a 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.8-b Same as Fig.8-a 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.(Hussain-Reynolds 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.12-a Velocity vectors on x2-x3 plane oo 00 40.00 90.00 t20.00 lS0.00' 200.00' 240.00 calculated by Adams-Bashforth method YE a =o. 03, t=9.5, x1=1.5, Re =500 + Exp.(Clark Re =1300) b)~- * Exp.(Hussain-Reynolds Re =1280) _- ~ Cal.(Moin-Kim Re =1280) 'I _ '2_ ·P ·, ·, 2 ~ - :~] t~ ~ ~ :1 Fig.12-b Same as Fig.12-a 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.12-c Velocity vectors on x2-x3 plane calculated by Euler method ~=0.15, t=9.5, x1=1.5, Re =500 294