Published online 31 December 2008 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/fld.1987
Roe scheme with preconditioning method for large eddy simulation
of compressible turbulent channel flow
Wu-Shung Fu
∗,†,‡, Chung-Gang Li, Wei-Fong Lin and Yu-Hsu Chen
Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta Hsueh Road,Hsinchu 30056, Taiwan
SUMMARY
Numerical investigations of fully compressible turbulent channel flows are conducted. Two different speeds of Mach numbers of 0.005 and 0.5 are separately taken into consideration to study mechanisms of compressible turbulent flows and increments of temperatures of working fluids. For computing the compressible low-speed turbulent flow efficiently, methods of Roe and preconditioning matching LES methods are adopted. In addition, usage of dual time stepping and curved linear coordinates transformation skills to stabilize transient situations and increase computing grids near the channel walls are performed. As for an assignation of initial conditions, random velocity fluctuations are imposed on the mean velocity distribution. The results show that the mean velocity distributions, root-mean-square values of fluctuating velocities and mean Reynolds stresses have good agreements with the existing results. The increments of temperatures of working fluids are remarkable for a Mach number of 0.5 and are also well consistent with the existing results. Copyrightq 2008 John Wiley & Sons, Ltd.
Received 8 July 2008; Revised 10 November 2008; Accepted 11 November 2008
KEY WORDS: compressible flow; LES; preconditioning method; Roe scheme; dual time stepping; mean Reynolds stresses
1. INTRODUCTION
It is well recognized that most real flow phenomena belong to turbulent flows.
In addition to high-speed turbulent flows, recently some phenomena such as compressible low-speed turbulent flows attract much attention as well, including acoustic noise induced by compressible turbulent flows, flows in combustion chamber, natural convection induced by high
∗Correspondence to: Wu-Shung Fu, Department of Mechanical Engineering, National Chiao Tung University, 1001
Ta Hsueh Road, Hsinchu 30056, Taiwan.
†E-mail: [email protected] ‡Professor.
temperature difference, etc. Fortunately with the promotion of computer capability and improve-ments in numerical schemes, the simulation of compressible turbulent flow inspite of high and low speeds could be expected to be solved more realistically.
Compared with the numerical turbulent solvers of the Reynolds averaged Navier–Stokes equa-tions and the direct numerical simulation (DNS), the method of large eddy simulation (LES), which has merits of practicability and efficiency is more adaptable to be used in simulating time-evolving turbulent flow. However, there is a critical problem of initial or inlet condition of flow field, which should be seriously considered during the LES solver being executed. Therefore, several methods of generation of initial condition used in the LES are proposed. In terms of the initial condition there are three ways to obtain: the results of the DNS method, random fluctuations imposed on theoretically mean velocities, and initial or inlet conditions matching with turbulent correlation. Tong [1] used the DNS results obtained by Gavrilakis [2] as the initial condition of LES method matching a modified Roe scheme to simulate a turbulent flow in a square duct. The results of Tong[1] were compared well with the results obtained by the DNS method. Stolz
et al.[3] simulated channel flows of Re=180 and Re=590 by the LES method with approximate
deconvolution model, and they had good results with DNS and LES methods despite the initial conditions obtained from the DNS method or the data from the LES at different parameters or synthetic data with random fluctuations. However, it is not easy to obtain the adaptable results with the DNS method as the initial condition of the LES method in practice except in special situations. Lund et al.[4] developed a new method called random fluctuation inflow generation method to generate reasonable inflow data used in the LES method as the initial condition. This method was based on turbulent statistical properties and the inflow data generated by the method could be conformed with turbulent correlations automatically. But the results obtained by the method were slightly different from those obtained by the DNS method as the flow conditions of inlet and outlet were assigned simultaneously. Klein et al.[5] used an autocorrelation function of turbulent flow to develop a new method of generation of inlet or initial condition of a turbulent flow. Further the new method matching the study investigated by Lund et al.[4] could be adapted for an inhomogeneous turbulent flow. The results had good agreements with experimental results.
Besides, consideration of compressibility of fluid under low-speed turbulent flows is another wrinkle. In this situation, the speed of the sound wave is much larger than that of the fluid flow. Therefore, in an explicit numerical method, the time step limited by Courant–Friedrichs–Levy condition should be extremely small and the convergence of calculation becomes difficult. In an implicit numerical method, a stiff phenomenon easily occurs in the matrix used in solving processes, which causes the situation to be solved only with great difficulty. Briley et al.[6] adopted a preconditioning method to add into an implicit method to improve the convergence of Navier– Stokes equations under a low Mach number situation. Turkel[7] developed a preconditioning matrix and discussed the applications of this method into incompressible and compressible flows. Choi and Merkel[8] added a preconditioning matrix into Euler equations and successfully improved the convergent problems caused by the stiff phenomenon and factorization error in an inviscid low Mach number flow (Ma no.=0.05) when an implicitly numerical method was used. Choi and Merkel [9] continuously investigated the former problems in[8] and proposed an adaptable preconditioning matrix to resolve the calculating convergence of low-speed viscous flows. For effectively resolving compressible flows troubles, Roe[10] proposed an ingenious method to define averaged variables that were used to calculate fluxes between interfaces of computational grids, and successfully resolved the problem of discontinuity happening at the interface of computational grids. This method has been widely applied to the problem of compressible flow recently. Weiss and Smith[11]
extended the research of Choi and Merkel[9], and applied the Roe and preconditioning methods into three-dimensional (3-D) Navier–Stokes equations. Besides, a computational process of dual time stepping was utilized to calculate transient flows under a low Mach number situation. Dennis
et al.[12] added multigrid method into OVERFLOW program to improve convergent efficiency,
and used Roe and preconditioning methods simultaneously to calculate compressible low-speed flows. Shishir et al.[13] continuously added the dual time-stepping process to compute transient flows. The results showed that the application of preconditioning into low-speed compressible flows remarkably improved both accuracy and convergent efficiency.
It is difficult to obtain the initial or inlet condition of LES method used in compressible turbulent flows under low-speed conditions, and there is comparatively little literature adopting the methods of LES and preconditioning simultaneously to investigate the compressible low-speed turbulent flows. Lessani et al.[14] applied the preconditioning method, with the multigrid skill into the LES method to simulate channel and cavity flows. The results showed that the convergence speed is increased 4–7 times compared with an explicit method. Xu et al.[15] used the finite volume and LES methods to calculate a pipe flow. The results compared with the results of experimental and DNS methods showed good agreements. Alkishriwi et al. [16] calculated channel and circular flows by the LES method matching the methods of preconditioning and multigrid under different Mach numbers and time intervals. This synthetic method had wonderful efficiency and improved the calculating convergence of low Mach numbers problem remarkably. The convergent speed of this method was quicker than that of the method of the fifth-order Runge–Kutta by about 4–60 times. Tong[1] also used Roe and LES methods, but without the preconditioning method to simulate duct flows, and improved a problem of excessively large dissipation which happens in the LES method. However, to the best of the authors’ knowledge, the simulation of compressible low-speed turbulent flows, which usually occur in many industrial problems, by the LES method matching Roe and preconditioning methods were preformed with difficulty.
Therefore, the aim of this study is to propose an appropriate method in which the combination of the Roe, preconditioning, and LES methods is performed to investigate compressible low (Ma no.=0.005) and high (Ma no.=0.5) speed turbulent channel flows. Random fluctuations imposed on the theoretical mean velocities are assumed as the initial condition. For stabilizing transient situations, a method of dual time stepping is conducted. Besides, in order to obtain accurate results of extremely near-wall region, a curved linear coordinates transformation is used, and five computation grids are distributed within the viscous sublayer region (x2+<10). The results show that the distributions of the mean velocity and the turbulent correlations are well consistent with the results of Kim et al. [17] under both low- and high-speed situations. Besides, in the high-speed situation the temperatures of fluids are apparently raised near the channel wall region, and the results have good agreement with the existing results of Lenormand et al.[18].
2. PHYSICAL MODEL
A turbulent channel flow in which the compressibility and viscosity of fluid are taken into consid-eration simultaneously is investigated . The LES method is adopted to describe the turbulent flow and the corresponding physical model is shown in Figure 1. The streamwise, vertical, and spanwise directions are x1, x2, and x3, respectively, and the corresponding velocities are u1, u2, and u3, respectively. The length, height, and width are l1, l2, and l3, respectively.
x1 x3 x2 l3 l1 l2 1atm 298.0592k P T = =
Figure 1. Physical model.
For facilitating the analysis, the following assumptions are made: (1) The fluid is air and ideal gas.
(2) The effect of gravity is neglected.
According to Tong[1], the governing equations are expressed as follows: *U *t + *F1 *x1 +*F2 *x2 +*F3 *x3 =0 (1)
The quantities included in U and Fi are separately shown in the following equations:
U= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ¯ ¯ ˜u1 ¯ ˜u2 ¯ ˜u3 ¯˜e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (2) and Fi= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ¯ ˜ui ¯ ˜ui˜u1+ ˜Pi 1−(+2¯t)Ai 1 ¯ ˜ui˜u2+ ˜Pi 2−(+2¯t)Ai 2 ¯ ˜ui˜u3+ ˜Pi 3−(+2¯t)Ai 3
(¯˜e+ ˜P) ˜ui−(¯+2¯t) Ai j˜uj− ¯+ ¯ t Prt * ˜T *xi ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ∀i =1,2,3 (3) where Ai j=* ˜uj *xi + * ˜uj *xi − 2 3(∇ · ˜u)i j and the filtered ideal gas is written as
Sutherland’s law is adopted to evaluate the viscosity shown as follows: =0 T T0 2/3 T0+110 T+110 (5)
where0and T0are 1.85×10−5N s/m2and 298.0592K, respectively. According to the Smagorinsky model,t is expressed as follows:
t=C2| ¯S| (6) where=(123)1/3,| ˜S|= 2 ˜Si j˜Si j ˜Si j= 1 2 * ˜u j *xi + * ˜uj *xi and Prt=0.71.
The parameter C in Equation (6) is determined by the Van Driest damping function of Smagorinsky model and expressed as follows:
C=0.01 1−exp − d+ 25 (7)
Where d+=ud/,u is a friction velocity and expressed in terms of √w/, and d is the distance from the wall.
3. NUMERICAL METHOD
A situation of fully developed compressible turbulent flow described by the LES model is inves-tigated. In order to validate the results of the present study with those of the previous study by Kim et al. [17], the magnitude of the height of channel is designed as a characteristic length and the same magnitude of Reynolds number Re used in the previous study of Kim et al.[17] is selected. In the streamwise direction, the length used for computation should be longer than the length of streak growing in the turbulent flow and is about 3.82 times the magnitude of the height. Furthermore, for economizing computational memory and time, the length of the channel is periodically used until a fully turbulent flow is formed. Based on Kazishima[19], for satisfying the characteristics of 3-D eddy motions of the turbulent flow, the width used for computation in the spanwise direction is about 1.92 times the magnitude of the height. Periodical conditions are held on the both sides of this width.
Except the usage of Roe and the preconditioning methods, two effective methods are used to avoid the occurrence of inefficiency and inaccuracy in computing processes. One is the method of 3-D curvilinear coordinates(, ,) to be used, then five grids are distributed within x2+=10 region. A dual time-stepping process is the other one. The original governing equation (1) is then transformed into the following equation:
* ¯Up * + * ¯U *t + * ¯F1 * + * ¯F2 * + * ¯F3 * =0 (8)
where is the preconditioning matrix derived by Weiss et al. [11], and ¯Upis the primitive form [ ˜P, ˜u1, ˜u2, ˜u3, ˜T ]/J in which J is the Jacobian matrix. and t are the artificial and physical times, respectively, ¯U is the conservative form of (¯, ¯ ˜u1, ¯ ˜u2, ¯ ˜u3, ¯˜e)/J. According to Dennis
et al.[12], the preconditioning parameter is chosen as =max(min(M2,1.0), min), where M
is the local Mach number and min≈3M∞2 .
Discrete Equation (8): The terms of* ¯Up/* and * ¯U/*t are differentiated by first-order forward difference and second-order backward difference, respectively, and the terms of* ¯F1/*, * ¯F2/* , and* ¯F3/* are differentiated by a central difference, the following equation can be obtained:
¯U k+1 p − ¯Upk + 3 ¯Uk+1−4 ¯Un+ ¯Un−1 2t + 1 ( ¯F1ki+1+1/2, j,k− ¯F k+1 1i−1/2, j,k) + 1 ( ¯F k+1 2i, j+1/2,k− ¯F k+1 2i, j−1/2,k)+ 1 ( ¯F k+1 3i, j,k+1/2− ¯F k+1 3i, j,k−1/2)=0 (9)
In which superscripts of k and n indicate iteration numbers of artificial time and real time, respectively. The real-time counts a proceeding step when the term of artificial time * ¯Up/* is convergent to zero, Equation (9) automatically transfers into the Navier–Stokes equation and the magnitudes at the iteration number of(k +1) of the artificial time in Equation (9) substantially become the magnitudes of the proceeding step of(n+1) of the real time.
Afterward the terms of ¯Uk+1 and ¯Fik+1 in Equation (9) are necessary to be linearized and expressed as follows, respectively
¯Uk+1= ¯Uk+ M ¯U
p (10)
where M=* ¯U/* ¯Upand ¯Up= ¯Upk+1− ¯Upk ¯Fk+1
1 = ¯F1k+ Ap ¯Up (11)
Where Ap=* ¯F1k/* ¯Up is the flux jacobian and the same methods of Bp=* ¯F2k/* ¯Up and
Cp=* ¯F3k/* ¯Up are used in linearization of ¯F2K+1 and ¯F3K+1, respectively.
Substituting Equations (10) and (11) into Equation (9), the following equation is obtained:
¯Up
+
3( ¯Uk+ M ¯Up)−4 ¯Un+ ¯Un−1 2t
+( ¯F1k+ Ap ¯Up)+ ( ¯F2k+ Bp ¯Up)+( ¯F3k+Cp ¯Up)=0 (12) Where, , and are central-difference operators. Equation (12) can be rearranged as follows:
I + M 3 2t+(Ap+ Bp+Cp) ¯Up= Rk (13) Where Rk=− 3 ¯Uk−4 ¯Un+ ¯Un−1 2t −(¯Fk 1+ ¯F2k+ ¯F3k)
Dividing in both sides, the following equation is obtained: I +−1M 3 2t+ −1( Akp+ Bpk+Cpk) ¯Up=−1Rk (14)
The solver of Equation (15) is the LUSGS implicit method proposed by Yoon et al.[20] which is used to solve Equation (14).
˜Ap= −1Akp ˜Bp= −1Bpk ˜Cp= −1Ckp
(15)
and divide ˜Ap, ˜Bp, and ˜Cp into two parts, respectively ˜Ap= ˜A+p+ ˜A−p ˜Bp= ˜Bp++ ˜Bp− ˜Cp= ˜Cp++ ˜Cp− (16) where ˜A± p = 12( ˜Ap±|˜A|I ) ˜B± p = 12( ˜Bp±|˜B|I ) ˜C± p = 12( ˜Cp±|˜C|I ) (17)
˜A,˜B, and ˜C are the largest eigenvalues of ˜Ap, ˜Bp, and ˜Cp, respectively.
Substituting Equations (15) and (17) into Equation (14), we obtain the following equation: I +−1M 3 2t+( ˜A + p + ˜A−p)+ ( ˜Bp++ ˜Bp−)+( ˜Cp++ ˜Cp−) ¯Up=−1Rk (18) The convective terms of ( ˜A+p+ ˜A−p) can be solved by a high-order scheme as the following equation:
( ˜A+p+ ˜A−p)=− ˜A+p++ ˜A−p = ˜A+ p,i− ˜A+p,i−1 + ˜A− p,i+1− ˜A−p,i (19)
where the backward difference approximation− is adopted for ˜A+p and the forward difference approximation+ is adopted for ˜A−p.
Substituting Equation (19) into (18), the following equation can be derived as: I +−1M 3 2t+ ˜A+ p,i− ˜A+p,i−1 + ˜A− p,i+1− ˜A−p,i + ˜B + p, j− ˜Bp+, j−1 + ˜B− p, j+1− ˜Bp−, j + ˜C+ p,k− ˜Cp+,k−1 + ˜C− p,k+1− ˜Cp−,k ¯Up=−1Rk (20)
Equation (20) can be arranged as follows: (L + D+U) ¯Up=−1Rk (21) Where L=− 1 ( ˜A+p)i−1, j,k+ 1 ( ˜Bp+)i, j−1,k+ 1 ( ˜Cp+)i, j,k−1 (22) D= I +−1M 3 2t+ 1 (( ˜A+p)i, j,k−( ˜A−p)i, j,k) + 1 (( ˜Bp+)i, j,k−( ˜Bp−)i, j,k)+ 1 (( ˜Cp+)i, j,k−( ˜Cp−)i, j,k) (23) U= 1 ( ˜A−p)i+1, j,k+ 1 ( ˜Bp−)i, j+1,k+ 1 ( ˜Cp−)i, j,k+1 (24)
As for the computation of Rk=−((3 ¯Uk−4 ¯Un+ ¯Un−1)/2t)−(¯F1k+ ¯F2k+¯F3k) in right-hand side (RHS) of Equation (21), the terms of Fi in Equation (3) based on Cartesian coordinates can be divided into two parts. One is the inviscid term Finviscid.
Finviscid= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ¯ ˜ui ¯ ˜ui˜u1+ ¯Pi 1 ¯ ˜ui˜u2+ ¯Pi 2 ¯ ˜ui˜u3+ ¯Pi 3 (¯˜e+ ¯P) ˜ui ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (25)
The other is viscous term Fviscous.
Fviscous=− ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 (+2¯t)Ai 1 (+2¯t)Ai 2 (+2¯t)Ai 3 (¯+2¯t) Ai j˜uj+ ¯+ ¯ t Prt * ˜T *xi ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (26)
The Roe upwind difference scheme[10] is employed in discretion of the term of Finviscidat the cells interface(i +12) and expressed as follows at a low Mach number situation:
Finviscid,i+1/2=12(FR+ FL)−12ε{|−1Ap|Up} (27) In whichε is less than 0.1 and appropriate for a turbulent flow [1].
The high-order MUSCL scheme proposed by Abalakin et al. [21] is used to compute Equation (27), and the related derivations are indicated as follows:
Up=uiL+1/2−uiR+1/2 (28)
uiL+1/2=ui+12uiL+1/2 (29)
uiR+1/2=ui−12uiR+1/2 (30) uL
i+1/2= (1− )(ui+1−ui)+ (ui−ui−1)+ c(−ui−1+3ui−3ui+1+ui+2) + d(−u
i−2+3ui−1−3ui+ui+1) (31)
uR
i+1/2= (1− )(ui+1−ui)+ (ui+2−ui+1) + c(−u
i−1+3ui−3ui+1+ui+2)+ d(−ui+3ui+1−3ui+2+ui+3) (32) Where , c, and d are tabulated in Table I. The fifth-order scheme is adopted to reduce the dissipation. When a high-order scheme is executed, a limiter function is usually used to suppress occurrence of numerical oscillations. However, to prevent the decay of the turbulent intensities of flow, the limiter function is then not adopted. The derivative terms of Ai j=* ˜uj/*xi+* ˜uj/*xi− 2
3(∇ · ˜u)i j in Equation (26) are computed by the fourth-order central difference * ˜u
*x=
˜ui−2−8 ˜ui−1+8 ˜ui+1− ˜ui+2
12x +o(x
4) (33)
In order to compute the terms of ¯F1k, ¯F2k, and ¯F3k of Rkshown in Equation (21) in the curved linear coordinates(, ,), the ¯F1k, ¯F2k, and ¯F3k can be expressed in terms of Fi indicated as follows:
¯Fk 1 = (x1F1+x2F2+x3F3)/J ¯Fk 2 = ( x1F1+ x2F2+ x3F3)/J ¯Fk 3 = ( x1F1+ x2F2+ x3F3)/J (34)
Table I. The parameters of the MUSCL scheme.
c d Order 1 3 0 0 3 1 3 − 1 6 0 4 1 3 0 − 1 6 4 1 3 −101 −151 5
then Rk= − 3 ¯Uk−4 ¯Un+ ¯Un−1 2t −{[(x1F1+x2F2+x3F3)/J] + [( x1F1+ x2F2+ x3F3)/J]+[( x1F1+ x2F2+ x3F3)/J]}
Substituting Equation (34) into Rk of Equation (21), Equation (21) can be expressed as follows:
L(i−1, j,k)Upk,(i−1, j,k)+L(i, j−1,k)Upk,(i, j−1,k)+L(i, j,k−1)Upk,(i, j,k−1)+D(i, j,k)Upk,(i, j,k)
+U(i+1, j,k)Upk,(i+1, j,k)+U(i, j+1,k)Upk,(i, j+1,k)+U(i, j,k+1)Upk,(i, j,k+1)=−1R(i, j,k)k (35) Based on Yoon et al.[20] and modification, Equation (21) could be expressed as follows:
(L + D)D−1(D+U) ¯Uk
p=−1Rk (36)
The solving procedures of Equation (36) are briefly shown in the following processes: 1. Derive the following equations to obtain forward sweep process
(L + D)Up∗=−1Rk (37) WhereUp∗= D−1(D+U) ¯Upk
Equation (37) can be rearranged as the following equations:
LUp∗+ DUp∗=−1Rk (38)
DUp∗=−1Rk− LUp∗ (39)
Equation (39) can be solved as follows:
Up∗= D−1(−1Rk− LUp∗) (40)
Equation (40) can be expanded as the following equation to execute the forward sweep process
Up∗,(i, j,k)= D(i, j,k)[−1R(i, j,k)k − L(i−1, j,k)Up∗,(i−1, j,k)
−L(i, j−1,k)Up∗,(i, j−1,k)− L(i, j,k−1)Up∗,(i, j,k−1)] (41) When i=1 or j =1 or k =1, the Up∗in RHS can be obtained from boundary conditions. 2. Derive the following equations to obtain backward sweep process.Upkin Equation (37) can
be expressed as the following equations:
(D+U) ¯Uk
p= DUp∗ (42)
Equation (42) can be rearranged as the following equations: ¯Uk
¯Uk
p can be expanded and solved as the following equations: Uk
p,(i, j,k)= Up∗,(i, j,k)− D−1(i, j,k)[U(i+1, j,k)Upk,(i+1, j,k)
+U(i, j+1,k)Upk,(i, j+1,k)+U(i, j,k+1)Upk,(i, j,k+1)] (44) When i=nx or j =ny or k =nz, the Upin RHS can be obtained from boundary conditions. 3. Then ¯Upk+1 can be solved as follows:
¯Uk+1
p = ¯Upk+ ¯Upk
4. Repeat step 1 to step 3 until ( ¯Upk+1− ¯Upk)/≈0. Then ¯Upn+1= ¯Upk+1, the next time step ¯Un+1
p can be obtained.
The adiabatic and no-slip conditions are adopted on the wall in case A tabulated in Table II. The equations are given as follows:
˜P(i,0,k) = ˜P(i,1,k) ˜u1(i,0,k) = − ˜u1(i,1,k) ˜u2(i,0,k) = − ˜u2(i,1,k) ˜u3(i,0,k) = − ˜u3(i,1,k)
˜T (i,0,k) = ˜T(i,1,k)
(45)
The isothermal and no-slip conditions are adopted on the wall in case B tabulated in Table III. The equations are given as follows:
˜P(i,0,k) = ˜P(i,1,k) ˜u1(i,0,k) = − ˜u1(i,1,k) ˜u2(i,0,k) = − ˜u2(i,1,k) ˜u3(i,0,k) = − ˜u3(i,1,k)
˜T (i,0,k) = 2Tfix− ˜T (i,1,k)
(46)
Table II. Computed parameters combinations (case A).
Parameters Magnitudes
Mean Mach number, c¯u 0.005
Mean streamwise velocity, ¯u 1.51m/s
Average friction velocity, u 0.11248m/s
Grids distribution, i× j ×k 64×64×64
l×m ×h 0.192×0.096×0.05m3
Eddy turnover time, te=0.5hu 0.222s
Initial temperature 298.05K
Table III. Computed parameters combinations (case B).
Parameters Magnitudes
Average Mach numbers, c¯u 0.5
Average streamwise velocity, ¯u 157.06m/s
Average friction velocity, u 11.7412m/s
Grid, i× j ×k 64×64×64
l×m ×h 0.00183936×0.000091968×0.000479m3
Eddy turnover time, te=0.5hu 2.04×10−5s
Initial temperature 298.05 K
Time interval,t 1×10−7s
1 0
wall
Figure 2. The grid distribution on the wall.
Where Tfix is wall temperature with a constant of 298.0592 K. 0 indicates the ghost cell and 1 indicates the cell nearest the wall. The relative positions are shown in Figure 2.
The computational length along the streamwise direction is periodically used, then the inlet and the outlet conditions can be expressed as the following equations, respectively.
˜P(0, j,k) = ˜P(nx, j,k) ˜u1(0, j,k) = ˜u1(nx, j,k) ˜u2(0, j,k) = ˜u2(nx, j,k) ˜u3(0, j,k) = ˜u3(nx, j,k) ˜T (0, j,k) = ˜T(nx, j,k) (47) ˜P(nx +1, j,k) = ˜P(1, j,k) ˜u1(nx +1, j,k) = ˜u1(1, j,k) ˜u2(nx +1, j,k) = ˜u2(1, j,k) ˜u3(nx +1, j,k) = ˜u3(1, j,k) ˜T (nx +1, j,k) = ˜T(1, j,k) (48)
0 indicates the cell at the inlet and nx+1 indicates the cell at the outlet. The relative positions are indicated in Figure 3.
The pressure gradient is the driving force that causes the channel flow to flow continuously. The local pressure gradient of * ˜P(i, j,k)/*x1 includes two parts. One is the mean pressure gradient of* ˜Pmean/*x1which mainly drives the flow, and the other is the fluctuating pressure gradient of
1 0 streamwise nx+1 nx Outlet Inlet
Figure 3. The grid distributions in the inlet and the outlet.
* ˜Pp/*x1 which is variable. The relationships of the above pressure gradients can be expressed as follows: * ˜P *x1= * ˜Pmean *x1 + * ˜Pp *x1= + * ˜Pp *x1 (49)
is a source term in computing processes and the periodic pressure conditions at the inlet and
outlet are indicated as follows:
˜Pp(0, j,k)= ˜Pp(nx, j,k) (50) ˜Pp(nx +1, j,k)= ˜Pp(1, j,k) (51) The mass flow rate in the channel flow is not easily kept conservatively due to the viscous friction along the walls and the numerical dissipation. Therefore, Equation (52) derived by Xu
et al. [15] is adopted to adjust the magnitude of to keep the mass flow rate constant during
computation. n+1= n− 1 t ˙m AC 0 −2 ˙m AC n + ˙m AC n−1 (52)
where AC is the cross section area of the channel andt is the physical time interval, and ˙m is the mass flow rate of the channel flow.
4. RESULTS AND DISCUSSION
In this work, air is the working fluid and Reynolds number Re=u/ is assigned and equal to 180. In order to investigate the flow mechanisms of fully compressible turbulent flow and the increments of temperatures of the working fluids, two different mean velocities of 0.005 (case A) and 0.5 (case B) Mach numbers of the flow are selected, but the Reynolds numbers of the both cases are the same.
The streamwise x1and spanwise x3directions are periodical situations to prevent the correlations in the above two directions from being interrupted, the lengths suggested by Kazishima[19] in the above two directions are about 1382.4+=(3.848l2) and 691.2+=(1.92l2), respectively.
Influences of the wall effect in a turbulent flow are apparent, then the density of grid distribution in the x2direction near the wall is increased by a curved linear coordinates transformation method and the related equation indicating the position of the grid is shown in Equation (53).
x2=h (2+ ) +1 −1 ( −)/(1−) +2− (2+1) 1+ +1 −1 ( −)/(1−) (53)
In which =0.5 means the densities of grid distribution near both top and bottom walls are equal. The more the magnitude of closes to 1, the more the difference of the intervals of grids is remarkable, and =1.1 is assigned. As a result, the smallest magnitude of grid x2+ which is closest to the channel wall is equal to 0.84, and the largest magnitude of gridx2+ is located at the center and equal to 9.39. As for the other two directions of x1 and x3, the uniform grid distributions are selected and the magnitudes of the two intervals arex1+=21.6 and x3+=10.8, respectively. The figure of grid distribution in x1and x2 directions is shown in Figure 4, and the grid distributions in x1, x2, and x3directions are 64×64×64. The results of the situation of grid distributions of 40×40×40 were executed, but the results of the situation of grid distributions of 64×64×64 are more accurate than those of 40×40×40 situation. Then the results of the grid distributions of 64×64×64 are indicated solely. According to numerical tests, in order to maintain the designed turbulent flow to flow continuously, computational time steps of an eddy turn over time need to be more than 400 steps approximately. Therefore, the computational time intervals t of 5×10−4s and 1×10−7s are selected for A and B cases, respectively.
The initial conditions of velocities of u1, u2 and u3 used in the computation are shown in Equation (54). u1= ¯u1+u1,max×× R u2= u1,max×× R u3= u1,max×× R (54) x2 x1
2 x
1 x
Figure 5. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=30%, te=0 of case A.
x2
x1
Figure 6. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=30%, te=2 of case A.
x2
x1
Figure 7. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=30%, te=4 of case A.
In which, a random magnitude R is selected by a random number with normal distribution, and
is an intensity parameter and varies from 0–100%. u1,maxis the maximum velocity of the mean velocity distribution. Other related parameters used in cases A and B are tabulated in Tables II and III, respectively.
Phenomena shown in Figures 5–7 are the developed history of case A and the intensity parameter of is 30%. The gray scale indicates the magnitude of the resultant velocity Vr=(u21+u22+u23)1/2. In Figure 5, the initial condition that causes the flow to be irregular is just imposed on the original mean flow. When two eddy turnover times pass away, the flow showed in Figure 6 affected by the initial condition becomes a wavy flow gradually. Figure 7 shows that, after four large eddy times pass, the wavy motions of the flow almost disappear and the flow changes into a well-arranged
x2
x1
Figure 8. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=40%, te=0 of case A.
x2
x1
Figure 9. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=40%, te=2 of case A.
x1 x2
Figure 10. The distribution of resultant velocity Vrat the central cross section of,
x1x2 plane for=40%, te=4 of case A.
flow in which the magnitudes of velocities of the streamwise direction vary smoothly from small to large in the x2direction. The flow is similar to a laminar flow. These phenomena indicate that the flow will decay to an orderly flow when the intensity parameter mentioned above is not large enough.
Figures 8–10 show that the intensity parameter is enlarged and equal to 40%. At the initial stage, the fluctuations of the flow shown in Figure 8 are more remarkable than those of Figure 5. In Figure 9, after completion of two eddy turnover times, the flow in the central part becomes a wavy flow, and near both channel walls turbulent structures are observed gradually. Figure 10 shows that, four eddy turnover times pass away, various magnitudes of velocities distributed in the whole region are found and the flow almost becomes a fully turbulent flow.
x1 x3
x2
x3+
u1(m / s)
Figure 11. The distributions of velocity intensities of u1 for=40% situation of case A.
In Figure 11, the distributions of velocity intensities of u1 in three different vertical planes are shown. On the x1x3 plane, the intervals of streaks are about 100+ which is consistent with the results of Kazishima[19] under Re=180 situation. On the x2x3 plane, several shapes of mushrooms are observed. On the x1x2 plane, eddies generated from the channel walls are qualitatively observed.
The theoretical[22] and present results of mean velocity profiles close to the channel wall are indicated in Figure 12. Within x2+<10, five computational grids are distributed, then the results calculated by this study have good agreements with the results of the DNS method conducted by Kim et al.[17]. The root-mean-square (rms) values of fluctuating velocities of u1,rms, u2,rms, and
u3,rms and the mean Reynolds stress<u1u2> of Kim et al. [17], and the present study are shown
in Figures 13 and 14, respectively. Except for slight differences between the present results and the results of Kim et al.[17] in the <u1> distribution, the other three results of the present study and Kim et al.[17] are well consistent.
In case B, the mean flow velocity is increased and equal to a Mach number of 0.5. The influences of compressibility of fluid need to be taken into consideration, however, according to the results of Lenormand et al.[18]. The results of the mean Reynolds stresses and turbulent intensities are close to those of the incompressible flow as the mean flow velocity is lower than a Mach number of 0.6. As a result, the rms values of fluctuating velocities of u1,rms, u2,rms, and u3,rms and the mean Reynolds stress of<u1u2> of the Kim et al. [17] are used to validate these of the present results
obtained from 0.5 Mach number situation, and both the results are shown in Figures 15 and 16. The results are well consistent with the corresponding result of Kim et al.[17]. The property of the compressibility of fluid may be a reason to cause a slight difference between both the results at the summit magnitude.
log x2+ u < u1 > 0 101 102 5 10 15 theoretical results present results of case A 20
Figure 12. The distributions of mean velocities of theoretical and present results of case A.
u ′1,rms ,u ′2,rms ,u ′3,rms u 2 0.5 x h Kim et al. [17] present Results of case A 0 0 0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1 u′1,rms u′2,rms u′3,rms
0.5 x h < u′1 u′2 > 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Kim et al. [17]
present results of case A
0.2 0.4 0.6 0.8 1
< u′1u′2 >
Figure 14. The distributions of mean Reynolds stress.
present Results of case B Kim et al. [17] u ′1 ,rms ,u ′2,rms ,u ′3,rms u u′1,rms u′2,rms u′3,rms 3 2.5 2 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 2 0.5 x h
Figure 15. The distributions of root-mean-square values of fluctuation velocities.
The compressible fluids are compressed due to being in a high-speed situation and the mean temperatures of fluids are raised and higher than the temperature of the wall, which is constant and equal to 298.05 k. In Figure 17, accompanying with the increment of x2/0.5h, the dimensionless
< u′1 u′2 > < u′1u′2 > 2 0.5 x h 0.9 Kim et al. [17]
present results of case B 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1
Figure 16. The distributions of mean Reynolds stress.
Tw < T > < T > 2 0.5 x h 1.1 1.05 1 0.95 0 0.2 0.4 0.6 0.8 1 Lenormand et al. [18] present results of case B
Figure 17. Mean temperature profile.
mean temperature< T >/Tw becomes larger gradually. As the magnitude of x2/0.5h is over 0.2, the distribution of < T >/Tw becomes flat, and this trend is reasonable because its position is far from the wall. The density of fluid assumed as an ideal gas is inversed to the temperature of
2 0.5 x h 1.1 1.05 1 0.95 0 0.2 0.4 Lenormand et al. [18] present results of case B
0.6 0.8 1
Figure 18. Mean density profile.
x2 x3 x1
Figure 19. The instant temperatures of case B.
the fluid. The longer the distance the fluid is from the channel wall, the smaller the density of the fluid becomes, which is shown in Figure 18. Both results shown in the above figures have good agreements with the results of Lenormand et al.[18]. Usage of a gray scale to indicate the variations of the fluid temperatures at a certain instant is shown in Figure 19. Corresponding to the appearances of streaks in the flow field indicated in Figure 11, the alternate arrangements of the low and high temperatures of fluids caused by the compressible turbulent flow are clearly
observed on the cross section very close to the channel wall of x1x3plane. Except for the near-wall regions, the fluids of which the temperatures are comparatively high occupy other regions of the channel.
5. CONCLUSIONS
Fully compressible turbulent channel flows obtained by the methods of Roe, preconditioning, and LES are investigated. In order to stabilize the transient situation and increase the number of computing grids near the channel walls, the methods of dual time stepping and curved linear coordinates transformation are adopted. A thorny problem of assignation of initial condition is successfully overcome by the imposition of random velocity fluctuations on the mean velocity distribution. The results show that the methods mentioned above effectively improve convergent efficiency of the above study. Achievements of economizing computing time and improving accu-racy are gained due to usage of a curved linear coordinates transformation. Under high-speed situations, the increments of the temperatures of working fluids are remarkable and well consistent with the previous results.
ACKNOWLEDGEMENTS
The authors gratefully acknowledge to the National Center for High-Performance Computing of Taiwan, R.O.C., Dr Trong T. Bui of Dryden Flight Research Center of Edwards, California, the Professor of National Chiao Tunge University, Wu-Ting Tsai and the Professor of National Taiwan University, Mei-Jiau Huang.
REFERENCES
1. Tong TB. A parallel finite-volume algorithm for large-eddy simulation of turbulence flows. Computers and Fluids 2000; 29:877–915.
2. Gavrilakis S. Numerical simulation of low-Reynolds-number turbulent flow through a straight square duct. Journal
of Fluid Mechanics 1992; 244:101–129.
3. Stolz S, Adams NA, Kleiser L. An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows. Physics of Fluids 2001; 13:997–1015.
4. Lund TS, Wu XH, Squires KD. Generation of turbulent inflow data for spatially developing boundary layer simulations. Journal of Computational Physics 1998; 140:233–258.
5. Klein M, Sadiki A, Janicka J. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulation. Journal of Computational Physics 2003; 186:652–665.
6. Briley WR, McDonald H, Shamroth SJ. At low Mach number Euler formulation and application to time iterative LBI schemes. AIAA Journal 1983; 21(10):1467–1469.
7. Turkel E. Preconditioned methods for solving the incompressible and low speed compressible equations. Journal
of Computational Physics 1987; 72:277–298.
8. Choi D, Merkel CL. Application of time-iterative schemes to incompressible flow. AIAA Journal 1985; 25(6): 1518–1524.
9. Choi D, Merkel CL. The application of preconditioning in viscous flows. Journal of Computational Physics 1993; 105(1193):207–223.
10. Roe PL. Approximation Riemann solver parameter vectors and difference schemes. Journal of Computational
Physics 1981; 43:357–372.
11. Weiss JM, Simth WA. Preconditioning applied to variable and constants density flows. AIAA Journal 1995; 33:2050–2506.
12. Dennis J, Thomas P, Pieter B. Recent enhancements to OVERFLOW. Thirty Fifth Aerospace Sciences Meeting
and Exhibit, Reno, NV, 1997.
13. Shishir AP, Sankaran V, Thomas HP. Implementation of preconditioned dual-time procedures in OVERFLOW.
Forty First Aerospace Sciences Meeting and Exhibit, Reno, NV, 2003.
14. Lessani B, Ramboer J, Lacor C. Efficient large-eddy simulations of low Mach number flows using preconditioning and multigrid. Journal of Computational Physics 2002; 18:221–233.
15. Xu XF, Lee JS, Pletcher RH. A compressible finite volume formulation for large eddy simulation of turbulent pipe flows at low Mach number in Cartesian coordinates. Journal of Computational Physics 2005; 203:22–48. 16. Alkishriwi N, Meinke M, Schroder W. A large-eddy simulation method for low Mach number flows using
precondition and multigrid. Computers and Fluids 2005; 35:1126–1136.
17. Kim J, Moin P, Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number. Journal
of Fluid Mechanics 1987; 177:133–166.
18. Lenormand E, Saguat P, Phuoc LT. Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number. International Journal for Numerical Methods in Fluids 2000; 32:369–406.
19. Kazishima T. Numerical Simulation of Turbulent Flows. Yokendo Co., Ltd., 1999; 149–155.
20. Yoon S, Jameson S. Lower-upper symmetric-Gauss–Seidel method for the Euler and Navier–Stokes equations.
AIAA Journal 1988; 26:1025–1026.
21. Abalakin I, Dervieux A, Kozubskaya T. A vertex centered high order MUSCL scheme applying to linearized Euler acoustics. INRIA No4459, 2002.
22. Durbin PA, Petterson Reif BA. Statistical Theory and Modeling for Turbulent Flows. Wiley: New York, 2001; 63.
