Chinese Journal of Aeronautics 25 (2012) 33-41

Contents lists available at ScienceDirect

Chinese Journal of Aeronautics

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

journal of

aeronautics

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 derivative 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 numerical 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; turbulent flow

1. Introduction

Accompanied with the development of computational fluid dynamics (CFD), numerically solving the so-called Reynolds-averaged Navier-Stokes (RANS) equations has been widely used for viscous flow simulation in many industry fields. Especially in the aeronautical engineering, a large number of high Reynolds number viscous flows always need to be simulated

*Corresponding author. Tel.: +86-29-88491342. E-mail address: wanggang@nwpu.edu.cn Foundation item: National Natural Science Foundation of China (10802067)

1000-9361/$ - see front matter © 2012 Elsevier Ltd. All rights reserved. doi: 10.1016/S1000-9361(11)60359-2

during the design process of an aircraft [1-2]. In order to efficiently simulate high Reynolds number viscous flow with current limited computer hardware, highly stretched grids with large aspect ratios should be used to calculate the viscous effect in the boundary layer region. Typically, the magnitudes of the aspect ratios for those "viscous" grids are more than 103. On such type of grid, explicit time integration scheme has very 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 scheme can hardly fulfill the engineering demands on high Reynolds number viscous flow simulation. Hence, an implicit time integration scheme is highly desired for the enhancement of efficiency since it has the advantage of no stability limitation in theory. Although the

theoretically unconditional stability cannot usually be achieved because some approximations must be introduced in practical implementation, the boundary of stability region of the implicit scheme can still be significantly extended. So larger time steps are allowed to be used in an implicit scheme and faster convergence speed can be obtained compared with the explicit scheme.

Many implicit schemes [3-10] have been developed and successfully applied to accelerating the convergence rate. Two sorts of them are widely used on the unstructured grids, which are the Gauss-Seidel type method represented by the well-known lower-upper symmetric Gauss-Seidel (LU-SGS) scheme [3-6] and the Krylov-subspace method typically represented by the generalized minimal residual (GMRES) scheme [7-8]. The strong point of GMRES scheme is that quadratic convergence speed can be achieved as the time step size is close to infinite. While the precondition for obtaining quadratic convergence is that the implicit system matrix should be constructed and inversed exactly. To fulfill this requirement, an accurate linearization of residuals and the inner iterations with preconditioning methodologies should be used. Consequently, the GMRES scheme is characterized with highly computational complexity and large memory consumption. Another popular implicit scheme is LU-SGS, which does not add so much extra storage requirement compared with explicit scheme and is free from direct inversion of large implicit system matrix. Compared with the GMRES scheme, LU-SGS is easier to be implemented and needs less computational cost for marching one time step. With respect to the computational efficiency, the LU-SGS scheme is not less competitive than the GMRES method in terms of CPU time. However, it should be pointed out that LU-SGS usually needs more iteration steps to reach the same convergence level compared with GMRES method. Such a drawback of LU-SGS scheme is mainly caused by several approximations introduced in the process of constructing the implicit system matrix and LDU decomposition [4-5]. If these approximations could be reduced by some kinds of modification, then the efficiency of LU-SGS scheme will be further enhanced. This untapped capacity of LU-SGS scheme gives rise to the basic motivation of the present modification work.

The principal ingredients of this paper are as follows. The basic algorithm of LU-SGS scheme is briefly described firstly. According to the analysis of the main approximations introduced by the original LU-SGS scheme [3], an improved LU-SGS scheme based on a new implicit system matrix construction method is presented in this paper. In order to investigate the convergence performance of the improved LU-SGS scheme, 2D and 3D viscous flows around the NACA-0012 airfoil, the RAE2822 airfoil and the LANN wing are simulated using hybrid unstructured grids. A comparison between the improved LU-SGS scheme and the original LU-SGS scheme is also presented.

2. Numerical Schemes

2.1. Governing equation

The integral form of RANS equations enclosed with Spalart-Allmaras (S-A) turbulence model [11] for a control volume Q with a surface element dS can be written as

where W = [p pu pv pw pE pu]T is the vector of conserved quantities with p, u, v, w, E and u denoting the density, the Cartesian velocity components, the specific total energy and the working variable of S-A turbulence model, respectively. n represents the surface normal. The source term Qsource = [0 0 0 0 0 QT]T only lies in turbulence model equation. In which QT represents the source term of S-A turbulence model. Fc and Fv comprise the convective and viscous flux vectors respectively; their detailed form can be found in Ref. [12]. Since the spatial discretization and time integration of turbulence model equation and mean flow equations are carried out in a loosely coupled way, this makes the turbulence model equation can be easily solved because only one unknown variable is involved. Hence the numerical treatments of RANS equations are mainly considered in this paper.

2.2. Spatial discretization

By using the finite volume method for each grid cell i, the spatial discretization of Eq. (1) can be expressed as

d(n,W,) ^ ( )

-dt--AQsource - ^ ((Fc)im - (Fv)im )Sim (2)

dt meN(i)

where Qt represents the volume of current grid cell i,

N(i) the set of face neighbor cells of cell i, Sim the normal vector area of the face, the subscript im denotes the current face which is shared by cell i and cell m. In order to capture the discontinuity accurately and suppress numerical oscillations, a series of schemes has been developed for the evaluation of convective flux (Fc)m , for example the central scheme of Jamerson [13], the Roe scheme [14] and the AUSM type schemes [15]. However, all these schemes can be rewritten to the following unified form:

(Fc)im = 1 ((WW.) + Fc(Wm,R))+ ^lW^R) (3)

where WmL and Wim,R denote the flow states on the left and the right side of the face shared by cell i and cell m. In first order scheme, Wim,L and Wim,R are directly specified as the solution quantities located at the center of cell i and cell m, namely Wi and Wm. In second order scheme, some reconstruction techniques should be

used to determine W^l and W^,r. The D(Wim,L, Wm,r) in Eq. (3) represents the dissipation term added to prevent non-physical oscillations, and the iconic features of different convective flux evaluation schemes mainly lie in the specific forms of this term. The evaluation of the viscous flux term Fv is relatively simple. The viscous flux is usually discretized by the standard central scheme because of the elliptic nature of this term. This means that the flow quantities and their first derivatives on the faces of the control volume, which exist in the formula of viscous flux term, can be obtained by the simple average, and then be used to calculate Fv straightforwardly.

To simplify the expression, the residual of Eq. (2) is denoted here by

R W) = £ ((FJm - (Fv),m ) Sm + — Q^ce (4)

meN (i)

When fixed-volume computational grid is used, Eq. (2) becomes

a -JT=-R (w ) dt

2.3. Analysis of LU-SGSscheme

where, in unstructured grid case, M is a large, sparse and generally non-symmetric matrix, which consists of DOFxDOF blocks. Here DOF represents the degrees of freedom, which equals the total number of control volumes in finite volume method.

The LU-SGS scheme is based on the LDU factorization of the implicit system matrix, which is given by

(— I+M) = L+D+U= (L+D)D- (U+D) - LD^U At

where L only consists of block terms in the strictly lower triangular matrix, U only consists of block terms in the strictly upper triangular matrix and D is a block diagonal matrix.

Suppose that the LD- U is a small term and its contribution to the implicit system matrix can be neglected, then Eq. (8) can be replaced by the following equation:

(L + D) D^(U + D)AWn =- Rn

Equation (10) can be inverted in the following for ward and backward sweep procedure:

[dAW * =-Rn - LAW *

I DAWn = DAW* - UAWn

Equation (5) is a system of coupled ordinary differential equations in time. By using the backward Euler scheme for the implicit time integration, we obtain

W"+1 - W.K

= - R (W n+1)

where the superscript n represents the number of time level. Since Wn+1 is unknown on current time level, the residual R(Wn+1) cannot be evaluated directly. However, it can be linearized by using first order Taylor expansion in the following way:

R (Wn+1) * R (Wn) + X dRfW"n ) (W

jeC (i)

1 - W ") (7)

where C(i) is the set of cell i and its neighbor cells. With the definitions of

3Rt (W n)

AWn = W n+1 - Wn

R(Wn+ ) in Eq. (6) is substituted by the linearization term in Eq. (7), then the following implicit system is obtained:

(a I + M )AW n = - R(W n) At

where AW is solution vector updated in the forward sweep.

According to the above analysis, we can draw a conclusion that the property of implicit system matrix M has strong influence on the performance of LU-SGS implicit scheme. Hence the construction method and properties of M need to be considered in detail. By using Eq. (4), 3Rt (Wn)/dWn is evaluated as

d R, (Wn)

d (FC)im 5 (Fv)„

m e N(i)

(FX d (Fv)

dWjn v J

where d—Qsource)/dW-n only lies in the turbulence model equations and its expression is the derivative of the source term in terms of the conservative variables.

The viscous flux Jacobians (A+ )jJ- =d(Fv)ijjdWand

(A-) j =3(Fv)j/dWjn in Eq. (12) have complex expressions even if the thin shear layer approximation is used. In order to reduce the computational cost, the viscous flux Jacobian can be approximately counted by using its spectral radius:

(A±)j *±(lv)j I =

3Pj ' P

^ Ml + Mt A v PrL

I (13)

where y, n and Pr denote the ratio of specific heat coefficient, dynamic viscosity coefficient and Prandtl number respectively, and ry is the distance between grid cell i and cell j. In original LU-SGS approach, the con-vective flux Jacobians d(Fc) y/dW? and d(Fc)y/dWyn

are not defined with the exact derivatives based on Eq. (3), but approximately evaluated by the splitting of convective flux Jacobian with its maximal eigenvalue

ew; d(F0) ,

( A ), = |( A + AnuK I )

!(A- ) , = "(A -AmaxI)

where A = (dFjdW).., and A„

is the maximal ei-■ n,, I + c,. Here, the

genvalue of A, namely2max = V „v

subscript ij denotes the face which is shared by grid cell i and cell j, c is the sonic speed, Vis the velocity.

The (Ac+) j in Eq. (14) only has non-negative eigenvalue, and(Ac-) j non-positive eigenvalue. According to the theories on propagation of perturbation, (Ac+) j reflects the variations of fluxes caused by the increment of flow variables on grid cell i, and (A-) j the variations of

fluxes received from the fluctuation of flow variables on neighboring grid cell j. Hence, there has certain rationality for using (A+)ijand (Ac-)j. instead of d(Fc)j/dW"

and d(Fe)j/dWj".

In fact, there are two kinds of iterative processes contained in the implicit time marching scheme of this paper. One is the non-linear time marching iteration for solving Eq. (6) with backward Euler method, the other one is an inner iteration on each time level for solving the linear equation system Eq. (8) with LU-SGS method. In constructing an efficient implicit scheme, convergence and stability performances of the above two iteration processes need to be considered on balance.

The construction of implicit system matrix with Eq. (14) is favorable to the inner iteration. In detail, two distinctive benefits can be identified. Firstly, the block diagonal matrix D in Eq. (10) degenerates into a pure diagonal matrix, which can be evaluated and inversed with very little CPU time. Secondly, it ensures the diagonal dominance of implicit system matrix for very large time stepping. This benefit not only guarantees the convergence of LU-SGS iteration, but also makes the neglected term LD- U really negligible compared to the entire matrix system.

However, we all know that convective flux Jacobian matrix A has five eigenvalues:

A =A2 =A3 -K 'nv., V

K •n, -c, A5:

3 3 ' 5

K'ny +c (15)

Their signs express the propagating directions of waves relative to the surface normal Sim, and the abso-

lute values of the above eigenvalues represent the propagating speed of corresponding waves. It is clear that (Ac+ )j and (A-)j given by Eq. (14) can reflect the propagating directions of perturbations correctly. But they cannot represent the propagation velocities of perturbations accurately, which leads to considerable distortions in evaluating the variations of fluxes caused by the fluctuations of flow states. If the values of

(IV-n-K)/Vnl and (V"ji + Cj)/(!VnI-Cj)

and are very large, serious distortions on the variations of fluxes may be introduced by using Eq. (14). And such kind of cases occurs not only in the low Mach flow simulations but also commonly lies in the high Reynolds number flow computations. The reason is that the absolute value of Vj • np could be very small on the cell faces which are parallel to the wall boundary faces or near the bottom region of boundary layer. Hence, too many distortions may be introduced into implicit system matrix with the above construction method based on the maximal eigenvalue splitting, and these distortions are harmful for the convergence behavior of non-linear time marching iteration.

In order to construct more accurate implicit system matrix, d(Fc)iJ/dW" and d(Fc)p/dW/ need to be evaluated by using the derivation of Eq. (3), namely

d(Fc) 3

dFc(Wp,L)

dW,n dFc(Wv,L)

5D(W,,L,Wj,R)

dwj ,l dFc (W,,R

dFc (W,,R )

) dD(W,,L, - + - j

w , ,r)

where dD(Wj l, Wj r) /W l and 3D(W„ l,Wj r) /Wr are defined according to the spatial discretization scheme. For example, if Roe scheme is used, and let A Roe represent Roe average Jacobian matrix, then we get

dD(W,,L,W,,R) 1 |A

— Au

SD(W,,l , W,,R )

=-2 Ar

In fact, the construction of the implicit system matrix with the above accurate convective flux Jacobian not only brings a notable increase in computational cost, but also causes the implicit system matrix M to lose its property of diagonal dominance as the time step size At takes a great value. It is well known that the convergence of the LU-SGS inner iteration cannot be guaranteed if the system matrix is not diagonally dominant. In order to achieve diagonal dominance, a

relatively small At should be used in the time marching process. In such case, the evolution speed of flow field is considerably decreased and more time marching steps are required to achieve the convergence state. Another way for overcoming the limitation of diagonal dominance is giving up the LU-SGS scheme but using some other inner iteration algorithms which have better stability, for example the GMRES scheme. The cost of this choice is that the computational complexity will be increased significantly.

2.4. Improved LU-SGS scheme

In the above analysis, we find out that the approximate convective flux Jacobian based on maximal eigenvalue splitting is helpful for the inner linear iteration while the accurate convective flux Jacobian based 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 LU-SGS implicit time marching scheme will be improved further. For this purpose, a hybrid method for constructing the implicit system matrix is proposed and described in the following way.

For the cell faces on which the contravariant speed has small value or close to sound speed, d(Fc)ijjdWln

and d(Fc)pJdW" are computed by using Eq. (16) so

that the excessive approximations on implicit matrix could be avoided; on the other cell faces, the classical maximal eigenvalue splitting method, namely Eq. (14) is used to evaluate the above convective Jacobians with the aim of reducing computational cost and strengthening the diagonal dominance of implicit system matrix. According to this principle, the hybrid method for determining d(Fc)../ dWln and d(Fc)../ dW" can be expressed as

If |Vy'n^cy <ei or 1 - ^ < V'jcy < ^2 ,

d(Fc)j

otherwise,

,dFc(W,L) + 'DW, -Wr)

dFc(Wj.R)

^ « 1( A + A„ I )

dWn 2 max

S( ^ 1( A-Im. I )

In Eq. (18), e1 and e2 represent non-negative empirical coefficients. The smaller values are given to e1 and e2, the more characteristics of the original LU-SGS will be assigned to the improved LU-SGS method. On the other hand, the larger values are given to ei and e2, the more characteristics of the block LU-SGS method [4, 6] will be represented. For example, if e1 = 1 is used in subsonic

flow simulation, the improved LU-SGS will be identified with the block LU-SGS method. The optimal values of s1 and s2 are hard to be known beforehand because they vary with different flow conditions and meshes. According to the practical experience, s1 = s2 ~ 0.25 is recommended in high Reynolds number flow simulations.

Finally, it should be pointed out that the existence of viscous flux Jacobian, which is computed by using Eq. (13), offers favorable term for the diagonal dominance of implicit system matrix. This means that the side effects of the usage of Eq. (18) on the stability of LU-SGS inner iteration can be eliminated partly in the boundary layer region of the high Reynolds number flow. Hence, the maximum allowable time step size will not be decreased obviously if the above improved LU-SGS scheme is used.

3. Test Cases and Results

Several typical test cases are selected here to illustrate the practicality and computational efficiency of the improved LU-SGS implicit time marching scheme. These test cases include the low Mach number and supersonic viscous flows around NACA0012 airfoil, transonic viscous flow around RAE2822 airfoil and 3D transonic flow around LANN wing. The S-A turbulence model, second order Roe scheme and the local time stepping [12] are chosen as common settings for all computational cases. The convergence performances of original LU-SGS and improved LU-SGS are analyzed and compared in each test case.

3.1. NACA0012 airfoil

As shown in Fig. 1, a hybrid unstructured mesh for viscous flow simulation on NACA0012 airfoil is generated with the method described in Ref. [16]. It consists of 200 wall boundary nodes, 7 991 computational field nodes and 11 552 mesh cells. Near the boundary layer region, 22 layers of quadrangle mesh are distributed and the spacing of the first layer normal to the wall is chosen as 3.0*10-6 times of chord length.

Firstly, the incompressible viscous flow is computed with an angle of attack a = 5°, a Mach number Ma = 0.15 and a Reynolds number Re = 3*106. The same Courant-Friedrichs-Lewy (CFL) number 30 is used for the computations with the improved and the original LU-SGS schemes. The convergence histories of the maximal residual in terms of the time iteration numbers are shown and compared in Fig. 2. As one can see, the convergence speed of the maximal residual computed by the original LU-SGS has slowed down quite markedly after it decreased five orders. While in the computation with the improved LU-SGS, the convergence rate is constant. In practice, we are more concerned about the convergence speeds of aerodynamic coefficients. In Fig. 3 the convergence curves of life coefficient are plotted. The dash-dotted lines in this

Fig. 1 Hybrid unstructured mesh around NACA0012 airfoil.

Fig. 2 Convergence histories of the maximal residual in terms of iterations for low speed flow past NA-CA0012 airfoil.

Fig. 3 Convergence histories of lift coefficient in terms of iterations for low speed flow past NACA0012 airfoil.

figure represent the upper and lower convergence tolerance of lift coefficients, which are defined

L,conveigence

±0.002 in this paper. The results shown in Fig. 3 indicate that with the same convergence tolerance, the improved LU-SGS only cost half number of iterations as much as that of original LU-SGS scheme. Because of the low Mach free stream condition of this test case, the implicit system equation is purely con-

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 a = 2°, Ma =1.5 and Re= 6x106 is also simulated with the improved and the original LU-SGS schemes. The computed pressure contours are shown in Fig. 4. It is easy to be observed that a detached shock has already occurred in the flow field. The convergence histories of the maximal residual obtained by using the original LU-SGS, the improved LU-SGS with S\ = s2 = 0.25 and the improved LU-SGS with £i=100 are presented in Fig. 5. The main purpose of choosing the setting of e1=100 is to investigate the convergence performance of the improved LU-SGS on condition 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 e1=100 is 15.5, while the corresponding CFL number for the original LU-SGS and

Fig. 4

Pressure contours for supersonic viscous flow past NACA0012 airfoil.

5 Convergence histories of the maximal residual in terms of iterations for supersonic viscous flow past NACA0012 airfoil.

the improved LU-SGS with e1 = s2 = 0.25 can be both 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 boundary nodes and 22 layers of quadrangle mesh in the boundary layer region. The spacing of the first layer

results agree very well with the experimental data.

normal to the wall is chosen as 5.0> length.

10 times of chord

Fig. 6 Hybrid unstructured mesh around RAE2822 airfoil.

Here, a transonic case at a = 2.79°, Ma = 0.734 and Re = 6.5x106 is selected to evaluate the performance of the original LU-SGS and the improved LU-SGS with different values of e1 and e2. The convergence behavior of the maximal residuals in terms of iteration numbers is shown in Fig. 7. As one can see from Fig. 7, the convergence speed can be greatly enhanced by specifying proper value for e1 and e2. In current test case, the main enhancement comes from the function of e1 in the range from 0 to 0.1 and e2 only plays a secondary role. In general, the recommended values of e1=e2=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 e1=e2=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 comparison 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 and chord length c. As one can see, the computational

Fig. 7 Convergence histories of the maximal residual in terms of iterations for transonic viscous flow past RAE2822 airfoil.

Fig. 8 Convergence histories of lift coefficient in terms of iterations for transonic viscous flow past RAE2822 airfoil.

Fig. 9 Comparison of computed surface pressure coefficients on RAE2822 airfoil with experimental data.

3.3. LANN wing

The LANN wing is frequently used as an unsteady flow demonstration case [17], while there are also some available experimental data for LANN wing in steady case. In this paper, a steady case of LANN wing at a=2.6°, Ma=0.82, and Re=7.3x106 is selected to evaluate the performance of the improved LU-SGS scheme in 3D 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.

(b) Sliced mesh distribution

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 £i=£2=0.25 individually. It needs to be mentioned that the maximal allowable integral CFL number for the improved LU-SGS with £i=100 is 12, while the corresponding CFL number for original LU-SGS and the improved LU-SGS with £i=£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 accelerating effect on the convergence rate of residual of the improved LU-SGS scheme is not significant if the 25% increment of CPU time is considered. But for the con-

vergence of aerodynamic coefficient, the efficiency of the improved LU-SGS scheme with e1=e2=0.25 is two times that of the original LU-SGS scheme. Figure 13 shows a comparison of the computed surface pressure coefficients at 32.5% span with the corresponding experimental data. As can be seen from Fig. 13, the computational result has a very good agreement with the experimental data. o

Fig. 11 Convergence histories of the maximal residual in terms of iterations for transonic viscous flow past LANN wing.

0.40 r

« 0.25 -i

1---Original LU-SGS

. - Improved LU-SGS with e,

0.20 ■ -----Improved LU-SGS with s,=100

------Convergence tolerance

Iteration/10'

Fig. 12 Convergence histories of lift coefficient in terms of iterations for transonic viscous flow past LANN wing.

-1.2 r

0.8 *—I—I-1-1 I I I I-1-i—l—:-1

0 0.2 0.4 0.6 0.8 1.0

Fig. 13 Comparison of computed surface pressure coefficients at 32.5% span section on LANN wing with experimental data.

Iteration/101

4. Conclusions

In this paper, the main approximations introduced in the primary theory and implementation process of LU-SGS scheme are analyzed in detail. It is found that the approximate convective flux Jacobian based on the maximal eigenvalue splitting is helpful for the inner linear iteration while the accurate convective flux Jacobian based on exact derivation is beneficial to the outer non-linear time marching. To make the advantages of them complement each other, an improved LU-SGS scheme with a hybrid construction method for the implicit system matrix is developed. The hybrid way is that the accurate derivative of the convective flux term is used to construct more accurate implicit system matrix on the cell faces which have small con-travariant velocity or contravariant velocity closing to sonic speed; while the original Jacobian splitting method remains on the other cell faces to reduce computational complexity and ensure the diagonally dominant property of the implicit system matrix.

A series of high Reynolds number viscous flows is simulated to demonstrate the performance of the improved LU-SGS scheme with hybrid unstructured mesh. The numerical results show that the improved LU-SGS scheme is significantly more efficient than the original LU-SGS scheme.

References

[1] Shang J S. Computational fluid dynamics application to aerospace science. The Aeronautical Journal 2009; 113: 619-632.

[2] Tinoco E N, Bogue D R, Kao T J, et al. Progress toward CFD for full flight envelope. The Aeronautical Journal 2005; 109: 451-460.

[3] Yoon S, Jameson A. An LU-SSOR scheme for the Euler and Navier-Stokes equations. AIAA Journal 1988; 26(9): 1025-1026.

[4] Chen R F, Wang Z J. Fast, block lower-upper symmetric Gauss Seidel scheme for arbitrary grids. AIAA Journal 2000; 38(12): 2238-2245.

[5] Kim J S, Kwon O J. Improvement on block LU-SGS scheme for unstructured mesh Navier-Stokes computations. AIAA-2002-1061, 2002.

[6] Zhang L P, Chang X H, Duan X P, et al. A block LU-SGS implicit unsteady incompressible flow solver on hybrid dynamic grids for 2D external bio-fluid simulations. Computers & Fluids 2009; 38(2): 290-308.

[7] Nejat A, Ollivier-Gooch C. Effect of discretization order on preconditioning and convergence of a high-order unstructured Newton-GMRES solver for the Euler equations. Journal of Computational Physics 2008; 227(4): 2366-2386.

[8] Nigro N, Storti M, Idelsohn S, et al. Physics based GMRES preconditioner for compressible and incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering 1998; 154(3-4): 203-228.

[9] Jameson A, Caughey D A. How many steps are required to solve the Euler equations of steady compressible flow: in search of a fast solution algorithm. AIAA-2001-2673, 2001.

[10] MacCormack R W. Iterative modified approximate factorization. Computers & Fluids 2001; 30(7-8): 917925.

[11] Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows. AIAA-1992-0439, 1992.

[12] Li C N, Ye Z Y, Wang G. Simulation of flow separation at the wing-body junction with different fairings. Journal of Aircraft 2008; 45(1): 258-266.

[13] Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes. AIAA-1981-1259, 1981.

[14] Roe P L. Characteristic based schemes for the Euler equations. Annual Review of Fluid Mechanics 1986; 18: 337-365.

[15] Liou M S. A sequel to AUSM, Part II: AUSM+-up for all speeds. Journal of Computational Physics 2006; 214(1): 137-170.

[16] Wang G, Ye Z Y. Generation of three dimensional mixed and unstructured grids and its application in solving Navier-Stokes equations. Acta Aeronautica et Astronautica Sinica 2003; 24(5): 385-390. [in Chinese]

[17] Zwaan I R J. LANN wing pitching oscillations, compendium of unsteady aerodynamic measurements. AGARD-R-702, 1982.

Biography:

WANG Gang received his M.S. and Ph.D. degrees of fluid dynamics from Northwestern Polytechnical University in 2001 and 2005 respectively. He is currently an associate professor at School of Aeronautics in Northwestern Polytechnical University. His main research interests are computational fluid dynamics and fluid-structure interaction. E-mail: wanggang@nwpu.edu.cn