A novel implicit difference algorithm of primitive variable flow equations with higher-order extra terms

Computers & Fluids

Vol. 26, 0

Pergamon PII: s0045-7930(%)ooo3&9

No. 2, pp. 163-182, 1997 1997 Elsevier Science Ltd

Printed in GreacBritain. All rights reserved 0045-7930/97 $17.00+ 0.00

A NOVEL IMPLICIT DIFFERENCE ALGORITHM OF PFlIMITIVE VARIABLE FLOW EQUATIONS WITH HIGHER-ORDER EXTRA TERMS E. Y. K. NG*’ and S. Z. LIU’ ‘School of Mechanical and Production Engineering, Nanyang Technological University, Nanyang Avenue, Singapore 639198 *Defence Science Organisation, 20 Science Park Drive, Singapore 118230 (Received 13 July 1995; in revised form 20 June 1996)

Abstract-A stable scheme for solving the Euler and Navier-Stokes compressible primitive variable equations by adding higher-order extra terms to the flux-vector splitting method has been developed. The approach eliminates approximate factorization and matrix operation. This gives an improvement in computational efficiency as compared to the standard Beam-Warming and Steger implicit factored schemes. The method is spatially second-order accurate, fully conservative and implemented with general coordinate transformations for treating the complex geometries. It shows good convergence rate and acceptable accuracy in capturing the shock waves. Results computed from the code include transonic flows through turbine stator, compressor rotor with large separate turbulent flow and supersonic compressor staggered wedge cascade. Comparisons with well-documented experimental data are included and the agreements are good. 0 1997 Elsevier Science Ltd. All rights reserved.

NOMENCLATURE static pressure = pe(y - 1) Jacobian Prandtl number = C&k molecular viscosity total internal energy = ~(e + 0.5~‘) specific internal energy e = c,T coefficient of thermal conductivity kinematic viscosity = p/p viscous stress tensor viscous and heat effect = t,.,u, + k/c, de/ax, Superscripts + stagnation condition

n i

time level index in marching direction

1. INTRODUCTION

Most of the popular time-marching algorithms use central-differencing or partial upwinding such as the Beam and Warming [l], MacCormack [2] and Thomas and Walters [3] schemes. Their common features; are the use of explicitly added artificial dissipation or smoothing terms to capture shocks and constrain neutrally stable decoupled modes for the steady-state solutions. Clearly, this artificially added viscosity could reduce the accuracy and robustness. Another class of algorithms gaining popularity are flux vector-splitting methods [46] and splitting coefficient matrix techniques [7,8]. Here, the flux vectors are split into forward and backward contributions based on an eigenvalu#e decomposition and differenced accordingly so as to gain the advantages of upwind-differencing. For example, the time-dependent Euler equations form a system of hyperbolic equations, and upwind-differencing models the characteristic nature of the equations in which information at each grid point from the directions is dictated by the theory of characteristics. Upwind schemes, have the benefits of being naturally dissipative [9]. Separate spatial dissipation terms, such as those generally needed in a central-difference methods to overcome oscillations or *Corresponding

author. 163

164

E. Y. K. Ng and S. Z. Liu

instabilities arising in regions of strongly varying gradients, need not be added [lo]. The splittings of Steger and Warming [6] and Van Leer [l l] are most popular. On the other hand, an implicit scheme which determines the unknowns at each new time-level through the solution of a linear algebraic system to avoid the CFL limitation for stability along the time direction has been implemented widely, especially in solving the Navier-Stokes equations for improving its convergence rate and computational efficiency. The development of an implicit scheme is classified into two approaches. (a) The implicit approximate factorization (AF) method, such as the schemes by Beam and Warming [l] and MacCormack [12]. These methods are widely used because they can solve the block-tridiagonal (or block-bidiagonal) algebra equations at each dimension, with less computational effort. The scheme has been modified and improved by Lyer and Von Lavante [ 131,Dawes [ 141and Jameson [ 151.They used the spectral radii technique instead of the matrix calculation. The characteristics of these methods are that they stress simplifying the implicit algebraic equations and enable the higher-dimension problem to be reduced to the lower-dimension one. But, the use of approximate factorization will inevitably increase the truncation error with large time step. Thus, the precision and stability of the scheme are affected, and this can cause an unstable situation when it is applied to the 3-D problem. (b) The method without any approximate factorization, solving the five diagonal (2-D) or seven diagonal (3-D) matrix equations directly by using some effective algebraic methods, such as the strongly implicit procedure [16] and Gauss-Seidel interaction method [17]. The method has good compatibility where higher CFL number can be used, however, its drawback is that a greater amount of computational effort is required. This paper will combine the advantages of the above two approaches (i.e. by using the spectral radii technique to simplify the calculation and by avoiding the approximate factorization to increase the CFL number and stability) to derive an improved implicit scheme characterized by exact factorization and not solving the block-diagonal system of equations. The first-order implicit scheme for Euler equations is developed initially and then extended to the second-order accuracy for Navier-Stokes equations. The method is based on the time-dependent higher-order extra terms in which the implicit smoothing operators are formed directly. It is a Jacobian spectral radius scheme and a non-approximate factorization scheme as compared to the conventional implicit approaches. 2. MATHEMATICAL

MODEL

2.1. The governing equations in body-fitted coordinate system The 2-D non-dimensional written as f18]:

Navier-Stokes

(N-S) equations, in conservative vector form can be

(1) where P

u=

PUY

;;

PUYUX PUYUY + p (E + VU,

F=

[ EY

Re

824, ziJ = 2P axj -

B, = rLj”j +

- -peLa* = Reynolds no., P _

23 ?j P axi when(i=j);ri,j=p(z

L = reference length.

+ z)when(i#j);

k ~(CuT) - Viscous effect and heat conduction. axi

C,

A novel .implicit difference algbithm

Similarly, the 2-D non-dimensional

165

Euler equations are written as:

~+cE+E=()

(2)

aY

where U, F and G are the same vectors as in equation (1). Here, a genera.1 transformation of the form

is applied to computational automatically, computational

transform the governing equations from the physical domain (x, y) to the domain (5, q) because this could make the physical boundary be the coordinate line and hence simplify the treatment of the boundary conditions and improve the accuracy at the boundaries. The inverse of the transformation is as follows:

x = x(54) = Y(t,rl)

Y

I

By using the Jacobian of the transformation J-

1 a0 _ 7 a(x,Y) Xt;Ys- X,Yl

we can derive

L = JY,,, 5, = - Jx,,

rlx =

-

JY,

qv = Jx,

If we apply the transformation to equations (1) and (2), together with the chain rule of partial differentiation, tlhe following 2-D transformed N-S (3) and Euler equation (4) are obtained: (3) and (4) where

Here

are the contravariant

velocities along < and 11coordinates,

respectively.

E. Y. K. Ng and S. Z. Liu

166

2.2. Explicit boundary conditions 2.2.1. Farfield boundary conditions. The solution at the farfield boundaries is not known a priori since it is affected by the flow within the computational domain. Thus, boundary conditions must evolve with calculation. Hedstrom [19] has shown that this can be done by performing a characteristic analysis of the 1-D equations normal to the boundaries. The number and type of conditions specified are determined by those characteristics that enter the domain as presented in Fig. 1. At the inlet, the right-moving characteristics enter the domain. Thus, four boundary conditions are specified for supersonic inflow, and three are specified for subsonic inflow.

(a) Supersonic inflow: the inlet solution is completely evaluated based on the inlet freestream conditions M>

1: u, = u,,,

In practical application, stagnation pressure P*, stagnation temperature T*, gas flow angle o! and inlet Mach number M, are specified. (b) Subsonic inflow: P*, T* and CIare specified, another condition corresponding to the characteristic u - a is evaluated based on information extrapolated from the interior, k?P/a< = 0 (this condition is also termed supplementary boundary condition). At the exit, the right-moving characteristics leave the domain. Thus, no boundary conditions are specified for supersonic flow, and only one is specified for subsonic flow. (c) Supersonic outflow: the exit solution is obtained completely by extrapolation

i.e. four supplementary

boundary conditions are given at the outlet. In practice, we use

ap/ag =

0,

afyag =

0,

au/a< = 0, avjag =0

(d) Subsonic outflow: static pressure is specified. Density p, velocities u and U, are extrapolated from the interior field, e.g.

ap/at = 0, au/at =

0,

au/at =0

2.2.2. Solid boundary conditions. This is a no-slip condition. Thus, for a stationary wall, we have

uies

Fig. 1. Inlet and exit farfield boundaries for a typical geometry.

A novel implicit difference algorithm

167

at wall u = 0, v = 0 and the wall pressure and density conditions

are used to determine these parameters. However, in the case of inviscid flow, it is necessary to satisfy the slip condition. This is expressed by the vanishing of the normal velocity u, = 0. Other parameters, such as P and p can be determined by the following relations:

[d(;

+ +:a,,)=0

where y is specific heat ratio ( = C,/C,) and R is curvature radius of the wall. In equation (S), the first relation considers the effect on the solid surface pressure P and density p with the radius of curvature R. The second and the third equations are state (by assuming T is approximately constant) and Bernoulli equations, respectively. 3. NUMERICAL

ALGORITHM

This section describes the development of the improved flux vector-splitting scheme. Various aspects of the algorithm including solution methodology, implicit boundary conditions, stability analysis, time-stelp size and convergence criteria are also included. 3.1. The basic implicit scheme From the Steger-Warming [6] technique of flux vector-splitting, for example, the generalized flux vectors F and G of every computational point can be split into two subvectors F * and G *, respectively, according to the positive and negative eigenvalues. The subvectors F + and G + which have positive eigenvalues use the backward-difference operator whereas the subvectors F- and Gthat have negative eigenvalues apply the forward-difference operator. The Euler equation (4) can be rewritten as

where A and B are the Jacobian matrices aF/aU and aG/aU , respectively. In general, they are homogeneous matrices which exist as a similarity transformation [19], such that Q-‘AQ = A(A) = diag[8,0,0

+ aL({),8

- aL(<)]

and Q-‘BQ = A(B) = diag[t,v,p+

aL(q),8-

aL(q)]

where A is the diagonal matrix, a is the speed of sound in the physical domain, and

A, = 1, = 8, Izs= 6 + aL(r) and A, = 8 - aL(r) for flux vector Fand R, = 2, = 8,1, = v + aL(q) and & = v - aL(q) for flux vector G. Based on these positive and negative eigenvalues, the expressions for calculating the subvectors

E. Y. K. Ng and S. Z. Liu

168

F * and G * can be simplified and summarized as follows

F+

=

C-1 Gu + C2aL

+

2Y

C,v CI(u2 + v2)/2 + GaLu

+

C2a& + Ga&.v + &a’/(7 - 1) 1

c, G’

=

Y!t! 2Y

0

+ C2arj,

C,v

+

C2afy

C,(u2 + v2)/2 + C,af.,u + C,at&v + C3a2/(y - 1)

1

where C, = 2(y - 1)1, + & + &, C, = A, - &, C, = ;1, + 1, and ?.Y= (,/L(5), & = &/L(5), &/L(q), ri, = q,,/L(q) are the constants on computational point (i,j). When calculating the F+ and G+ terms, the positive eigenvalues are used with any negative eigenvalue set to zero. Conversely, when calculating F- and G- terms, the negative eigenvalues are used with the positive one set to zero. The final equations for solving inviscid flow in the flux vector-splitting form are thus $+-

8F’ a<

+ -L’F- + -dG+ a< alt

I dGaq

_.

(6)

For the N-S equation (3), the right-hand side terms, i.e. viscous flux vectors Q and R mainly consist of the first-order derivative of variables U, v and T. Firstly, the direction of the information-dependent region for N-S solutions corresponding to Q and R does not change with the velocity, therefore it is not necessary to choose the difference direction according to the eigenvalues. Secondly, the first-order derivatives pc, z+, v, and I; in the Q and R only bring the upward influence for the solution while the second-order derivatives per etc. induce both the upward and backward influences for the solution. Therefore, flux vectors Q and R should use the two-side differencing and need not be split into the subvectors. The final equations with flux vector-splitting form for viscous flow, are therefore obtained accordingly as follows:

cc!+- aF+ + -aF- + -aG+ +

at

at

all

aG_+@+~)

all

In brief, the problem that exists in all of these schemes is that the approximation factorization techniques are all adopted in their smoothing for implicit calculation. These cause the difficulties in the stability for the 3-D calculations. 3.2. The formation of the present algorithm In order to form the implicit part, combined with vector flux characterization direction splitting, the following higher-order extra small terms:

of different

can be added to the right-hand side of equations (6) and (7) directly to form the implicit schemes for the Euler and N-S equations, respectively. However, the last two second-order terms are zero for the Euler equations. The higher-order extra terms - At/2 alag (&+&!Jlat ), At/2 a/@ (2; tXJ/at ) are the implicit construction items for the convection fluxes F respectively, along the 5 direction, where J: are the scalars corresponding to the eigenvalues of the matrices F Likewise, - At/2 al&l (A,,+aUl& ), At/2 a/@ (A; atif/& ) are the implicit construction items for the fluxes G respectively, along the

*,

*.

*,

169

A novel implicit difference algorithm

q direction, where 1: are the scalars corresponding to the eigenvalues of matrices G * . This newly formed implicit scheme is stable if we could select the correct values for the scalars (i.e. nt and n,,+ are relative to the spectral radii of the Jacobians of the convection and dissipation terms of equations (6) and (7). These four extra terms have some order with the second-order time derivative terms. They are the unsteady items corresponding to the bU/at , therefore, we can say that as far as the theory is concerned, they have no influence on the solution once it has converged to the steady state. Depending on the scheme used, these terms may adapt minus or plus to allow the genuine coefficient to be minus or plus to form dominant main diagonal elements. The implicit smoothing operator can be established directly with these extra terms. By differencing the Euler equation (6), we obtain

U?+‘I.1

-At-

*VP+

IY;~= -At

+ *AF-

+ (VG’ A?

+ A,G Arl

1

V&++“U,l;‘) - Vx(&+“U;,) + At A.&“U;,?‘) - A&“q,) 265 2A<

_,&(J;“U;;‘)

- VY(l,,YJ;j) + At A#,“U;,!‘)

2A1

- AY(n;“U;j) 2Aq

Letting

AU;, = - At

VF+

$q-

AF-

+ %

VG+

+ +-

+

A,G A,,

1

and SU;/ ’ = Uj’j’ ’ - U;’

the above equations can thus be rearranged as follows: VF+ AQj = - At *

AF+ *

+ VxG+ A?

I

AGx Atl

1

where V, = backiward operator = 0.5(3Uj - 4uj _ I + Uj- *) and, A, = forward operator = 0.5( - 3uj + 411j+,- u,+*). In the above scheme, equation 8b can be considered as an explicit step that describes the physical conservation law while equation 8c can be considered as an implicit step to ensure the calculation stability for a large time step. For inviscid flow, the explicit step AUtj can be differenced with either first-order accuracy by using two-point backward and forward methods or second-order accuracy by using three-point backward and forward methods. In scheme (8a), as the terms nt and n: are all scalars, the computational efforts are thus reduced greatly as compared to the other coefficient matrix implicit schemes for all time steps, such as those in [7] and [8]. The major characteristic of the present scheme is the use of exact factorization approach. Furthermore, the above scheme is a one-step method which has the added advantage for calculation wi.th a large time-step where the same difference scheme should be used to obtain the convergent solution without oscillations. Similarly, we could obtain the difference scheme for N-S equation (7) as follows: ASUl?,f

+ BSU;,+_L, + CSU;j”

+ DSU,l/‘,‘l + ESU~~/j = AU!’

(9)

E. Y. K. Ng and S. Z. Liu

170

where AU;’ = - At

--aF+

ar

+

aF-

ag

+

aG+

all

+

and A, B, C, D and E are coefficients of difference equations that are related to eigenvalues. It is important to note that, to achieve high accuracy results, higher-order discretization schemes can be employed only in the explicit part, Au”, which do not increase computing effort and complicity in the implicit part. Here, a single-step three-node upwind-difference method for vector flux and a central-difference approximation to viscous terms are adapted. Thus, the novel scheme has no second-order numerical dissipation or artificial dissipation and it has advantages in aspects of robustness, accuracy and ease of computation. 3.3. Solving methodology In brief, the established implicit vector flux-splitting scheme consists of EXPLICIT

Au” = - At

IMPLICIT

B6lJ[‘_‘, •l- CSU$

UPDATING

Vi:; ’ = U; j + 6 U:; ’

“L;f-+

+

A<_F-

+ A,,+G+ + A,-G+

At

•l- DSU;T+‘l = AU:‘-

Art AGU:T/j-

A? EdUF+lj (10)

where A,lA< , A,/Aq , mean the central-difference operators while AC*/At , A,,* /Aq , stand for the second-order backward/forward-difference operators along <- and q-directions, respectively. The equation (10) can be solved easily by the Thomas algorithm [18]. Just two sweeps are required to solve the equations. Calculation begins with I = 1 station until I = II along the r-direction. At each I station, the SU:.?& term is just obtained and SU;, ,j is available with the value of the last time step, so the right-hand sides of equation (10) are all known. In this way, the unknown W;,?‘,, &!I;~ ’ and SU;/‘,‘, can be obtained directly by the Thomas algorithm along the q-direction. Hence, an approximate factorization and coefficient block-diagonal system of equations can be avoided successfully by solving equation (10) and the algebraic equations derived from implicit operator are simple and easy to solve. Notice that from equation (lo), this implicit flux vector-splitting scheme consists of an explicit

171

A novel implicit difference algorithm

part AU,?, and an implicit part SU,,j. Explicit boundary conditions have been given in Section 2.2. Here, we discuss only the implicit boundary conditions 6U(,,. We propose 6U = 0 for the inlet, outlet and solid implicit boundary conditions. For the steady-state calculation, when it approaches a convergent solution, XJjar will tend towards zero, i.e. 6U -+ 0. Therefore, at least in theory, this hypothesis is correct. For the periodical boundary condition that exists in the cascade calculation, the more realistic approach is to let SUi,l = SUi,,. This relation can be realized by using the Thomas algorithm [ 181 which includes the periodical boundary conditions to solve equation (IO). 3.4, Stability

study

The 2-D N-S equations in flux vector-splitting

form (7) can be written as

+

N

t?!! (11)

‘I afj*

where A*, B*, MC, Ad,,, N, and N, are all matrices with

N,=

L

and N, = ~

aR

It is very difficult and tedious to perform a stability analysis for an implicit scheme directly. However, according to the I-D model equations’ [20] and the McCormack scheme’s [12] stability analysis together with the characteristic theory of equations, we propose the following stability conditions for the current implicit scheme (10):

A? = max P(A*) + m

A: = max p(B*)

+ &

2

A5

]P(M,) +

PWJ - z

4

MN,)

PWJ - 2

4

+

where p(A *), p(B*), p(M,), p(M,,), p(N,) and p(N,,) are the spectral radii of the matrices A *, B*, M,, M,,, N, and N,,, respectively. By solving these spectral radii, and putting them back into equation (12), the stability conditions can be obtained as follows:

172

If ii - CJ,/~

E. Y. K. Ng and S. Z. Liu

< 0, then:

else, 1; = 0.

If p - a,/qr + qf < 0, then:

1; =max

1

IP+aJml+

z

( $+r]i 1[ ,)-

$o};

(13)

else, A,,- = 0; where o = max(p/Re ,p’ + 2plRe ,,uylPrRe ), a is the speed of sound, $ is the second coefficient of p, and c = 2. Likewise, we can also obtain the stability conditions for an inviscid flow as follows:

A,+ =max If

0 - a,/‘=

I

l0+aJFLjj-

s,O

IV+aJml-

g,O

I

-c 0, then:

else, A; = 0. &+ =max If v - aJm

1

I

< 0, then:

(14) else, 1; = 0. In principle, when Izs are chosen from equations (10) or (1 l), the implicit method (7) or (6) can be granted as a stable scheme at any time step. Therefore, we could say that the schemes (6) and (7) are unconditionally stable from this point of view. 3.5. Time-step length analysis In practice, the stability conditions suggested in equations (13) and (14) cannot guarantee that the vector-splitting schemes (7) and (6) are unconditionally stable. The selection of time-step length is still constricted by various factors. The main factor is that the smoothing for implicit simulation is approximate, and this will certainly cause some errors. Other factors, such as local linearization in differencing the non-linear equations, the complex geometries and its grid distribution etc., have not been included in the stability analysis conditions (13) and (14). Therefore, they are not the sufficient conditions for a specific problem to conform to unconditional stability in which the time-step must be related with a CFL number and is determined by At = min{At,,At,}

A novel implicit difference algorithm

N (time step)

N (time step) Fig. 2. Typical convergence history with different CFL numbers (inviscid flow-convergence-divergence nozzle).

For

Fig. 3. Typical convergence history with different CFL numbers (viscous flow-convergence-divergence nozzle).

invicid flow: At, I

At,, I

CFL.A<

CFL*Aq

For viscous flow: At, I

At,, I

CFL*A<

CFL.Aq

It is obvious that when the CFL number is smaller than or equal to 1.0, At is the time-step length for an explicit scheme, and nt and A,,+of the whole field will go to zero automatically. Therefore, the implicit scheme will reduce to the explicit scheme in this case. In the current calculation, we take a large time step and the CFL number for the implicit scheme is chosen to be 2 30 (i.e. CFL can be raised by 1 to 2 orders of magnitude). 3.6. Convergence consideration There are various methods to check if the result has converged to the steady-state solution. In the present work, the following two convergence criteria are used:

where U is the same term as used in the schemes (6) and (7), AU is the difference between the two iteration steps, cp is any other gas dynamics parameter and (Acp/cp),,,., is the maximum relative difference of thal. gas dynamics parameter at any grid point.

174

E. Y. K. Ng and S. Z. Liu

4. COMPUTATIONAL

RESULTS

AND

SAMPLE

APPLICATIONS

Based on the efficient flux vector-splitting implicit scheme, the algebraic H-grid generator and Baldwin-Lomax turbulent model [21], a software package for solving the 2-D Euler and time-averaged N-S equations is developed. Code verifications have been conducted for various inviscid and viscous cases, with separated and turbulent flow, that include transonic turbine stator, compressor rotor and supersonic compressor cascades. 4.1. Study of the typical convergence histories with d@rent

CFL numbers

Typical convergence histories with both inviscid and viscous flows are presented in Figs 2 and 3 for the computation of convergence-divergence nozzle flows [22] with their grid distributions of (31 x 25) and (41 x 25), respectively, where the operating mode of supersonic flows at exit are being considered. Figure 2 reveals the convergence history with CFL = 0.8, 5, 10,25, 100, respectively, for inviscid flow. From this figure, we can see that the convergence rate is improved significantly when larger CFL numbers are used. For example, the norm residual reduces four orders of magnitudes with 300 iteration time steps for CFL = 25, as compared to 1000 iteration steps for the case of CFL = 5 to reduce the same order of magnitude. Also, for the case of CFL = 0.8, which means the implicit scheme reduces to the explicit scheme, it only drops 1.5 orders of magnitudes with 2000 iteration steps. However, there is no obvious effect on improving the convergence rate if the CFL number is continually increased. For example, there are little differences for the convergent rate among the cases with CFL = 100, 50 and 25. Figure 3 indicates the convergence history with CFL = 0.8, 2, 5, 8, 10, 12, respectively, for viscous flow. Because the flux vectors associated with viscous terms Q and R are not split into the upwind forms (central-differencing is used here) as highlighted in Section 3.2, the convergence ability is influenced and thus the CFL number cannot be given as large as in inviscid flow calculation. Nevertheless, the convergence rates are improved significantly when larger CFL numbers are used. For example, the norm residual reduces four orders of magnitudes with 600 iteration time steps for CFL = 10, as compared to the case of CFL = 0.8, which represents the explicit scheme, and only drops one order of magnitude with 2000 iteration steps. In general, it is better to choose the CFL number in the range of 20-30 for inviscid calculation and between 5 and 10 for viscous flow computation when using the present efficient flux vector-splitting implicit scheme. In such a case, only about 300-500 iteration time steps are needed to obtain the satisfactory convergent result for inviscid calculation and about 600 iteration time steps for viscous prediction. The code runs at about 1.042 x 10 -4 s per point per time step for inviscid flow and about 1.902 x 10m4s per point per time step for viscous flow on a VAX-9000 machine. However, this is still a research state with a larger amount of overhead computation. It is believed that the run time can be reduced with a simple code restructuring. 4.2. Grid refinement studies The effects of mesh refinement with high-grid-aspect ratios are well-known to impact convergence adversely for most types of computational algorithms (e.g. relaxation methods). These issues are considered here by means of vector stability theory and computational experiments but not as an end in itself. The results indicate that the current novel ‘intermediate’ scheme experiences slight convergence deterioration for the same CFL number due to non-optimum local time-stepping procedures despite the fact that the increased need for inviscid and viscous preconditioning is not required here since approximate factorization and matrix operation are avoided. An improved selection of the local time and perhaps proper implementation of the boundary conditions [25] are needed to provide uniformly efficient convergence at all aspect ratios (see the Appendix for more details). A systematic grid refinement study is also conducted to provide quantitative estimates of grid requirement for accurate resolution for typical N-S calculations on a transonic compressor rotor [25]. It suggests that the mesh of 86 x 45 is sufficient for this flow. Figure 4 shows the impact on the accuracy with the level of grid refinement. The computed error in mass flux (defined as I[(& - riz,xP)/riz,,IIx 100%) against level of grid refinement (i.e. 2 x , 4 x , 6 x and 8 x of the base mesh 43 x 23) tends to converge/reduce to a certain level gradually (i.e. expected asymptotic improvement in accuracy).

175

A novel implicit difference algorithm

Error in mass flux

Gradually “level-off” -__w_ *‘(TsyTnptotic improvement)

Ok

,

4x

I

6x

I

1

f

I

8x

NX

Level of Grid refinement Fig. 4. Grid refinement studies.

4.3.

Fig. 5. Grid

used for transonic turbine (inviscid and viscous).

flows

Transonic twrbine cascade [23]

The second test case is the inviscid and viscous calculations through a transonic turbine cascade. The code used a simple H-grid system with (41 x 25) points for both viscous and inviscid flows as shown in Fig. 5. The flow conditions are as follows: Inlet stagnation pressure Inlet stagnation temperature Reynolds number {Outlet static pressure Inlet flow angle Mass flux Cascade pitch

PO* = 199 555 NmV2 TJ’ = 288’K

= 5.102e +6 Pb = 102,318.3 Nm-*

CI= - 21.2” = 22.12 kg s’ t = 0.09163 m

Comparisons between the inviscid and viscous calculations are conducted. Figure 6 presents the surface pressure distributions. The triangle and square symbols indicate the measured data obtained from Ref. [23]. The short and long dash lines illustrate the current inviscid prediction on suction and pressure surfaces, respectively, while the solid and dash-dot-dot lines represent viscous prediction results on suction and pressure walls. The shock wave is captured accurately in both cases as compared to the experimental result, except a slightly higher prediction occurs in the suction surface ‘pressure distribution with the inviscid flow calculation (due to the zero blade boundary layer as expected). Finally, Fig. 7 illustrates the Mach contour distributions for the inviscid and viscous flow predictions; their differences can be clearly discerned. About 600 iteration time steps are needed for inviscid flow with CFL = 18, whereas 900 time steps are required in the viscous prediction with CFL = 9. -Predicted (suction, NS) _____ Predicted (suction, Euler) .-..- Predicted (pressure, NS) --Predicted (pressure, Euler) Experimental data (suction) Experimental data (pressure)

Fig. 6. Surface static pressure distributions of transonic turbine cascades.

176

E. Y. K. Ng and S Z. Liu

Level MACH D c B A 9 8 7 6 5 4 3 2 I

(a) lnviscid

I.0118 0.950194 0.888588 0.826982 0.765375 0.703769 0.642163 0580557 0518951 0.457345 0 395738 0.334132 0.272526

Level MACH

(b) Viscous

D C B A 9 8 7 6 5 4 3 2 I

flow prediction

Fig. 7. Mach

contours

distributions

(a) Pressure surface 1.1 ,.,,I.““,..,I....I,. 0 I 2 3 Distance

(inviscid

4

1.116140 1.048449 0.946912 0.879221 0811530 0743838 0676147 0.608456 0.540765 0.473074 0.405382 0.337691 0.270000

and viscous).

Fig. 8. Grid (71 x 17) for the supersonic wedge cascade.

staggered

.a.. 5

along axial chord

Level MACH 7 1.61507 6 1.54256 5 1.47005 4 1.39753 3 1.32502 2 1.2525 I 1 1.18

1.6

I .5

1.2 (b) suction 1.1

Distance Fig.

surface

“‘.‘..“I”..“‘..‘.” 0 1 2

9. Isentropic

math

3

4

I”

5

along axial chord

number cascade.

distributions

for wedge

Fig. 10. Predicted

isentropic

isoMach

lines (inviscid

only).

A novel implicit

difference

algorithm

177

2D Inviscid

0 (a) Inviscid

0.02

0.04

(86 x 45 uniform)

Fig. 11. H-Mesh

(b) N-S (86 x 45 uniform7

for DFVLR’s

compressor

rotor

flow.

4.4. Supersonic compressor staggered wedge cascade [24]

This is a fully supersonic flow in a compressor cascade. We compared the inviscid prediction only with the exact solution given in [24] which was derived from characteristic theory and oblique shock relations. A simple H-grid (71 x 17) with 45 points on the wedge was used. The quality of the mesh is rather poor due to the high skewing of the cells as shown in Fig. 8. The inlet boundary condition is supersonic. For the outlet, uniform static pressure is imposed. This is an approximation, because the outlet is not uniform due to the shocks and expansion waves leaving the domain. The outlet has been defined at three times the wedge length, in order to allow the shocks be attenuated by the expansions.

--+---

Expt-SS

-LL

Expt-PS

L

NS-SS

0 0 _..-..-.A -__

Fraction Fig. 12. Mach

number

NS-PS _ _._ Euler-SS

A- -.-.- Eulert-PS

of axial chord distributions

around

blade.

178

E. Y. K. Ng and S. Z. Liu

m: 17, Streamsurface:

4

Blade height:

(a) Measurement

0

0.02

0.04

Fig. 13. Relative

Level

Mach

K J I H G F E D c B A 9 8 I 6 5 4 3 2 1

1.28 1.22 1.15 1.09 1.03 0.96 0.90 0.83 0.77 0.71 0.64 0.58 0.51 0.45 0.38 0.32 0.26 0.19 0.13 0.06

45%

from [24]

Level

Mach

0.06 isoMach

contours

between

prediction

and measurement.

The computed isentropic Mach number distributions of the pressure and suction surfaces are illustrated in Fig. 9. The dashed lines are the exact distribution as calculated using oblique shock and Prandtl-Meyer formulas [24] whereas the circles show the present inviscid prediction. The strong shocks do not produce any significant under and over-shoots. This result is reasonable, taking into account that the mesh is highly skewed, and that the shocks are not aligned with the grid lines. The major defect in the calculated distribution is the varying pressure (which is equivalent to isentropic Mach number aft of the expansion fan). An improvement of the shock prediction could be achieved by using orthogonal and adaptive unstructured meshes. Furthermore, Mach number contours are presented in Fig. 10, where the shocks are reasonably resolved, with almost uniform regions in between. The reflected bow shock is nearly cancelled by an expansion at the corner on which it impinges. In the above calculation, 1500 time steps are used with CFL = 15 to obtain this result.

A novel implicit difference algorithm

179

Level PO 151116.00 140675.09 130234.20 119793.30 109352.40 98911.50 88470.60 78029.70

Level PoRel

0.00

0.02

0.04

0.06 Fig. 14. Stagnation pressure contours.

4.5. DFVLR’s transonic compressor rotor flow [24] The next test case computes the flow field through the transonic axial compressor rotor containing areas of turbulent flow with separation or large recirculation (N-S) and is compared to measurements made with an advanced laser velocimeter at DFVLR [24]. The comparison is made at design speed at a pressure ratio corresponding to peak efficiency. At this measured condition the fan rotor efficiency is 89.6%. The information needed includes blade geometry and rotor rotational speed. To initialize the flow field, the upstream relative Mach number has been obtained from the measurement using velocity triangle analysis. In the absence of such information, approximate guesses may be made for this variable. The calculation was performed with a grid network formed by the intersections of 86 axial and 45 tangential faces. Typical inviscid and viscous computing grids are shown in Figs 1l(a) and (b), respectively.

Vector with stream ilines

(a) Inviscid Fig. 15. Velocity vector with streamlines plot.

(b) N-S

180

E. Y. K. Ng and S. 2. Liu

A comparison between rotor local Mach number measurements (at 45% span) and computations around the blade surface (50% span for both Euler and N-S) is shown in Fig. 12. The general trend of the predicted values is very similar to the measurement. However, the inviscid result shows an over-prediction of the shock wave on the suction surface. Figures 13(a)-(c) compare the measured and computed relative Mach number flow field at 45% (near mid-span) and 50% blade span, respectively. The incoming flow is slightly supersonic. The detached shock wave can be identified at the blade channel entrance. The shock wave is also standing off the blade leading edge implying spilling of the flow. This feature of the flow field is also observed by the present calculations. The shock decelerates from Mach 1.15 to 0.7 at the suction surface. The match between the measurement and prediction is very close in the entrance region. In the covered passage the agreement is less precise. A trailing edge Mach number of 0.65 is experimentally observed, the N-S calculation giving 0.66. This is probably due to the excessive viscous blockage effects in the calculation (which produce less effective area and slightly higher average Mach number) as shown in the stagnation pressure contours plot (Fig. 14(b)). Figure 14(a) [Euler] and (b) [N-S] present the contour of entropy with inlet entropy as zero. There is a much greater rise of entropy in the boundary layer of the blade, wake and near shock; for the other region it is nearly zero. This shows the scheme has a good accuracy. Finally, the velocity vector with streamlines plots is presented in Fig. 15(a) [Euler] and (b) [N-S]. In all, comparisons of the predicted and experimentally determined Mach number contours indicate close agreement in the entrance region where the viscous blockage effects are small. In the above calculations, only about 400 time steps with CFL = 25 and about 700 time steps with CFL = 10 are required for the inviscid and viscous calculations, respectively. 5. CONCLUSIONS The programme with sound structural and versatile ability including pre-treatment, time-marching and post-treatment modules is developed using FORTRAN-77 for both inviscid and viscous 2-D flow calculations. The code is spatially second-order accurate, fully conservative and with general coordinate transformations for treating the complex geometries. Code verifications include transonic flows through convergence-divergence nozzle, turbine stator, compressor rotor and supersonic compressor staggered wedge cascade. Besides possessing the general features of the implicit scheme, it also has the following advantageous characteristics: 1. A significant improvement in computational efficiency is obtained as compared to the well-known Beam and Warming, and Steger implicit factored schemes because it formulates higher-order numerical implicit terms. 2. No approximate factorization occurs in the time-marching process, hence there is no numerical error caused by an approximate factorization, and it helps to increase the CFL number for calculation. 3. Matrix operation is not necessary in the implicit part and the computing effort in each time step can be significantly reduced. 4. The spectral radii technique is used to replace the coefficient matrix calculation, therefore, the 3 x 3 systems of equations are being transformed to the three scalar equations in the 2-D calculation with comparatively less computational effort. 5. The scheme preserves the strong points of the flux vector-splitting scheme with reasonable performance on shock capturing. 6. The present scheme is a one-step scheme which helps to improve the convergent ability especially for large time-step computations, because it can converge to the numerical solution with zero computer error output, in principle. 7. The flux vector-splitting implicit scheme is an upwind scheme, the flux vector is split according to the physical characteristics and corresponding difference schemes are adopted to restrain the numerical errors. Therefore, it is a stable scheme without the need of adding any artificial viscosity (i.e. it has superiority in aspects of robustness, accuracy, calculation and convergence). In summary, the concepts involved in the construction of an implicit scheme with an exact factorization, and avoiding the coefficient block-diagonal system of equations together with their successful applications, provide a prospect for further development into the 3-D problems in the future.

A novel implicit difference algorithm

181

REFERENCES 1. Beam, R. M. and Warming, R. F., An implicit finite-difference algorithm for hyperbolic systems in conservation-law form. J. Compufa~‘. Phys., 1976, 22, 87-110. 2. MacCormack, R. W., An efficient numerical method for solving the time-dependent compressible Navier-Stokes equations at high Reynolds number. NASA TM X-73, July 1976; p. 129. 3. Thomas. J. L. and Walters. R. W.. Unwind relaxation algorithms for the Navier-Stokes eauations. Proc. AIAA 7th Computkonal Fluid Dynamics Coif., Paper 85-1501, June 1985. 4. Beam, R. M. and ‘Warming, R. F., An implicit factored scheme for the compressible Navier-Strokes equations. AIAA J., 1978, 16, 4, 393402. 5. Steger, J. L., Coefficient matrices for implicit finite difference solution of the inviscid fluid conservation law equations. Computer Meth. 61 Appl. Mech. and Engng, 1978, 13, 175. 6. Steger, J. L. and Warming, R. F., Flux vector splitting of the inviscid gasdynamics equation with application to finite difference methods. J. Computat. Phys., 1981, 40, 2, 263-293. 7. Chakravarthy, S. R., Anderson, D. A. and Salas, M. D., The split coefficient matrix method for hyperbolic systems of gasdynamics equations. AIAA Paper 80-0268, 1980. 8. Moretti, G., An old integration scheme for compressible flows revisited. Refurbished and put to work. Compurat. Flu&, 1979, 7, 191. 9. Anderson, W. K., Thomas, J. L. and David, L. W., Multigrid acceleration of the flux split Euler equations. AIAA .I., 1988, 26, 6, 649-654.

10. Taylor, A. C., Ng, W. and Walters, R., An improved upwind finite volume relaxation method for high speed viscous flows. AIAA Paper 89-0549, 1989. 11. Van Leer, B., Flux-vector splitting for the Euler equations. Lecture Notes in Physics, 1982. 170, Springer-Verlag, p. 507. 12. MacCormack, R. ‘W.,A numerical method for solving the equations of compressible viscous flow. AIAA Paper 81-0110, 1981. 13. Von Lavante, E. and Lyer, V. S. V., Simplified implicit block-bidiagonal finite difference method for solving the Navier-Strokes equations. AIAA .I., 1985, 23, 7, 113&l 132. 14. Dawes, W. N., A pre-processed implicit algorithm for 3D viscous compressible flow. In Notes on Numerical Fluid Mechanics. Vol. 12. Vieweg, 1986, pp. 7&77. 15. Jameson, A., Successes and challenges in computational aerodynamics. Proc. AIAA 8th CFD Conf., 1987. AIAA, l-35. 16. Sankar, N. L., Malone, J. B. and Tassa, Y., An implicit conservative algorithm for steady and unsteady three-dimensional potential flows. AIAA Paper 81-1016, 1981. 17. Ames, W. F., Nwnerical Methods for Partial Differential Equations, 2nd edn. Academic, New York, 1977. 18. Anderson, D. A., Tannehill J. C. and Pletcher R. H., Compurarional Fluid Mechanics and Heat Transfer. Hemisphere Publishing Corporation, New York, 1984. 19. Hedstrom, G. W.. Nonreflecting boundary conditions for nonlinear hyperbolic systems. J. Computar. Phys., 1979, 32. 20. Liu, J. M., Liu, B. and Xiang, Y. M., An explicit-implicit numerical method for 3-D compressible N-S equations. Proc. 4th Int. Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Honolulu, 1992, pp. 544-552. 21. Baldwin, B. S. and Lomax, H., Thin layer approximation and algebraic model for separated turbulent flows. AIAA Paper 78-257, 1978. 22. Zheng, Z., Transonic flow calculation for two-dimensional nozzle. J. Aerodynamics, 1989, 7, 4. 23. Zhang, Y., Gong, Z. and Sheng, M., On the solution of transonic flow of plane cascade by a time-marching method. J. Thermophy. Engng, 1980, 1, 4, 337-340. P.R. China. 24. Fottner, i., .Test cases for computation of internal flows in aero engine components. AGARD-AR-275, 1990. 25. Buelow, P. E. O., Venkateswaran and Merkle, C. L., Effect of grid aspect ratio on convergence. AIAA J., 1994, 32, 12, 2401-2408.

APPENDIX TIME-STEP

DEFINITION

AND THE DIFFICULTY

AT HIGH-ASPECT

RATIOS

The present work applies local time stepping to specify a different physical time step (At) at every grid point (i.e. by setting the dominant or maximum CFL at any grid location to a specified ‘optimum’ number). In the physical domain, this max-CFL definition is calculated mathematically by evaluating the minimum time step in the two directions

where n, and zv are the acoustic eigenvalues in the respective coordinate directions. The CFL numbers are typically chosen in the ranges of 20-30 for inviscid and 5-10 for viscous flow optimum convergence. The above time-step selection shows the problem experienced with high-aspect-ratio meshes. If the aspect ratio (AR = AxlAy) is ve&much greater than one {a million is often used locall; in G-S calculations of turbulent‘flow), the CF& (i.e. CFL,) is maintained at the oDtimum value, but CFL, is verv small (i.e. CFL,. = CFL and CFL. z CFLIAR. assuming that n,‘&d or, are of the same drder). This disparity results in-very slow propagation of waves in the streadwise direction, giving poor convergence. The converse is true for ARs that are less than one (i.e. CFL, = CFL and CFL, z CFL.AR, and errors propagate very slowly in the y-direction). As suggested in Ref. [25], for efficient convergence at all aspect ratios, the local time step must be selected based on the min-CFL definition as follows: Ar=max

CFLAx nx

-

CFLA Y

’

-

KY

182

E. Y. K. Ng and S. Z. Liu

At the inlet, the right-moving characteristics enter the domain. Thus, four boundary conditions are specified for supersonic inflow, and three are specified for subsonic inflow. (a) Supersonic inflow: the inlet solution is completely evaluated based on the inlet freestream conditions M> 1: u, = u,, In practical application, stagnation pressure P*,stagnation temperature T*, gas flow angle 01and inlet Mach number M, are specified. (b) Subsonic inflow: P*,T* and a are specified, another condition which corresponding to the characteristic u - (I is evaluated based on information extrapolated from the interior, d P/at = 0 (this condition is also termed the supplementary boundary condition).At the exit, the right-moving characteristics leave the domain. Thus, no boundary conditions are specified for supersonic flow, and only one is specified for subsonic flow. (c) Supersonic outflow: the exit solution is obtained completely by extrapolation A4 > 1: u, = u,,,,a,