- Email: [email protected]
Computational study of incompressible turbulent flows with method of characteristics
Journal of Computational and Applied Mathematics (
)
–
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Computational study of incompressible turbulent flows with method of characteristics Ali Atashbar Orang a,∗ , Seyed Esmail Razavi b , Hamoon Pourmirzaagha c a
Department of Mechanical and Aerospace Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
b
Faculty of Mechanical Engineering, Tabriz University, Tabriz, Iran
c
Department of Mechanical and Aerospace Engineering, Ramsar Branch, Islamic Azad University, Ramsar, Iran
highlights • • • • •
We developed a new characteristic-based scheme for incompressible turbulent flows. The Navier–Stokes partial differential equations were solved numerically. Numerical tests were conducted to flow past a circular cylinder. High stability range was achieved which led to fast convergence. Accurate results were obtained by applying the new CB method.
article
info
Article history: Received 26 December 2012 Received in revised form 5 May 2013 Keywords: Incompressible turbulent flow Navier–Stokes equations RANS Spalart–Allmaras turbulence model Method of characteristics Upwind scheme
abstract The current study presents numerical investigation of incompressible turbulent flow. Momentum and continuity equations are coupled by the Artificial Compressibility Method (ACM). The one-equation Spalart–Allmaras turbulence model along with RANS (Reynolds Averaged Navier–Stokes) equations is used to calculate Reynolds stresses in 2-D problems. For convective fluxes a characteristic-based (CB) finite-volume solution has been developed. Also, a revised Jameson averaging (AVR) method was implemented. In comparison, the proposed CB upwind scheme, demonstrated very accurate results, high stability and faster convergence. For time discretization, the fifth-order Runge–Kutta scheme is used. As the convergence acceleration techniques, the local time stepping and implicit residual smoothing were applied. Considering flow behavior, suitable boundary conditions have been implemented. Cross flow around a horizontal circular cylinder at high Reynolds numbers has been studied as a benchmark. The results obtained using the new proposed scheme, are in good agreement with the standard available solutions in the literature. © 2013 Elsevier B.V. All rights reserved.
1. Introduction The most encountered incompressible flows in many industrial cases include fully-developed turbulence. Majority of numerical methods have been developed for compressible flows with decreasing Mach number which because of stiffness of the convective terms in incompressible flows cause convergence and accuracy problems, as mentioned by Volpe [1]. This increases sensitivity to boundary conditions as reported in the works of Gresho [2]. The ACM of Chorin [3] as a pre-conditioner enables the coupling of continuity and momentum equations thus allowing time-marching compressible schemes to be applied for incompressible ones. It has been shown that the ACM requires a fast-converging scheme to satisfy the incompressibility conditions [4]. The CB method for ACM, was proposed by Drikakis et al. [5] based on the definition of primitive
∗
Corresponding author. E-mail address: [email protected] (A. Atashbar Orang).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.05.019
2
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Nomenclature Cb1 , Cb2 , Cv 1 , Cv 2 , σ , κ, Cw1 , Cw2 , Cw3 Constants of the standard Spalart–Allmaras model Ct1 , Ct2 , Ct3 , Ct4
CFL FN p Re R S
Ω′ u
v
U∞ uN uP W
Courant number Normal convective flux vector along the x direction Pressure Reynolds number Viscous flux vector along the x direction Viscous flux vector along the y direction Area of secondary cells Velocity along the x direction Velocity along the y direction Reference velocity Normal to face velocity Parallel to face velocity Initial variable vector
Greek symbols
α λ β ν˜ ∆t νT ρ φ Ω
Pseudo-sound velocity Wave propagation speed Artificial compressibility Eddy-viscosity variable Time step Turbulent viscosity Density Symbol for initial stream variables Initial cells’ area
variables as functions of their values on the characteristics and formation of compatibility equations. Recently Neofytou [6] proposed revisions in the derivation of the compatibility equations which are used to derive the CB formulas for the calculations of the primitive flow variables within the iterative process which is directly related to the convergence rate. The above references show that the CB scheme along with AC was successfully used in a large number of applications [7]. Recently Razavi and Zamzamian [8] have proposed a multidimensional characteristic-based scheme coupled with ACM. This scheme has been successful in reducing the convergence steps and giving accurate results. Other CB method is the Roe scheme which was originally developed for simulation of compressible Euler equations. Razavi and Atashbar extended the Roe scheme to estimate convective fluxes of incompressible viscous flows past a hydrofoil. This method is based on the propagations of virtual acoustic waves which enables damping the numerical oscillations [9]. Nithiarasu and Liu [10] implemented artificial compressibility to solve incompressible turbulent flows at medium Reynolds numbers with characteristic-based scheme. It has been used for solutions of steady and unsteady flows. They used three different turbulence models namely Wolfenstein l, Spalart–Allmaras and standard k −ε . They showed that the artificial compressibility based on CB is suitable for solutions of both steady and unsteady incompressible turbulent flows at medium Reynolds numbers. The finite-volume method along with hybrid grid and k − ε turbulence model was proposed by Feng and Shanhong [11]. Hybrid grid can be used for accelerating convergence in the unsteady time step. Yang and Chang [12] used 2-D numerical simulation of turbulent flow fields for the NREL S809 airfoil. They assumed the flow as steady, turbulent and incompressible. The turbulent flow by the revised k −ε model was adopted to calculate effects of low Reynolds numbers near the wall. Marx [13] implemented the implicit MUSCL scheme to calculate 3-D compressible Navier–Stokes equations. They also described turbulence models for separated flows with emphasis on the non-equilibrium numerical model. Ng et al. [14] developed a finite volume method based on vortex flux flow for compressible flows. In their method, the second-order central finite difference and Runge–Kutta schemes are used for flux and time discretizations, respectively. Pentaris et al. [15] presented predictions for 2-D, steady incompressible flows in laminar and turbulent regimes by the standard k − ε model. Mass conservation is coupled to other equations by artificial compressibility and pressure correction methods. They also used the second- and the fourth-order artificial dissipation to obtain good convergence and efficiency of k − ε equations. Erturk et al. [16] used stream function–vorticity equations for steady incompressible Navier–Stokes equations. They solved steady flow inside a driven cavity up to Reynolds numbers 21 000 using a 601 × 601 fine grid. Their formulation proved to be stable and effective at very high Reynolds numbers and non-orthogonal grids. In addition LES1 modeling was done
1 Large Eddy Simulation
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
3
by Catalano et al. [17] to simulate flows past a circular cylinder in the supercritical regime by the finite difference method. Elmiligui et al. [18] implemented two multiscale-type turbulence models in the PAB3D solver. The models are based on modifying RANS. Recently, Atashbar and Razavi [19] developed the Roe-like scheme for incompressible turbulent flows which is advantageous in terms of accuracy and convergence. In this study, a finite volume model is presented for incompressible turbulent flows past a circular cylinder. Considering Spalart–Allmaras as a turbulence model the CB method is used in formulations. The behavior and performance of the Spalart–Allmaras turbulence model are investigated. The numerical results are validated by other available ones. 2. Turbulence model equation Utilizing the Spalart–Allmaras model, turbulent flows are successfully calculated over a wide range of regimes using RANS equations. The Spalart–Allmaras model is consisted of transport equations for an eddy-viscosity variable ν˜ following [20]:
∂ ν˜ ∂ ν˜ ∂ ∂ ∂ ν˜ ∂ ν˜ 1 + + Cb2 ν˜ vj = Cb1 (1 − ft2 )S˜ ν˜ + (νL + ν˜ ) ∂t ∂ xj σ ∂ xj ∂ xj ∂ xj ∂ xj 2 Cb1 ν˜ − Cw1 fw − 2 ft2 + ft1 ∥∆ν˜ ∥22 (1) κ d where νL denotes the kinematic viscosity and d is the distance to the closest wall. In the above transport equation, the first term on the right hand side is the eddy-viscosity production term. The production term S˜ is evaluated as follows: χ3 χ −3 ˜S = fv3 S + ν˜ fv2 , fv 1 = , , fv 2 = 1 + κ 2 d2 Cv 2 χ 3 + Cv31 (2) ν˜ (1 + χ fv1 ) (1 − fv2 ) fv 3 = , χ= . Max (χ , 0.001) νL In Eq. (2) S is the magnitude of the mean rotation rate: S=
2ωij ωij
ωij =
1
2
∂vi ∂vj − ∂ xj ∂ xi
.
In Eq. (1) the second and third terms are conservative diffusion and non-conservative diffusion respectively. Diffusion terms 1
σ
∂ ∂ ν˜ ∂ ν˜ ∂ ν˜ + Cb2 (νL + ν˜ ) ∂ xj ∂ xj ∂ xj ∂ xj
can be expressed as follows: 1 + Cb2 ∂
σ
∂ xj
∂ ν˜ Cb2 − (νL + ν˜ ) (νL + ν˜ ) ∇ 2 ν˜ . ∂ xj σ
Then difficulties with discretization of nonlinear term
(3)
∂ ν˜ ∂ ν˜ ∂ xj ∂ xj
are circumvented [20]. The fourth term in Eq. (1) denotes
near-wall turbulence destruction which reads
1 + Cw6 3
16
g6
6
+ Cw 3
,
g = r + Cw 2 r 6 − r , r =
ν˜
. (4) S˜ κ 2 d2 Finally the fifth and sixth terms are transition damping the production and transition source of turbulence. Those transition functions are given by fw = g
ωt2 2 2 2 ft1 = Ct1 gt exp −Ct2 d + gt dt , ∆U 2 gt = Min (0.1, ∥∆v ⃗ ∥2 / (ωt ∆xt ))
ft2 = Ct3 exp −Ct4 χ 2
(5)
where dt is the distance from the nearest trip point, ωt is the vorticity at the wall at the trip point, ∆U is the 2-norm of the difference between the velocity at the trip (zero if the wall stationary) and ∆xt stands for the spacing along the wall at the trip point. Since the present study assumes fully turbulent flow, this term is not used. Therefore, the values of ft1 and ft2 which are associated with these terms are taken as zero. Also, the turbulent kinematic eddy viscosity is obtained via
νT = fv1 ν˜ .
(6)
The various constants of this model have the following values: Cb1 = 0.1355,
σ = 2/3, Cw2 = 0.3,
Cb2 = 0.622,
κ = 0.41, Cw3 = 2,
Cv 1 = 7.1,
Cv 2 = 5,
Cw1 = Cb1 /κ + (1 + Cb2 ) /σ , 2
Ct1 = 1,
Ct2 = 2,
Ct3 = 1.3,
(7) Ct4 = 0.5.
The above model is valid through the assumption of homogeneous turbulent flow.
4
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 1. Schematic representation of the primary cell to evaluate normal convective fluxes.
3. Governing equations The non-dimensional form of Navier–Stokes equations plus the eddy-viscosity transport equation for incompressible turbulent flow in the finite-volume form with the ACM read:
∂ ∂t
WdA +
FN dL =
∂Ω
Ω
∂Ω
(Rdy − Sdx) + Z
(8) 0
p u
W = , v
ν˜
0 0 0
Z=
Cb1 S˜ ν˜ +
Cb2
σ
∂ ν˜ ∂ ν˜ ∂ xj ∂ xj
∂u (1 + νT ) ∂x 1 ∂v R= , ReD (1 + νT ) ∂x 1 ∂ ν˜ (1 + ν˜ ) σ ∂x
β uN uu + p cos ζ FN = N , v uN + p sin ζ ν˜ uN
/ReD − Cw1 fw
2 ν˜ d
/ReD
0
∂ u (1 + νT ) ∂y 1 S= (1 + ν ) ∂v , T ReD ∂y ∂ ν˜ 1 (1 + ν˜ ) σ ∂y
(9)
where, x=
X D
,
y=
Y D
,
u=
U U∞
,
v=
V U∞
,
∆y ∆x , cos ζ = , sin ζ = − , ν ∆L ∆L − → − → − → ∆L = ∆x2 + ∆y2 , N = cos ζ i + sin ζ j , up = v cos ζ − u sin ζ , a = u2N + β ReD =
t =
τ U∞ D
,
p=
P − P∞ 2 ρ∞ U∞
,
U∞ D
(10) uN = u cos ζ + v sin ζ ,
− →
where, β, uN , up , D, ∆L, a, N are the artificial compressibility parameter, the normal and parallel components of the velocity, the diameter of the cylinder, the length, the pseudo-sound speed and the normal unit vector or the direction of virtual wave propagation respectively as shown in Fig. 1. 3.1. Convective fluxes Two schemes were used to evaluate convective fluxes at cell interfaces. Considering Fig. 1, one has
∂Ω
FN dL =
4
(K )
FN ∆L
(11)
K =1
where, K is attributed to each cell side. Two schemes were applied to convective fluxes (FN ) namely Averaging (AVR) and characteristic-based (CB) schemes to compare and show the superiority of the introduced CB scheme. Considering Fig. 1 at primary cell interfaces like AB in Fig. 1, these schemes are expressed as follows.
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
5
Fig. 2. τ − ξ co-ordinate.
3.1.1. Averaging scheme (AVR) In this scheme with pressure-based modification we have: W = ηWL + (1 − η)WR ;
η=
pL
(12)
pL + pR
where, R and L indicate right and left cells respectively by p as the driving force [9]. The first and higher accuracy can be achieved in Eq. (12). To overcome numerical oscillations, artificial dissipation was implemented. 3.1.2. CB scheme The original CB method, defines primitive variables as functions of corresponding values along the virtual characteristics. Considering Fig. 2, where τ is the pseudo time and ξ is the direction of the outward normal vector of the surface of the cell, the flow variable vector at pseudo time n + 1 can be calculated along a characteristic, k, using Taylor series expansion and the initial value at pseudo time level n: W = Wk + Wξ ξτ 1τ + Wτ1τ and Wτ =
W − Wk
1τ
− Wξ ξτ .
(13)
Substituting Wτ into the inviscid form of Eq. (8), after some algebraic manipulations, leads to the following compatibility equations along characteristics:
λ(1) = uN ,
λ(2) = uN + a,
λ(3) = uN − a
(u − u1 ) sin ζ − (v − v1 ) cos ζ = 0 for λ = λ1
(14)
p − p2 = −λ2 [(u − u2 ) cos ζ − (v − v2 ) sin ζ ] for λ = λ2
(15)
p − p3 = −λ3 [(u − u3 ) cos ζ − (v − v3 ) sin ζ ] for λ = λ3 .
(16)
Then u, v and p at cell interfaces are determined using preceding equations as [5]: u = Z cos ζ + sin ζ (u1 sin ζ − v1 cos ζ )
(17)
v = Z sin ζ + cos ζ (v1 cos ζ − u1 sin ζ )
(18)
p=
λ2 [p3 + λ3 (u3 cos ζ + v3 sin ζ )] − λ3 [p2 + λ2 (u2 cos ζ + v2 sin ζ )] 2a
(19)
where Z =
1 2a
[(p2 − p3 ) + cos ζ (λ2 u2 − λ3 u3 ) + sin ζ (λ2 v2 − λ3 v3 )].
The variables at n pseudo time level are approximately evaluated based on the signs of the characteristics: Wj =
1 2
[(1 + sign(λj ))WL + (1 − sign(λj ))WR ]
WLi+ 1 = Wi , 2
WRi+ 1 = Wi+1 . 2
(20) (21)
6
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 3. Schematic representation of primary and secondary cells to evaluate viscous terms.
For the eddy-viscosity of the Spalart–Allmaras model considering Eq. (12), we have:
ν˜ = ην˜ L + (1 − η)˜νR .
(22)
3.2. Viscous fluxes Calculation of viscous fluxes on the cell boundary requires the second-order derivatives of flow parameters within cell sides. The calculation method is illustrated in Fig. 3. Calculation of viscous flux vectors is similar to estimation of convective vectors as follows:
∂Ω
(Rdy − Sdx) =
∂Ω
∂φ ∂x
dy −
∂φ ∂y
dx =
4 ∂φ K =1
∂x
∆yK −
∂φ ∂y
∆xK
(23)
where, ϕ stands for u, v, ν˜ and K is attributed to each cell side, namely AB, BC, CD and DA. For example for face AB one has
∂φ ∂x
= AB
1
Ω′
∂Ω′
φ dy =
4 1
Ω′
φK ∆yK =
K =1
1
Ω′
[φMB ∆MB + φBN ∆BN + φNA ∆NA + φAM ∆AM ].
(24)
For example one can take the following averaged values:
ϕM =
1 4
(ϕA + ϕB + ϕC + ϕD ),
ϕB =
1 4
(ϕN + ϕNE + ϕE + ϕM ) and ϕMB =
1 2
(ϕM + ϕB ).
(25)
4. Time discretization and boundary conditions 4.1. Time discretization The fifth-order Runge–Kutta (RK5) schemes are considered because of wider stability range. (CFL ≤ 3)
∂U + Q = 0, ∂t
U(n+1) = U(n) − εm ∆tQ(n) (U),
εm = 1/4, 1/6, 3/8, 1/2, 1 · (m = 1, . . . .5)
(26)
where Q includes the discretized form of integral convective and viscous terms. The number of iterations to achieve the steady state is reduced by using a variable time step that depends on the local flow changes and the grid spacing as [9]
∆x2i.j + ∆y2i.j CFL ∆t = Max , Max (γ ) 2 2 n n Max a + ui,j + vi,j γn = Max λ(1) , λ(2) , λ(3) , λ(4) ∆x2i.j + ∆y2i.j , CFL × Min
(27) n = 1, 2, 3, 4.
To accelerate convergence the local time stepping and implicit residual smoothing were utilized. 4.2. Boundary conditions The explicit boundary conditions are used throughout the computational domain. The consistent boundary treatment ensures the disturbance dissipation in the discretized domain without reflection. At the inlet boundary, all of the flow
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
7
Fig. 4. A part of grid around a circular cylinder.
Fig. 5. Mean pressure distribution on the cylinder: (—) CB, characteristic-based at ReD = 6.7E5; (•) experiment by Falchsbart [17] (in Zdravkovich, 1997) at ReD = 6.7E5; (H) experiment by Warschauer and Leene [21] (1971) at ReD = 1.2E6; (- - -} LES by Catalano et al. [17] at ReD = 1E6.
parameters except pressure are set to free stream values, whereas the pressure is extrapolated from the interior domain. The free stream value of the turbulent eddy-viscosity is usually taken as ν˜ = 0.1 [20]. At the outlet boundary, only the pressure equals its reference value in free stream. Other flow parameters are extrapolated from the interior domain. At the solid boundary considering the no-slip condition zero mass flux through the surface is prescribed. The pressure on the solid surface is found by extrapolation. It was set ν˜ = 0 at the solid boundary and hence νT = 0 [20]. 5. Numerical results and discussion Numerical tests were performed with fully turbulent flows for cross flow around a horizontal circular cylinder by two schemes: averaging (AVR) and characteristic-based (CB). The grid consisted of 100 × 100 cells along both radial and circumferential directions and the outer boundary is located 15 diameters away from the cylinder center as shown in Fig. 4. In Fig. 5 the wall pressure distribution at ReD = 6.7E5 by the CB scheme is compared to experimental data of Falchsbart [17] at ReD = 6.7E5 and Warschauer and Leene [21] at ReD = 1.2E6 and also to the LES model by Catalano et al. at ReD = 1E6. The maximum pressure is located at the stagnation point where θ = 0°. It can be observed that near 140° where separation occurs, the wall pressure falls down.
8
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 6. Comparison of convergence histories: (a) ReD = 6.7E5; (b) ReD = 1.2E6. CB, characteristic-based; AVR, averaging.
Fig. 7. Mean streamwise velocity distribution at ReD = 1E6 predicted: (a) AVR, averaging; (b) CB, characteristic-based; (c) LES and URANS by Catalano et al. [17].
In Fig. 6, convergence histories, for ReD = 6.7E5 and ReD = 1.2E6 have been shown. Because of numerical oscillations in AVR, artificial dissipation is unavoidable. It is shown that the CB scheme resulted in a faster convergence which demonstrates its appropriate capability to handle incompressible flows. In Fig. 7, u-velocity contours for ReD = 1E6 are presented and compared with that of Catalano et al. by LES and URANS.2 As shown, the CB method has smooth iso-velocity lines but in the AVR method, the results were polluted by wiggles. In Fig. 8 considering other RANS solutions, the wall pressure distribution at ReD = 1.4E5 is compared to that of Elmiligui et al. [18] in which the k − ε turbulence model has been used and validated with experimental data of Roshko’s [22]. Skin friction distribution for ReD = 6.7E5 utilizing the CB scheme is compared to that of experimental data by Achenbach [22] and LES by Catalano et al. in Fig. 9. At the separation area the skin friction decreases because of the low velocity gradient versus the vertical distance from the wall. Substantial agreement between the results of CB and experimental data can be noticed. In Fig. 10, the drag coefficient versus ReD is plotted. It can be seen that at higher Reynolds numbers, the regime is changed to turbulence. Hence, the drag falls down because of separation delay and less pressure difference between upstream and downstream. However, in this situation both the drag and skin friction grow as a result of high velocity gradient.
2 Unsteady Reynolds Averaged Navier–Stokes
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
9
Fig. 8. Mean pressure distribution on the cylinder at ReD = 1.4E5: (—) CB, characteristic-based; (N) k − ε by Elmiligui et al. [18]; (◦) experiment by Roshko [22].
Fig. 9. Skin friction distribution on the cylinder: (—) CB, characteristic-based at ReD = 6.7E5; (◦) experiments by Achenbach [23] (1968) at ReD = 3.62E6; (- - -) LES by Catalano et al. [17] at ReD = 1E6.
In Fig. 11 streamlines for ReD = 6.7E5 and ReD = 1.2E6 with the CB method are presented. As shown, in the wake region vortices are appeared and by increasing ReD the separation point moves to downstream. The vector plot of downstream by the CB method at ReD = 6.7E5 is demonstrated in Fig. 12. 6. Conclusion and outlook This paper presents 2-D incompressible turbulent flow by the finite volume method and artificial compressibility. Utilizing RANS, the one equation Spalart–Allmaras turbulence model was implemented. Numerical tests were conducted to flow past a circular cylinder as a case study. Time discretization was done by the fifth-order Runge–Kutta scheme. Local time stepping and residual smoothing were utilized to accelerate the convergence process. In the convective flux treatment, the propagation of information along characteristics was taken into account. From the mathematical point of view,
10
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 10. Drag coefficient as a function of the ReD : (—) Achenbach [23] (1968); (•) CB, characteristic-based; () LES, (△) URANS by Catalano et al. [17].
Fig. 11. Stream lines: (a) ReD = 6.7E5; (b) ReD = 1.2E6. CB, characteristic-based.
Fig. 12. Vector plot: ReD = 6.7E5. CB, characteristic-based.
the characteristic-based scheme was developed which shows the true wave paths. This scheme which is based on the propagation of information along characteristics, allows relatively high CFL numbers and has high convergence rate with less numerical oscillations without any requirement to damping. The substantial consensus was achieved utilizing the CB method in comparison with that of others. Application to three-dimensional and unsteady turbulent flows is recommended as future extensions of the present scheme.
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
11
References [1] G. Volpe, Performance of compressible flow codes at low Mach numbers, AIAA Journal 31 (1) (1993) 46–56. [2] P.M. Gresho, Some interesting issues in incompressible fluid dynamics, both in the continuum and in numerical simulation, Journal of Advances in Applied Mechanics 28 (1992) 45–140. [3] A.J. Chorin, A numerical method for solving incompressible viscous problems, Journal of Computational Physics 2 (1967) 12–26. [4] D. Kwak, C. Kiris, C.S. Kim, Computational challenges of viscous incompressible flows, Journal of Computational Fluids 34 (2004) 283–299. [5] D. Drikakis, P.A. Govatsos, D.E. Papatonis, A characteristic based method for incompressible flows, International Journal of Numerical Methods in Fluids 19 (1994) 667–685. [6] P. Neofytou, Revision of the characteristic-based scheme for incompressible flows, Journal of Computational Physics 222 (2007) 475–484. [7] X. Su, Y. Zhao, X. Huang, On the characteristics-based ACM for incompressible flows, Journal of Computational Physics 227 (2007) 1–11. [8] S.E. Razavi, K. Zamzamian, A. Farzadi, Genuinely multidimensional characteristic-based scheme for incompressible flows, International Journal of Numerical Methods in Fluids 57 (2008) 929–949. [9] S.E. Razavi, A. Atashbar Orang, A comparative investigation of hydrofoil at angles of attack, International Journal of Numerical Methods in Fluids 68 (2012) 1087–1101. [10] P. Nithiarasu, C.B. Liu, An artificial compressibility based characteristic based split (CBS) scheme for steady and unsteady turbulent incompressible flows, Computational Methods in Applied Mechanics and Engineering 195 (2005) 2961–2982. [11] L. Feng, J. Shanhong, Unsteady flow calculations with a multigrid Navier–Stokes method, AIAA Journal 34 (1996) 2047–2053. [12] S.L. Yang, Y.L. Chang, O. Arici, Navier–Stokes computations of the NREL airfoil using k − ω turbulent model at high angles of attack, ASME 117 (1995) 304–310. [13] Y.P. Marx, A practical implementation of turbulence models for the computation of three-dimensional separated flows, International Journal of Numerical Methods in Fluids 13 (1991) 775–796. [14] K.C. Ng, M.Z. Yusoff, T.F. Yusaf, Simulations of two-dimensional high speed turbulent compressible flow in a diffuser and a nozzle blade cascade, American Journal of Applied Sciences 9 (2005) 1325–1330. [15] A. Pentaris, K. Nikolados, S. Tsangaris, Development of projection and artificial compressibility methodologies using the approximate factorization technique, International Journal of Numerical Methods in Fluids 19 (1994) 1013–1018. [16] E. Erturk, T.C. Corke, C. Gokcol, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, International Journal of Numerical Methods in Fluids 48 (2005) 747–774. [17] P. Catalano, M. Wang, G. Iaccarino, P. Moin, Numerical simulation of the flow around a circular cylinder at high Reynolds numbers, International Journal of Heat and Fluid Flow 24 (2003) 463–469. [18] A. Elmiligui, Kh.S. Abdol-Hamid, S.J. Massey, S.P. Pao, Numerical study of flow past circular cylinder using hybrid turbulence formulations, Journal of Aircraft 47 (2010) 434–440. [19] A. Atashbar Orang, S.E. Razavi, Incompressible Roe-like scheme for turbulent flows, International Journal for Numerical Methods in Fluids 72 (4) (2013) 414–426. [20] J. Blazek, Computational Fluid Dynamics: Principles and Applications, 2005. [21] K.A. Warschauer, J.A. Leene, Experiments on mean and fluctuating pressures of circular cylinders at cross flow at very high Reynolds numbers, in: Proc. International Conference on Wind Effects on Buildings and Structures, Tokyo, Japan, 1971, pp. 305–315. [22] A. Roshko, Experimental on flow past a circular cylinder at very high Reynolds number, Journal of Fluid Mechanics 10 (1961) 345–356. [23] E. Achenbach, Distribution of local pressure and skin friction around a circular cylinder in cross-flow up to Re = 5 × 106 , Journal of Fluid Mechanics 34 (1968) 625–639.
)
–
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Computational study of incompressible turbulent flows with method of characteristics Ali Atashbar Orang a,∗ , Seyed Esmail Razavi b , Hamoon Pourmirzaagha c a
Department of Mechanical and Aerospace Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
b
Faculty of Mechanical Engineering, Tabriz University, Tabriz, Iran
c
Department of Mechanical and Aerospace Engineering, Ramsar Branch, Islamic Azad University, Ramsar, Iran
highlights • • • • •
We developed a new characteristic-based scheme for incompressible turbulent flows. The Navier–Stokes partial differential equations were solved numerically. Numerical tests were conducted to flow past a circular cylinder. High stability range was achieved which led to fast convergence. Accurate results were obtained by applying the new CB method.
article
info
Article history: Received 26 December 2012 Received in revised form 5 May 2013 Keywords: Incompressible turbulent flow Navier–Stokes equations RANS Spalart–Allmaras turbulence model Method of characteristics Upwind scheme
abstract The current study presents numerical investigation of incompressible turbulent flow. Momentum and continuity equations are coupled by the Artificial Compressibility Method (ACM). The one-equation Spalart–Allmaras turbulence model along with RANS (Reynolds Averaged Navier–Stokes) equations is used to calculate Reynolds stresses in 2-D problems. For convective fluxes a characteristic-based (CB) finite-volume solution has been developed. Also, a revised Jameson averaging (AVR) method was implemented. In comparison, the proposed CB upwind scheme, demonstrated very accurate results, high stability and faster convergence. For time discretization, the fifth-order Runge–Kutta scheme is used. As the convergence acceleration techniques, the local time stepping and implicit residual smoothing were applied. Considering flow behavior, suitable boundary conditions have been implemented. Cross flow around a horizontal circular cylinder at high Reynolds numbers has been studied as a benchmark. The results obtained using the new proposed scheme, are in good agreement with the standard available solutions in the literature. © 2013 Elsevier B.V. All rights reserved.
1. Introduction The most encountered incompressible flows in many industrial cases include fully-developed turbulence. Majority of numerical methods have been developed for compressible flows with decreasing Mach number which because of stiffness of the convective terms in incompressible flows cause convergence and accuracy problems, as mentioned by Volpe [1]. This increases sensitivity to boundary conditions as reported in the works of Gresho [2]. The ACM of Chorin [3] as a pre-conditioner enables the coupling of continuity and momentum equations thus allowing time-marching compressible schemes to be applied for incompressible ones. It has been shown that the ACM requires a fast-converging scheme to satisfy the incompressibility conditions [4]. The CB method for ACM, was proposed by Drikakis et al. [5] based on the definition of primitive
∗
Corresponding author. E-mail address: [email protected] (A. Atashbar Orang).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.05.019
2
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Nomenclature Cb1 , Cb2 , Cv 1 , Cv 2 , σ , κ, Cw1 , Cw2 , Cw3 Constants of the standard Spalart–Allmaras model Ct1 , Ct2 , Ct3 , Ct4
CFL FN p Re R S
Ω′ u
v
U∞ uN uP W
Courant number Normal convective flux vector along the x direction Pressure Reynolds number Viscous flux vector along the x direction Viscous flux vector along the y direction Area of secondary cells Velocity along the x direction Velocity along the y direction Reference velocity Normal to face velocity Parallel to face velocity Initial variable vector
Greek symbols
α λ β ν˜ ∆t νT ρ φ Ω
Pseudo-sound velocity Wave propagation speed Artificial compressibility Eddy-viscosity variable Time step Turbulent viscosity Density Symbol for initial stream variables Initial cells’ area
variables as functions of their values on the characteristics and formation of compatibility equations. Recently Neofytou [6] proposed revisions in the derivation of the compatibility equations which are used to derive the CB formulas for the calculations of the primitive flow variables within the iterative process which is directly related to the convergence rate. The above references show that the CB scheme along with AC was successfully used in a large number of applications [7]. Recently Razavi and Zamzamian [8] have proposed a multidimensional characteristic-based scheme coupled with ACM. This scheme has been successful in reducing the convergence steps and giving accurate results. Other CB method is the Roe scheme which was originally developed for simulation of compressible Euler equations. Razavi and Atashbar extended the Roe scheme to estimate convective fluxes of incompressible viscous flows past a hydrofoil. This method is based on the propagations of virtual acoustic waves which enables damping the numerical oscillations [9]. Nithiarasu and Liu [10] implemented artificial compressibility to solve incompressible turbulent flows at medium Reynolds numbers with characteristic-based scheme. It has been used for solutions of steady and unsteady flows. They used three different turbulence models namely Wolfenstein l, Spalart–Allmaras and standard k −ε . They showed that the artificial compressibility based on CB is suitable for solutions of both steady and unsteady incompressible turbulent flows at medium Reynolds numbers. The finite-volume method along with hybrid grid and k − ε turbulence model was proposed by Feng and Shanhong [11]. Hybrid grid can be used for accelerating convergence in the unsteady time step. Yang and Chang [12] used 2-D numerical simulation of turbulent flow fields for the NREL S809 airfoil. They assumed the flow as steady, turbulent and incompressible. The turbulent flow by the revised k −ε model was adopted to calculate effects of low Reynolds numbers near the wall. Marx [13] implemented the implicit MUSCL scheme to calculate 3-D compressible Navier–Stokes equations. They also described turbulence models for separated flows with emphasis on the non-equilibrium numerical model. Ng et al. [14] developed a finite volume method based on vortex flux flow for compressible flows. In their method, the second-order central finite difference and Runge–Kutta schemes are used for flux and time discretizations, respectively. Pentaris et al. [15] presented predictions for 2-D, steady incompressible flows in laminar and turbulent regimes by the standard k − ε model. Mass conservation is coupled to other equations by artificial compressibility and pressure correction methods. They also used the second- and the fourth-order artificial dissipation to obtain good convergence and efficiency of k − ε equations. Erturk et al. [16] used stream function–vorticity equations for steady incompressible Navier–Stokes equations. They solved steady flow inside a driven cavity up to Reynolds numbers 21 000 using a 601 × 601 fine grid. Their formulation proved to be stable and effective at very high Reynolds numbers and non-orthogonal grids. In addition LES1 modeling was done
1 Large Eddy Simulation
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
3
by Catalano et al. [17] to simulate flows past a circular cylinder in the supercritical regime by the finite difference method. Elmiligui et al. [18] implemented two multiscale-type turbulence models in the PAB3D solver. The models are based on modifying RANS. Recently, Atashbar and Razavi [19] developed the Roe-like scheme for incompressible turbulent flows which is advantageous in terms of accuracy and convergence. In this study, a finite volume model is presented for incompressible turbulent flows past a circular cylinder. Considering Spalart–Allmaras as a turbulence model the CB method is used in formulations. The behavior and performance of the Spalart–Allmaras turbulence model are investigated. The numerical results are validated by other available ones. 2. Turbulence model equation Utilizing the Spalart–Allmaras model, turbulent flows are successfully calculated over a wide range of regimes using RANS equations. The Spalart–Allmaras model is consisted of transport equations for an eddy-viscosity variable ν˜ following [20]:
∂ ν˜ ∂ ν˜ ∂ ∂ ∂ ν˜ ∂ ν˜ 1 + + Cb2 ν˜ vj = Cb1 (1 − ft2 )S˜ ν˜ + (νL + ν˜ ) ∂t ∂ xj σ ∂ xj ∂ xj ∂ xj ∂ xj 2 Cb1 ν˜ − Cw1 fw − 2 ft2 + ft1 ∥∆ν˜ ∥22 (1) κ d where νL denotes the kinematic viscosity and d is the distance to the closest wall. In the above transport equation, the first term on the right hand side is the eddy-viscosity production term. The production term S˜ is evaluated as follows: χ3 χ −3 ˜S = fv3 S + ν˜ fv2 , fv 1 = , , fv 2 = 1 + κ 2 d2 Cv 2 χ 3 + Cv31 (2) ν˜ (1 + χ fv1 ) (1 − fv2 ) fv 3 = , χ= . Max (χ , 0.001) νL In Eq. (2) S is the magnitude of the mean rotation rate: S=
2ωij ωij
ωij =
1
2
∂vi ∂vj − ∂ xj ∂ xi
.
In Eq. (1) the second and third terms are conservative diffusion and non-conservative diffusion respectively. Diffusion terms 1
σ
∂ ∂ ν˜ ∂ ν˜ ∂ ν˜ + Cb2 (νL + ν˜ ) ∂ xj ∂ xj ∂ xj ∂ xj
can be expressed as follows: 1 + Cb2 ∂
σ
∂ xj
∂ ν˜ Cb2 − (νL + ν˜ ) (νL + ν˜ ) ∇ 2 ν˜ . ∂ xj σ
Then difficulties with discretization of nonlinear term
(3)
∂ ν˜ ∂ ν˜ ∂ xj ∂ xj
are circumvented [20]. The fourth term in Eq. (1) denotes
near-wall turbulence destruction which reads
1 + Cw6 3
16
g6
6
+ Cw 3
,
g = r + Cw 2 r 6 − r , r =
ν˜
. (4) S˜ κ 2 d2 Finally the fifth and sixth terms are transition damping the production and transition source of turbulence. Those transition functions are given by fw = g
ωt2 2 2 2 ft1 = Ct1 gt exp −Ct2 d + gt dt , ∆U 2 gt = Min (0.1, ∥∆v ⃗ ∥2 / (ωt ∆xt ))
ft2 = Ct3 exp −Ct4 χ 2
(5)
where dt is the distance from the nearest trip point, ωt is the vorticity at the wall at the trip point, ∆U is the 2-norm of the difference between the velocity at the trip (zero if the wall stationary) and ∆xt stands for the spacing along the wall at the trip point. Since the present study assumes fully turbulent flow, this term is not used. Therefore, the values of ft1 and ft2 which are associated with these terms are taken as zero. Also, the turbulent kinematic eddy viscosity is obtained via
νT = fv1 ν˜ .
(6)
The various constants of this model have the following values: Cb1 = 0.1355,
σ = 2/3, Cw2 = 0.3,
Cb2 = 0.622,
κ = 0.41, Cw3 = 2,
Cv 1 = 7.1,
Cv 2 = 5,
Cw1 = Cb1 /κ + (1 + Cb2 ) /σ , 2
Ct1 = 1,
Ct2 = 2,
Ct3 = 1.3,
(7) Ct4 = 0.5.
The above model is valid through the assumption of homogeneous turbulent flow.
4
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 1. Schematic representation of the primary cell to evaluate normal convective fluxes.
3. Governing equations The non-dimensional form of Navier–Stokes equations plus the eddy-viscosity transport equation for incompressible turbulent flow in the finite-volume form with the ACM read:
∂ ∂t
WdA +
FN dL =
∂Ω
Ω
∂Ω
(Rdy − Sdx) + Z
(8) 0
p u
W = , v
ν˜
0 0 0
Z=
Cb1 S˜ ν˜ +
Cb2
σ
∂ ν˜ ∂ ν˜ ∂ xj ∂ xj
∂u (1 + νT ) ∂x 1 ∂v R= , ReD (1 + νT ) ∂x 1 ∂ ν˜ (1 + ν˜ ) σ ∂x
β uN uu + p cos ζ FN = N , v uN + p sin ζ ν˜ uN
/ReD − Cw1 fw
2 ν˜ d
/ReD
0
∂ u (1 + νT ) ∂y 1 S= (1 + ν ) ∂v , T ReD ∂y ∂ ν˜ 1 (1 + ν˜ ) σ ∂y
(9)
where, x=
X D
,
y=
Y D
,
u=
U U∞
,
v=
V U∞
,
∆y ∆x , cos ζ = , sin ζ = − , ν ∆L ∆L − → − → − → ∆L = ∆x2 + ∆y2 , N = cos ζ i + sin ζ j , up = v cos ζ − u sin ζ , a = u2N + β ReD =
t =
τ U∞ D
,
p=
P − P∞ 2 ρ∞ U∞
,
U∞ D
(10) uN = u cos ζ + v sin ζ ,
− →
where, β, uN , up , D, ∆L, a, N are the artificial compressibility parameter, the normal and parallel components of the velocity, the diameter of the cylinder, the length, the pseudo-sound speed and the normal unit vector or the direction of virtual wave propagation respectively as shown in Fig. 1. 3.1. Convective fluxes Two schemes were used to evaluate convective fluxes at cell interfaces. Considering Fig. 1, one has
∂Ω
FN dL =
4
(K )
FN ∆L
(11)
K =1
where, K is attributed to each cell side. Two schemes were applied to convective fluxes (FN ) namely Averaging (AVR) and characteristic-based (CB) schemes to compare and show the superiority of the introduced CB scheme. Considering Fig. 1 at primary cell interfaces like AB in Fig. 1, these schemes are expressed as follows.
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
5
Fig. 2. τ − ξ co-ordinate.
3.1.1. Averaging scheme (AVR) In this scheme with pressure-based modification we have: W = ηWL + (1 − η)WR ;
η=
pL
(12)
pL + pR
where, R and L indicate right and left cells respectively by p as the driving force [9]. The first and higher accuracy can be achieved in Eq. (12). To overcome numerical oscillations, artificial dissipation was implemented. 3.1.2. CB scheme The original CB method, defines primitive variables as functions of corresponding values along the virtual characteristics. Considering Fig. 2, where τ is the pseudo time and ξ is the direction of the outward normal vector of the surface of the cell, the flow variable vector at pseudo time n + 1 can be calculated along a characteristic, k, using Taylor series expansion and the initial value at pseudo time level n: W = Wk + Wξ ξτ 1τ + Wτ1τ and Wτ =
W − Wk
1τ
− Wξ ξτ .
(13)
Substituting Wτ into the inviscid form of Eq. (8), after some algebraic manipulations, leads to the following compatibility equations along characteristics:
λ(1) = uN ,
λ(2) = uN + a,
λ(3) = uN − a
(u − u1 ) sin ζ − (v − v1 ) cos ζ = 0 for λ = λ1
(14)
p − p2 = −λ2 [(u − u2 ) cos ζ − (v − v2 ) sin ζ ] for λ = λ2
(15)
p − p3 = −λ3 [(u − u3 ) cos ζ − (v − v3 ) sin ζ ] for λ = λ3 .
(16)
Then u, v and p at cell interfaces are determined using preceding equations as [5]: u = Z cos ζ + sin ζ (u1 sin ζ − v1 cos ζ )
(17)
v = Z sin ζ + cos ζ (v1 cos ζ − u1 sin ζ )
(18)
p=
λ2 [p3 + λ3 (u3 cos ζ + v3 sin ζ )] − λ3 [p2 + λ2 (u2 cos ζ + v2 sin ζ )] 2a
(19)
where Z =
1 2a
[(p2 − p3 ) + cos ζ (λ2 u2 − λ3 u3 ) + sin ζ (λ2 v2 − λ3 v3 )].
The variables at n pseudo time level are approximately evaluated based on the signs of the characteristics: Wj =
1 2
[(1 + sign(λj ))WL + (1 − sign(λj ))WR ]
WLi+ 1 = Wi , 2
WRi+ 1 = Wi+1 . 2
(20) (21)
6
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 3. Schematic representation of primary and secondary cells to evaluate viscous terms.
For the eddy-viscosity of the Spalart–Allmaras model considering Eq. (12), we have:
ν˜ = ην˜ L + (1 − η)˜νR .
(22)
3.2. Viscous fluxes Calculation of viscous fluxes on the cell boundary requires the second-order derivatives of flow parameters within cell sides. The calculation method is illustrated in Fig. 3. Calculation of viscous flux vectors is similar to estimation of convective vectors as follows:
∂Ω
(Rdy − Sdx) =
∂Ω
∂φ ∂x
dy −
∂φ ∂y
dx =
4 ∂φ K =1
∂x
∆yK −
∂φ ∂y
∆xK
(23)
where, ϕ stands for u, v, ν˜ and K is attributed to each cell side, namely AB, BC, CD and DA. For example for face AB one has
∂φ ∂x
= AB
1
Ω′
∂Ω′
φ dy =
4 1
Ω′
φK ∆yK =
K =1
1
Ω′
[φMB ∆MB + φBN ∆BN + φNA ∆NA + φAM ∆AM ].
(24)
For example one can take the following averaged values:
ϕM =
1 4
(ϕA + ϕB + ϕC + ϕD ),
ϕB =
1 4
(ϕN + ϕNE + ϕE + ϕM ) and ϕMB =
1 2
(ϕM + ϕB ).
(25)
4. Time discretization and boundary conditions 4.1. Time discretization The fifth-order Runge–Kutta (RK5) schemes are considered because of wider stability range. (CFL ≤ 3)
∂U + Q = 0, ∂t
U(n+1) = U(n) − εm ∆tQ(n) (U),
εm = 1/4, 1/6, 3/8, 1/2, 1 · (m = 1, . . . .5)
(26)
where Q includes the discretized form of integral convective and viscous terms. The number of iterations to achieve the steady state is reduced by using a variable time step that depends on the local flow changes and the grid spacing as [9]
∆x2i.j + ∆y2i.j CFL ∆t = Max , Max (γ ) 2 2 n n Max a + ui,j + vi,j γn = Max λ(1) , λ(2) , λ(3) , λ(4) ∆x2i.j + ∆y2i.j , CFL × Min
(27) n = 1, 2, 3, 4.
To accelerate convergence the local time stepping and implicit residual smoothing were utilized. 4.2. Boundary conditions The explicit boundary conditions are used throughout the computational domain. The consistent boundary treatment ensures the disturbance dissipation in the discretized domain without reflection. At the inlet boundary, all of the flow
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
7
Fig. 4. A part of grid around a circular cylinder.
Fig. 5. Mean pressure distribution on the cylinder: (—) CB, characteristic-based at ReD = 6.7E5; (•) experiment by Falchsbart [17] (in Zdravkovich, 1997) at ReD = 6.7E5; (H) experiment by Warschauer and Leene [21] (1971) at ReD = 1.2E6; (- - -} LES by Catalano et al. [17] at ReD = 1E6.
parameters except pressure are set to free stream values, whereas the pressure is extrapolated from the interior domain. The free stream value of the turbulent eddy-viscosity is usually taken as ν˜ = 0.1 [20]. At the outlet boundary, only the pressure equals its reference value in free stream. Other flow parameters are extrapolated from the interior domain. At the solid boundary considering the no-slip condition zero mass flux through the surface is prescribed. The pressure on the solid surface is found by extrapolation. It was set ν˜ = 0 at the solid boundary and hence νT = 0 [20]. 5. Numerical results and discussion Numerical tests were performed with fully turbulent flows for cross flow around a horizontal circular cylinder by two schemes: averaging (AVR) and characteristic-based (CB). The grid consisted of 100 × 100 cells along both radial and circumferential directions and the outer boundary is located 15 diameters away from the cylinder center as shown in Fig. 4. In Fig. 5 the wall pressure distribution at ReD = 6.7E5 by the CB scheme is compared to experimental data of Falchsbart [17] at ReD = 6.7E5 and Warschauer and Leene [21] at ReD = 1.2E6 and also to the LES model by Catalano et al. at ReD = 1E6. The maximum pressure is located at the stagnation point where θ = 0°. It can be observed that near 140° where separation occurs, the wall pressure falls down.
8
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 6. Comparison of convergence histories: (a) ReD = 6.7E5; (b) ReD = 1.2E6. CB, characteristic-based; AVR, averaging.
Fig. 7. Mean streamwise velocity distribution at ReD = 1E6 predicted: (a) AVR, averaging; (b) CB, characteristic-based; (c) LES and URANS by Catalano et al. [17].
In Fig. 6, convergence histories, for ReD = 6.7E5 and ReD = 1.2E6 have been shown. Because of numerical oscillations in AVR, artificial dissipation is unavoidable. It is shown that the CB scheme resulted in a faster convergence which demonstrates its appropriate capability to handle incompressible flows. In Fig. 7, u-velocity contours for ReD = 1E6 are presented and compared with that of Catalano et al. by LES and URANS.2 As shown, the CB method has smooth iso-velocity lines but in the AVR method, the results were polluted by wiggles. In Fig. 8 considering other RANS solutions, the wall pressure distribution at ReD = 1.4E5 is compared to that of Elmiligui et al. [18] in which the k − ε turbulence model has been used and validated with experimental data of Roshko’s [22]. Skin friction distribution for ReD = 6.7E5 utilizing the CB scheme is compared to that of experimental data by Achenbach [22] and LES by Catalano et al. in Fig. 9. At the separation area the skin friction decreases because of the low velocity gradient versus the vertical distance from the wall. Substantial agreement between the results of CB and experimental data can be noticed. In Fig. 10, the drag coefficient versus ReD is plotted. It can be seen that at higher Reynolds numbers, the regime is changed to turbulence. Hence, the drag falls down because of separation delay and less pressure difference between upstream and downstream. However, in this situation both the drag and skin friction grow as a result of high velocity gradient.
2 Unsteady Reynolds Averaged Navier–Stokes
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
9
Fig. 8. Mean pressure distribution on the cylinder at ReD = 1.4E5: (—) CB, characteristic-based; (N) k − ε by Elmiligui et al. [18]; (◦) experiment by Roshko [22].
Fig. 9. Skin friction distribution on the cylinder: (—) CB, characteristic-based at ReD = 6.7E5; (◦) experiments by Achenbach [23] (1968) at ReD = 3.62E6; (- - -) LES by Catalano et al. [17] at ReD = 1E6.
In Fig. 11 streamlines for ReD = 6.7E5 and ReD = 1.2E6 with the CB method are presented. As shown, in the wake region vortices are appeared and by increasing ReD the separation point moves to downstream. The vector plot of downstream by the CB method at ReD = 6.7E5 is demonstrated in Fig. 12. 6. Conclusion and outlook This paper presents 2-D incompressible turbulent flow by the finite volume method and artificial compressibility. Utilizing RANS, the one equation Spalart–Allmaras turbulence model was implemented. Numerical tests were conducted to flow past a circular cylinder as a case study. Time discretization was done by the fifth-order Runge–Kutta scheme. Local time stepping and residual smoothing were utilized to accelerate the convergence process. In the convective flux treatment, the propagation of information along characteristics was taken into account. From the mathematical point of view,
10
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 10. Drag coefficient as a function of the ReD : (—) Achenbach [23] (1968); (•) CB, characteristic-based; () LES, (△) URANS by Catalano et al. [17].
Fig. 11. Stream lines: (a) ReD = 6.7E5; (b) ReD = 1.2E6. CB, characteristic-based.
Fig. 12. Vector plot: ReD = 6.7E5. CB, characteristic-based.
the characteristic-based scheme was developed which shows the true wave paths. This scheme which is based on the propagation of information along characteristics, allows relatively high CFL numbers and has high convergence rate with less numerical oscillations without any requirement to damping. The substantial consensus was achieved utilizing the CB method in comparison with that of others. Application to three-dimensional and unsteady turbulent flows is recommended as future extensions of the present scheme.
A. Atashbar Orang et al. / Journal of Computational and Applied Mathematics (
)
–
11
References [1] G. Volpe, Performance of compressible flow codes at low Mach numbers, AIAA Journal 31 (1) (1993) 46–56. [2] P.M. Gresho, Some interesting issues in incompressible fluid dynamics, both in the continuum and in numerical simulation, Journal of Advances in Applied Mechanics 28 (1992) 45–140. [3] A.J. Chorin, A numerical method for solving incompressible viscous problems, Journal of Computational Physics 2 (1967) 12–26. [4] D. Kwak, C. Kiris, C.S. Kim, Computational challenges of viscous incompressible flows, Journal of Computational Fluids 34 (2004) 283–299. [5] D. Drikakis, P.A. Govatsos, D.E. Papatonis, A characteristic based method for incompressible flows, International Journal of Numerical Methods in Fluids 19 (1994) 667–685. [6] P. Neofytou, Revision of the characteristic-based scheme for incompressible flows, Journal of Computational Physics 222 (2007) 475–484. [7] X. Su, Y. Zhao, X. Huang, On the characteristics-based ACM for incompressible flows, Journal of Computational Physics 227 (2007) 1–11. [8] S.E. Razavi, K. Zamzamian, A. Farzadi, Genuinely multidimensional characteristic-based scheme for incompressible flows, International Journal of Numerical Methods in Fluids 57 (2008) 929–949. [9] S.E. Razavi, A. Atashbar Orang, A comparative investigation of hydrofoil at angles of attack, International Journal of Numerical Methods in Fluids 68 (2012) 1087–1101. [10] P. Nithiarasu, C.B. Liu, An artificial compressibility based characteristic based split (CBS) scheme for steady and unsteady turbulent incompressible flows, Computational Methods in Applied Mechanics and Engineering 195 (2005) 2961–2982. [11] L. Feng, J. Shanhong, Unsteady flow calculations with a multigrid Navier–Stokes method, AIAA Journal 34 (1996) 2047–2053. [12] S.L. Yang, Y.L. Chang, O. Arici, Navier–Stokes computations of the NREL airfoil using k − ω turbulent model at high angles of attack, ASME 117 (1995) 304–310. [13] Y.P. Marx, A practical implementation of turbulence models for the computation of three-dimensional separated flows, International Journal of Numerical Methods in Fluids 13 (1991) 775–796. [14] K.C. Ng, M.Z. Yusoff, T.F. Yusaf, Simulations of two-dimensional high speed turbulent compressible flow in a diffuser and a nozzle blade cascade, American Journal of Applied Sciences 9 (2005) 1325–1330. [15] A. Pentaris, K. Nikolados, S. Tsangaris, Development of projection and artificial compressibility methodologies using the approximate factorization technique, International Journal of Numerical Methods in Fluids 19 (1994) 1013–1018. [16] E. Erturk, T.C. Corke, C. Gokcol, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, International Journal of Numerical Methods in Fluids 48 (2005) 747–774. [17] P. Catalano, M. Wang, G. Iaccarino, P. Moin, Numerical simulation of the flow around a circular cylinder at high Reynolds numbers, International Journal of Heat and Fluid Flow 24 (2003) 463–469. [18] A. Elmiligui, Kh.S. Abdol-Hamid, S.J. Massey, S.P. Pao, Numerical study of flow past circular cylinder using hybrid turbulence formulations, Journal of Aircraft 47 (2010) 434–440. [19] A. Atashbar Orang, S.E. Razavi, Incompressible Roe-like scheme for turbulent flows, International Journal for Numerical Methods in Fluids 72 (4) (2013) 414–426. [20] J. Blazek, Computational Fluid Dynamics: Principles and Applications, 2005. [21] K.A. Warschauer, J.A. Leene, Experiments on mean and fluctuating pressures of circular cylinders at cross flow at very high Reynolds numbers, in: Proc. International Conference on Wind Effects on Buildings and Structures, Tokyo, Japan, 1971, pp. 305–315. [22] A. Roshko, Experimental on flow past a circular cylinder at very high Reynolds number, Journal of Fluid Mechanics 10 (1961) 345–356. [23] E. Achenbach, Distribution of local pressure and skin friction around a circular cylinder in cross-flow up to Re = 5 × 106 , Journal of Fluid Mechanics 34 (1968) 625–639.