JOURNAL OF COMPUTATIONAL PHYSICS146, 464–487 (1998) ARTICLE NO. CP986062
Implicit Weighted ENO Schemes for the
Three-Dimensional Incompressible
Navier–Stokes Equations
Jaw-Yen Yang,∗ Shih-Chang Yang,† Yih-Nan Chen,† and Chiang-An Hsu‡
∗Institute of Applied Mechanics and † Department of Mechanical Engineering, National Taiwan University, Taipei, Taiwan, Republic of China; and†Sinotech Engineering Consultants, Inc.,
Taipei, Taiwan, Republic of China
Received January 22, 1998; revised July 20, 1998
A class of lower–upper approximate-factorization implicit weighted essentially nonoscillatory (ENO) schemes for solving the three-dimensional incompressible Navier–Stokes equations in a generalized coordinate system is presented. The algo-rithm is based on the artificial compressibility formulation, and symmetric Gauss– Seidel relaxation is used for computing steady-state solutions. Weighted essentially nonoscillatory spatial operators are employed for inviscid fluxes and fourth-order central differencing for viscous fluxes. Two viscous flow test problems, laminar en-try flow through a 90◦bent square duct and three-dimensional driven square cavity flow, are presented to verify the numerical schemes. The use of the weighted ENO spatial operator not only enhances the accuracy of solutions but also improves the convergence rate for steady-state computation as compared with that using the ENO counterpart. It is found that the present solutions compare well with experimental data and other numerical results. ° 1998 Academic Pressc
1. INTRODUCTION
The design and construction of the WENO (weighted ENO) schemes for hyperbolic conservation laws are based on ENO (essentially nonoscillatory) schemes which were first introduced by Harten et al. [1] in the form of cell averages. Later Shu and Osher [2,3] devised a class of flux-based efficient ENO schemes. The main concept of ENO schemes is to use the “smoothest” stencil (in the asymptotic sense) among several candidates to approximate the fluxes at cell boundaries to a high-order accuracy and at the same time to avoid oscillations near discontinuities. ENO schemes are uniformly high-order accurate right up to the shock and are very robust to use. However, they also have certain drawbacks, as Jiang and Shu [4] have pointed out. One problem is that the freely adaptive stencil could change even by a round-off perturbation near zeroes of the solution and its derivatives.
464 0021-9991/98 $25.00
Copyright c° 1998 by Academic Press All rights of reproduction in any form reserved.
IMPLICIT WEIGHTED ENO SCHEMES 465
This free adaptation of the stencil is also not necessary in regions where the solution is smooth. The convergence rate for the implicit ENO scheme is generally less efficient. Another problem is that ENO schemes are not effective on vector supercomputers because the stencil-choosing step involves heavy usage of logical statements which perform poorly on such machines. The WENO schemes introduced by Liu et al. [5] and extended by Jiang and Shu [4] can overcome these drawbacks while keeping the robustness and high-order accuracy of ENO schemes. The concept of WENO schemes is the following: instead of approximating the numerical flux using only one of the candidate stencils, one uses a convex combination of all the candidate stencils. Each of the candidate stencils is assigned a weight which determines the contribution of this stencil to the final approximation of the numerical flux. The weights are defined in such a way that in smooth regions the stencil approaches certain optimal weights to achieve a higher order of accuracy, while in regions near discontinuities, the stencils which contain the discontinuities are assigned a nearly zero weight. Thus the essentially nonoscillatory property is achieved by emulating ENO schemes around discontinuities and a higher order of accuracy is obtained by emulating upstream central schemes with the optimal weights away from the discontinuities. Both efficient ENO and weighted ENO schemes have been extensively tested and applied to the compressible Euler/Navier–Stokes equations.
The solution methodology for viscous incompressible flows is rather different from that for compressible flows, due to the fact that there exists no time derivative in the continuity equation for incompressible flows. In order to apply compressible flow solution algorithms to incompressible flow problems, the continuity equation needs to be modified to couple with the momentum equation so that the whole system of equations can be put into the same formulation and solved efficiently. To achieve this goal, artificial compressibility may be introduced by adding the time derivative of pressure to the continuity equation, as was first proposed by Chorin [6]. The modified continuity equation, together with the unsteady momentum equations, yields a hyperbolic–parabolic-type time-dependent system of equations. Thus, fast implicit schemes developed for compressible flows, such as the approximate-factorization scheme by Beam and Warming [7], can be implemented. Various applications which evolved from this artificial compressibility concept have been reported for obtaining steady-state solution [8–16]. Merkle and Athavale [17], Rogers and Kwak [18], Rogers et al. [19], and Rosenfeld et al. [20] have reported successful computations using the pseudo-time-iteration approach for the time-dependent flow problems. The preconditioning methods for solving the incompressible flow problems were reviewed by Turkel [21]. Further developments of numerical methods for incompressible viscous flows can be found in the work by Anderson et al. [22] and by Briley et al. [23].
In this work, the WENO schemes of Jiang and Shu [4] are adopted to solve incompressible flow problems. An implicit code of WENO schemes is developed for the artificial com-pressibility formulation of the three-dimensional incompressible Navier–Stokes equations. The lower–upper symmetric Gauss–Seidel (LU-SGS) implicit algorithm [16] is adopted to solve the steady-state flow problems. This algorithm is not only unconditionally stable but also completely vectorizable in any dimensions. We apply the resulting schemes to com-pute several standard laminar flow problems including the entry flow through a 90◦bent square duct and a three-dimensional driven square cavity flow. It is found that the present solutions are in good agreement with available experimental results and other numerical results. Meanwhile, the convergence rate to a steady-state solution using implicit weighted ENO schemes is found to be much superior to that using the implicit ENO counterpart.
466 YANG ET AL.
2. GOVERNING EQUATIONS
The Navier–Stokes equations in the integral conservation law form for an incompressible, three-dimensional viscous flow with artificial compressibility can be written as
∂ ∂t µ 1 V Z V Q dV ¶ + 1 V I S EF · dES = 0, (1)
where V is the volume of an arbitrary control volume, S is the area of an arbitrary control surface, the direction of d ES is outward, Q is the conservative variables, and EF= (E − Ev)Ei+
(F − Fv) Ej+ (G − Gv)Ek is the flux vector. In Cartesian coordinates system, Eq. (1) can be
expressed as ∂Q ∂t + ∂(E − Ev) ∂x + ∂(F − Fv) ∂y + ∂(G − Gv) ∂z = 0, (2) with Q= p u v w , E= βu u2+ p uv uw , F= βv vu v2+ p vw , G= βw wu wv w2+ p , Ev= Re−1 0 2ux uy+ vx uz+ wx , Fv = Re−1 0 vx+ uy 2vy vz+ wy , Gv= Re−1 0 wx+ uz wy+ vz 2wz ,
whereβ is the artificial compressibility parameter and Re = ρV∞L/µ is the Reynolds
number. The Cartesian velocity components u, v, w are scaled with the freestream velocity
V∞and the Cartesian coordinates x, y, z are normalized with the characteristic length L. The nondimensional pressure is defined as p= (P − P∞)/ρV2
∞, and the densityρ and
dynamic viscosityµ are assumed to be constant.
Conventionally, Eq. (2) is transformed into the generalized coordinates(ξ, η, ζ) as
∂ ˆQ ∂t + ∂( ˆE − ˆEv) ∂ξ + ∂( ˆF − ˆFv) ∂η + ∂( ˆG − ˆGv) ∂ζ = 0, (3) where ˆ Q= h p u v w , Eˆ = h βU uU+ ξxp vU + ξyp wU + ξzp , Fˆ = h βV uV+ ηxp vV + ηyp wV + ηzp , ˆ G= h βW uW+ ζxp vW + ζyp wW + ζzp ,
IMPLICIT WEIGHTED ENO SCHEMES 467
ˆ
Ev= h[ξxEv+ ξyFv+ ξzGv], Fˆv= h[ηxEv+ ηyFv+ ηzGv], ˆ
Gv= h[ζxEv+ ζyFv+ ζzGv],
U = ξxu+ ξyv + ξzw, V = ηxu+ ηyv + ηzw, W = ζxu+ ζyv + ζzw,
and h is the Jacobian of the coordinate transformation (the cell volume) given by
h= ¯¯ ¯¯ ¯¯ xξ xη xζ yξ yη yζ zξ zη zζ ¯¯ ¯¯ ¯¯ =xξyηzζ + xηyζzξ + xζyξzη− xζyηzξ− xηyξzζ − xξyζzη.
The Jacobians of the inviscid fluxes ˆE, ˆF, ˆG are needed for the flux-difference splitting
and for the implicit algorithm. Let the Jacobian matrices ˆA, ˆB, ˆC ( ˆA =∂ ˆE
∂ ˆQ, ˆB =∂ ˆQ∂ ˆF, ˆC =∂ ˆG∂ ˆQ) be represented by ˆ Ai= 0 kxβ kyβ kzβ kx 2 + kxu kyu kzu ky kxv 2 + kyv kzv kz kxw kyw 2 + kzw , (4)
where ˆAi= ˆA, ˆB, ˆC for i = 1, 2, 3, respectively, and 2 = kxu+ kyv + kzw
kx= (ξi)x, ky= (ξi)y, kz= (ξi)z, ξi = (ξ, η or ζ) for i = 1, 2, 3.
A similarity transform for the Jacobian matrix is introduced, ˆ Ai = Ri3iR−1i , (5) with 3i = 2 0 0 0 0 2 0 0 0 0 2 + c 0 0 0 0 2 − c , (6)
where c is the scaled artificial speed of sound given by
c=p22+ β. (7)
The matrix of the right eigenvectors is given by
Ri= 0 0 −λ4c λ3c x2 x1 u− λ4kx u+ λ3kx y2 y1 v − λ4ky v + λ3ky z2 z1 w − λ4kz w + λ3kz , (8)
468 YANG ET AL.
and its inverse is given by
R−1i = 1 2c2 2(x1a2+ y1a3+ z1a1) 2(z1d2− y1d3) 2(x1d3− z1d1) 2(y1d1− x1d2) −2(x2a2+ y2a3+ z2a1) 2(y2d3− z2d2) 2(z2d1− x2d3) 2(x2d2− y2d1) 1 λ3kx λ3ky λ3kz 1 λ4kx λ4ky λ4kz , (9) where x1= ∂x ∂ξi+1, y1= ∂y ∂ξi+1, z1= ∂z
∂ξi+1, and ξi+1= η, ζ, or ξ for i = 1, 2, and 3, respectively
x2= ∂x ∂ξi+2, y2= ∂y ∂ξi+2, z2= ∂z
∂ξi+2, and ξi+2= ζ, ξ, or η for i = 1, 2, and 3, respectively λ3= 2 + c, λ4= 2 − c
a1= kxv − kyu, a2= kyw − kzv, a3= kzu− kxw
d1= kxβ + 2u, d2= kyβ + 2v, d3= kzβ + 2w.
3. NUMERICAL METHOD 3.1. Spatial Discretization
A semidiscrete finite volume method is used to solve Eq. (3) to ensure that the final con-verged solution is independent of the integration procedure and to avoid metric singularity problems. The finite volume method is based on the local flux balance of each mesh cell. The semidiscrete form of Eq. (3) can be written as
∂ ˆQ
∂t == −
1
Vi, j,k
{[( ˜E − ˜Ev)S]i+1/2, j,k− [( ˜E − ˜Ev)S]i−1/2, j,k]}
− 1 Vi, j,k {[( ˜F − ˜Fv)S]i, j+1/2,k− [( ˜F − ˜Fv)S]i, j−1/2,k]} − 1 Vi, j,k {[( ˜G − ˜Gv)S]i, j,k+1/2− [( ˜G − ˜Gv)S]i, j,k−1/2]}, (10) where(i, j, k) is the (i, j, k)th computational cell with volume Vi, j,k, and S is the area of each control surface and the direction is outward. The spatial differencing of numerical fluxes adopts fifth-order accurate (r= 3) weighted ENO scheme (WENO3) [4] for the inviscid convective fluxes( ˜E, ˜F, ˜G) and fourth-order central differencing for the viscous fluxes( ˜Ev, ˜Fv, ˜Gv).
By adopting WENO3 schemes, we split the physical fluxes (say, ˆF) locally into positive
and negative parts as
ˆ
F( ˆQ) = ˆF+( ˆQ) + ˆF−( ˆQ), (11)
IMPLICIT WEIGHTED ENO SCHEMES 469
chosen. In this paper, we use the local Lax–Friedrichs flux splitting method, i.e.,
ˆ
F±( ˆQ) = 1
2( ˆF( ˆQ) ± |3| ˆQ), (12)
where|3| = diag(|λ1|, |λ2|, |λ3|, |λ4|) and λ1, λ2, λ3, λ4 are the local eigenvalues2, 2,
2+c, and 2−c, respectively. For easy understanding, we first consider the one-dimensional
scalar conservation laws. For example,
ut+ f (u)x= 0. (13)
Let us discretize the space into uniform intervals of size1x and denote xj= j1x. Various
quantities at xj will be identified by the subscript j . The spatial operator of the WENO3
schemes which approximates− f (u)xat xjwill take the conservative form
L
= − 1
1x( ˜fj+1/2− ˜fj−1/2), (14)
where ˜fj+1/2and ˜fj−1/2are the numerical fluxes. Designate ˜f+j+1/2and ˜f−j+1/2respectively
the numerical fluxes obtained from the positive and negative parts of f(u); then we have ˜
fj+1/2= ˜f+j+1/2+ ˜f−j+1/2. (15)
Here we first describe the approximation of the numerical flux ˜fj+1/2 in the
one-dimensional scalar conservation law. The WENO3 numerical flux for the positive part of f(u) is ˜ f+j+1/2= ω+0 µ 2 6 f + j−2− 7 6 f + j−1+ 11 6 f + j ¶ + ω+ 1 µ −1 6f + j−1+ 5 6 f + j + 2 6f + j+1 ¶ + ω+ 2 µ 2 6 f + j + 5 6 f + j+1− 1 6 f + j+2 ¶ , (16) where ω+ k = α+ k α+ 0 + α1++ α2+ , k= 0, 1, 2 α+0 = 1 10(² + I S + 0)−2, α1+= 6 10(² + I S + 1)−2, α2+= 3 10(² + I S + 2)−2, ² = 10−6 and I S0+= 13 12( f + j−2− 2 fj+−1+ fj+) 2+1 4( f + j−2− 4 fj+−1+ 3 fj+) 2 I S1+= 13 12( f + j−1− 2 fj++ f+j+1) 2+1 4( f + j−1− fj++1) 2 I S2+= 13 12( f + j − 2 fj++1+ f+j+2) 2+1 4(3 f + j − 4 f+j+1+ f+j+2) 2.
470 YANG ET AL.
Similarly, the WENO3 numerical flux for the negative part of f(u) is
˜ f−j+1/2= ω−0 µ −1 6 f − j−1+ 5 6 f − j + 2 6 f − j+1 ¶ + ω− 1 µ 2 6 f − j + 5 6 f − j+1− 1 6 f − j+2 ¶ + ω− 2 µ 11 6 f − j+1− 7 6 f − j+2+ 2 6 f − j+3 ¶ , (17) where ω− k = α− k α− 0 + α1−+ α−2 , k= 0, 1, 2 α0−= 3 10(² + I S − 0)−2, α1−= 6 10(² + I S − 1)−2, α2−= 1 10(² + I S − 2)−2, ² = 10−6 and I S0−= 13 12( f − j−1− 2 fj−+ f−j+1) 2+1 4( f − j−1− 4 fj−+ 3 f−j+1) 2 I S1−= 13 12( f − j − 2 fj−+1+ f−j+2) 2+1 4( f − j − f−j+2) 2 I S2−= 13 12( f − j+1− 2 fj−+2+ fj−+3) 2+1 4(3 f − j+1− 4 fj−+2+ fj−+3) 2.
Next we consider the system of three-dimensional incompressible Navier–Stokes equa-tions; the numerical flux at a cell surface m+ 1/2 in direction m is usually approximated in the local characteristic fields.
Now, we denote by rs(column vector) and ls(row vector) the sth right and left
eigenvec-tors of ˆAm+1/2 (the average Jacobian at ξm+1/2), respectively. Then the scalar WENO3
scheme can be applied to each of the characteristic fields, i.e.,
¯ Fm+1/2,s = 2 X k=0 ωk,sqk(ls· ˆFm+k−2, . . . , ls· ˆFm+k), (18) which gives the numerical flux in the sth characteristic field. Hereωk,s, k = 0, 1, 2, are the weights in the sth characteristic field,
ωk,s= ωk(ls· ˆFm−2, . . . , ls· ˆFm+2), (19)
which is a nonlinear function (ωkis defined previously), and qkare the stencils as in Eqs. (16)
and (17). The numerical fluxes obtained in each characteristic field can then be projected back to the physical space by
˜ Fm+1/2= 4 X s=1 ˆ Fm+1/2,srs. (20)
IMPLICIT WEIGHTED ENO SCHEMES 471
3.2. Time Discretization
The lower–upper (LU) factored implicit scheme which was developed by Jameson and Yoon [24] is unconditionally stable in any number of space dimensions. In the framework of aproximate factorization implicit scheme, the flux vectors can be linearized by setting
ˆ
En+1 = ˆEn+ ˆAn1 ˆQ + O(k1 ˆQk2)
ˆ
Fn+1 = ˆFn+ ˆBn1 ˆQ + O(k1 ˆQk2)
ˆ
Gn+1 = ˆGn+ ˆCn1 ˆQ + O(k1 ˆQk2)
ˆ
Env+1 = ˆEnv+ ˆAnv1 ˆQ + O(k1 ˆQk2)
ˆ
Fnv+1 = ˆFnv+ ˆBnv1 ˆQ + O(k1 ˆQk2)
ˆ
Gnv+1 = ˆGnv+ ˆCnv1 ˆQ + O(k1 ˆQk2),
where ˆA, ˆB, ˆC, ˆAv, ˆBv, ˆCv are the Jacobian matrices of inviscid fluxes ˆE, ˆF, ˆG and
vis-cous fluxes ˆEv, ˆFv, ˆGv, respectively, and1 ˆQ = ˆQn+1− ˆQnis the increment of conservative
variables.
The inviscid Jacobians can be split according to the sign of the eigenvalues, ˆ
Ai = ˆA+i + ˆA−i = Ri3+i R−1i + Ri3−i R−1i . (21)
Here3+i is formed by the nonnegative part of the3imatrix and3−i by the nonpositive part. An Euler implicit time discretization of Eq. (10) can be written as
Vi, j,k ¡ˆ
Qni, j,k+1 − ˆQni, j,k¢
1t = −
©
[( ˜E − ˜Ev)S]ni+1/2, j,k+1 −£( ˜E − ˜Ev)S]ni−1/2, j,k+1 ¤ª
−©[( ˜F − ˜Fv)S]ni, j+1/2,k+1 −£( ˜F − ˜Fv)S]ni, j−1/2,k+1 ¤ª
−©[( ˜G − ˜Gv)S]ni, j,k+1/2+1 −£( ˜G − ˜Gv)S]ni, j,k−1/2+1 ¤ª, (22) where n is the time level. An unfactored implicit scheme can be obtained by substituting the above relations into Eq. (22) and dropping terms of second and higher orders. This results in the governing equation in diagonally dominant form
Vi, j,k
1t I1 ˆQi, j,k+ {[( ˆA+− ˆAv)S]i+1/2, j,k1 ˆQi, j,k− [( ˆA+− ˆAv)S]i−1/2, j,k1 ˆQi−1, j,k + [( ˆA−+ ˆAv)S]i+1/2, j,k1 ˆQi
+1, j,k− [( ˆA−+ ˆAv)S]i−1/2, j,k1 ˆQi, j,k + [( ˆB+− ˆBv)S]i, j+1/2,k1 ˆQi
, j,k− [( ˆB+− ˆBv)S]i, j−1/2,k1 ˆQi, j−1,k + [( ˆB−+ ˆBv)S]i, j+1/2,k1 ˆQi
, j+1,k− [( ˆB−+ ˆBv)S]i, j−1/2,k1 ˆQi, j,k + [( ˆC+− ˆC
v)S]i, j,k+1/21 ˆQi, j,k− [( ˆC+− ˆCv)S]i, j,k−1/21 ˆQi, j,k−1 + [( ˆC−+ ˆC
v)S]i, j,k+1/21 ˆQi, j,k+1− [( ˆC−+ ˆCv)S]i, j,k−1/21 ˆQi, j,k}n = −{[( ˜E − ˜Ev)S]i+1/2, j,k− [( ˜E − ˜Ev)S]i−1/2, j,k]}n
− {[( ˜F − ˜Fv)S]i, j+1/2,k− [( ˜F − ˜Fv)S]i, j−1/2,k]}n − {[( ˜G − ˜Gv)S]i, j,k+1/2− [( ˜G − ˜Gv)S]i, j,k−1/2]}n
≡ RHS, (23)
472 YANG ET AL.
The implicit viscous Jacobians are also considered here to enhance the convergence rate, especially for high-Reynolds-number flows in which grid systems with high aspect ratio near the walls are used to resolve the boundary layer.
In order to maximize the efficiency, Jacobian matrices of the flux vectors are approx-imately constructed to give diagonal dominance. ˆA+, ˆA−, ˆB+, ˆB−, ˆC+, and ˆC− are con-structed so that the eigenvalues of “+” matrices are nonnegative and those of “−” matrices are nonpositive, i.e.,
ˆ Ai±=1 2 £ˆ Ai± ρAˆiI ¤ , (24)
with the spectral radius of Jacobians
ρAˆi = κ max[|λ( ˆAi)|], (25)
whereλ( ˆAi) represent eigenvalues of Jacobian matrix ˆAiandκ is a constant that is greater than or equal to 1 to ensure the splitting of flux Jacobians diagonally dominant.
The unfactored implicit scheme, Eq. (23), produces a large block banded matrix that is very costly to invert and requires large amounts of storage. This difficulty can be solved by adopting the LU factored implicit scheme. The lower–upper symmetric successive over-relaxation (LU-SSOR) scheme of Yoon and Jameson [25] has the advantages of LU factor-ization and SSOR relaxation. In this paper, we adopt the LU-SSOR implicit factorfactor-ization scheme to solve the flow problems.
Equation (23) can be simplified if all the Jacobians that should be evaluated at the indicated cell faces are calculated at the local cell centers, and this can be achieved if two-point, one-sided differences are used. In addition, if we assume that the adjacent cell faces on the diagonal are approximately equal, say in i direction,
Si+1/2, j,k' Si−1/2, j,k = SI = 0.5(Si+1/2, j,k+ Si−1/2, j,k) (26)
and recognize that
ˆ
A+− ˆA−= ρAˆ (27)
and replace all viscous Jacobians with their spectral radius approximation ˆ
Av' ρAˆ
v= νSI
V I, (28)
then, using the above relations, the LU-SSOR scheme can be written as
[LD−1U]n1 ˆQ = RHSn, (29) where L= Vi, j,k 1t I+ ©£¡ ρAˆ + 2ρAˆv ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆv ¢ SK ¤ i, j,k −£¡Aˆ++ ρˆ Av ¢ i−1, j,kSi−1/2, j,k+ ¡ˆ B++ ρˆB v ¢ i, j−1,kSi, j−1/2,k +¡Cˆ++ ρˆ Cv ¢ i, j,k−1Si, j,k−1/2 ¤ª
IMPLICIT WEIGHTED ENO SCHEMES 473
D= Vi, j,k 1t I+ £¡ ρAˆ + 2ρAˆv ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆv ¢ SK ¤ i, j,k U= Vi, j,k 1t I+ ©£¡ ρAˆ + 2ρAˆv ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆv ¢ SK ¤ i, j,k +£¡Aˆ−− ρˆ Av ¢ i+1, j,kSi+1/2, j,k+ ¡ˆ B−− ρˆBv¢i, j+1,kSi, j+1/2,k +¡Cˆ−− ρˆ Cv ¢ i, j,k+1Si, j,k+1/2 ¤ª . (30)
The LU-SSOR implicit scheme reduces to the LU-SGS implicit algorithm [16] in the limit1t → ∞. Thus, Eq. (30) reduces to
L=£¡ρAˆ + 2ρAˆv ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆv ¢ SK ¤ i, j,k −£¡Aˆ++ ρˆ Av ¢ i−1, j,kSi−1/2, j,k+ ¡ˆ B++ ρˆB v ¢ i, j−1,kSi, j−1/2,k +¡Cˆ++ ρˆ Cv ¢ i, j,k−1Si, j,k−1/2 ¤ D=£¡ρAˆ + 2ρAˆv ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆv ¢ SK ¤ i, j,k U=£¡ρAˆ + 2ρAˆ v ¢ SI + ¡ ρˆB+ 2ρˆBv ¢ SJ+ ¡ ρCˆ + 2ρCˆ v ¢ SK ¤ i, j,k +£¡Aˆ−− ρˆ Av ¢ i+1, j,kSi+1/2, j,k+ ¡ˆ B−− ρˆB v ¢ i, j+1,kSi, j+1/2,k +¡Cˆ−− ρˆ Cv ¢ i, j,k+1Si, j,k+1/2 ¤ . (31)
It is interesting to note that the present implicit algorithm (LU-SGS) permits scalar diagonal inversion.
Equation (31) is solved in the following three steps:
Step 1: L1 ˆQ∗ = RHSn
Step 2: U1 ˆQn = D1 ˆQ∗
Step 3: Qˆn+1 = ˆQn+ 1 ˆQn.
(32)
The LU-SGS algorithm employs a series of corner-to-corner sweeps through the flow-fields and uses the latest available data for the off-diagonal terms to solve Eq. (32). This algorithm is completely vectorizable on i+ j + k = constant oblique planes of sweep.
3.3. Boundary Conditions
The boundary conditions imposed on the solid surface are the no-slip conditions. A zero normal pressure gradient on the wall is applied. In the far field, a locally one-dimensional characteristic type of boundary condition is used. The procedures employed here are similar to those usually used for the compressible flows. The Riemann invariants for the present system of equations are now given by
R±= p + 1 2u 2 n± 1 2[unc+ β ln(un+ c)], (33)
where un is the component of the velocity normal to the boundary. In all calculations, the
474 YANG ET AL.
4. RESULTS AND DISCUSSION
Presented here are the results of two different three-dimensional laminar flow computa-tions. These are the flow through a 90◦bending square duct and lid-driven cavity flow. 4.1. Flow through a 90◦Bending Square Duct
Ducts with rectangular/square cross sections are very frequently used in many engineer-ing applications, such as aircraft intakes, turbomachinery blade passages, diffusers, and heat exchangers. A distinguished characteristic of the flow in ducts with strong curvature is the generation of streamwise vorticity caused by the centrifugal forces which gener-ate substantial secondary flow and redistribution of the streamwise velocity in the radial direction.
The experiment of Humphrey et al. [26], in which the flow through a strongly curved 90◦ square bend duct was measured, is selected as a test case in the present study. The measurements were carried out at Reynolds number Re= 790, based on the inflow bulk velocity and the hydraulic diameter, with the corresponding Dean’s number De= 368; i.e., the problem was nondimensionalized using the side of the square cross section as the unit length and the average inflow velocity as the unit velocity. In the present work, three different grid systems with mesh sizes of 25× 17 × 17, 49 × 33 × 33, and 73 × 49 × 49 with the same artificial compressibility parameterβ = 1.0 were used to solve this problem.
The geometry and the grid system of 25× 17 × 17 are shown in Fig. 1. The straight inflow section before the bend was set to a length of 5.0 and the outflow section downstream of
IMPLICIT WEIGHTED ENO SCHEMES 475
FIG. 2. The convergence history of ENO and WENO schemes for the flow through a 90◦bending square duct at the grid systems of 49× 33 × 33.
the bend was also set to a length of 5.0. The radius curvature of the inner wall(ri) in the curved section was 1.8, while that of the outer wall(ro) was 2.8. A fully developed inflow velocity profile is prescribed at the inlet boundary and the Neumann boundary conditions (zero normal derivatives for all velocity components) are imposed at the outflow boundary.
The convergence history of the ENO2 scheme(r = 2) and the WENO scheme (r = 3) for this problem at the grid system of 49× 33 × 33 is shown in Fig. 2. It can be seen that the rapid convergence rate and monotonous curve were obtained for the WENO3 scheme. Meanwhile, the corresponding convergence rate of the ENO2 scheme [27] is very poor. It is clear to see the effectiveness in applying the WENO3 scheme in this three-dimensional case.
For the grid independence consideration, three different grid systems of 25× 17 × 17, 49× 33 × 33, and 73 × 49 × 49 are investigated. Figure 3 shows the comparison of com-puted results of the WENO3 scheme at those grid systems with the experimental data of Humphrey et al. [26]. The streamwise velocity(Vθ) profile is presented in this figure for six different cross sections along the duct. The location of these cross flow planes is shown in Fig. 1. In Fig. 3, the x-axis is the normalized radial distance and the y-axis is in the azimuthal direction. Except for the results of coarse grid system, the computed results com-pare well with the experimental results, particularly at the first four streamwise stations. However, some discrepancy is found between the numerical and experiment results at the two downstream planes. This deviation also can be found for the other numerical calcula-tions of Rogers et al. [19] and Rosenfeld et al. [20]. Nevertheless, the peaks of streamwise velocity near the outside wall at those stations are very well captured.
476 YANG ET AL.
FIG. 3. Comparison of streamwise velocity(Vθ) profiles at different streamwise locations (midspan) on three
different grids with the experimental results.
Figure 4 shows the comparison of computed results of ENO2 and WENO3 schemes at the middle grid system (49× 33 × 33) with the experimental data. It can be seen that even thougth the convergences rate of the ENO2 scheme is poor, the accuracy is as good as that of the WENO3 scheme.
Figure 5 shows the cross-sectional velocity vector fields at the plane ofθ = 0◦, 30◦, 60◦, and 90◦. The figures show how a pair of secondary vortices are generated. The centers of these vortices seen to move toward the inner wall between theθ = 30◦ station and the
θ = 60◦ position, and then tend to center again further downstream (θ = 90◦), and at the
same time a secondary pair of vortices near the outer corners is established. This agrees qualitatively with the observations of the experiment of Humphrey et al. [26].
4.2. Driven Cavity Flow
The lid-driven cavity flow, a classic recirculating flow, is an idealization of many envi-ronmental, geophysical, and industrial flows. It is a typical benchmark problem for solvers
IMPLICIT WEIGHTED ENO SCHEMES 477
FIG. 4. Comparison of the computed streamwise velocity(Vθ) profiles of ENO and WENO schemes at
different streamwise locations (midspan) with the experimental results.
of the incompressible Navier–Stokes equations. Since the cavity flow problem in either its two- or three-dimensional case is an ideal configuration for studying complex flow physics in a simple geometry, this problem has been extensively studied for more than three decades and draws continuous attention.
This problem choice is prompted by numerous experimental observations of Koseff and Street [28–30] and Aidun et al. [31]. Three-dimensional calculations have been performed by Ku et al. [32], Guj and Stella [33], Jiang et al. [34], Fujima et al. [35], and Ho and Lin [36] for the spanwise aspect ratio (SAR)= 1.0, and by Freitas et al. [37], Freitas and Street [38], and Chiang et al. [39] for SAR= 3.0. For the code validation, numerical simulations using the WENO3 scheme have been conducted first for the upper-lid-driven flow in a cubic cavity (SAR= 1.0) at three different Reynolds numbers, Re = 100, 400 and 1000, and then for the case of SAR= 3.0 over a wide range of Reynolds numbers from Re = 100 to Re= 3200.
The geometry and grid systems of 33× 33 × 33 (SAR = 1.0) is shown in Fig. 6. In Fig. 7, the computed velocity profiles of u on the vertical centerline andv on the horizontal
478 YANG ET AL.
FIG. 5. The cross-sectional velocity vector fields at the plane ofθ = 0◦, 30◦, 60◦, and 90◦for 73× 49 × 49 grid.
IMPLICIT WEIGHTED ENO SCHEMES 479
FIG. 7. The computed velocity profiles of u on the vertical centerline andv on the horizontal centerline of the symmetry plane(z = 0.5) at Re = 100, 400, and 1000 for the driven square cavity flow (SAR = 1.0).
centerline of the symmetry plane(z = 0.5) at Re = 100, 400, and 1000 are compared with the other calculations by Jiang et al. [34]. It is shown that our numerical results compare very well with the results of Jiang et al.
Figures 8 and 9 show the steady velocity vectors plots on three midplanes (a) z= 0.5, (b) y= 0.5, and (c) x = 0.5 for Re = 400 and 1000, respectively. We can observe on the symmetric plane (z= 0.5), in Figs. 8a and 9a that the secondary vortices appear in the two lower corners and the primary vortex moves toward the center of the cube as the Reynolds number increases. This phenomenon is similar to that in the two-dimensional lid-driven square cavity, but there does not exist a secondary vortex near the left upper corner. Figures 8b and 8c illustrate a pair of primary contrarotating vortices near the upstream wall and near the bottom wall, respectively. Meanwhile, another pair of secondary vortices appears near the upper corners on the plane x= 0.5. Those pairs of primary and sec-ondary vortices strengthen with increasing Reynolds number and become more distinctive at Re= 1000, as shown in Figs. 9b and 9c. Those characteristics have also been observed in other numerical studies [32–36].
The other test case is the lid-driven cavity flow for SAR= 3.0. The geometry and grid systems of 33× 33 × 91 are shown in Fig. 10. This flow problem was calculated for a series of Reynolds numbers on a fixed nonuniform grid system of 33× 33 × 91. Figure 11 shows the convergence history of lid-driven cavity flow for different Reynolds numbers. We were able to obtain converged solutions at Reynolds number up to 1200. For higher Reynolds numbers, attempts to obtain the converged solutions failed. For the Reynolds number around 1200 (i.e., Re= 1000, 1200, 1250, 1300, and 1500), the downstream secondary eddy (DSE)
480 YANG ET AL.
FIG. 8. Velocity vectors for Re= 400 on the midplanes: (a) z = 0.5, (b) y = 0.5, and (c) x = 0.5.
size with iteration numbers are plotted in Fig. 12. Figure 12 shows the steady fixed DSE for Reynolds numbers 1000 and 1200. For Reynolds numbers beyond 1250, the fluctuation of DSE becomes more distinctive when the Reynolds number is increasing. From Figs. 11 and 12, we can see that flow patterns remain steady up to Re= 1200. With increasing Reynolds numbers, the flow unsteadiness becomes appreciable at Re= 1250 approximately. As Re takes on values larger than the critical Reynolds number, the Taylor–G¨ortler-like (TGL) vortices appear.
Figure 13 shows the comparison of the steady flow separation length DSE of predicted and the experimental results of Aidun et al. [31]. It shows good agreement with the experimental
IMPLICIT WEIGHTED ENO SCHEMES 481
FIG. 9. Velocity vectors for Re= 1000 on the midplanes: (a) z = 0.5, (b) y = 0.5, and (c) x = 0.5.
results. For the cases of Re larger than the critical Reynolds number, the flow patterns are unsteady. It is difficult to compare quantitatively with experimental results. In Fig. 14, we try to compare the normalized mean u andv velocity profiles at symmetry plane (z = 1.5) for Re= 3200. It shows that predicted results compare well with the experimental data of Koseff and Street [30].
Figures 15 and 16 show the velocity vector plots on three midplanes (a) z= 1.5, (b) y = 0.5 and (c) x = 0.5 for Re = 1000 and 3200, respectively. In Fig. 15, we can see that the flow characteristic of Re= 1000 is still steady and similar to the case of SAR = 1.0. The velocity vectors also have a stationary pair of primary contrarotating vortices near the upstream
FIG. 10. The geometry and grid systems of 33× 33 × 91 of lid-driven cavity flow for SAR = 3.0.
FIG. 11. The convergence history of lid-driven cavity flow for different Reynolds numbers.
FIG. 12. The downstream secondary eddy (DSE) size plotted with iteration numbers for different Reynolds numbers.
IMPLICIT WEIGHTED ENO SCHEMES 483
FIG. 13. The comparison of the steady flow separation length DSE of the predicted and experimental results.
wall for the y= 0.5 plane and near the bottom wall for the x = 0.5 plane and still have a stationary pair of secondary vortices near the upper corners on the plane x= 0.5. For Re= 3200 (Fig. 16), the TGL vortices which were first predicted by Freitas et al. [37] were observed. For the different iteration numbers, the structure of TGL vortices is different and is no longer stationary.
FIG. 14. The normalized mean u-velocity component along vertical centerline andv-velocity component along horizontal centerline for Re= 3200 at symmetry plane (z = 1.5).
484 YANG ET AL.
FIG. 15. The velocity vector plots on three midplanes: (a) z= 1.5, (b) y = 0.5, and (c) x = 0.5 for Re = 1000.
5. CONCLUSIONS
An efficient three-dimensional incompressible Navier–Stokes code based on the artificial compressibility formulation of Chorin has been developed using the implicit LU-SGS and LU-SSOR time stepping and the weighted essentially nonoscillatory spatial operator. Ap-plications to several three-dimensional steady viscous incompressible flow problems have been carried out to validate and illustrate the code. For the flow problems considered, the flow through a 90◦bending square duct and the lid-driven cavity flow, the LU-SGS implicit algorithm is employed. The use of a weighted ENO spatial operator for the inviscid fluxes not only enhances the accuracy but also improves the convergence rate for steady-state computation as compared with using the ENO counterpart. It is found that the solutions of the present algorithm compare well with experimental data and other numerical results.
IMPLICIT WEIGHTED ENO SCHEMES 485
FIG. 16. The velocity vector plots on three midplanes: (a) z= 1.5, (b) y = 0.5, and (c) x = 0.5 for Re = 3200. ACKNOWLEDGMENT
This work was sponsored by the National Science Council, ROC, through Grant NSC 84-0210-D002-019. Dr. Wun-Wen Lin of Chung-Shan Institute of Science and Technology is the technical monitor.
REFERENCES
1. A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, Uniformly high-order accurate nonoscillatory schemes,
J. Comput. Phys. 71, 231 (1987).
2. C.-W. Shu and S. Osher, Efficient implementation of nonoscillatory shock capturing schemes, J. Comput.
Phys. 77, 439 (1988).
3. C.-W. Shu and S. Osher, Efficient implementation of nonoscillatory shock capturing schemes, II, J. Comput.
Phys. 83, 32 (1989).
4. G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126, 202 (1996).
486 YANG ET AL.
5. X.-D. Liu, S. Osher, and T. Chan, Weighted essentially nonoscillatory schemes, J. Comput. Phys. 115, 200 (1994).
6. A. J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2, 12 (1967).
7. R. M. Beam and R. F. Warming, An implicit factored scheme for the compressible Navier–Stokes equations,
AIAA J. 16, 393 (1978).
8. J. L. Steger and P. Kutler, Implicit finite-difference procedures for the computation of vortex wakes, AIAA J. 15, 581 (1977).
9. D. Choi and C. L. Merkle, Application of time-iterative schemes to incompressible flow, AIAA J. 23, 1518 (1985).
10. D. Kwak, J. L. C. Chang, S. P. Shanks, and S. R. Chakravarthy, A three-dimensional incompressible Navier– Stokes flow solver using primitive variables, AIAA J. 24, 390 (1986).
11. J. L. C. Chang, D. Kwak, and S. C. Dao, A Three Dimensional Incompressible Flow Simulation Method and
Its Application to the Space Shuttle Main Engine. I. Laminar Flow, AIAA Paper 85-0175 (1985).
12. J. L. C. Chang, D. Kwak, S. C. Dao, and R. Rosen, A Three Dimensional Incompressible Flow Simulation
Method and Its Application to the Space Shuttle Main Engine. II. Turbulent Flow, AIAA Paper 85-1670
(1985).
13. J. L. C. Chang and D. Kwak, Numerical Study of Turbulent Internal Shear Layer Flow in an Axisymmetric
U-Duct, AIAA Paper 88-0596 (1988).
14. J. L. C. Chang, D. Kwak, S. E. Rogers, and R. J. Yang, Numerical simulation methods of incompressible flows and an application to the space shuttle main engine, Int. J. Numer. Methods Fluids 8, 1241 (1988). 15. P. M. Hartwich, C. H. Hsu, and C. H. Liu, Vectorizable implicit algorithms for the flux-difference split,
three-dimensional Navier–Stokes equations, ASME J. Fluids Eng. 110, 297 (1988).
16. S. Yoon, D. Kwak, and L. Chang, LU-SGS Implicit Algorithm for Three-Dimensional Incompressible Navier–
Stokes Equations with Source Term, AIAA-89-1964-cp (1989).
17. C. L. Merkle and M. Athavale, Time-Accurate Unsteady Incompressible Flow Algorithms Based on Artificial
Compressibility, AIAA Paper 87-1137 (1987).
18. S. E. Rogers and D. Kwak, Upwind differencing schemes for the time-accurate incompressible Navier–Stokes equations, AIAA J. 28, 253 (1990).
19. S. E. Rogers, D. Kwak, and C. Kiris, Steady and unsteady solutions of the incompressible Navier–Stokes equations, AIAA J. 29, 603 (1991).
20. M. Rosenfeld, D. Kwak, and M. Vinokur, A fractional step solution method for the unsteady incompressible Navier–Stokes equations in generalized coordinate systems, J. Comput. Phys. 94, 102 (1991).
21. E. Turkel, Review of preconditioning methods for fluid dynamics, Appl. Numer. Math. 12, 257 (1993). 22. W. K. Anderson, R. D. Rausch, and D. L. Bonhaus, Implicit/Multigrid Algorithms for Incompressible Turbulent
Flows on Unstructured Grids, AIAA-95-1740-CP (1995).
23. W. R. Briley, S. S. Neerarambam, and D. L. Whitfield, Implicit Lower-Upper/Approximate-Factorization
Algorithms for Viscous Incompressible Flows, AIAA-95-1742-CP (1995).
24. A. Jameson and S. Yoon, Lower-upper implicit schemes with multiple grids for the Euler equations, AIAA
J. 25, 929 (1987).
25. S. Yoon and A. Jameson, An LU-SSOR Scheme for the Euler and Navier–Stokes Equations, AIAA Paper 87-0600 (1987).
26. J. A. C. Humphrey, A. M. K. Taylor, and J. H. Whitelaw, Laminar flow in a square duct of strong curvature,
J. Fluid Mech. 83, 509 (1977).
27. J. Y. Yang and C. A. Hsu, High-resolution, nonoscillatory schemes for unsteady compressible flows, AIAA J. 30, 1570 (1992).
28. J. R. Koseff and R. L. Street, Visualization studies of a shear driven three-dimensional recirculating flow,
ASME J. Fluids Eng. 106, 21 (1984).
29. J. R. Koseff and R. L. Street, On end wall effects in a lid-driven cavity flow, ASME J. Fluids Eng. 106, 385 (1984).
IMPLICIT WEIGHTED ENO SCHEMES 487
30. J. R. Koseff and R. L. Street, The lid-driven cavity flow: A synthesis of qualitative and quantitative observations,
ASME J. Fluids Eng. 106, 390 (1984).
31. C. K. Aidun, N. G. Triantafillopoulos, and J. D. Benson, Global stability of a lid-driven cavity with throughflow: Flow visualization studies, Phys. Fluids A 3, 2081 (1991).
32. H. C. Ku, R. S. Hirsh, and T. D. Taylor, A pseudospectral method for solution of the three-dimensional incompressible Navier–Stokes equations, J. Comput. Phys. 70, 439 (1987).
33. G. Guj and F. Stella, A vorticity–velocity method for the numerical solution of 3D incompressible flows,
J. Comput. Phys. 106, 286 (1993).
34. B. N. Jiang, T. L. Lin, and L. A. Povinelli, Large-scale computation of incompressible viscous flow by least-squares finite element method, Comput. Methods Appl. Mech. Eng. 114, 213 (1994).
35. S. Fujima, M. Tabata, and Y. Fukasawa, Extension to three-dimensional problems of the upwind finite ele-ment scheme based on the choice of up- and downwind points, Comput. Methods Appl. Mech. Eng. 112, 109 (1994).
36. C. J. Ho and F. H. Lin, Numerical simulation of three-dimensional incompressible flow by a new formulation,
Int. J. Numer. Methods Fluids 23, 1073 (1996).
37. C. J. Freitas, R. L. Street, A. N. Findikakis, and J. R. Koseff, Numerical simulation of three-dimensional flow in a cavity, Int. J. Numer. Methods Fluids 5, 561 (1985).
38. C. J. Freitas and R. L. Street, Non-linear transient phenomena in a complex recirculating flow: A numerical investigation, Int. J. Numer. Methods Fluids 8, 769 (1988).
39. T. P. Chiang, R. R. Hwang, and W. H. Sheu, Finite volume analysis of spiral motion in a rectangular lid-driven cavity, Int. J. Numer. Methods Fluids 23, 325 (1996).
