Euler pressure relaxation and flux vector splitting

Computers & Fluids Vol. 16, No. 4, pp. 485-490, 1988 Printed in Great Britain. All rights reserved

0045-7930/88 $3.00+ 0.00 Copyright © 1988 Pergamon Press plc

TECHNICAL RNS/EULER

PRESSURE VECTOR

NOTE

RELAXATION SPLITTING

AND

FLUX

S. G. RUBIN Department of Aerospace Engineering and Engineering Mechanics, University of Cincinnati, Cincinnati, OH 45221, U.S.A.

(Received 10 February 1988)

1. I N T R O D U C T I O N

The present author and co-workers have previously presented composite primitive variable [1-4] and velocity/potential [5-8] formulations for the computation of viscous/inviscid interacting flows. These procedures have been designed to emulate large Reynolds number (Re) asymptotic behavior, but, at the same time, to allow for the efficient solution of a single composite system of equations. This methodology is intermediate to that of full time-dependent Navier-Stokes solvers and to that of interacting or inverse matched boundary layer-inviscid solvers. The resulting system of reduced Navier-Stokes or RNS equations is a composite of the full Euler and second-order boundary layer systems. The neglected Navier-Stokes or diffusive terms are higher-order in Re for appropriate "streamline" coordinates. Although, these terms can be retained through an explicit deferredcorrector, it must be emphasized that the boundary conditions and the discrete approximation to the differential equations are dictated solely by the form of the lowest-order implicit RNS operator. The RNS model, which is valid throughout the entire Mach number (M) range, represents an enhancement of the hypersonic merged layer [9, 10] and thin layer approximations [1 1], wherein the axial (x) or streamwise pressure gradient (p~) is prescribed, or of the PNS model [12, 13], where Px is approximated with the truncated form copx. The parameter o~(M) is a function of the local Math number and is such that 0 ~
486

Technical Note

This also allows for more efficient application of far-field boundary conditions and more accurately represents the behavior of thin viscous layers, see Refs [ 1 4 , 14] for further details. Other investigators concerned with Euler, thin-layer Navier-Stokes or full Navier-Stokes equations have also presented algorithms that are designed to relieve the "artificial-viscosity" difficulties associated with central-difference discretizations. These are "upwind" approximations that are associated with the movement of the physical forward and backward moving waves, i.e. flux-difference splitting, or with the movement of the discrete forward and backward moving particles, i.e. flux-vector splitting. In previous PNS or RNS investigations, the choice of pressure/convective upwinding has resulted from characteristic or stability considerations. However, these results should also be derivable from a form of flux-vector splitting, since the characteristic analysis is simply a reflection of the eigenvalue or flux behavior. Therefore, in the present note, the pressure-velocity decomposition, that has been applied for the initial value PNS or boundary value RNS methodology, is re-examined from the point of view of flux vector splitting. The results represent a variation of the methods previously presented by Steger-Warming [16] and Van Leer [17].

2. ANALYSIS For the purpose of the present analysis we consider the 1-D Euler system given by: Q,+{E(Q)}x=Q,+{E+(Q)}~+{E

(Q)}x=0,

(la)

where Q is the solution vector Q = [p, pu, pe] v, e is the specific internal energy, p the density and u the velocity. E(Q)= E+(Q)+ E-(Q) and is such that all eigenvalues 2 + of dE+/dQ, dE /dQ are non-negative and non-positive, respectively. In order to reflect the initial value characteristic behavior for supersonic flows, we also require that E+(Q)= E(Q), E-(Q)= 0 for M = u/a >~ 1, where a is the sound speed (a2= 7P/P for a perfect gas), 7 is the ratio of specific heats and p is the pressure (p = (7 - 1)pc). Other desirable properties [17] are that E+-(Q), dE++-(Q)/dQand 2 + be continuous for all M, or for all (U/Urer)for incompressible flow, and that dE+-/dQ have at least one zero eigenvalue for IM[ < 1. This latter condition allows for a steady state shock structure with a maximum of two interior zones [17] and therefore results in more accurate shock capturing. In the present analysis, we have abandoned both the uniqueness condition [17] that E(Q) is the lowest degree polynomial in M and the symmetry condition for M ~ - M [17]. The objective herein is to identify the flux-splitting as strongly as possible with the respective convective (axial) and acoustic disturbances in order to maximize the "marching" or relaxation bias of the RNS formulation. This is in keeping with the methodology of interacting boundary layer calculations, in which the axial pressure gradient provides the only boundary value contribution for attached flows. This is also true for the Euler system, where the acoustic upstream influence is also pressure based. For regions of reversed flow, the convective fluxes must be appropriately upwinded in order to reflect the direction of the particle flow. The specification of the acoustic and mass flux-splitting is non-unique [17] and the form presented herein is particularly suited to relaxation or marching procedures with confined regions of reversed flow. In addition, at separation or reattachment points, this form of upwinding provides less numerical damping than previous upwind methods [16, 17] that are Euler based. Two new flux vector splittings that satisfy many of the properties on E +(Q) delineated previously are introduced here. In one model, the continuity condition on E +(Q) at separation, reattachment or stagnation points is no longer satisfied. For a second more satisfactory model, all of the properties including continuity are satisfied; however, the additional requirement that E-(Q) = E(Q), E+(Q) = 0 for M ~< - 1 is not enforced. Since we are not concerned here with reversed supersonic regions, this additional condition does not appear to be a severe limitation. Moreover, it is possible to "blend" the two forms of flux splitting in order to move smoothly from one to the other and thereby satisfy all conditions. This would lead to a "central" difference approximation over a limited range where - liftl[ <~ M ~< 0 or - [ ~ [ ~< u <~O. This is important primarily in the vicinity of stagnation points.

Technical Note

487

In order to evaluate the pressure/convective form of flux-vector splitting, we consider the simplified, partially conservative, Euler system [1]. The differential system is given by:

Q, + (dE+/dO)O; + (dE-/dO)O+~ = 0

(lb)

p, + (puL = 0

(lc)

( p u ) , - uZpx + 2u(pu)x + co(7 - 1)(pe)x + (1 - oJ)( 7 - 1)(pe)x = 0

(ld)

up (pe), + u(pe)x + e(pu)x - eup~ + p- (pu)~ - - - P x = 0 P P

(le)

or in scalar form

where co = co (M). In the RNS methodology for attached flows, all spatial gradients, save the pressure gradient (1 - c o ) ( 7 - 1)(pe)x term in eqn (ld), are associated with E+(Q). Therefore: dE+(Q) 0 1 - - - - u 2 2u

dQ

- uh

0 co(7-1),

h

dE-(Q) d~-

0 0

0 0

u

0 (1-o9)(7-

0

1) (2)

0

where h = e + p/p and a 2 = (7 - - 1)h for a perfect gas. The values of E+(Q), E - ( Q ) can be approximated by the following generalized expressions:

E+(Q)=

E-(Q)=

[

pu, p u ~ + ( 7 - 1 )

O, p u r l + ( 7 - 1 )

[

cope-

r

]

ped~ , p a e +

do

(1-co)pe+

rx,

pd¢+

do

ped~ , pile+

do

hd~

(3a)

do

do

pd¢-

do

hd~

(3b)

where 8 = (u, 0)maxand fi = - ( - u, 0)max. Note that the model equations are not in full conservation form and this leads to the integral flux contribution in E+-(Q). Also note that the contribution in E+-(Q) from the pressure gradient Px or (7 - 1)(pe)x in eqn (ld) is not given by pu 2 + cop and (1 - co)p, respectively. The parameter co is functionally dependent on M and has a discontinuity in slope at M = 1. However, with the definitions eqn (3) for E+-(Q), cox contributions will not appear in dE+-(Q)/dQ. These terms are not present in dE(Q)/dQ. Therefore, for the RNS splitting given by eqns (2), (3), the continuity of dE +-(Q)/dQ is maintained at sonic points to ensure smooth shock capturing. The present form of E+-(Q) leads directly to the RNS pressure decomposition that has been applied in earlier analysis. In particular, the discretization for the pressure gradient terms in eqn (ld) are given by:

[

COiPi--COi IP i l - -

L ,]/, p d

i

Xi--Xi

I)

(4a)

I

and (1

-coi+l)Pi+l-(1 - c o i ) P i -

pd~

(xi+l-xi)

i

for the positive and negative fluxes, respectively; coiP~ are the values at the mesh point x,. Since

~w'i p d ¢ ~(coi--coi_l) ~ i-I

where p ~[p~_ i, Pi], we obtain for eqn (4a): co,

,_,. (p~-p~ ,)

for

/~=

~+ ,_,

(4b)

488

Technical Note

and similarly for eqn (4b):

1

o~;+~o,+~. ( P i + l - P i )

•

2

for

/~=/p,+p;+t/.

,

\

2

/

These are precisely the discrete a p p r o x i m a t i o n s specified in previous R N S c o m p u t a t i o n s [1--4, 14, 15]. The eigenvalue (2 +--) analysis for the matrices eqn (2) leads to the following for u >f 0; 2 ~+2,3= (u + a ~ 1'2, u - ao9 j2, u)

(5a)

all of which are non-negative for u ~> 0, and:

2i2.3 = (0, 0, O)

(5b)

all of which are zero; therefore, in regions o f attached flow u p s t r e a m influence is of purely diffusive character. These results require that the function ~o is given by 0 ~< ~o ~< com = M 2. F o r ~o = c%:

2~2,3 = (2u, 0, u)

(5c)

and therefore dE+(Q)/dQ has one zero eigenvalue for all }MI < 1. F o r separated flows (u < 0), one choice of flux splitting is to reverse the roles of all fluxes so that E+-(Q)~ET-(Q) and hence 2 +-~ 2 ;. As noted previously, this satisfies all of the positivity (negativity) conditions. Although, the eigenvalues o f OE+-(Q)/dQ are continuous, b o t h E+(Q) and dE+(Q)/dQ are discontinuous when u changes sign. Although this splitting has been applied successfully for several time dependent c o m p u t a t i o n s [14], it is sensitive to the temporal increment and marginally stable for N e w t o n linearized steady state relaxation. An alternate and more desirable f o r m for m a r c h i n g m e t h o d s is obtained with the fluxes eqns (2), (3), but with ~o = 0 for all u ~< 0. The flux expressions eqn (3) for u < 0 become:

E+(Q)=

pu, O,

h d~

(6a)

do

[

O, p u 2 + ( 7 - 1 ) p e , p u e +

E (Q)=

;o i: ] pd~-

hd~

(6b)

so that for u < 0 dE+(Q) dQ

0

-- 00

1 0

0

00

,

dE-(Q) - -dQ

0

--

h

/22

0 2u

0 7 - 1

-- uh

0

u

-

(6c)

In this f o r m E+(Q) and dE+-(Q)/dQ are continuous at u = 0 , as well as at M = 1, as are the eigenvalues 2 ±, which for u < 0 are all zero or non-positive, respectively: )-,~2.3 = (0, 0, 0)

(7a)

Zi.3.3 = (0, 2u, u).

(7b)

Identical results for the eigenvalues 2 + are obtained with eqn (1) given in full conservation form in terms o f entropy, or in terms o f total energy (er = e + u 2/2), if the upx work term in the resulting energy equation is split in the same m a n n e r as is p.,. in the m o m e n t u m equation. If the up~ term in the total energy equation is unsplit, i.e. constant stagnation enthalpy or fully b a c k w a r d differenced, an alternate form of flux splitting results and the eigenvalues 2 ± for u >/0 become: ~,.,3 =

]

u , ~ [2 + (7 -- 1)(1 -- co)] +__a({(;~ -- 1)M(1 - e))} 2 + 4o)) 1/2

(8a)

and 2i~2, 3 = [0, (1 -- co)u(T -- 1), O]

(Sb)

Technical Note

489

where again 0 ~
2u).

(9a) (9b)

When the total energy is one of the dependent variables, an additional condition is required to maintain the full upwind form. Since p = (7 - l)pe r -[(7 - 1)/2] pu 2, the (pU2)x flUX arising from Px also changes sign with u; therefore, these terms which appear in the momentum and energy equations must also be given in ~, fi form. The results of eqn (9) are then recovered. This additional condition does not arise when the temperature or the specific internal energy is one of the dependent variables; the results of eqn (7) are then recovered.

3. S T A B I L I T Y

It can be shown [18] that both the full shift and co = 0 shift flux-splittings are unconditionally stable for steady state line relaxation. These methods are sensitive, however, to the assumed initial conditions for pressure, etc. that are required for steady-state calculations, for terms entering from E (Q). It has been shown [18] that this sensitivity can be controlled by (i) an improved initial guess, e.g. several inviscid relaxation steps, (ii) local non-linear iteration at each axial location or (iii) underrelaxation of the velocities in separated regions, by the inclusion of temporal damping, i.e. finite At, in the early stages of the relaxation process. After a few steps e.g. 5-20, full relaxation without local iteration can be restored. It can be shown that the discretization for the elliptic pressure term arising from E-(Q) corresponds to complete overrelaxation. While this is desirable close to convergence, it is undesirable in the initial stages of the relaxation process. This difficulty can be tempered by a pressure correction term developed by Israeli [19]. The correction, which vanishes in the steady state introduces pressure underrelaxation and also residual smoothing for multi-grid transfers. A combination of the pressure smoothing with full overrelaxation has provided a useful approach for the semi-coarsening or uni-directional multi-grid accelerator [15, 18]. 4. C O N C L U S I O N

Global pressure relaxation for a pressure split form of the Euler or reduced Navier-Stokes (RNS) equations has been shown to be equivalent to a form of flux-vector splitting that satisfies the major eigenvalue and continuity constraints on the fluxes and flux derivatives. The upwinding is applied only in the axial or "streamline" direction. In keeping with the asymptotic form of the RNS system, two-point or trapezoidal discretization is used for appropriate normal gradients in order to allow for the accurate evaluation of shear layers and a consistent specification of far field boundary conditions. If the pressure gradient parameter 09(M) is given by its maximum value in subsonic regions, one eigenvalue of 2 ÷ is always zero and therefore shock resolution is improved. For regions of reversed flow convective upwinding is combined with the condition co = 0. This ensures that the fluxes, flux derivatives and eigenvalues remain continuous throughout the flow. This form

490

Technical Note

o f flUX v e c t o r s p l i t t i n g is specifically d e s i g n e d to m a i n t a i n a bias in the d i r e c t i o n o f t h e c o n v e c t i v e fluxes a n d t h e r e f o r e a b a n d o n s the ( _ ) s y m m e t r y [17] o f the e a r l i e r f o r m s o f flux v e c t o r splitting. T h e u p w i n d i n g l e a d s to a r e l a x a t i o n m e t h o d t h a t is solely a c o u s t i c d r i v e n t h r o u g h o u t s u b s o n i c r e g i o n s , b u t also i n c l u d e s c o n v e c t i v e r e l a x a t i o n in r e g i o n s o f r e v e r s e d flow w i t h o u t i n t r o d u c i n g large a m o u n t s o f n u m e r i c a l d i s s i p a t i o n . Acknowledgements--This work was supported by the Air Force Office of Scientific Research under Contract F49620-85-C-0027 (Dr J. Wilson, Technical Monitor) and under an AFOSR-University Resident Research Agreement, while the author was in residence at the Air Force Wright Aeronautical Research Laboratories, WPAFB. The author wishes to thank Dr Joseph Shang of the Flight Dynamics Laboratory, AFWAL for his considerable support and discussions during this period. REFERENCES 1. S.G. Rubin, A Review of Marching ProceduresJor PNS Equations Numerical and Physical Aspects o/Aerodynamic [.'.lows (Edited by T. Cebeci), pp. 171-186. Springer, Berlin (1981). 2. S. G. Rubin and D. R. Reddy, Global solution procedures for laminar flow with strong pressure interaction and separation. Comput. Fluids 11, 281-306 (1983). 3. S. G. Rubin and P. K. Khosla, Global relaxation procedures for a reduced form of the Navier Stokes. Lecture Notes in Physics, Vol. 218, pp, 62-71. Springer, Berlin (1984). 4. P. K. Khosla and H. T. Lai, Global relaxation procedures for compressible solutions of the steady-state Euler equations. Comput. Fluids 15, 215-228 (1987). 5. S. G. Rubin and P. K. Khosla, A composite solution procedure for the incompressible Navier-Stokes equations. Lecture Notes in Physics, Vol. 170. Springer, Berlin (1982). 6. P. K. Khosla and S. G. Rubin, A composite solution procedure for the compressible Navier-Stokes equations. A1AA J. 21, 1546-q551. 7. R. C. Swanson, S. G. Rubin and P. K. Khosla, Calculation of afterbody flows with a composite velocity formulation. AIAA Paper No. 83-1736 (1983). 8. R. Gordnier and S. G. Rubin, Transonic flow solutions using a composite velocity procedures for potential, Euler and RNS equations. Comput. Fluids. Accepted for publication. 9. S. Rudman and S. G. Rubin, Hypersonic viscous flow over slender bodies having sharp leading edge. AIAA J, 6, 1883-90 (1968). 10. T. C. Lin and S. G. Rubin, Viscous flow over a cone at moderate incidence I, II. Comput. Fluids 1, 37-57 (1973): J. Fluid Mech. 59, 593-620 (1973). 11. L. B. Schiff and J. L. Steger, Numerical simulation of steady supersonic viscous flow. AIAA J. 18, (1980). 12. Y. C. Vigneron, J. V. Rakich and J. C. Tannehill, Calculation of supersonic viscous flow over delta wings with sharp leading edges. AIAA Paper No. 78-1137 (1978). 13. R.T. Davis, M. Barnett and J. V. Rakich, The calculation of supersonic viscous flows using the PNS equations. Comput. Fluids 14, 197-224 (1986). 14. S. V. Ramakrishnan and S. G. Rubin, Time consistent pressure relaxation procedure for compressible RNS equations. AIAA J. 25, 905-913 (1987). 15. A. Himansu and S. G. Rubin, Multigrid acceleration of a relaxation procedure for RNS equations. AIAA Paper 87-1145C8 (1987). 16. J. L. Steger and R. F. Warming, Flux vector splitting of the inviscid gas dynamic equations, with application to finite-difference methods. J. Comp. Phys. 40, 263 293 (1981). 17. B. Van Leer, Flux-vector splitting for the Euler equations. Lecture Notes Phys. 170, 507-572 (1982). 18. S. G. Rubin and A. Himansu, Rapid convergence of laminar and turbulent strong interaction flows. Comput. Fluids Accepted for publication. 19. M. Israeli and M. Rosenfield, Numerical solution Of incompressible flows by a marching multigrid nonlinear method. AIAA Paper No. 85-1000 (1985).