Chinese Journal of Aeronautics 25 (2012) 33-41

Contents lists available at ScienceDirect

Chinese Journal of Aeronautics

journal homepage: www.elsevier.com/locate/cja

An Improved LU-SGS Implicit Scheme for High Reynolds Number

Flow Computations on Hybrid Unstructured Mesh
WANG Gang*, JIANG Yuewen, YE Zhengyin
National Key Laboratory of Science and Technology on Aerodynamic Design and Research, Northwestern Polytechnical
University, Xi’an 710072, China

Received 28 February 2011; revised 7 April 2011; accepted 8 June 2011

Abstract
The lower-upper symmetric Gauss-Seidel (LU-SGS) implicit relaxation has been widely used because it has the merits of less
dependency on grid topology, low numerical complexity and modest memory requirements. In original LU-SGS scheme, the
implicit system matrix is constructed based on the splitting of convective flux Jacobian according to its spectral radius. Although
this treatment has the merit of reducing computational complexity and helps to ensure the diagonally dominant property of the
implicit system matrix, it can also cause serious distortions on the implicit system matrix because too many approximations are
introduced by this splitting method if the contravariant velocity is small or close to sonic speed. To overcome this shortcoming,
an improved LU-SGS scheme with a hybrid construction method for the implicit system matrix is developed in this paper. The
hybrid way is that: on the cell faces having small contravariant velocity or transonic contravariant velocity, the accurate deriva-
tive of the convective flux term is used to construct more accurate implicit system matrix, while the original Jacobian splitting
method is adopted on the other cell faces to reduce computational complexity and ensure the diagonally dominant property of the
implicit system matrix. To investigate the convergence performance of the improved LU-SGS scheme, 2D and 3D turbulent
flows around the NACA0012 airfoil, RAE2822 airfoil and LANN wing are simulated on hybrid unstructured meshes. The nu-
merical results show that the improved LU-SGS scheme is significantly more efficient than the original LU-SGS scheme.

Keywords: LU-SGS scheme; hybrid unstructured mesh; Navier-Stokes equations; flux Jacobian; convergence performance; tur-
bulent flow

1. Introduction1 during the design process of an aircraft [1-2]. In order to

efficiently simulate high Reynolds number viscous
Accompanied with the development of computa- flow with current limited computer hardware, highly
tional fluid dynamics (CFD), numerically solving the stretched grids with large aspect ratios should be used
so-called Reynolds-averaged Navier-Stokes (RANS) to calculate the viscous effect in the boundary layer
equations has been widely used for viscous flow simu- region. Typically, the magnitudes of the aspect ratios
lation in many industry fields. Especially in the aero- for those “viscous” grids are more than 103. On such
nautical engineering, a large number of high Reynolds type of grid, explicit time integration scheme has very
number viscous flows always need to be simulated low efficiency because very small time stepping size
should be used with the strict limitation of the stability
condition. And this means that explicit time integration
*Corresponding author. Tel.: +86-29-88491342. scheme can hardly fulfill the engineering demands on
E-mail address: [email protected] high Reynolds number viscous flow simulation. Hence,
Foundation item: National Natural Science Foundation of China
(10802067) an implicit time integration scheme is highly desired
for the enhancement of efficiency since it has the adva-
1000-9361 © 2012 Elsevier Ltd Open access under CC BY-NC-ND license. ntage of no stability limitation in theory. Although the
doi: 10.1016/S1000-9361(11)60359-2
· 34 · WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 No.1

theoretically unconditional stability cannot usually be 2. Numerical Schemes

achieved because some approximations must be intro-
duced in practical implementation, the boundary of sta- 2.1. Governing equation
bility region of the implicit scheme can still be signifi-
cantly extended. So larger time steps are allowed to be The integral form of RANS equations enclosed with
used in an implicit scheme and faster convergence speed Spalart-Allmaras (S-A) turbulence model [11] for a con-
can be obtained compared with the explicit scheme. trol volume Ω with a surface element dS can be written
Many implicit schemes [3-10] have been developed as
and successfully applied to accelerating the conver-
gence rate. Two sorts of them are widely used on the (1)
unstructured grids, which are the Gauss-Seidel type
method represented by the well-known lower-upper where W = [ ρ ρ u ρ v ρ w ρ E ρυ% ]T is the vector of
symmetric Gauss-Seidel (LU-SGS) scheme [3-6] and the
conserved quantities with ρ, u, v, w, E and υ~ denoting
Krylov-subspace method typically represented by the
generalized minimal residual (GMRES) scheme [7-8]. the density, the Cartesian velocity components, the
The strong point of GMRES scheme is that quadratic specific total energy and the working variable of S-A
convergence speed can be achieved as the time step turbulence model, respectively. n represents the sur-
size is close to infinite. While the precondition for ob- face normal. The source term Qsource = [0 0 0 0 0
taining quadratic convergence is that the implicit sys- QT]T only lies in turbulence model equation. In which
tem matrix should be constructed and inversed exactly. QT represents the source term of S-A turbulence model.
To fulfill this requirement, an accurate linearization of Fc and Fv comprise the convective and viscous flux
residuals and the inner iterations with preconditioning vectors respectively; their detailed form can be found
methodologies should be used. Consequently, the in Ref. [12]. Since the spatial discretization and time
GMRES scheme is characterized with highly computa- integration of turbulence model equation and mean
tional complexity and large memory consumption. An- flow equations are carried out in a loosely coupled way,
other popular implicit scheme is LU-SGS, which does this makes the turbulence model equation can be easily
not add so much extra storage requirement compared solved because only one unknown variable is involved.
with explicit scheme and is free from direct inversion Hence the numerical treatments of RANS equations are
of large implicit system matrix. Compared with the mainly considered in this paper.
GMRES scheme, LU-SGS is easier to be implemented
and needs less computational cost for marching one 2.2. Spatial discretization
time step. With respect to the computational efficiency,
the LU-SGS scheme is not less competitive than the By using the finite volume method for each grid cell
GMRES method in terms of CPU time. However, it i, the spatial discretization of Eq. (1) can be expressed
should be pointed out that LU-SGS usually needs more as
iteration steps to reach the same convergence level d(Ω iWi )
compared with GMRES method. Such a drawback of
dt
= Ω i Qsource − ∑ ( (Fc )im − (Fv )im ) Sim (2)
m∈ N ( i )
LU-SGS scheme is mainly caused by several approxi-
mations introduced in the process of constructing the where Ω i represents the volume of current grid cell i,
implicit system matrix and LDU decomposition [4-5]. If
these approximations could be reduced by some kinds N(i) the set of face neighbor cells of cell i, Sim the nor-
of modification, then the efficiency of LU-SGS scheme mal vector area of the face, the subscript im denotes the
will be further enhanced. This untapped capacity of current face which is shared by cell i and cell m. In
LU-SGS scheme gives rise to the basic motivation of order to capture the discontinuity accurately and sup-
the present modification work. press numerical oscillations, a series of schemes has
The principal ingredients of this paper are as follows. been developed for the evaluation of convective flux
The basic algorithm of LU-SGS scheme is briefly de- (Fc)im , for example the central scheme of Jamerson [13],
scribed firstly. According to the analysis of the main the Roe scheme [14] and the AUSM type schemes [15].
approximations introduced by the original LU-SGS However, all these schemes can be rewritten to the
scheme [3], an improved LU-SGS scheme based on a following unified form:
new implicit system matrix construction method is 1
presented in this paper. In order to investigate the con- (Fc )im =
2
( Fc (Wim,L ) + Fc (Wim,R )) + D(Wim,L ,Wim,R ) (3)
vergence performance of the improved LU-SGS
scheme, 2D and 3D viscous flows around the NACA- where Wim,L and Wim,R denote the flow states on the left
0012 airfoil, the RAE2822 airfoil and the LANN wing and the right side of the face shared by cell i and cell m.
are simulated using hybrid unstructured grids. A com- In first order scheme, Wim,L and Wim,R are directly
parison between the improved LU-SGS scheme and the specified as the solution quantities located at the center
original LU-SGS scheme is also presented. of cell i and cell m, namely Wi and Wm. In second or-
der scheme, some reconstruction techniques should be
No.1 WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 · 35 ·

used to determine Wim,L and Wim,R. The D(Wim,L, Wim,R) where, in unstructured grid case, M is a large, sparse
in Eq. (3) represents the dissipation term added to pre- and generally non-symmetric matrix, which consists of
vent non-physical oscillations, and the iconic features DOF×DOF blocks. Here DOF represents the degrees
of different convective flux evaluation schemes mainly of freedom, which equals the total number of control
lie in the specific forms of this term. The evaluation of volumes in finite volume method.
the viscous flux term Fv is relatively simple. The vis- The LU-SGS scheme is based on the LDU factoriza-
cous flux is usually discretized by the standard central tion of the implicit system matrix, which is given by
scheme because of the elliptic nature of this term. This Ω
means that the flow quantities and their first deriva- ( I + M ) = L+ D +U = ( L+ D ) D −1 (U + D) − LD −1U
Δt
tives on the faces of the control volume, which exist in
the formula of viscous flux term, can be obtained by (9)
the simple average, and then be used to calculate Fv where L only consists of block terms in the strictly
straightforwardly. lower triangular matrix, U only consists of block terms
To simplify the expression, the residual of Eq. (2) is in the strictly upper triangular matrix and D is a block
denoted here by diagonal matrix.
Suppose that the LD−1U is a small term and its con-
Ri (W ) = ∑ ( (Fc )im − (Fv )im ) Sim + Ωi Qsource (4)
tribution to the implicit system matrix can be neglected,
m∈N ( i )
then Eq. (8) can be replaced by the following equation:
When fixed-volume computational grid is used,
Eq. (2) becomes ( L + D) D −1 (U + D)ΔW n = − Rn (10)
dWi Equation (10) can be inverted in the following for-
Ωi = − Ri (W ) (5) ward and backward sweep procedure:
dt
⎪⎧ DΔW = − R − LΔW
* n *

2.3. Analysis of LU-SGS scheme ⎨ n * n

(11)
⎪⎩ DΔW = DΔW − U ΔW
Equation (5) is a system of coupled ordinary differ- where ΔW * is solution vector updated in the forward
ential equations in time. By using the backward Euler sweep.
scheme for the implicit time integration, we obtain According to the above analysis, we can draw a con-
Wi n +1 − Wi n clusion that the property of implicit system matrix M
Ωi = − Ri (W n +1 ) (6) has strong influence on the performance of LU-SGS
Δt
implicit scheme. Hence the construction method and
where the superscript n represents the number of time properties of M need to be considered in detail. By
level. Since W n+1 is unknown on current time level, the
using Eq. (4), ∂Ri (W n ) ∂Wj n is evaluated as
residual Ri(W n+1) cannot be evaluated directly. How-
ever, it can be linearized by using first order Taylor ex- ∂ R i (W n )
pansion in the following way: =
∂Wj n
∂Ri (W n )
Ri (W n +1 ) ≈ Ri (W n ) + ∑ ∂W j n
(W j n +1 − W j n ) (7) ⎧
⎪ ∑ ⎜
⎛ ∂ (Fc )im ∂ (Fv)im ⎞ (Ω i Qsource)
∂　
j∈C ( i ) − ⎟ Sim + i= j
⎪⎪ m∈ N(i) ⎝⎜ ∂Wi ∂Wi ⎠⎟
n n
∂Win
where C(i) is the set of cell i and its neighbor cells. ⎨
⎪ ⎛ ∂ (Fc )　 ⎞ (12)
With the definitions of ij ∂ (Fv )　
ij
⎪ ⎜⎜ − ⎟ S 　 　i ≠ j
∂Wjn ⎟⎠
ij
⎧ n

⎪ Mi, j =
∂R i (W n ) ⎩⎪ ⎝ ∂Wj
⎪ ∂W jn where ∂ (Ωi Qsource ) ∂Wi n only lies in the turbulence
⎪
⎪ ⎡ M1,1 M1,2 L M1,DOF ⎤ model equations and its expression is the derivative of
⎨ ⎢ M
⎪ ⎢ 2,1 M 2,2 L M 2,DOF ⎥⎥ the source term in terms of the conservative variables.
⎪ M = The viscous flux Jacobians ( Av+ )ij = ∂(Fv )ij ∂Wi n and
⎢ M M M ⎥
⎪ ⎢ ⎥
⎪⎩ ⎢⎣ M DOF,1 M DOF,2 L M DOF,DOF ⎥⎦ ( Av− )ij = ∂(Fv )ij ∂W j n in Eq. (12) have complex ex-
and pressions even if the thin shear layer approximation is
used. In order to reduce the computational cost, the
ΔW n = W n +1 − W n viscous flux Jacobian can be approximately counted by
n+1
Ri(W ) in Eq. (6) is substituted by the linearization using its spectral radius:
term in Eq. (7), then the following implicit system is ( Av±)ij ≈ ±(λ v)ij I =
obtained:
Ω Sij ⎡ ⎛ 4 γ ⎞ ⎛ μL μ ⎞ ⎤
( I + M )ΔW n = − R(W n ) (8) ± ⎢ max ⎜ , ⎟⎜ + T ⎟ ⎥I (13)
Δt rij ⎢ ⎜ 3ρij ρij ⎟ ⎝ PrL PrT ⎠ ⎥
⎣ ⎝ ⎠ ij ⎦
· 36 · WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 No.1

where γ, μ and Pr denote the ratio of specific heat coef- lute values of the above eigenvalues represent the pro-
ficient, dynamic viscosity coefficient and Prandtl num- pagating speed of corresponding waves. It is clear that
ber respectively, and rij is the distance between grid ( Ac+ )ij and ( Ac− )ij given by Eq. (14) can reflect the pro-
cell i and cell j. In original LU-SGS approach, the con-
pagating directions of perturbations correctly. But they
vective flux Jacobians ∂(Fc )ij ∂Wi n and ∂(Fc )ij ∂W j n cannot represent the propagation velocities of pertur-
are not defined with the exact derivatives based on bations accurately, which leads to considerable distor-
Eq. (3), but approximately evaluated by the splitting of tions in evaluating the variations of fluxes caused by
convective flux Jacobian with its maximal eigenvalue the fluctuations of flow states. If the values of
⎧ ∂ ( Fc )ij (|Vij・nij | + cij ) |Vij・nij | and (|Vij・nij | + cij ) (|Vij・nij | −cij )
1
⎪ n
≈ ( Ac+ )ij = ( A + λ max I ) and are very large, serious distortions on the variations
⎪ ∂ W i 2
⎨ (14) of fluxes may be introduced by using Eq. (14). And
⎪ ∂ ( Fc )ij ≈ ( A− ) = 1 ( A − λ I ) such kind of cases occurs not only in the low Mach
⎪ ∂W n c ij
2
max
flow simulations but also commonly lies in the high
⎩ j
Reynolds number flow computations. The reason is
where A = ( ∂ Fc ∂W )ij , and λmax is the maximal ei- that the absolute value of Vij・nij could be very small on
the cell faces which are parallel to the wall boundary
genvalue of A, namely λmax = Vij・nij + cij . Here, the faces or near the bottom region of boundary layer.
subscript ij denotes the face which is shared by grid Hence, too many distortions may be introduced into
cell i and cell j, c is the sonic speed, V is the velocity. implicit system matrix with the above construction
The ( Ac+ )ij in Eq. (14) only has non-negative eigenval- method based on the maximal eigenvalue splitting, and
these distortions are harmful for the convergence be-
ue, and ( Ac− )ij non-positive eigenvalue. According to the havior of non-linear time marching iteration.
theories on propagation of perturbation, ( Ac+ )ij reflects In order to construct more accurate implicit system
the variations of fluxes caused by the increment of flow matrix, ∂( Fc )ij ∂Wi n and ∂(Fc )ij ∂W j n need to be
variables on grid cell i, and ( Ac− )ij the variations of evaluated by using the derivation of Eq. (3), namely
fluxes received from the fluctuation of flow variables on ⎧ ∂ ( Fc )ij 1 ∂Fc (Wij ,L ) ∂D(Wij ,L ,Wij ,R )
⎪ = ・ + ≈
neighboring grid cell j. Hence, there has certain rational- ⎪ ∂Wi
n
2 ∂Wi n ∂Wi n
ity for using ( Ac+)ij and ( Ac− )ij instead of ∂( Fc )ij ∂Wi n ⎪ 1 ∂Fc (Wij ,L ) ∂D (Wij ,L , Wij ,R )
⎪ ・ +
and ∂(Fc )ij ∂W j n . ⎪⎪ 2 ∂Wij ,L ∂Wij ,L
In fact, there are two kinds of iterative processes ⎨ (16)
contained in the implicit time marching scheme of this ⎪ ∂ ( Fc )ij = 1・∂Fc (Wij ,R ) + ∂D(Wij ,L , Wij ,R ) ≈
⎪ ∂W n 2 ∂W j n ∂W j n
paper. One is the non-linear time marching iteration for ⎪ j

solving Eq. (6) with backward Euler method, the other ⎪ 1 ∂Fc (Wij ,R ) ∂D (Wij ,L , Wij ,R )
one is an inner iteration on each time level for solving ⎪ ・ +
⎪⎩ 2 ∂Wij ,R ∂Wij ,R
the linear equation system Eq. (8) with LU-SGS
method. In constructing an efficient implicit scheme, where ∂D(Wij, L,Wij, R) /∂Wij, L and ∂D(Wij, L,Wij, R) /∂Wij, R
convergence and stability performances of the above are defined according to the spatial discretization
two iteration processes need to be considered on bal- scheme. For example, if Roe scheme is used, and let
ance. A Roe represent Roe average Jacobian matrix, then we
The construction of implicit system matrix with Eq.
get
(14) is favorable to the inner iteration. In detail, two
distinctive benefits can be identified. Firstly, the block ⎧ ∂D(Wij ,L , Wij ,R ) 1
diagonal matrix D in Eq. (10) degenerates into a pure ⎪ = ARoe
⎪ ∂Wij ,L 2 ij
diagonal matrix, which can be evaluated and inversed ⎨ (17)
with very little CPU time. Secondly, it ensures the dia- ⎪ ∂D(Wij ,L , Wij ,R ) = − 1 A
gonal dominance of implicit system matrix for very ⎪ ∂Wij ,R 2
Roe ij
⎩
large time stepping. This benefit not only guarantees
the convergence of LU-SGS iteration, but also makes In fact, the construction of the implicit system ma-
the neglected term LD−1U really negligible compared trix with the above accurate convective flux Jacobian
to the entire matrix system. not only brings a notable increase in computational
However, we all know that convective flux Jacobian cost, but also causes the implicit system matrix M to
matrix A has five eigenvalues: lose its property of diagonal dominance as the time
λ1 = λ 2 = λ3 =Vij・nij , λ 4 =Vij・nij -c, λ 5 =Vij・nij +c (15) step size Δt takes a great value. It is well known that
the convergence of the LU-SGS inner iteration cannot
Their signs express the propagating directions of be guaranteed if the system matrix is not diagonally
waves relative to the surface normal Sim, and the abso- dominant. In order to achieve diagonal dominance, a
No.1 WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 · 37 ·

relatively small Δt should be used in the time marching flow simulation, the improved LU-SGS will be identi-
process. In such case, the evolution speed of flow field fied with the block LU-SGS method. The optimal values
is considerably decreased and more time marching of ε1 and ε2 are hard to be known beforehand because
steps are required to achieve the convergence state. they vary with different flow conditions and meshes.
Another way for overcoming the limitation of diagonal According to the practical experience, ε1 = ε2 ≈ 0.25 is
dominance is giving up the LU-SGS scheme but using recommended in high Reynolds number flow simula-
some other inner iteration algorithms which have better tions.
stability, for example the GMRES scheme. The cost of Finally, it should be pointed out that the existence of
this choice is that the computational complexity will be viscous flux Jacobian, which is computed by using
increased significantly. Eq. (13), offers favorable term for the diagonal domi-
nance of implicit system matrix. This means that the side
2.4. Improved LU-SGS scheme effects of the usage of Eq. (18) on the stability of
LU-SGS inner iteration can be eliminated partly in the
In the above analysis, we find out that the approxi-
boundary layer region of the high Reynolds number flow.
mate convective flux Jacobian based on maximal ei-
Hence, the maximum allowable time step size will not be
genvalue splitting is helpful for the inner linear itera-
decreased obviously if the above improved LU-SGS
tion while the accurate convective flux Jacobian based
scheme is used.
on exact derivation is beneficial to the outer non-linear
time stepping. If we can make the advantages of them
complement each other, then the convergence speed of 3. Test Cases and Results
LU-SGS implicit time marching scheme will be im-
proved further. For this purpose, a hybrid method for Several typical test cases are selected here to illus-
constructing the implicit system matrix is proposed and trate the practicality and computational efficiency of
described in the following way. the improved LU-SGS implicit time marching scheme.
For the cell faces on which the contravariant speed These test cases include the low Mach number and
supersonic viscous flows around NACA0012 airfoil,
has small value or close to sound speed, ∂( Fc )ij ∂Wi n transonic viscous flow around RAE2822 airfoil and 3D
and ∂(Fc )ij ∂W j n are computed by using Eq. (16) so transonic flow around LANN wing. The S-A turbu-
lence model, second order Roe scheme and the local
that the excessive approximations on implicit matrix
time stepping [12] are chosen as common settings for all
could be avoided; on the other cell faces, the classical
computational cases. The convergence performances of
maximal eigenvalue splitting method, namely Eq. (14)
original LU-SGS and improved LU-SGS are analyzed
is used to evaluate the above convective Jacobians with
and compared in each test case.
the aim of reducing computational cost and strength-
ening the diagonal dominance of implicit system ma-
trix. According to this principle, the hybrid method 3.1. NACA0012 airfoil
for determining ∂( Fc )ij ∂Wi n and ∂ ( Fc )ij ∂W j n can be As shown in Fig. 1, a hybrid unstructured mesh for
expressed as viscous flow simulation on NACA0012 airfoil is gen-
If Vij・nij cij ≤ ε1 or 1 - ε 2 ≤ Vij・nij cij ≤ 1+ε 2 , erated with the method described in Ref. [16]. It con-
sists of 200 wall boundary nodes, 7 991 computational
⎧ ∂ ( Fc )ij 1 ∂Fc (Wij ,L ) ∂D(Wij ,L , Wij ,R ) field nodes and 11 552 mesh cells. Near the boundary
⎪ ≈ ・ + layer region, 22 layers of quadrangle mesh are distrib-
⎪ ∂Wi
n
2 ∂Wij ,L ∂Wij ,L
(18) uted and the spacing of the first layer normal to the
⎨
⎪ ∂ ( Fc )ij ≈ 1・∂Fc (Wij ,R ) + ∂D(Wij ,L , Wij ,R ) wall is chosen as 3.0×10−6 times of chord length.
⎪ ∂W n 2 ∂Wij ,R ∂Wij ,R Firstly, the incompressible viscous flow is computed
⎩ j
with an angle of attack α = 5º, a Mach number Ma =
otherwise, 0.15 and a Reynolds number Re = 3×106. The same
⎧ ∂ ( Fc )ij 1 Courant-Friedrichs-Lewy (CFL) number 30 is used for
⎪ n
≈ ( A + λ max I ) the computations with the improved and the original
⎪ ∂Wi 2 LU-SGS schemes. The convergence histories of the
⎨ (19)
⎪ ∂ ( Fc )ij ≈ 1 ( A − λ I )
maximal residual in terms of the time iteration numbers
⎪ ∂W n max are shown and compared in Fig. 2. As one can see, the
⎩ j 2
convergence speed of the maximal residual computed
In Eq. (18), ε1 and ε2 represent non-negative empirical by the original LU-SGS has slowed down quite mark-
coefficients. The smaller values are given to ε1 and ε2, edly after it decreased five orders. While in the com-
the more characteristics of the original LU-SGS will be putation with the improved LU-SGS, the convergence
assigned to the improved LU-SGS method. On the other rate is constant. In practice, we are more concerned
hand, the larger values are given to ε1 and ε2, the more about the convergence speeds of aerodynamic coeffi-
characteristics of the block LU-SGS method [4, 6] will be cients. In Fig. 3 the convergence curves of life coeffi-
represented. For example, if ε1 = 1 is used in subsonic cient are plotted. The dash-dotted lines in this
· 38 · WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 No.1

structed by using accurate convective flux Jabobians

(Eq. (16)) for the improved LU-SGS scheme. This leads
to the improved LU-SGS costs 125% of CPU time for
marching one time step compared with the original
LU-SGS method. Nonetheless, the convergence speed
is nearly accelerated with a factor of 160% in terms of
the CPU time by using the improved LU-SGS scheme.
On the same mesh, the supersonic viscous flow past
the NACA0012 airfoil at α = 2º, Ma =1.5 and Re=
6×106 is also simulated with the improved and the origi-
nal LU-SGS schemes. The computed pressure contours
are shown in Fig. 4. It is easy to be observed that a de-
tached shock has already occurred in the flow field. The
Fig. 1 Hybrid unstructured mesh around NACA0012 airfoil. convergence histories of the maximal residual obtained
by using the original LU-SGS, the improved LU-SGS
with ε1 = ε2 = 0.25 and the improved LU-SGS with
ε1=100 are presented in Fig. 5. The main purpose of
choosing the setting of ε1=100 is to investigate the con-
vergence performance of the improved LU-SGS on con-
dition that the implicit system matrix is constructed
completely with accurate convective flux Jacobians. In
this test case, the maximal allowable CFL number for the
improved LU-SGS with ε1=100 is 15.5, while the corre-
sponding CFL number for the original LU-SGS and

Fig. 2 Convergence histories of the maximal residual in

terms of iterations for low speed flow past NA-
CA0012 airfoil.

Fig. 4 Pressure contours for supersonic viscous flow past

NACA0012 airfoil.

Fig. 3 Convergence histories of lift coefficient in terms of

iterations for low speed flow past NACA0012 air-
foil.

figure represent the upper and lower convergence

tolerance of lift coefficients, which are defined
as CL,convergence ±0.002 in this paper. The results shown
in Fig. 3 indicate that with the same convergence tol-
erance, the improved LU-SGS only cost half number of
iterations as much as that of original LU-SGS scheme. Fig. 5 Convergence histories of the maximal residual in
Because of the low Mach free stream condition of this terms of iterations for supersonic viscous flow past
test case, the implicit system equation is purely con- NACA0012 airfoil.
No.1 WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 · 39 ·

the improved LU-SGS with ε1 = ε2 = 0.25 can be both results agree very well with the experimental data.
given as 50. The convergence histories in Fig. 5 show
that due to the stability limitation on CFL number, the
improved LU-SGS scheme based on pure accurate con-
vective flux Jacobians even has less efficiency compared
with the same scheme based on the hybrid definition of
convective flux Jacobians, its convergence rate is 1.5
times faster than that of the original LU-SGS based on
the maximal eigenvalue flux Jacobian splitting.

3.2. RAE2822 airfoil

The hybrid unstructured mesh generated around the

RAE2822 airfoil is shown in Fig. 6. It consists of 10 041
nodes and 14 420 mesh cells, including 256 wall bound-
ary nodes and 22 layers of quadrangle mesh in the
boundary layer region. The spacing of the first layer
normal to the wall is chosen as 5.0×10−6 times of chord
length.
Fig. 7 Convergence histories of the maximal residual in
terms of iterations for transonic viscous flow past
RAE2822 airfoil.

Fig. 6 Hybrid unstructured mesh around RAE2822 airfoil.

Here, a transonic case at α = 2.79º, Ma = 0.734 and

Re = 6.5×106 is selected to evaluate the performance of
Fig. 8 Convergence histories of lift coefficient in terms of
the original LU-SGS and the improved LU-SGS with iterations for transonic viscous flow past RAE2822
different values of ε1 and ε2. The convergence behavior airfoil.
of the maximal residuals in terms of iteration numbers is
shown in Fig. 7. As one can see from Fig. 7, the con-
vergence speed can be greatly enhanced by specifying
proper value for ε1 and ε2. In current test case, the main
enhancement comes from the function of ε1 in the range
from 0 to 0.1 and ε2 only plays a secondary role. In gen-
eral, the recommended values of ε1=ε2=0.25 are close to
the optimal combination. The convergence histories of
lift coefficient computed by the original LU-SGS and the
improved LU-SGS with ε1=ε2=0.25 are compared in
Fig. 8. Clearly there is a significant improvement in
convergence rate with the improved LU-SGS, as the
number of time iterations to achieve the convergence
standard is decreased by half. Figure 9 shows a com-
parison of the computed surface pressure coefficients
with the experimental data, where Cp is the surface
pressure coefficient, x/c means a ratio of the x-coordinate Fig. 9 Comparison of computed surface pressure
and chord length c. As one can see, the computational coefficients on RAE2822 airfoil with experimental
data.
· 40 · WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 No.1

3.3. LANN wing vergence of aerodynamic coefficient, the efficiency of

the improved LU-SGS scheme with ε1=ε2=0.25 is two
The LANN wing is frequently used as an unsteady times that of the original LU-SGS scheme. Figure 13
flow demonstration case [17], while there are also some shows a comparison of the computed surface pressure
available experimental data for LANN wing in steady coefficients at 32.5% span with the corresponding ex-
case. In this paper, a steady case of LANN wing at perimental data. As can be seen from Fig. 13, the com-
α=2.6°, Ma=0.82, and Re=7.3×106 is selected to evaluate putational result has a very good agreement with the
the performance of the improved LU-SGS scheme in 3D experimental data.
computational mesh with large number of grid cells. A
hybrid unstructured mesh for LANN wing, which has
979 844 nodes and 2 840 153 cells, is generated with the
mesh generation methodology described in Ref. [16].
Figure 10(a) shows the distribution of the wing surface
mesh, which consists of 81 855 triangle cells and 40 983
vortices. Figure 10(b) presents a sliced distribution of the
spatial field mesh.

Fig. 11 Convergence histories of the maximal residual in

terms of iterations for transonic viscous flow past
LANN wing.

Fig. 12 Convergence histories of lift coefficient in terms of

iterations for transonic viscous flow past LANN wing.
Fig. 10 Hybrid unstructured mesh around the LANN wing.

Based on the above mesh, three computations are

conducted by using the original LU-SGS, the improved
LU-SGS with ε1=100 and the improved LU-SGS with
ε1=ε2=0.25 individually. It needs to be mentioned that
the maximal allowable integral CFL number for the
improved LU-SGS with ε1=100 is 12, while the corre-
sponding CFL number for original LU-SGS and the
improved LU-SGS with ε1=ε2=0.25 can be both given
as 30. The obtained convergence rates of the maximal
residual in terms of iteration numbers are compared in
Fig. 11 and the convergence histories of lift coefficient
are shown in Fig. 12. As shown in Fig. 11, the acceler-
ating effect on the convergence rate of residual of the Fig. 13 Comparison of computed surface pressure coeffi-
improved LU-SGS scheme is not significant if the 25% cients at 32.5% span section on LANN wing with
increment of CPU time is considered. But for the con- experimental data.
No.1 WANG Gang et al. / Chinese Journal of Aeronautics 25(2012) 33-41 · 41 ·

4. Conclusions [7] Nejat A, Ollivier-Gooch C. Effect of discretization order

on preconditioning and convergence of a high-order un-
In this paper, the main approximations introduced structured Newton-GMRES solver for the Euler equa-
tions. Journal of Computational Physics 2008; 227(4):
in the primary theory and implementation process of
2366-2386.
LU-SGS scheme are analyzed in detail. It is found that [8] Nigro N, Storti M, Idelsohn S, et al. Physics based
the approximate convective flux Jacobian based on the GMRES preconditioner for compressible and incom-
maximal eigenvalue splitting is helpful for the inner pressible Navier-Stokes equations. Computer Methods
linear iteration while the accurate convective flux in Applied Mechanics and Engineering 1998; 154(3-4):
Jacobian based on exact derivation is beneficial to the 203-228.
outer non-linear time marching. To make the advan- [9] Jameson A, Caughey D A. How many steps are re-
tages of them complement each other, an improved quired to solve the Euler equations of steady com-
LU-SGS scheme with a hybrid construction method for pressible flow: in search of a fast solution algorithm.
AIAA-2001-2673, 2001.
the implicit system matrix is developed. The hybrid
[10] MacCormack R W. Iterative modified approximate
way is that the accurate derivative of the convective factorization. Computers & Fluids 2001; 30(7-8): 917-
flux term is used to construct more accurate implicit 925.
system matrix on the cell faces which have small con- [11] Spalart P R, Allmaras S R. A one-equation turbulence
travariant velocity or contravariant velocity closing to model for aerodynamic flows. AIAA-1992-0439, 1992.
sonic speed; while the original Jacobian splitting [12] Li C N, Ye Z Y, Wang G. Simulation of flow separation
method remains on the other cell faces to reduce com- at the wing-body junction with different fairings. Jour-
putational complexity and ensure the diagonally domi- nal of Aircraft 2008; 45(1): 258-266.
nant property of the implicit system matrix. [13] Jameson A, Schmidt W, Turkel E. Numerical solution
A series of high Reynolds number viscous flows is of the Euler equations by finite volume methods using
Runge Kutta time stepping schemes. AIAA-1981-1259,
simulated to demonstrate the performance of the im-
1981.
proved LU-SGS scheme with hybrid unstructured [14] Roe P L. Characteristic based schemes for the Euler
mesh. The numerical results show that the improved equations. Annual Review of Fluid Mechanics 1986; 18:
LU-SGS scheme is significantly more efficient than the 337-365.
original LU-SGS scheme. [15] Liou M S. A sequel to AUSM, Part II: AUSM+-up for
all speeds. Journal of Computational Physics 2006;
References 214(1): 137-170.
[16] Wang G, Ye Z Y. Generation of three dimensional
[1] Shang J S. Computational fluid dynamics application to mixed and unstructured grids and its application in
aerospace science. The Aeronautical Journal 2009; 113: solving Navier-Stokes equations. Acta Aeronautica et
619-632. Astronautica Sinica 2003; 24(5): 385-390. [in Chinese]
[2] Tinoco E N, Bogue D R, Kao T J, et al. Progress to- [17] Zwaan I R J. LANN wing pitching oscillations, com-
ward CFD for full flight envelope. The Aeronautical pendium of unsteady aerodynamic measurements.
Journal 2005; 109: 451-460. AGARD-R-702, 1982.
[3] Yoon S, Jameson A. An LU-SSOR scheme for the Euler
and Navier-Stokes equations. AIAA Journal 1988; Biography:
26(9): 1025-1026.
[4] Chen R F, Wang Z J. Fast, block lower-upper symmet- WANG Gang received his M.S. and Ph.D. degrees of fluid
ric Gauss Seidel scheme for arbitrary grids. AIAA dynamics from Northwestern Polytechnical University in
Journal 2000; 38(12): 2238-2245. 2001 and 2005 respectively. He is currently an associate pro-
[5] Kim J S, Kwon O J. Improvement on block LU-SGS fessor at School of Aeronautics in Northwestern Polytechni-
scheme for unstructured mesh Navier-Stokes computa- cal University. His main research interests are computational
tions. AIAA-2002-1061, 2002. fluid dynamics and fluid-structure interaction.
[6] Zhang L P, Chang X H, Duan X P, et al. A block E-mail: [email protected]
LU-SGS implicit unsteady incompressible flow solver
on hybrid dynamic grids for 2D external bio-fluid
simulations. Computers & Fluids 2009; 38(2): 290-308.