- Email: [email protected]
Implicit gas-kinetic BGK scheme with multigrid for 3D stationary transonic high-Reynolds number flows
Computers & Fluids 66 (2012) 21–28
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Implicit gas-kinetic BGK scheme with multigrid for 3D stationary transonic high-Reynolds number flows Jin Jiang, Yuehong Qian ⇑ Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, PR China
a r t i c l e
i n f o
Article history: Received 1 August 2011 Received in revised form 4 March 2012 Accepted 29 April 2012 Available online 16 May 2012 Keywords: Gas-kinetic BGK scheme Computational efficiency Acceleration techniques 3D steady transonic viscous flows
a b s t r a c t Instead of solving the Euler and Navier–Stokes equations directly, the gas-kinetic BGK scheme based on the Boltzmann equation has been developed and attracted more and more attentions since the early 1990s. It shows high accuracy and robustness for a wide range of flow regimes. But an obvious disadvantage of the BGK scheme is the low computational efficiency, in particular for multidimensional problems. Till now it has not been widely used as a practical tool for science and engineering applications. To overcome this drawback, in this paper some acceleration techniques, including local time stepping, implicit LU-SGS method, and multigrid strategy, are implemented into the original BGK approach and the new scheme is applied to study 3D steady transonic viscous flows. Numerical results show the significant speed up of the scheme to capture the steady state solution. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Transonic aerodynamics plays an important role in the operation of long range aircrafts [1]. The cruise speed of most civil airplanes is in the transonic regime. The physical natures associated with transonic flow, including shock wave, discontinuity, transition from laminar to turbulent flow, etc., are still not fully understood. Mathematically, the governing equations are nonlinear and hybrid, which show elliptic feature in the subsonic region and hyperbolic property in the supersonic parts. Thus it is hard to find analytical solutions to predict the flow patterns and then to guide the airplane designing. Following the pioneering work of Murman and Cole [2] in the early 1970s, many numerical methods for transonic aerodynamics have been developed, from the small disturbance equation to the full potential equation, the Euler equations, and the Navier–Stokes equations. Some brief reviews can be found in [3,4]. With the ever-increasing high performance of super-computer and sophisticated software, Computational Fluid Dynamics (CFD) has been considered as a convenient and useful tool for aerodynamic analysis and designing. However, to simulate high-Reynolds number flows or around complex configurations flows is still a challenge [5]. The main difficulty is not the vast need of computing resources, but the physical modeling and numerical approaches. Gas-kinetic BGK scheme (GKS) [6,7] has been developed for compressible flows since the early 1990s. Since then, many ⇑ Corresponding author. Tel.: +86 021 56331974 (O), mobile: +86 13636443941. E-mail address: [email protected] (Y. Qian). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.04.029
progress has been obtained on the validation, analysis and improvement of the GKS during the last two decades. There are numerous literatures on this topic. Xu [8] compared the gas-kinetic scheme with the Godunov method on the accuracy, efficiency, robustness, and claimed that the BGK method is more physical than the latter. Ohwada [9] explained how to construct the kinetic scheme for evolutionary equations and investigated the theoretical background and numerical errors of the scheme. Xu et al. [10] simulated hypersonic viscous flows by a multidimensional BGK flow solver with implicit LU-SGS method. Su et al. [11] and Xu et al. [12] extended the GKS to low Mach number flows. May et al. [13] improved the efficiency and convergence of the BGK scheme and presented a simplified scheme to reduce the CPU time. Issues with respect to the construction of high-order scheme [14], the extension to rarefied flows [15] and MHD [16] have also been published. Based on its constructive modeling, the GKS has a sound physical basis and can be employed as a robust and reliable flow solver to simulate problems from low Mach number flows to hypersonic flows. Nevertheless the GKS has an obvious drawback in terms of its weak computational efficiency and expensive computational cost, which limits its practical applications. Kim et al. [17] stated that in comparison with the Riemann solver, the ratio of CPU time between convective upwind and split pressure (CUSP) scheme and standard BGK scheme per grid point per step is CUSP:BGK = 1:1.908. To overcome this disadvantage, in this paper we will present an implicit multigrid BGK scheme and extend it to solve 3D high-Reynolds number flow problems.
22
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
and can be determined through the relationships
2. Numerical method
(
2.1. Gas-kinetic BGK solver
Ql ¼ r
Q ¼ After introducing BGK [18] model into 3D Boltzmann equation, in the absence of external force term, we get:
@f @f @f @f @f @f f g þu þv þw ¼ þ ub ¼ ; @t @x @y @z @t @xb s
ð1Þ
where f is the gas distribution function and g is the equilibrium state. s is the relaxation time. Both f and g are functions of time t, particle velocities ub = (u, v, w) (b = 1, 2, 3), space xb = (x, y, z) and internal variable n. The equilibrium state g is Maxwellian
Kþ3 k 2
g¼q
p
e
2
2
2
kððuUÞ þðv VÞ þðwWÞ þn
2
Þ;
q
1
C B B qU C Z C B C Q ¼B B qV C ¼ f wa dN; C B @ qW A qE 0
Fq
a ¼ 1; 2; 3; 4; 5;
ð2Þ
1 ð3Þ
b
ð4Þ
To simulate multidimensional problems, the directionalsplitting method has been adopted. Here, the x-directional equation, for example, is taken:
@f @f f g : þu ¼ @t @x s
ð5Þ
Following Xu’s [7] suggestion, the initial gas distribution f0 at control volume (CV) interface xi+1/2 can be expressed as
f0 ¼
8 l < @Q l ¼ Q iþ1=2 Q i ¼ R g l al w dN; a @x Dxl : @Q r Q iþ1 Q riþ1=2 R r r ¼ ¼ g a wa dN; @x Dxr
ð9Þ
where Dx is the distance from the CV center to the interface. To avoid the numerical oscillation close to the interface, the conservation flow variables Q liþ1=2 and Q riþ1=2 are reconstructed by MUSCL interpolation [19]:
(
Q liþ1=2 ¼ Q i þ 12 uðQ iþ1 Q i ; Q i Q i1 Þ;
where u is a limiter function uðy1 ; y2 Þ ¼
ð10Þ
y1 ðy22 þÞþy2 ðy21 þÞ y21 þy22 þ2
. The param-
eter is a small constant with the order being proportional to the local grid scale. On the fact that the non-equilibrium terms of Eq. (6) have no contributions to the conservative variables, we know:
(R R
g l ðal u þ Al Þwa dN ¼ 0;
ð11Þ
g r ðar u þ Ar Þwa dN ¼ 0:
By solving Eq. (11) the temporal slope A will be obtained. The equilibrium state g is assumed to have two slopes too [7]
r ðx xiþ1=2 Þ þ At; þ H½x xiþ1=2 a
H½x ¼
ðg f Þwa dN ¼ 0:
(
ð8Þ
g r wa dN;
g l ½1 þ al ðx xiþ1=2 Þ sðal u þ Al Þ;
x < xiþ1=2 ;
g r ½1 þ ar ðx xiþ1=2 Þ sðar u þ Ar Þ; x > xiþ1=2 ;
u2 þ v 2 þ w2 þ n2 ¼ aa wa ; 2
1; x P 0; 0; x < 0:
ð13Þ
has the similar form of a: a ¼a a wa and can be determined by the a relation of
8 0 < Q iþ1=2 Q i ¼ R g a l 0 wa dN; Dxl 0 : Q iþ1 Q iþ1=2 R r wa dN; ¼ g0 a Dxr
ð14Þ
where Q 0iþ1=2 is the vector of the initial conservative values at CV interface xi+1/2, which can be calculated by
Q 0iþ1=2 ¼
Z
Z Z
g 0 wa dN ¼
u>0
g l wa dN þ
Z Z
u<0
g r wa dN:
ð15Þ
So far, the only unknown coefficient is A. The general solution f of BGK Eq. (1) can be expressed as:
1 f xiþ1=2 ; t; ub ; n ¼
s
Z
t
0 0 g xiþ1=2 uðt t 0 Þ; t 0 ; ub ; n eðtt Þ=s dt
0
þ et=s f0 ðxiþ1=2 ut; 0; ub ; nÞ;
ð16Þ
Substituting Eqs. (6) and (12) into the above formula, after integrating, we get the gas distribution function f:
f ðxiþ1=2 ; t; ub ; nÞ ¼ ð1 et=s Þg 0 þ sðt=s 1 þ et=s ÞAg 0 l H½u þ a r ð1 þ ðsð1 þ et=s Þ þ tet=s Þða
ð6Þ
here the indexes for the other two directions (j, k) are omitted. The superscripts 0 l0 and 0 r0 denote the left- and right-side of the interface xi+1/2. a is the spatial slope of f, which has the form as:
a ¼ a1 þ a2 u þ a3 v þ a4 w þ a5
ð12Þ
where H[x] is the Heaviside function
where dN = du dv dw dn with dn = dn1dn2 dnK, Twa are the collision invariants wa ¼ 1; u; v ; w; 12 ðu2 þ v 2 þ w2 þ n2 Þ , and E is the total energy per unit mass. According to the conservation laws, the right hand side of Eq. (1) at any point in space and time should satisfy the compatibility conditions:
Z
g l wa dN;
l ðx xiþ1=2 Þ g ¼ g 0 ½1 þ ð1 H½x xiþ1=2 Þa
BF C B qU C Z C B C Fb ¼ B B F qV C ¼ ub f wa dN; C B @ F qW A F qE
R
Q riþ1=2 ¼ Q iþ1 12 uðQ iþ1 Q i ; Q iþ2 Q iþ1 Þ;
where q is the density, k = 1/2RT with the gas constant R and the temperature T. K denotes the number of internal freedom of n and equals to (5 3c)/(c 1), here c is the ratio of specific heats. (U, V, W) are the macroscopic velocities along x-, y- and z-directions. And n2 = n1n1 + n2n2 + + nKnK. Details on derivation from Eq. (1) to the Euler and Navier– Stokes equations can be found in [8]. Once the gas distribution is known, the macroscopic flow variables and their fluxes (including convection term and viscous term) can be obtained by the following integrals:
0
R
ð7Þ
H½uÞÞug 0 set=s ðAl H½ug l þ Ar ð1 H½uÞg r Þ þ et=s ðð1 uðt þ sÞal ÞH½ug l þ ð1 uðt þ sÞar Þð1 H½uÞg r Þ:
ð17Þ
Then substituting Eqs. (12) and (17) into the compatibility conditions (4), a simplified algorithm for A is given as below [20,21]
23
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
Z
Ag 0 wa dN ¼
Z
ðal g l H½u þ ar g r ð1 H½uÞÞuwa dN:
ð18Þ
Till now, all parameters are obtained. Taking the moments of BGK equation
@Q @Fb þ ¼ 0; @t @xb
ð19Þ
where,
4 c l lt DS2i t1 ¼ ðjV ni j þ cs ÞDSi þ C max ; ; þ 3q q Pr Prt X 4 c l lt DS2j ; ; þ t2 ¼ ðjV nj j þ cs ÞDSj þ C max 3q q Pr Prt X 2 4 c l l t D Sk ; ; þ t3 ¼ ðjV nk j þ cs ÞDSk þ C max 3q q Pr Prt X
ð25Þ
we could update the macroscopic flow variables by explicit or implicit time marching method. The Prandtl number Pr of the BGK scheme has a unit value and does not match the physical observation. There are many attempts have been proposed to modify it. In this paper, Xu’s treatment [7] is employed.
where CFL is the Courant–Friedrichs–Lewy number, X denotes the volume of CV (i, j, k), (DSi, DSj, DSk) stand for areas along i-, j- and k-directions, (ni, nj, nk) are unit normal vectors, V = (U, V, W) defines the components of macroscopic velocity, cs is the local sound speed, Prt is the turbulent Prandtl number, and C is equal to 4 for the Navier–Stokes equations while 0 for the Euler equations.
2.2. Relaxation time
3.2. Implicit LU-SGS method
Based on the Chapman–Enskog expansion of the BGK equation, the connection among the relaxation time s, the viscosity coefficient l, and the pressure p is obtained:
Implicit scheme is utilized in order to use large CFL number. In [10,26,27], Xu et al. argued the advantages of available implicit time integration method. Here we follow the way in traditional CFD and show how to implement LU-SGS method [28] briefly as below. Defining the residual as
l ¼ sp:
ð20Þ
Thus the relaxation time is calculated as
s ¼ l=p:
ð21Þ
RðQ Þ ¼
Z X
@Fb dX ¼ @xb
Z
ðr HÞdX ¼
X
I
H dS;
ð26Þ
@X
In order to keep the stability of current scheme, numerical dissipation is added through
then we get the semi-discrete finite volume scheme of Eq. (19)
l p pr s ¼ þ Dt l r ; p p þp
X
l
ð22Þ
where Dt is the time step for time advancing method. 2.3. Turbulence model As the increasing of Reynolds number, turbulent effect must be taken into account. Generally, there are two ways to solve this problem [22]: one is to develop one-equation or equations to describe the turbulence field based on the kinetic theory; while another one is to incorporate engineering turbulence model into the BGK method directly. Here we follow the second approach. Three different models, including Baldwin–Lomax (BL) model [23], Spalart–Allmaras (SA) model [24] and Menter’s SST (SST) model [25], have been implemented and tested. The effect of turbulent viscosity lt is included in the determination of the relaxation time of the BGK scheme:
s¼
l þ lt p
l p pr : þ Dt l r p þp
ð23Þ
In multigrid formulation, the SA model or SST model will be solved only on the finest grid, and then the turbulent viscosity is transferred to the coarser grids, where it is frozen. Furthermore, the BL model will be calculated on all grid levels. 3. Convergence acceleration techniques for steady flows
dQ þ RðQ Þ ¼ 0; dt
where H is the vector of time-averaged numerical fluxes on x-, yH and z-directions. In discrete space, the integral @ X H dS equals the sum of fluxes crossing the bounded boundaries of CV. Integrating implicitly and linearizing the residual, we obtain:
X
n dQ nþ1 @RðQ Þ DQ ¼ RðQ n Þ; þ @Q dt
8 @ HdS > > > A ¼ ; > @Q > > i < @ HdS B ¼ @Q ; > j > > > > > @ HdS :C ¼ ; @Q k
ðL þ D þ UÞDQ ¼ RðQ n Þ;
8 þ þ þ > L ¼ A þ B þ C ; > i1 j1 k1 < þ þ þ X > D ¼ Dt þ Ai Ai þ Bj Bj þ Ck Ck ; > : U ¼ Aiþ1 þ Bjþ1 þ Ckþ1 ; 8 1 > < A ¼ 2 ðA rA IÞ; r A P rA ; B ¼ 12 ðB r B IÞ; r B P rB ; > : 1 C ¼ 2 ðC r C IÞ; r C P rC ;
;
ð24Þ
ð30Þ
with
For stationary problems, the possible largest time step for each CV can be used to accelerate the solution. Here we employ the maximum time step for each CV (i, j, k) as
X
ð29Þ
expanding the second term in the bracket on the left hand side of Eq. (28), and splitting the Jacobian using upwind-bias method, the final discretization form is gotten
and
t1 þ t2 þ t3
ð28Þ
where DQ = Qn+1 Qn, superscript 0 n0 denotes the nth time level. After introducing the Jacobian matrices,
3.1. Local time stepping
Dt ¼ CFL
ð27Þ
ð31Þ
ð32Þ
where I is the unit matrix, ðrA ; rB ; rC Þ are the spectral radii of the Jacobian matrices.
24
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28 Table 1 Comparisons on CPU time, work units and factor of speedup of different numerical approaches for ONERA M6. The left column under 0 CPU time0 is the total cost and the right one is the time per step. 0 V20 and 0 V30 stand for 2-level and 3-level V-cycle multigrid. SST model is employed.
Y X Z
Numerical scheme
CPU time (s)
LU-SGS JST LU-SGS JST/V3 Explicit GKS LU-SGS GKS LU-SGS GKS/V2 LU-SGS GKS/V3
43,418 11,357 93,874 41,258 18,633 12,115
6.76 9.26 6.04 6.94 11.79 11.99
Work units
Speedup
6420 1226.875 15,535 5945 1581 1010.3125
1 5.23 – 1 3.76 5.88
2. Transfer the flow variables to coarse grid by volume-weighted method
Q 02h ¼ Fig. 1. Computational mesh for ONERA M6 wing in detail. The red grid is the symmetry plane of k = 1 and the green one is the wing surface grid. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
X
Rh ðQ h Þ R2h Q 02h ;
Rnew 2h ðQ 2h Þ ¼ R2h ðQ 2h Þ þ P 2h :
Explicit/single grid LU-SGS/single grid LU-SGS/2-level multigrid LU-SGS/3-level multigrid
-1
-2
0 dQ Nh ¼ Q m Nh Q Nh :
-3
ð39Þ
Then, get the new solution on the finest grid by interpolating the coarse-grid correction successively
Q new ðN1Þh ¼ Q ðN1Þh þ I Nh
ðN1Þh
0
4000
8000
12000
16000
Fig. 2. Convergence history of continuity equation in GKS framework by different time marching methods with single grid or multigrid for transonic viscous flow over ONERA M6 wing: C–O grid; CFL = 0.8 for explicit scheme and CFL = 5 for LU-SGS scheme; SST model.
Then, factorize the left hand side of Eq. (30) approximately into
ðL þ DÞD1 ðD þ UÞDQ ¼ RðQ n Þ;
ð33Þ
and use two-step sweeping way to get the solution DQ:
ðL þ DÞDQ ¼ RðQ n Þ; ðD þ UÞDQ ¼ DDQ :
¼ Q n þ DQ :
ð40Þ
where INh is the prolongation operator. 5. Repeat the work done from step 1 to step 4 until the solution is convergent. 4. Numerical cases To validate and analyze the accuracy and robustness of the current implicit gas-kinetic BGK approach with multigrid, two test cases have been carried out. 4.1. ONERA M6 wing
ð34Þ
Sequently, the macroscopic flow variables are updated by
Q
dQ Nh ;
ðN1Þh
Work Units
nþ1
ð38Þ
Finally, calculate the correction and update the solution on the coarse grid in the same way as the fine grid. 3. Do what have been done in step 2 on the coarser grid. 4. Obtain the coarse-grid correction after m-iteration (in this paper, m is set to 2N1)
-4
ð37Þ
and get the residual on the coarse grid
0
-5
ð36Þ
then collect the residual over the eight CVs (for 3D Cartesian grid) on the fine grid for the coarse grid. Define a forcing function as
P 2h ¼
Log10 (Residual)
.X X ðQ h Xh Þ Xh ;
ð35Þ
3.3. Multigrid Multigrid has been proved as a powerful acceleration tool for stationary flows. The following gives a brief procedure of V-cycle multigrid [29,30], where the subscript 0 Nh0 denotes the Nhth level grid. 1. Calculate the approximate solution Qh and residual Rh(Qh) on the fine grid.
The flow conditions are Ma = 0.8395, angle of attack (AOA) of 3.06° and mean-chord based Reynolds number of Re = 11.72 106. Computational grid is C–O type of i j k = 193 49 33 (seen in Fig. 1). Firstly, we compare the convergence rate between implicit LU-SGS method and explicit time stepping technique and check the effect of multigrid formulation. Fig. 2 illustrates the time history of residual of continuity equation. It shows that (1) the implicit LU-SGS method is faster to reach the convergence state than explicit one; (2) convergence is accelerated significantly by applying multigrid to current BGK scheme. Table 1 lists the comparisons on CPU cost, work units and speedup of various approaches. The needed CPU time per step of implicit GKS is nearly the same of that of implicit JST scheme [31]. When compared with the single grid calculation, even more time is spent for each step, the total CPU
25
1
1
0.5
0.5
0.5
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.20
-0.5
-1
0
0.2
0.4
0.6
0.8
-Cp
1
-Cp
-Cp
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.44
-0.5
-1
1
0
0.2
0.4
0.6
0.8
-0.5
-1
1
1
0.5
0.5
0.5
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.80
0.2
0.4
0.6
0.8
-Cp
1
-1 0
0
0.2
0.4
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.90
-0.5
-1 0
1
0.2
0.4
x/c
0.6
0.8
1
x/c
1
-0.5
Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.65
x/c
-Cp
-Cp
x/c
0
0.6
0.8
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.95
-0.5
-1
1
0
0.2
0.4
x/c
0.6
0.8
1
x/c
Fig. 3. Comparison of pressure distributions at selected span-wise locations of ONERA M6 from implicit GKS and implicit JST scheme with various turbulence models and experimental data [32]. c is the chord length, Z is the distance to the symmetry plane and b is the semi-span. V3: 3-level V-cycle multigrid.
higher resolution, better capability on shock capturing than JST scheme and (2) both SA model and SST model are more suitable for this case than BL model. We also present some primitive results with respect to grid independence. Three meshes are tested: one is that described above and marked as CO-I; the second is a finer C–O type grid CO-II of dimensions i j k = 289 64 49; and the last one is C–H type of i j k = 161 51 55 with 31 grids imposed on the wing surface along span direction. Fig. 4 shows the pressure distributions at some selected span-wise locations. Here BL model is used. With different grids the solutions are quite the same. The structure of two branches of the shock at 80% span location is verified clearly only by GKS on mesh C–H. This illustrates the solution dependence on the computational grid. On C–O mesh k-planes (x–y
1
1
0.5
0.5
0
Exp. GKS, mesh C-H JST, mesh C-H GKS, mesh CO-I JST, mesh CO-I GKS, mesh CO-II JST, mesh CO-II Z/b = 0.65
-0.5
-1
0
0.2
0.4
0.6
x/c
0.8
1
-Cp
-Cp
cost of multigrid method is much lower. According to this case, the more levels used the higher speedup is obtained. With 3-level V-cycle multigrid, the factor of speedup is up to about 6 for GKS, in contrast to above 5 for JST scheme. Secondly, we focus on the influence of turbulence models and analyze the resolution of shocks. Fig. 3 plots the pressure distributions at selected span-wise locations. All solutions are obtained with CFL = 5. The overall numerical results are in good agreement with the experimental data [32]. But we can see GKS captures the shocks more sharply than JST scheme. Beside 20% span location, the positions of shock wave given by GKS are closer to the experimental data [32] than results of JST scheme. In the GKS framework, SA model and SST model can predict more accurate shock positions than BL model. After all we can say (1) GKS shows
0
Exp. GKS, mesh C-H JST, mesh C-H GKS, mesh CO-I JST, mesh CO-I GKS, mesh CO-II JST, mesh CO-II Z/b = 0.80
-0.5
-1
0
0.2
0.4
0.6
0.8
1
x/c
Fig. 4. Comparison of pressure distributions at selected span-wise locations of ONERA M6 from implicit GKS and implicit JST scheme on mesh CO-I, mesh CO-II and mesh C–H with experimental data [32]. c is the chord length, Z is the distance to the symmetry plane and b is the semi-span.
26
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28 Y
Y X
X
Z
Z
Fig. 5. Schematic surface grid for ARA M100 wing/body configuration and the symmetry plane (red). The body grid (green) covers a region of i j k = (33–289) 1 (1–17), while the wing surface grid (blue) of i j k = (65–257) 1 (17–49). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
schematic surface grid and symmetry plane are illustrated in Fig. 5, which is C–O type of dimensions of i j k = 321 57 49. Flow field is solved by GKS with BL model and SST model. Comparison of pressure distributions on wing surface at various span-wise locations from GKS, CFL3D [33] and experimental data [34] is plotted in Fig. 6. The turbulence model is the key point for this case. Current GKS approach coupled with SST turbulence model gives solutions in good match to experimental data [34]. The shocks gotten with BL model situate too far aft than results with the other two turbulence models, which can be confirmed from the flow patterns on the upper wing surface (see Fig. 7). When SST model was employed, the separated streamline is located
planes) are not totally parallel to the free stream and thus the solutions selected at these special wing sections are not the precise ones we wanted to observe. By projecting or interpolating the neighboring data, the results may be improved and become better. Similarly, the good shock capturing capability of GKS is validated again. 4.2. ARA M100 wing/body Flow around a complex configuration of ARA M100 wing/body at free stream condition Ma = 0.8027, AOA = 2.873° and chord based Reynolds number of Re = 13.1 106 is simulated. The
1
1
1
0.5
0.5
0.5
-Cp
1.5
-Cp
1.5
-Cp
1.5
0
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5 Z/b = 0.123
-1
0
0.2
0.4
0.6
0.8
Z/b = 0.231
-1
1
0
0.2
0.4
x/c
0.6
0.8
Z/b = 0.325
-1
1
1
0.5
0.5
0.5
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
0.6
x/c
0.8
1
0.8
1
Exp. GKS, BL model GKS, SST model CFL3D, SA model
Exp.
GKS, BL model GKS, SST model CFL3D, SA model
-0.5 Z/b = 0.455
0.4
0.6
0
0
0.2
0.4
-Cp
1
-Cp
1
-Cp
1.5
0
0.2
x/c
1.5
-1
0
x/c
1.5
-0.5
Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5
-0.5 Z/b = 0.633
-1
0
0.2
0.4
0.6
x/c
0.8
1
Z/b = 0.817
-1
0
0.2
0.4
0.6
0.8
1
x/c
Fig. 6. Comparison of pressure distributions at selected span-wise locations of ARA M100 wing/body from GKS, CFL3D [33] and experimental data [34]. c is the chord length, Z is the distance to the symmetry plane and b is the wing span.
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
27
Fig. 7. Extreme streamlines on the upper wing surface of ARA M100 wing/body in the context of pressure contours. Red dash-dot-dot line is a reference of i = 209. Pink solid lines identify the various stations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 2 Comparison of lift CL and drag CD for ARA M100 wing/body among different turbulence models.
CL CD
BL model/GKS
SST model/GKS
SA model/CFL3D [33]
0.7358 0.04920
0.6768 0.04430
0.6632 0.04855
upstream while the reattached line is further downstream, that forms a wider reverse flow regime than that with BL model and leads to great variances in lift and drag (listed in Table 2). Those differences may be caused by the distinct natures of the turbulence models: (1) BL model is an algebraic two-layer model developed from boundary layer theory. Therefore it cannot trace the time evolution of the turbulent viscosity, and is only valid for attached flows or small separated flows. (2) SST model is developed on the transportation theory of the turbulence field. It is suitable for flows involving adverse pressure gradients and has ability to predict pressure-induced separation [25]. The SST model is the preferable one for future applications.
1. multigrid optimization, such as strategy (V-cycle or W-cycle), number of iterations, restriction/prolongation operators, and so on, 2. comparison and option of turbulence models, 3. effectiveness of the acceleration techniques for unsteady problems.
Acknowledgments The present work is supported partly by Program for Changjiang Scholars and Innovative Research Team in University, China (No. IRT0844) and Shanghai Science and Technique Committee (No. 11XD1402300). The work of Jiang is supported by Innovative Foundation for Graduate Students of Shanghai University too. The authors are grateful to Dr. Aiming Yang at Fudan University for his fruitful discussion and helpful suggestions on the implementation of multigrid. We appreciate Prof. Kun Xu, at HongKong University of Science and Technology, for his kindness and communication on BGK scheme. We thank Dr. Jianping Meng for discussion. In addition, the authors would also like to thank the reviewers for their constructive comments and to refine the current manuscript carefully.
5. Conclusions In this paper, we develop and extend the GKS into implicit multigrid formulation for high-Reynolds number steady flows. Some convergence acceleration techniques, including local time stepping, implicit LU-SGS method, and multigrid formulation, are proposed to the original scheme presented by Prendergast and Xu [6]. It is encouraged to improve the computational efficiency by applying those techniques into GKS for its significant acceleration to achieve convergence steady state and dramatic reduction of CPU time. At the same time, we incorporate three turbulence models (BL model, SA model and SST model) into current GKS code and check the validity of those models in combination with GKS. Based on the results of the two 3D test cases, SST model is identified as a best choice for further investigation. Additionally, there are still some problems need to be solved in the future:
References [1] Jameson A. Computational transonics. Commun Pure Appl Math 1988;41: 507–49. [2] Murman EM, Cole JD. Calculation of plane steady transonic flows. AIAA J 1971;9:114–21. [3] Jameson A, potential Full. Euler and Navier–Stokes schemes in applied computational aerodynamics. Prog Astronaut Aeronaut 1990;125:39–88. [4] Caughey DA, Jameson A. Development of computational techniques for transonic flows: an historical perspective; 2002. [5] Jameson A. Essential elements of computational algorithms for aerodynamic analysis and design. NASA/CR-97-206268. ICASE report no. 97-68. [6] Prendergast KH, Xu K. Numerical hydrodynamics from gas-kinetic theory. J Comput Phys 1993;109:53–66. [7] Xu K. A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method. J Comput Phys 2001;171:289–335. [8] Xu K. Gas-kinetic schemes for unsteady compressible flow simulations. VKI report 1998-03. [9] Ohwada T. On the construction of kinetic schemes. J Comput Phys 2002;177:156–75.
28
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
[10] Xu K, Mao ML, Tang L. A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow. J Comput Phys 2005;203:405–21. [11] Su MD, Xu K, Ghidaoui MS. Low-speed flow simulation by the gas-kinetic scheme. J Comput Phys 1999;150:17–39. [12] Xu K, He XY. Lattice Boltzmann method and gas-kinetic BGK scheme in the low Mach number viscous flow simulations. J Comput Phys 2003;190:100–17. [13] May G, Srinivasan B, Jameson A. An improved gas-kinetic BGK finite-volume method for three-dimensional transonic flow. J Comput Phys 2007;220: 856–78. [14] Li QB, Xu K, Fu S. A high-order gas-kinetic Navier–Stokes flow solver. J Comput Phys 2010;229:6715–31. [15] Xu K, Huang JC. A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 2010;229:7747–64. [16] Tang HZ, Xu K. A high-order gas-kinetic method for multidimensional ideal magnetohydrodynamics. J Comput Phys 2000;165:69–88. [17] Kim C, Jameson A. A robust and accurate LED-BGK solver on unstructured adaptive meshes. J Comput Phys 1998;143:598–627. [18] Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev 1954;94:511–25. [19] Blazek J. Computational fluid dynamics: principles and applications. Elsevier; 2001. [20] Li QB. Numerical study of compressible mixing layer with BGK scheme. Ph.D. thesis. Tsinghua University; 2002. [21] Li QB, Fu S, Xu K. A compressible Navier–Stokes flow solver with scalar transport. J Comput Phys 2005;204:692–714.
[22] May G, Jameson A. Improved gaskinetic multigrid method for threedimensional computation of viscous flow. AIAA Paper 2005-5106. [23] Baldwin BS, Lomax H. Thin-layer approximation and algebraic model for separated turbulent flows. AIAA 78-257. [24] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. AIAA 92-0439. [25] Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32:1598–605. [26] Chae D, Kim C, Rho OH. Development of an improved gas-kinetic BGK scheme for inviscid and viscous flows. J Comput Phys 2000;158:1–27. [27] Li QB, Fu S. Applications of implicit BGK scheme in near-continuum flow. Int J Comput Fluid Dyn 2006;20:453–61. [28] Yoon S, Jameson A. An LU-SSOR scheme for the Euler and Navier–Stokes equations. NASA CR-179556, AIAA-87-0600. [29] Yoon S, Jameson A. A multigrid LU-SSOR scheme for approximate Newton iteration applied to the Euler equations. NASA CR-179524; 1986. [30] Yang AM, Weng PF, Qiao ZD. Euler solutions of transonic flow for a helicopter rotor in hover using multigrid method. Acta Aerodyn Sin 2004;22:313–8. [31] Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes. AIAA 19811259. [32] Schmitt V, Charpin F. Pressure distributions on the ONERA M6 wing at transonic Mach numbers. AGARD AR-138; 1979. [33] CFL3D Version 6 home page.
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Implicit gas-kinetic BGK scheme with multigrid for 3D stationary transonic high-Reynolds number flows Jin Jiang, Yuehong Qian ⇑ Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, PR China
a r t i c l e
i n f o
Article history: Received 1 August 2011 Received in revised form 4 March 2012 Accepted 29 April 2012 Available online 16 May 2012 Keywords: Gas-kinetic BGK scheme Computational efficiency Acceleration techniques 3D steady transonic viscous flows
a b s t r a c t Instead of solving the Euler and Navier–Stokes equations directly, the gas-kinetic BGK scheme based on the Boltzmann equation has been developed and attracted more and more attentions since the early 1990s. It shows high accuracy and robustness for a wide range of flow regimes. But an obvious disadvantage of the BGK scheme is the low computational efficiency, in particular for multidimensional problems. Till now it has not been widely used as a practical tool for science and engineering applications. To overcome this drawback, in this paper some acceleration techniques, including local time stepping, implicit LU-SGS method, and multigrid strategy, are implemented into the original BGK approach and the new scheme is applied to study 3D steady transonic viscous flows. Numerical results show the significant speed up of the scheme to capture the steady state solution. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Transonic aerodynamics plays an important role in the operation of long range aircrafts [1]. The cruise speed of most civil airplanes is in the transonic regime. The physical natures associated with transonic flow, including shock wave, discontinuity, transition from laminar to turbulent flow, etc., are still not fully understood. Mathematically, the governing equations are nonlinear and hybrid, which show elliptic feature in the subsonic region and hyperbolic property in the supersonic parts. Thus it is hard to find analytical solutions to predict the flow patterns and then to guide the airplane designing. Following the pioneering work of Murman and Cole [2] in the early 1970s, many numerical methods for transonic aerodynamics have been developed, from the small disturbance equation to the full potential equation, the Euler equations, and the Navier–Stokes equations. Some brief reviews can be found in [3,4]. With the ever-increasing high performance of super-computer and sophisticated software, Computational Fluid Dynamics (CFD) has been considered as a convenient and useful tool for aerodynamic analysis and designing. However, to simulate high-Reynolds number flows or around complex configurations flows is still a challenge [5]. The main difficulty is not the vast need of computing resources, but the physical modeling and numerical approaches. Gas-kinetic BGK scheme (GKS) [6,7] has been developed for compressible flows since the early 1990s. Since then, many ⇑ Corresponding author. Tel.: +86 021 56331974 (O), mobile: +86 13636443941. E-mail address: [email protected] (Y. Qian). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.04.029
progress has been obtained on the validation, analysis and improvement of the GKS during the last two decades. There are numerous literatures on this topic. Xu [8] compared the gas-kinetic scheme with the Godunov method on the accuracy, efficiency, robustness, and claimed that the BGK method is more physical than the latter. Ohwada [9] explained how to construct the kinetic scheme for evolutionary equations and investigated the theoretical background and numerical errors of the scheme. Xu et al. [10] simulated hypersonic viscous flows by a multidimensional BGK flow solver with implicit LU-SGS method. Su et al. [11] and Xu et al. [12] extended the GKS to low Mach number flows. May et al. [13] improved the efficiency and convergence of the BGK scheme and presented a simplified scheme to reduce the CPU time. Issues with respect to the construction of high-order scheme [14], the extension to rarefied flows [15] and MHD [16] have also been published. Based on its constructive modeling, the GKS has a sound physical basis and can be employed as a robust and reliable flow solver to simulate problems from low Mach number flows to hypersonic flows. Nevertheless the GKS has an obvious drawback in terms of its weak computational efficiency and expensive computational cost, which limits its practical applications. Kim et al. [17] stated that in comparison with the Riemann solver, the ratio of CPU time between convective upwind and split pressure (CUSP) scheme and standard BGK scheme per grid point per step is CUSP:BGK = 1:1.908. To overcome this disadvantage, in this paper we will present an implicit multigrid BGK scheme and extend it to solve 3D high-Reynolds number flow problems.
22
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
and can be determined through the relationships
2. Numerical method
(
2.1. Gas-kinetic BGK solver
Ql ¼ r
Q ¼ After introducing BGK [18] model into 3D Boltzmann equation, in the absence of external force term, we get:
@f @f @f @f @f @f f g þu þv þw ¼ þ ub ¼ ; @t @x @y @z @t @xb s
ð1Þ
where f is the gas distribution function and g is the equilibrium state. s is the relaxation time. Both f and g are functions of time t, particle velocities ub = (u, v, w) (b = 1, 2, 3), space xb = (x, y, z) and internal variable n. The equilibrium state g is Maxwellian
Kþ3 k 2
g¼q
p
e
2
2
2
kððuUÞ þðv VÞ þðwWÞ þn
2
Þ;
q
1
C B B qU C Z C B C Q ¼B B qV C ¼ f wa dN; C B @ qW A qE 0
Fq
a ¼ 1; 2; 3; 4; 5;
ð2Þ
1 ð3Þ
b
ð4Þ
To simulate multidimensional problems, the directionalsplitting method has been adopted. Here, the x-directional equation, for example, is taken:
@f @f f g : þu ¼ @t @x s
ð5Þ
Following Xu’s [7] suggestion, the initial gas distribution f0 at control volume (CV) interface xi+1/2 can be expressed as
f0 ¼
8 l < @Q l ¼ Q iþ1=2 Q i ¼ R g l al w dN; a @x Dxl : @Q r Q iþ1 Q riþ1=2 R r r ¼ ¼ g a wa dN; @x Dxr
ð9Þ
where Dx is the distance from the CV center to the interface. To avoid the numerical oscillation close to the interface, the conservation flow variables Q liþ1=2 and Q riþ1=2 are reconstructed by MUSCL interpolation [19]:
(
Q liþ1=2 ¼ Q i þ 12 uðQ iþ1 Q i ; Q i Q i1 Þ;
where u is a limiter function uðy1 ; y2 Þ ¼
ð10Þ
y1 ðy22 þÞþy2 ðy21 þÞ y21 þy22 þ2
. The param-
eter is a small constant with the order being proportional to the local grid scale. On the fact that the non-equilibrium terms of Eq. (6) have no contributions to the conservative variables, we know:
(R R
g l ðal u þ Al Þwa dN ¼ 0;
ð11Þ
g r ðar u þ Ar Þwa dN ¼ 0:
By solving Eq. (11) the temporal slope A will be obtained. The equilibrium state g is assumed to have two slopes too [7]
r ðx xiþ1=2 Þ þ At; þ H½x xiþ1=2 a
H½x ¼
ðg f Þwa dN ¼ 0:
(
ð8Þ
g r wa dN;
g l ½1 þ al ðx xiþ1=2 Þ sðal u þ Al Þ;
x < xiþ1=2 ;
g r ½1 þ ar ðx xiþ1=2 Þ sðar u þ Ar Þ; x > xiþ1=2 ;
u2 þ v 2 þ w2 þ n2 ¼ aa wa ; 2
1; x P 0; 0; x < 0:
ð13Þ
has the similar form of a: a ¼a a wa and can be determined by the a relation of
8 0 < Q iþ1=2 Q i ¼ R g a l 0 wa dN; Dxl 0 : Q iþ1 Q iþ1=2 R r wa dN; ¼ g0 a Dxr
ð14Þ
where Q 0iþ1=2 is the vector of the initial conservative values at CV interface xi+1/2, which can be calculated by
Q 0iþ1=2 ¼
Z
Z Z
g 0 wa dN ¼
u>0
g l wa dN þ
Z Z
u<0
g r wa dN:
ð15Þ
So far, the only unknown coefficient is A. The general solution f of BGK Eq. (1) can be expressed as:
1 f xiþ1=2 ; t; ub ; n ¼
s
Z
t
0 0 g xiþ1=2 uðt t 0 Þ; t 0 ; ub ; n eðtt Þ=s dt
0
þ et=s f0 ðxiþ1=2 ut; 0; ub ; nÞ;
ð16Þ
Substituting Eqs. (6) and (12) into the above formula, after integrating, we get the gas distribution function f:
f ðxiþ1=2 ; t; ub ; nÞ ¼ ð1 et=s Þg 0 þ sðt=s 1 þ et=s ÞAg 0 l H½u þ a r ð1 þ ðsð1 þ et=s Þ þ tet=s Þða
ð6Þ
here the indexes for the other two directions (j, k) are omitted. The superscripts 0 l0 and 0 r0 denote the left- and right-side of the interface xi+1/2. a is the spatial slope of f, which has the form as:
a ¼ a1 þ a2 u þ a3 v þ a4 w þ a5
ð12Þ
where H[x] is the Heaviside function
where dN = du dv dw dn with dn = dn1dn2 dnK, Twa are the collision invariants wa ¼ 1; u; v ; w; 12 ðu2 þ v 2 þ w2 þ n2 Þ , and E is the total energy per unit mass. According to the conservation laws, the right hand side of Eq. (1) at any point in space and time should satisfy the compatibility conditions:
Z
g l wa dN;
l ðx xiþ1=2 Þ g ¼ g 0 ½1 þ ð1 H½x xiþ1=2 Þa
BF C B qU C Z C B C Fb ¼ B B F qV C ¼ ub f wa dN; C B @ F qW A F qE
R
Q riþ1=2 ¼ Q iþ1 12 uðQ iþ1 Q i ; Q iþ2 Q iþ1 Þ;
where q is the density, k = 1/2RT with the gas constant R and the temperature T. K denotes the number of internal freedom of n and equals to (5 3c)/(c 1), here c is the ratio of specific heats. (U, V, W) are the macroscopic velocities along x-, y- and z-directions. And n2 = n1n1 + n2n2 + + nKnK. Details on derivation from Eq. (1) to the Euler and Navier– Stokes equations can be found in [8]. Once the gas distribution is known, the macroscopic flow variables and their fluxes (including convection term and viscous term) can be obtained by the following integrals:
0
R
ð7Þ
H½uÞÞug 0 set=s ðAl H½ug l þ Ar ð1 H½uÞg r Þ þ et=s ðð1 uðt þ sÞal ÞH½ug l þ ð1 uðt þ sÞar Þð1 H½uÞg r Þ:
ð17Þ
Then substituting Eqs. (12) and (17) into the compatibility conditions (4), a simplified algorithm for A is given as below [20,21]
23
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
Z
Ag 0 wa dN ¼
Z
ðal g l H½u þ ar g r ð1 H½uÞÞuwa dN:
ð18Þ
Till now, all parameters are obtained. Taking the moments of BGK equation
@Q @Fb þ ¼ 0; @t @xb
ð19Þ
where,
4 c l lt DS2i t1 ¼ ðjV ni j þ cs ÞDSi þ C max ; ; þ 3q q Pr Prt X 4 c l lt DS2j ; ; þ t2 ¼ ðjV nj j þ cs ÞDSj þ C max 3q q Pr Prt X 2 4 c l l t D Sk ; ; þ t3 ¼ ðjV nk j þ cs ÞDSk þ C max 3q q Pr Prt X
ð25Þ
we could update the macroscopic flow variables by explicit or implicit time marching method. The Prandtl number Pr of the BGK scheme has a unit value and does not match the physical observation. There are many attempts have been proposed to modify it. In this paper, Xu’s treatment [7] is employed.
where CFL is the Courant–Friedrichs–Lewy number, X denotes the volume of CV (i, j, k), (DSi, DSj, DSk) stand for areas along i-, j- and k-directions, (ni, nj, nk) are unit normal vectors, V = (U, V, W) defines the components of macroscopic velocity, cs is the local sound speed, Prt is the turbulent Prandtl number, and C is equal to 4 for the Navier–Stokes equations while 0 for the Euler equations.
2.2. Relaxation time
3.2. Implicit LU-SGS method
Based on the Chapman–Enskog expansion of the BGK equation, the connection among the relaxation time s, the viscosity coefficient l, and the pressure p is obtained:
Implicit scheme is utilized in order to use large CFL number. In [10,26,27], Xu et al. argued the advantages of available implicit time integration method. Here we follow the way in traditional CFD and show how to implement LU-SGS method [28] briefly as below. Defining the residual as
l ¼ sp:
ð20Þ
Thus the relaxation time is calculated as
s ¼ l=p:
ð21Þ
RðQ Þ ¼
Z X
@Fb dX ¼ @xb
Z
ðr HÞdX ¼
X
I
H dS;
ð26Þ
@X
In order to keep the stability of current scheme, numerical dissipation is added through
then we get the semi-discrete finite volume scheme of Eq. (19)
l p pr s ¼ þ Dt l r ; p p þp
X
l
ð22Þ
where Dt is the time step for time advancing method. 2.3. Turbulence model As the increasing of Reynolds number, turbulent effect must be taken into account. Generally, there are two ways to solve this problem [22]: one is to develop one-equation or equations to describe the turbulence field based on the kinetic theory; while another one is to incorporate engineering turbulence model into the BGK method directly. Here we follow the second approach. Three different models, including Baldwin–Lomax (BL) model [23], Spalart–Allmaras (SA) model [24] and Menter’s SST (SST) model [25], have been implemented and tested. The effect of turbulent viscosity lt is included in the determination of the relaxation time of the BGK scheme:
s¼
l þ lt p
l p pr : þ Dt l r p þp
ð23Þ
In multigrid formulation, the SA model or SST model will be solved only on the finest grid, and then the turbulent viscosity is transferred to the coarser grids, where it is frozen. Furthermore, the BL model will be calculated on all grid levels. 3. Convergence acceleration techniques for steady flows
dQ þ RðQ Þ ¼ 0; dt
where H is the vector of time-averaged numerical fluxes on x-, yH and z-directions. In discrete space, the integral @ X H dS equals the sum of fluxes crossing the bounded boundaries of CV. Integrating implicitly and linearizing the residual, we obtain:
X
n dQ nþ1 @RðQ Þ DQ ¼ RðQ n Þ; þ @Q dt
8 @ HdS > > > A ¼ ; > @Q > > i < @ HdS B ¼ @Q ; > j > > > > > @ HdS :C ¼ ; @Q k
ðL þ D þ UÞDQ ¼ RðQ n Þ;
8 þ þ þ > L ¼ A þ B þ C ; > i1 j1 k1 < þ þ þ X > D ¼ Dt þ Ai Ai þ Bj Bj þ Ck Ck ; > : U ¼ Aiþ1 þ Bjþ1 þ Ckþ1 ; 8 1 > < A ¼ 2 ðA rA IÞ; r A P rA ; B ¼ 12 ðB r B IÞ; r B P rB ; > : 1 C ¼ 2 ðC r C IÞ; r C P rC ;
;
ð24Þ
ð30Þ
with
For stationary problems, the possible largest time step for each CV can be used to accelerate the solution. Here we employ the maximum time step for each CV (i, j, k) as
X
ð29Þ
expanding the second term in the bracket on the left hand side of Eq. (28), and splitting the Jacobian using upwind-bias method, the final discretization form is gotten
and
t1 þ t2 þ t3
ð28Þ
where DQ = Qn+1 Qn, superscript 0 n0 denotes the nth time level. After introducing the Jacobian matrices,
3.1. Local time stepping
Dt ¼ CFL
ð27Þ
ð31Þ
ð32Þ
where I is the unit matrix, ðrA ; rB ; rC Þ are the spectral radii of the Jacobian matrices.
24
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28 Table 1 Comparisons on CPU time, work units and factor of speedup of different numerical approaches for ONERA M6. The left column under 0 CPU time0 is the total cost and the right one is the time per step. 0 V20 and 0 V30 stand for 2-level and 3-level V-cycle multigrid. SST model is employed.
Y X Z
Numerical scheme
CPU time (s)
LU-SGS JST LU-SGS JST/V3 Explicit GKS LU-SGS GKS LU-SGS GKS/V2 LU-SGS GKS/V3
43,418 11,357 93,874 41,258 18,633 12,115
6.76 9.26 6.04 6.94 11.79 11.99
Work units
Speedup
6420 1226.875 15,535 5945 1581 1010.3125
1 5.23 – 1 3.76 5.88
2. Transfer the flow variables to coarse grid by volume-weighted method
Q 02h ¼ Fig. 1. Computational mesh for ONERA M6 wing in detail. The red grid is the symmetry plane of k = 1 and the green one is the wing surface grid. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
X
Rh ðQ h Þ R2h Q 02h ;
Rnew 2h ðQ 2h Þ ¼ R2h ðQ 2h Þ þ P 2h :
Explicit/single grid LU-SGS/single grid LU-SGS/2-level multigrid LU-SGS/3-level multigrid
-1
-2
0 dQ Nh ¼ Q m Nh Q Nh :
-3
ð39Þ
Then, get the new solution on the finest grid by interpolating the coarse-grid correction successively
Q new ðN1Þh ¼ Q ðN1Þh þ I Nh
ðN1Þh
0
4000
8000
12000
16000
Fig. 2. Convergence history of continuity equation in GKS framework by different time marching methods with single grid or multigrid for transonic viscous flow over ONERA M6 wing: C–O grid; CFL = 0.8 for explicit scheme and CFL = 5 for LU-SGS scheme; SST model.
Then, factorize the left hand side of Eq. (30) approximately into
ðL þ DÞD1 ðD þ UÞDQ ¼ RðQ n Þ;
ð33Þ
and use two-step sweeping way to get the solution DQ:
ðL þ DÞDQ ¼ RðQ n Þ; ðD þ UÞDQ ¼ DDQ :
¼ Q n þ DQ :
ð40Þ
where INh is the prolongation operator. 5. Repeat the work done from step 1 to step 4 until the solution is convergent. 4. Numerical cases To validate and analyze the accuracy and robustness of the current implicit gas-kinetic BGK approach with multigrid, two test cases have been carried out. 4.1. ONERA M6 wing
ð34Þ
Sequently, the macroscopic flow variables are updated by
Q
dQ Nh ;
ðN1Þh
Work Units
nþ1
ð38Þ
Finally, calculate the correction and update the solution on the coarse grid in the same way as the fine grid. 3. Do what have been done in step 2 on the coarser grid. 4. Obtain the coarse-grid correction after m-iteration (in this paper, m is set to 2N1)
-4
ð37Þ
and get the residual on the coarse grid
0
-5
ð36Þ
then collect the residual over the eight CVs (for 3D Cartesian grid) on the fine grid for the coarse grid. Define a forcing function as
P 2h ¼
Log10 (Residual)
.X X ðQ h Xh Þ Xh ;
ð35Þ
3.3. Multigrid Multigrid has been proved as a powerful acceleration tool for stationary flows. The following gives a brief procedure of V-cycle multigrid [29,30], where the subscript 0 Nh0 denotes the Nhth level grid. 1. Calculate the approximate solution Qh and residual Rh(Qh) on the fine grid.
The flow conditions are Ma = 0.8395, angle of attack (AOA) of 3.06° and mean-chord based Reynolds number of Re = 11.72 106. Computational grid is C–O type of i j k = 193 49 33 (seen in Fig. 1). Firstly, we compare the convergence rate between implicit LU-SGS method and explicit time stepping technique and check the effect of multigrid formulation. Fig. 2 illustrates the time history of residual of continuity equation. It shows that (1) the implicit LU-SGS method is faster to reach the convergence state than explicit one; (2) convergence is accelerated significantly by applying multigrid to current BGK scheme. Table 1 lists the comparisons on CPU cost, work units and speedup of various approaches. The needed CPU time per step of implicit GKS is nearly the same of that of implicit JST scheme [31]. When compared with the single grid calculation, even more time is spent for each step, the total CPU
25
1
1
0.5
0.5
0.5
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.20
-0.5
-1
0
0.2
0.4
0.6
0.8
-Cp
1
-Cp
-Cp
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.44
-0.5
-1
1
0
0.2
0.4
0.6
0.8
-0.5
-1
1
1
0.5
0.5
0.5
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.80
0.2
0.4
0.6
0.8
-Cp
1
-1 0
0
0.2
0.4
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.90
-0.5
-1 0
1
0.2
0.4
x/c
0.6
0.8
1
x/c
1
-0.5
Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.65
x/c
-Cp
-Cp
x/c
0
0.6
0.8
0 Exp. GKS/V3, BL model GKS/V3, SA model GKS/V3, SST model JST, SST model Z/b = 0.95
-0.5
-1
1
0
0.2
0.4
x/c
0.6
0.8
1
x/c
Fig. 3. Comparison of pressure distributions at selected span-wise locations of ONERA M6 from implicit GKS and implicit JST scheme with various turbulence models and experimental data [32]. c is the chord length, Z is the distance to the symmetry plane and b is the semi-span. V3: 3-level V-cycle multigrid.
higher resolution, better capability on shock capturing than JST scheme and (2) both SA model and SST model are more suitable for this case than BL model. We also present some primitive results with respect to grid independence. Three meshes are tested: one is that described above and marked as CO-I; the second is a finer C–O type grid CO-II of dimensions i j k = 289 64 49; and the last one is C–H type of i j k = 161 51 55 with 31 grids imposed on the wing surface along span direction. Fig. 4 shows the pressure distributions at some selected span-wise locations. Here BL model is used. With different grids the solutions are quite the same. The structure of two branches of the shock at 80% span location is verified clearly only by GKS on mesh C–H. This illustrates the solution dependence on the computational grid. On C–O mesh k-planes (x–y
1
1
0.5
0.5
0
Exp. GKS, mesh C-H JST, mesh C-H GKS, mesh CO-I JST, mesh CO-I GKS, mesh CO-II JST, mesh CO-II Z/b = 0.65
-0.5
-1
0
0.2
0.4
0.6
x/c
0.8
1
-Cp
-Cp
cost of multigrid method is much lower. According to this case, the more levels used the higher speedup is obtained. With 3-level V-cycle multigrid, the factor of speedup is up to about 6 for GKS, in contrast to above 5 for JST scheme. Secondly, we focus on the influence of turbulence models and analyze the resolution of shocks. Fig. 3 plots the pressure distributions at selected span-wise locations. All solutions are obtained with CFL = 5. The overall numerical results are in good agreement with the experimental data [32]. But we can see GKS captures the shocks more sharply than JST scheme. Beside 20% span location, the positions of shock wave given by GKS are closer to the experimental data [32] than results of JST scheme. In the GKS framework, SA model and SST model can predict more accurate shock positions than BL model. After all we can say (1) GKS shows
0
Exp. GKS, mesh C-H JST, mesh C-H GKS, mesh CO-I JST, mesh CO-I GKS, mesh CO-II JST, mesh CO-II Z/b = 0.80
-0.5
-1
0
0.2
0.4
0.6
0.8
1
x/c
Fig. 4. Comparison of pressure distributions at selected span-wise locations of ONERA M6 from implicit GKS and implicit JST scheme on mesh CO-I, mesh CO-II and mesh C–H with experimental data [32]. c is the chord length, Z is the distance to the symmetry plane and b is the semi-span.
26
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28 Y
Y X
X
Z
Z
Fig. 5. Schematic surface grid for ARA M100 wing/body configuration and the symmetry plane (red). The body grid (green) covers a region of i j k = (33–289) 1 (1–17), while the wing surface grid (blue) of i j k = (65–257) 1 (17–49). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
schematic surface grid and symmetry plane are illustrated in Fig. 5, which is C–O type of dimensions of i j k = 321 57 49. Flow field is solved by GKS with BL model and SST model. Comparison of pressure distributions on wing surface at various span-wise locations from GKS, CFL3D [33] and experimental data [34] is plotted in Fig. 6. The turbulence model is the key point for this case. Current GKS approach coupled with SST turbulence model gives solutions in good match to experimental data [34]. The shocks gotten with BL model situate too far aft than results with the other two turbulence models, which can be confirmed from the flow patterns on the upper wing surface (see Fig. 7). When SST model was employed, the separated streamline is located
planes) are not totally parallel to the free stream and thus the solutions selected at these special wing sections are not the precise ones we wanted to observe. By projecting or interpolating the neighboring data, the results may be improved and become better. Similarly, the good shock capturing capability of GKS is validated again. 4.2. ARA M100 wing/body Flow around a complex configuration of ARA M100 wing/body at free stream condition Ma = 0.8027, AOA = 2.873° and chord based Reynolds number of Re = 13.1 106 is simulated. The
1
1
1
0.5
0.5
0.5
-Cp
1.5
-Cp
1.5
-Cp
1.5
0
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5 Z/b = 0.123
-1
0
0.2
0.4
0.6
0.8
Z/b = 0.231
-1
1
0
0.2
0.4
x/c
0.6
0.8
Z/b = 0.325
-1
1
1
0.5
0.5
0.5
0 Exp. GKS, BL model GKS, SST model CFL3D, SA model
0.6
x/c
0.8
1
0.8
1
Exp. GKS, BL model GKS, SST model CFL3D, SA model
Exp.
GKS, BL model GKS, SST model CFL3D, SA model
-0.5 Z/b = 0.455
0.4
0.6
0
0
0.2
0.4
-Cp
1
-Cp
1
-Cp
1.5
0
0.2
x/c
1.5
-1
0
x/c
1.5
-0.5
Exp. GKS, BL model GKS, SST model CFL3D, SA model
-0.5
-0.5 Z/b = 0.633
-1
0
0.2
0.4
0.6
x/c
0.8
1
Z/b = 0.817
-1
0
0.2
0.4
0.6
0.8
1
x/c
Fig. 6. Comparison of pressure distributions at selected span-wise locations of ARA M100 wing/body from GKS, CFL3D [33] and experimental data [34]. c is the chord length, Z is the distance to the symmetry plane and b is the wing span.
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
27
Fig. 7. Extreme streamlines on the upper wing surface of ARA M100 wing/body in the context of pressure contours. Red dash-dot-dot line is a reference of i = 209. Pink solid lines identify the various stations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 2 Comparison of lift CL and drag CD for ARA M100 wing/body among different turbulence models.
CL CD
BL model/GKS
SST model/GKS
SA model/CFL3D [33]
0.7358 0.04920
0.6768 0.04430
0.6632 0.04855
upstream while the reattached line is further downstream, that forms a wider reverse flow regime than that with BL model and leads to great variances in lift and drag (listed in Table 2). Those differences may be caused by the distinct natures of the turbulence models: (1) BL model is an algebraic two-layer model developed from boundary layer theory. Therefore it cannot trace the time evolution of the turbulent viscosity, and is only valid for attached flows or small separated flows. (2) SST model is developed on the transportation theory of the turbulence field. It is suitable for flows involving adverse pressure gradients and has ability to predict pressure-induced separation [25]. The SST model is the preferable one for future applications.
1. multigrid optimization, such as strategy (V-cycle or W-cycle), number of iterations, restriction/prolongation operators, and so on, 2. comparison and option of turbulence models, 3. effectiveness of the acceleration techniques for unsteady problems.
Acknowledgments The present work is supported partly by Program for Changjiang Scholars and Innovative Research Team in University, China (No. IRT0844) and Shanghai Science and Technique Committee (No. 11XD1402300). The work of Jiang is supported by Innovative Foundation for Graduate Students of Shanghai University too. The authors are grateful to Dr. Aiming Yang at Fudan University for his fruitful discussion and helpful suggestions on the implementation of multigrid. We appreciate Prof. Kun Xu, at HongKong University of Science and Technology, for his kindness and communication on BGK scheme. We thank Dr. Jianping Meng for discussion. In addition, the authors would also like to thank the reviewers for their constructive comments and to refine the current manuscript carefully.
5. Conclusions In this paper, we develop and extend the GKS into implicit multigrid formulation for high-Reynolds number steady flows. Some convergence acceleration techniques, including local time stepping, implicit LU-SGS method, and multigrid formulation, are proposed to the original scheme presented by Prendergast and Xu [6]. It is encouraged to improve the computational efficiency by applying those techniques into GKS for its significant acceleration to achieve convergence steady state and dramatic reduction of CPU time. At the same time, we incorporate three turbulence models (BL model, SA model and SST model) into current GKS code and check the validity of those models in combination with GKS. Based on the results of the two 3D test cases, SST model is identified as a best choice for further investigation. Additionally, there are still some problems need to be solved in the future:
References [1] Jameson A. Computational transonics. Commun Pure Appl Math 1988;41: 507–49. [2] Murman EM, Cole JD. Calculation of plane steady transonic flows. AIAA J 1971;9:114–21. [3] Jameson A, potential Full. Euler and Navier–Stokes schemes in applied computational aerodynamics. Prog Astronaut Aeronaut 1990;125:39–88. [4] Caughey DA, Jameson A. Development of computational techniques for transonic flows: an historical perspective; 2002. [5] Jameson A. Essential elements of computational algorithms for aerodynamic analysis and design. NASA/CR-97-206268. ICASE report no. 97-68. [6] Prendergast KH, Xu K. Numerical hydrodynamics from gas-kinetic theory. J Comput Phys 1993;109:53–66. [7] Xu K. A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method. J Comput Phys 2001;171:289–335. [8] Xu K. Gas-kinetic schemes for unsteady compressible flow simulations. VKI report 1998-03. [9] Ohwada T. On the construction of kinetic schemes. J Comput Phys 2002;177:156–75.
28
J. Jiang, Y. Qian / Computers & Fluids 66 (2012) 21–28
[10] Xu K, Mao ML, Tang L. A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow. J Comput Phys 2005;203:405–21. [11] Su MD, Xu K, Ghidaoui MS. Low-speed flow simulation by the gas-kinetic scheme. J Comput Phys 1999;150:17–39. [12] Xu K, He XY. Lattice Boltzmann method and gas-kinetic BGK scheme in the low Mach number viscous flow simulations. J Comput Phys 2003;190:100–17. [13] May G, Srinivasan B, Jameson A. An improved gas-kinetic BGK finite-volume method for three-dimensional transonic flow. J Comput Phys 2007;220: 856–78. [14] Li QB, Xu K, Fu S. A high-order gas-kinetic Navier–Stokes flow solver. J Comput Phys 2010;229:6715–31. [15] Xu K, Huang JC. A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 2010;229:7747–64. [16] Tang HZ, Xu K. A high-order gas-kinetic method for multidimensional ideal magnetohydrodynamics. J Comput Phys 2000;165:69–88. [17] Kim C, Jameson A. A robust and accurate LED-BGK solver on unstructured adaptive meshes. J Comput Phys 1998;143:598–627. [18] Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev 1954;94:511–25. [19] Blazek J. Computational fluid dynamics: principles and applications. Elsevier; 2001. [20] Li QB. Numerical study of compressible mixing layer with BGK scheme. Ph.D. thesis. Tsinghua University; 2002. [21] Li QB, Fu S, Xu K. A compressible Navier–Stokes flow solver with scalar transport. J Comput Phys 2005;204:692–714.
[22] May G, Jameson A. Improved gaskinetic multigrid method for threedimensional computation of viscous flow. AIAA Paper 2005-5106. [23] Baldwin BS, Lomax H. Thin-layer approximation and algebraic model for separated turbulent flows. AIAA 78-257. [24] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. AIAA 92-0439. [25] Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32:1598–605. [26] Chae D, Kim C, Rho OH. Development of an improved gas-kinetic BGK scheme for inviscid and viscous flows. J Comput Phys 2000;158:1–27. [27] Li QB, Fu S. Applications of implicit BGK scheme in near-continuum flow. Int J Comput Fluid Dyn 2006;20:453–61. [28] Yoon S, Jameson A. An LU-SSOR scheme for the Euler and Navier–Stokes equations. NASA CR-179556, AIAA-87-0600. [29] Yoon S, Jameson A. A multigrid LU-SSOR scheme for approximate Newton iteration applied to the Euler equations. NASA CR-179524; 1986. [30] Yang AM, Weng PF, Qiao ZD. Euler solutions of transonic flow for a helicopter rotor in hover using multigrid method. Acta Aerodyn Sin 2004;22:313–8. [31] Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes. AIAA 19811259. [32] Schmitt V, Charpin F. Pressure distributions on the ONERA M6 wing at transonic Mach numbers. AGARD AR-138; 1979. [33] CFL3D Version 6 home page.