Computation of complex turbulent flow using matrix-free implicit dual time-stepping scheme and LRN turbulence model on unstructured grids

Computers & Fluids 33 (2004) 119–136 www.elsevier.com/locate/compfluid

Computation of complex turbulent flow using matrix-free implicit dual time-stepping scheme and LRN turbulence model on unstructured grids Yong Zhao

*

School of Mechanical and Production Engineering, Nanyang Technological University, Nanyang Avenue 639798, Singapore Received 28 January 2002; received in revised form 11 September 2002; accepted 23 January 2003

Abstract In this study, a matrix-free implicit dual time-stepping method has been developed. It is implemented, together with a low-Reynolds-number q–x turbulence model, in a high-order upwind finite-volume solver on unstructured grids. Semi-implicit treatment of the source terms of the q and x equations is also introduced to further stabilize the numerical solution. It has been found that these techniques provide strong stabilization in the computation of a supersonic flow with complex shock–boundary-layer interactions in a channel with a backward-facing step. The proposed method has a low-memory overhead, similar to an explicit scheme, while it shows good stability and computational efficiency as an implicit scheme. The method developed has been validated by comparing the computed results with the corresponding experimental measurements and other calculated results, which shows good agreement. Research is being done to extend the method to calculate unsteady turbulent flows. Ó 2003 Elsevier Ltd. All rights reserved. Keywords: Matrix-free method; Implicit scheme; Turbulence modelling; Dual time-stepping scheme; Unstructured grids

1. Introduction Transport type two-equation turbulence models, especially the low-Reynolds-number (LRN) versions, are suitable candidates for tackling complex high-speed turbulent flows with massive flow separations because they can provide the history effects in turbulence and they do not need

*

Tel.: +65-790-4545; fax: +65-791-1859. E-mail address: [email protected] (Y. Zhao).

0045-7930/04/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0045-7930(03)00031-8

120

Y. Zhao / Computers & Fluids 33 (2004) 119–136

ad hoc near wall treatment. However, it has been found that their use often leads to unstable solution, or even incorrect solution, in both the near wall and far field regions. This is especially true in the near wall regions. There are sharp gradients of turbulence quantities in these regions which make the k–e equations difficult to solve due to their very stiff source terms. This difficulty is further complicated by the use of highly stretched grids near the solid wall. Gerolymos [1] introduced limits to the values of k and e in order to ensure positivity of these values. Special correlations were also used to specify the initial distribution of the two quantities. Others [2,3] imposed a growth rate for the turbulence quantities and used the solution with the Baldwin– Lomax model [4] as an initial condition for the k–e computation. However, it was found by the above authors that although such treatment might work for certain types of flows, the convergence rate is still very low. In addition, they are not generally applicable and sometimes these methods lead to completely incorrect solution. Mavriplis and Martinelli [3] used a point-implicit method with source term linearization to solve the transformed k and e equations with a multigrid scheme. The transformed k equation is obtained by dividing the k by a damping function which will vanish at the wall in order to avoid the singularity problem there. Results obtained show that even with the use of an uniform initial condition, rapid convergence rate can be obtained. On the other hand, according to Vandromme [5], the inclusion of the positive source terms in an implicit scheme is unnecessary and even destabilizing. Zhao and Ding [6] adopted a strategy of using a k–e model in the high-Reynolds region while using an one-equation model in the near wall region, in order to achieve a stable solution in the computation of 3D channel flow with shock/boundarylayer interaction. The source terms were treated by a point-implicit method in both regions and the time step size was also limited to ensure positivity of turbulence quantities. A very sophisticated method was developed to ensure the stability of the turbulence quantities based on asymptotic development of these quantities in a uniform turbulence field [7]. Similar approach was also used by Mohammadi [8]. The q–x turbulence model was first proposed by Coakley [9] to overcome the stiffness problems associated with other LRN turbulence models mentioned above. The model was successfully used in an implicit LU-SGS multigrid solver on structured grids [10]. Recently, Luo et al. [11] have developed a matrix-free implicit LU-SGS method with GMRES on unstructured grids for unsteady inviscid flow simulation. Low storage was achieved by using an approximate flux function and replacing the Jacobian-state vector product by a flux difference between two successive time levels in the SGS framework. The dual time-stepping scheme has been used for unsteady flow simulation, see [12–15] for details. A multigrid scheme can be used together with implicit physical time-stepping and explicit pseudo-time-stepping to achieve high-computational efficiency. However, it has been pointed out in [12] that the pseudo-time-step size is limited by 2/3 of the physical time step size found in a stability analysis. This will result in a significant reduction in convergence rate if the physical time step size is very small due to certain physical requirement. Thus it would be desirable to perform the pseudo-time-stepping implicitly in order to decouple the two time steps. Even for steady flows, an implicit pseudo-time-stepping method coupled with a proper matrix-free method will improve the computational efficiency without excessive memory requirements. This paper describes the development and validation of a matrix-free fully implicit dual timestepping scheme which is used in conjunction with a LRN q–x model of Coakley [9] in a highorder upwind finite-volume (FV) solver on unstructured grids.

Y. Zhao / Computers & Fluids 33 (2004) 119–136

121

In the following sections, the governing equations are introduced, and this is followed by the description of the numerical methods. And finally the results and discussion are presented and concluding remarks are then provided.

2. Governing equations The governing equations are the Reynolds averaged Navier–Stokes (RANS) equations in integral form (see [16] for detailed derivation) including the q–x transport equations [9]: Z Z Z Z o ~ ~ Fv  ~ W dV þ Fc  ~ ns dS ¼ ns dS þ St dV ; ð1Þ ot X s s X where W ¼ ðq; qu; qv; qE; qq; qxÞT ; 2 3 ~ qU 7 6 quU 6 ~ þ p~i 7 7 6 ~ þ p~ 7 6 qvU j 7; ~ Fc ¼ 6 6 qH U ~ 7 7 6 7 6 ~ 5 4 qqU ~ qxU 2

3 0 7 6 ~ sx 7 6 7 6 7 6 ~ s y 7 6 6 ~ 7 ~Q 7 6 sU  7; ~ Fv ¼ 6 7 6 l ~q 7 6 ll þ t r 7 6 rq 7 6 7 6 lt ~ 7 6 4 ll þ rx 5 rx

 S 2D St ¼ 0; 0; 0; 0; Cq Cl Dq 2   1 qxq; x 3x



 Cx1

S D Cl 2  Cx3 x x





2

 Cx2 qx

:

The terms in the above equations are as follows: qE ¼ qe þ 12qðu2 þ v2 Þ; ~ si ¼ six~i þ siy~ j; le ¼ ll þ lt ; Sij ¼

sy~ j; s ¼~ sx~i þ~

sij ¼ le Sij ; ~ ¼ cp Q



ll lt ~ þ rT ; Prl Prt

oui ouj 2 ouk þ  dij ; oxj oxi 3 oxk

D¼

ou ov þ ; ox oy

The low-Reynolds damping terms introduced in [9] are employed: Dq ¼ 1  expð0:022Req Þ;

Cx1 ¼ 0:5Dq þ 0:055;

here Req ¼ ðqqlw =lm Þ is the turbulent Reynolds number. The eddy viscosity is calculated by

122

Y. Zhao / Computers & Fluids 33 (2004) 119–136

qq2 : x The above system of equations are closed by the state equation of perfect gas: l t ¼ Cl Dq

p ¼ qðc  1Þe:

ð2Þ

The various constants in the q–x model are listed as follows: Cq 0.5

Cx2 0.833

rq 0.8

Cx3 2/3

rx 2.0

Cl 0.09

The above equations are non-dimensionalized using a reference length and free stream flow conditions. Only the non-dimensional form of the equations are used in the numerical calculations. 3. Numerical methods The above equations in integral form are discretized on unstructured triangular grids. Here a cell-vertex scheme is used, and this means that all computed variables in vector W are stored at vertices of the cells. In two-dimensional cases, for every cell-vertex P , a control volume is constructed by using the median dual of the triangular mesh as shown in Fig. 1. Integrating Eq. (1) over the control volume surrounding node P results in the following: ZZ ZZ ZZ ZZ o W dS þ r~ Fc dS ¼ r~ FV dS þ St dS ð3Þ ot Scv Scv Scv Scv The convection term is ZZ Z ~ Fc  ~ r~ Fc dS ¼ n dL; Scv

Lcv

where Scv is the area of a particular control volume, and Lcv is the control volume surface length. N3 Triangular element

Wj

2

N2 Control Volume (Shaded Part)

Cij Wi

N4

P

1

N1 N5

Fig. 1. A control volume.

Y. Zhao / Computers & Fluids 33 (2004) 119–136

123

All the first derivatives of fluid properties at the centres of cells are evaluated by using the GreenÕs theorem, the first derivatives at nodes are calculated by volume averaging those at the centres of the surrounding cells [17]: PN r/j  Scelli ; r/jnodeP ¼ i¼1PN celli i¼1 Scelli where N is the number of cells associated with the node P . Finally for a certain node p, Eqs. (3) are discretized as o Wp DScv ¼ RðWp Þ ot

ð4Þ

and RðWp Þ ¼ 

nbseg X

ð~ Fc  ~ Fv Þij  ~ nLn þ St DScv ;

ð5Þ

n¼1

where nbseg is the number of the edges associated with the node P . 3.1. Edged-based flux calculation The control volume surface enclosing a node is broken down into segments with each segment connected to an edge, whose one end is the node under consideration. The total inviscid flux entering the control volume surface surrounding a node will then be gathered by calculating the fluxes for the two ends of every edge, which can be performed in a loop over all edges. The inviscid flux at the control volume surface associated with an edge ij (shown in Fig. 1 as edge pN2 ) Fc Þij  ~ nij is evaluated based on RoeÕs approximate Riemann solver [18]: ðFc Þij ¼ ð~ ðFc Þij ¼ 12½FcL þ FcR  jAjðW R  W L Þ :

ð6Þ

The left and right state vectors W L and W R are evaluated using a nominally third-order upwind biased interpolation scheme, the MUSCL (monotonic upstream schemes for conservation laws) scheme [19] (if the left and right state vectors are set to Wi and Wj , it is a first-order upwind scheme): þ WL ¼ Wi þ 14½ð1  jÞD i þ ð1 þ jÞDi ;  WR ¼ Wj  14½ð1  jÞDþ j þ ð1 þ jÞDj :

Here Dþ i ¼ Wj  Wi ; ! ! þ D i ¼ Wi  Wi1 ¼ 2 ij  rWi  Wj  Wi ¼ 2 ij  rWi  Di ; þ D j ¼ Wj  Wi ¼ Di ;

! ! þ Dþ j ¼ Wjþ1  Wj ¼ 2 ij  rWj  Wj  Wi ¼ 2 ij  rWj  Dj :

124

Y. Zhao / Computers & Fluids 33 (2004) 119–136

Here j is set to 1=3 to obtain nominally third-order accuracy. It has been shown by Zhao and Zhang [20] that this interpolation is more accurate than the central difference scheme with fourthorder artificial viscosity in the computation of incompressible flows. If the solution exhibits strong flow gradients, a slope limiting procedure can be introduced in the process of interpolating the state vector W at L and R. In the present study, the Minmod limiter [21] is employed. 3.2. Matrix-free fully implicit dual time-stepping scheme A pseudo-time-dependent term is added to the discretized equations (4) and the physical timedependent term can be calculated by a three-point second-order implicit scheme: Wi nþ1;mþ1  Wi nþ1;m 1:5Wi nþ1;mþ1  2:0Wi n þ 0:5Wi n1 e ðWi Þ: DScv ¼ Rnþ1;mþ1 ðWi Þ  DScv ¼ R Ds Dt

ð7Þ

Here n indicates the n time level in physical time and m is the m pseudo-time level. It is clear that in the above equation both time integration in physical time and pseudo-time are implicit, thus the term fully implicit dual time-stepping scheme is adopted. The matrix-free implicit pseudo-timestepping is proposed as follows. Similar to [11], in order to obtain the Taylor expansion of the residual term in pseudo-time, the total numerical flux (including convection and viscous fluxes) associated with edge ij is approximated as: Fij 1=2½Fc ðWi Þ þ Fc ðWj Þ  jkij jðWi  Wj Þ ; here kij is the spectral radius associated with edge ij. Thus  oFij oFc ðWi Þ  jkij j ; ¼ 1=2 oWi oWi  oFij oFc ðWj Þ þ jkij j : ¼ 1=2 oWj oWj The residual term can be expanded in pseudo-time at a node i as: Wi nþ1;mþ1 ¼ Wi nþ1;m þ DWi ; X oRi oRi DWi þ DWj : oWi oWj j¼1 jm

Rnþ1;mþ1 ¼ Rnþ1;m þ i i

For the whole flow field, Rnþ1;mþ1 ¼ Rnþ1;m þ ADW ; here



oRi A¼ oWj

¼ ½Aij :

Y. Zhao / Computers & Fluids 33 (2004) 119–136

125

With the above Taylor expansion, the whole-field equivalent of Eq. (7) can be rewritten as  Dt þ 1:5Ds ADs DW 1:5W nþ1;m  2:0W n þ 0:5W n1 ¼ Rnþ1;m  DScv ; DScv  Dt DScv Ds Dt i.e. e DScv A

DW e nþ1;m ; ¼R Ds

e ¼ Dt þ 1:5Ds  ADs ; A Dt Scv e nþ1;m ¼ Rnþ1;m  1:5W R

nþ1;m

 2:0W n þ 0:5W n1 DScv : Dt

Further approximation can be introduced in order to achieve a matrix-free computation. If we employ a point implicit treatment to the above equations, only diagonal terms in Ae are used in the pseudo-time-stepping. As a result, the equation for every node can now be written as: e ii DWi ¼ R e nþ1;m : DScv A i Ds Thus DScv

DWi e 1 e nþ1;m e e nþ1;m ; ¼R ¼ A ii R i i Ds

and e nþ1;m e 1 e nþ1;m e R ¼ A ii R i ; e 1 ¼ A ii



Dt þ 1:5Ds Aii Ds I Dt Scv

1 :

Further stability enhancement is introduced by decoupling the equations when evaluating the matrix Aii . This is obtained by assuming that the inviscid flux of a particular equation is a function of the variable of that equation while other variables are locally frozen. This has been found to be effective and stable. It also makes the solution method truly matrix-free. Pseudo-time-stepping is then performed on the above equation. For a five-stage scheme, the stage coefficients are a1 ¼ 14;

a2 ¼ 16;

a3 ¼ 38;

a4 ¼ 12;

a5 ¼ 1:

To further enhance solution convergence, implicit residual smoothing is also employed. The idea behind the implicit residual smoothing is to replace the residual at one point of the flow field with a smoothed or weighted average of the residuals at the neighboring points. The averaged residuals are calculated implicitly in order to increase the maximum CFL number, thus increasing the convergence rate. Normally this procedure allows the CFL number to be increased by a factor of 2 or 3. The smoothing equation for a node k can be expressed as

126

Y. Zhao / Computers & Fluids 33 (2004) 119–136

Rk ¼ Rk þ r2 Rk ;

ð8Þ

where R is the smoothed residual, R is the original residual and e is the smoothing coefficient: ( " # ) 2 1 CFL  1 ;0 :  ¼ max 4 CFL Here CFL is the maximum CFL number of the basic scheme. The solution to the above equations can be obtained on a unstructured grid by a Gauss–Seidel iterative method: ðmÞ

ðoÞ

Rk ¼ Rk þ 

numnodðkÞ X

ðmÞ

½Ri

ðmÞ

 Rk ;

i¼1

i.e., ðmÞ Rk

P ðm1;mÞ ðoÞ Rk þ  numnodðkÞ Ri i¼1 ¼ : 1 þ   numnodðkÞ

Here the numnodðkÞ is number of neighboring nodes of node k.

4. Treatment of turbulence equations The q–x equations are numerically solved together with the mean flow equations by the same matrix-free implicit dual time scheme. A simple upwind scheme is used for the turbulence equations. To enhance stability and remove the stiffness associated with the source terms, it is also necessary to adopt a point implicit method in the calculation of these terms, which are given as follows for node i:  n oSit nþ1 n ðWitnþ1  Witn Þ: ð9Þ Sit ¼ Sit þ oWit In order to simplify the calculation and avoid inversion of matrices, the Jacobian of the source terms is diagonalized by introducing implicit treatment for the negative parts of the source terms only, which can be considered to be the major factors producing non-positivity in turbulence quantities:   n oSt nþ1 n DWt ; St St þ oWt 3 2 oStq oStq 6 oðqqÞ oðqxÞ 7 oSt 7 6 ð10Þ ¼6 7   5 4 oStx oWt oStx ¼

oðqqÞ oðqxÞ 2 ðD þ xÞCq 3 0

0

Cx1 Cx3 D þ Cx2x

:

ð11Þ

Y. Zhao / Computers & Fluids 33 (2004) 119–136

127

The implicit source terms will then be added to the residual implicit terms and they become part of the Aii matrix. And the turbulence equations are integrated together with the mean flow equations in exactly the same manners. 5. Boundary conditions At the solid wall, non-slip, no-injection and adiabatic boundary conditions are imposed. The inlet and outlet boundary conditions for the flow equations are based on the flux vector splitting technique for the inviscid fluxes as proposed by Stegger and Warming [22]. For the supersonic channel flow considered here, the flow conditions upstream of the inlet edges are given, such as non-dimensional velocity profile, non-dimensional density, pressure and temperature profiles and turbulence quantities. Pressure downstream of the outlet edges may be prescribed if the flow at the outlet edge is subsonic. If it is supersonic no conditions are specified and all the boundary conditions are extrapolated from the flow field. These conditions are used as parts of the flux vector splitting calculations.

Fig. 2. The geometry and in flow conditions.

Fig. 3. The fine mesh (no. nodes ¼ 38,962 and no. elements ¼ 76,402).

128

Y. Zhao / Computers & Fluids 33 (2004) 119–136

6. Results and discussion The proposed method has been tested with a case of turbulent supersonic flow over a backward facing step and the calculated results are compared with the corresponding experimental measurements of McDaniels et al. [23] and Eklund et al. [24]. The flow conditions and flow field 0

10

CFL=0.5 No. of Sub-iterations: 50

-1

Log10(RMS)

10

-2

10

-3

10

-4

10

400

(a)

600

800

1000

No. of Iterations 0

10

CFL = 0.5 No. of Sub-iterations: 50

-1

Log10(RMS)

10

-2

10

-3

10

-4

10

(b)

400

600

800

1000

1200

1400

No. of Iterations 0

10

CFL= 10 No. sub-iter= 15

-1

10

Log10 (RMS)

-2

10

-3

10

-4

10

-5

10

10

(c)

20 30 Number of physical iterations

40

50

Fig. 4. Convergence history: (a) explicit pseudo-time-stepping scheme with the medium grid; (b) explicit pseudo-timestepping scheme with the fine grid; (c) the implicit dual time-stepping scheme with the fine grid.

Y. Zhao / Computers & Fluids 33 (2004) 119–136 PLIIF Step4 Step2 Step1

0

0

-1

-1

-1

-1

-2

-2

-2

-2

-3

-3

-3

-3

Y/ L

0

Y/L

0

Y/L

Y/L

PLIIF Step4 Step2 Step1

PLIIF Step4 Step2 Step1

PLIIF Step4 Step2 Step1

129

-4

-4

-4

-4

-5

-5

-5

-5

-6

-6

-6

-6

0

0.5

u/Uin

1

0

(a)

0.5

1

u/Uin

0

(b)

0.5

1

u/Uin

0

(c)

0.5

1

u/Uin

(d)

Fig. 5. U-velocity profiles calculated with the three meshes (denoted by lines, step1 ¼ fine, step2 ¼ medium and step4 ¼ coarse) and their comparisons with experimental measurements (denoted by symbols) at (a) X =L ¼ 0:97; (b) X =L ¼ 2:76; (c) X =L ¼ 3:99; (d) X =L ¼ 7:67. PLIIF Step4 Step2 Step1

PLIIF Step4 Step2 Step1

PLIIF Step4 Step2 Step1

0

0

0

-1

-1

-1

-1

-2

-2

-2

-2

-3

-3

-3

-3

Y/L

Y/L

0

Y/L

Y/L

PLIIF Step4 Step2 Step1

-4

-4

-4

-4

-5

-5

-5

-5

-6

-6

-6

-6

-0.2 0 v/Uin (a)

-0.2 0 v/Uin (b)

-0.25 0 v/Uin (c)

-0.25 0 v/Uin (d)

Fig. 6. V-velocity profiles calculated with the three meshes (denoted by lines, step1 ¼ fine, step2 ¼ medium and step4 ¼ coarse) and their comparisons with experimental measurements (denoted by symbols) at (a) X =L ¼ 0:97; (b) X =L ¼ 2:76; (c) X =L ¼ 3:99; (d) X =L ¼ 7:67.

130

Y. Zhao / Computers & Fluids 33 (2004) 119–136

dimensions are given in Fig. 2. A total of three grids are used in the calculations to ensure grid convergence. The finest grid shown in Fig. 3 has 38,962 nodes and 76,402 cells and the distance of nodes adjacent to the solid wall in the form of y þ is kept within unit one. The medium grid has 19,495 nodes and 37,537 cells while the coarse grid has 14,182 nodes and 26,943 cells. Fig. 4 shows the convergence histories in the form of RMS of all the flow quantities in the whole field versus the numbers of physical time steps. Fig. 4(a) shows the convergence history with the explicit pseudotime-stepping scheme (implicit physical time-stepping) on the medium grid, while Fig. 4(b) represents the convergence history with the same scheme on the fine grid. Fig. 4(c) presents the convergence history with the implicit dual time-stepping scheme on the fine mesh. It is noted that for the explicit pseudo-time-stepping, the CFL is 0.5 and number of subiterations per physical time step is 50. The horizontal axis represents the number of physical time steps used in the calculations. It is obvious that the fine mesh calculation requires about 40% more steps or computing time for the same amount of residual reduction. And it is interesting to note that for the implicit dual time-stepping scheme, the CFL is 10 and the number of subiterations is 15. It takes 50 physical time steps to reduce the residual by almost four order of magnitudes, which is substantially more efficient than the explicit pseudo-time-stepping scheme. The computer used in this study is a digital AlphaServer 8400 5/440 running Digital UNIX. The CPU time for the explicit pseudo-time scheme using the fine grid is 5800 min, while the matrix-free implicit scheme uses 190 min for four orders of magnitude in residual reduction. And both methods have almost the same memory requirement, which is 511 Bytes per node.

calculated PLIF

calculated PLIF

calculated PLIF

calculated PLIF

0

0

-1

-1

-1

-2

-2

-2

-3

-3

-3

0

-1

Y/L

-3

Y/L

Y/L

Y/L

-2

-4

-4

-4

-5

-5

-5

-6

-6

-6

-4

-5

-6 0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

0.6

0.8

P/Pin

P/Pin

P/Pin

P/Pin

(a)

(b)

(c)

(d)

1

Fig. 7. Non-dimensional pressure profiles calculated with the fine mesh (denoted by lines) and their comparisons with experimental measurements (denoted by symbols) at (a) X =L ¼ 0:97; (b) X =L ¼ 2:76; (c) X =L ¼ 3:99; (d) X =L ¼ 7:67.

Y. Zhao / Computers & Fluids 33 (2004) 119–136

131

Shown in Figs. 5 and 6 are the calculated velocity profiles in u and v directions with the three sets of grids and their comparisons with the corresponding measurements at X =L ¼ 0:97, 2.76, 3.99 and 7.67, where L is the height of the step. It is clear that the three grids do not produce significantly different results and the fine grid was later adopted for all the calculations. Good agreement between the calculated u and v profiles and their measurements are observed. Fig. 7 shows a comparison of the calculated pressure profiles using the fine grid at different x locations and their corresponding measurements. The expansion fan can be found in all the figures and the shock wave can be clearly seen in the lower parts of the last two figures. And it is noted that the agreement between calculations and measurements is good. Fig. 8 shows the calculated pressure contours together with those calculated (using the Cross-MacCormack scheme [25] and the Baldwin–Lomax turbulence model [4]) and measured by Eklund et al. [24], Although a direct, quantitative comparison is impossible due to the use of different scales, they all show very similar

Fig. 8. Pressure contours: (a) calculated Eklund et al. [24]; (b) measured by Eklund et al. [24]; (c) calculated by the present method.

132

Y. Zhao / Computers & Fluids 33 (2004) 119–136

flow structures qualitatively. The expansion fans behind the step and the reattachment shock after the expansion fan have all been correctly captured. Similarly, the calculated and measured temperature contours as well as velocity vectors and streamlines, in the vicinity of the step, are compared in Figs. 9 and 10 and qualitative agreement is noted in the comparisons. The calculated reattachment length is 3.96 step heights which agrees with the value of 4 step heights given by Eklund et al. [24]. Steep temperature gradients are observed in the shear layer and reattached boundary layer in both the results. The effects of expansion fan and the reattachment shock on the velocity field can be clearly found, which turn the flow toward and away from the solid wall behind the step. Finally, the calculated density, Mach number and q contours are presented in Figs. 11–13. In the density contours, there is a strong density gradient across the free shear layer, and the expansion fan as well as the oblique shock are also clearly shown. From the Mach contours, the recirculation zone can be identified and the point of reattachment is found to agree

Fig. 9. Temperature contours: (a) calculated Eklund et al. [24]; (b) measured by Eklund et al. [24]; (c) calculated by the present method.

Y. Zhao / Computers & Fluids 33 (2004) 119–136

133

Fig. 10. Streamlines and velocity vectors: (a) calculated Eklund et al. [24]; (b) measured by Eklund et al. [24]; (c) calculated by the present method.

well with experimental measurement. As expected, high-production rate of turbulent kinetic energy along the shear layer results in high q values in this narrow stripe with maximum q just above the reattachment point.

7. Concluding remarks A robust matrix-free implicit dual time-stepping scheme has been developed and implemented with a LRN q–x turbulence model in a edge-based high-order upwind FV solver on unstructured grids. The developed method is as memory efficient as an explicit scheme and it can also use large CFL number like an implicit scheme even for complex turbulent flow calculations with good

134

Y. Zhao / Computers & Fluids 33 (2004) 119–136

Fig. 11. Calculated density contours.

Fig. 12. Calculated Mach contours.

numerical stability. In addition, it can be used as an unified approach for calculating both steady and unsteady flows with equally large CFL number in pseudo-time-stepping while the physical time step is controlled for numerical accuracy. The method has been validated by comparing the calculated results of a supersonic turbulent flow over a backward facing step with experimental measurements. These characteristics will make it a good candidate for hybrid RANS/LES simulation. Research is being done to evaluate this method for the calculation of unsteady turbulent flows.

Y. Zhao / Computers & Fluids 33 (2004) 119–136

135

Fig. 13. Calculated q contours.

Acknowledgement The author wishes to thank Caroline Tay for providing some of the data.

References [1] Gerolymos GA. Implicit multiple-grid solution of the compressible Navier–Stokes equations using k–e turbulence closure. AIAA J 1990;28(10):1707–17. [2] VUB/FFA. Turbulence models in EURANUS and the 3D Delery bump. Technical Report SNWP3.3/01, VUB, Pleinlaan 2, 1050 Brussels, Belgium and FFA, P.O. Box 11021, S-161 11 Bromma, Sweden, 1993. [3] Mavriplis DJ, Martinelli L. Multigrid solution of compressible turbulent flow on unstructured meshes using a twoequation model. Int J Numer Meth Fluids 1994;18:887–914. [4] Baldwin B, Lomax H. Thin layer approximation and algebraic model for separated turbulent flows. AIAA paper AIAA-78-257, AIAA, 1978. [5] Vandromme D. Introduction to the modelling of turbulence. Lecture series, Von Karman Institute for Fluid Dynamics, 1992–02. [6] Zhao Y, Ding Z. Computation of shock/boundary-layer interactions in bump channels with transport-type of turbulence models. Comput Meth Appl Mech Eng 1999;173(1–2):55–69. [7] Zhao Y. Stable computation of turbulent flows with a LRN two-equation turbulence model and explicit solver. Adv Eng Software 1997;28:487–99. [8] Mohammadi B. Complex turbulent compressible flow computation using a two-layer approach. Int J Numer Meth Fluids 1992;15:747–71. [9] Coakley TJ. Turbulence modeling methods for compressible Navier–Stokes equations. Paper 83-1693, AIAA, 1983. [10] Gerlinger P, Bruggemann D. An implicit multigrid scheme for the compressible Navier–Stokes equations with lowReynolds-number turbulence closure. ASME J Fluids Eng 1998;120:257–62. [11] Luo H, Baum JD, Lohner R. An accurate, fast, matrix-free implicit method for computing unsteady flows on unstructured grids. Comput Fluids 2001;30:137–59.

136

Y. Zhao / Computers & Fluids 33 (2004) 119–136

[12] Lin PT. Implicit time dependent calculations for compressible and incompressible flows on unstructured meshes. MasterÕs thesis, Department of Mechanical and Aerospace Engineering, Princeton University, 1994. [13] Ahmed F, Tai CH, Zhang BL, Zhao Y. Simulation of unsteady incompressible flows with moving boundaries. In: Wu JH, Zhu ZJ, editors. ACFD 2000 Beijing; 2000. p. 108–17. [14] Venkatakrishnan V, Mavriplis DJ. Implicit method for the computation of unsteady flows on unstructured grids. J Computat Phys 1996;127:380–97. [15] Rogers SE, Kwak D. Steady and unsteady solutions of the incompressible Navier–Stokes equations. AIAA J 1991;29(April):603–10. [16] Hirsch C. Numerical computation of internal and external flows, vols. 1–2. John Wiley & Sons; 1988. [17] Zhao Y, Winterbone DE. The finite volume FLIC method and its stability analysis. Int J Mech Sci 1995;37: 1147–60. [18] Roe PL. Approximate Riemman solvers, parameter vectors and difference schemes. J Comput Phys 1981;43: 357–72. [19] van Leer B. Towards the ultimate conservation difference scheme V, a second-order sequel to GodunovÕs method. J Comput Phys 1979;32:101–36. [20] Zhao Y, Zhang B. A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids. Comput Meth Appl Mech Eng 2000;190(5–7):733–56. [21] Yee HC. Upwind and Symmetric Shock-Capturing Schemes. Technical Memorandum 89464, NASA, Ames Research Center, Moffett Field, California 94035, USA, 1989. [22] Steger JL, Warming RF. Flux vector splitting of the inviscid gas dynamic equations with application to finite difference methods. J Computat Phys 1981;40:263–93. [23] Mcdaniels J, Fletcher D, Hatfield Jr R, Hallo S. Staged transverse injection into Mach 2 flow behind a reawardfacing step: a 3-D compressible test case for hypersonic combustion code validation. Paper 91-5071, AIAA, 1991. [24] Eklund DR, Fletcher DG, Hartfield RJ, Northam GB, Dancey CL. A comparative computational/experimental investigation of Mach 2 flow over a rearward-facing step. Comput Fluids 1995;24(5):593–608. [25] Carpenter MH. A high order numerical algorithm for supersonic flows. In: 12th Int Conf Numer Meth Fluid Dynamics. Springer, 1990.