The calculation of supersonic viscous flows using the parabolized Navier-Stokes equations

Computers & Fluids Vol. 14, No. 3, pp. 197-224, 1986 Printed in Great Britain

THE

0045-7930/86 $3.00 + 0.00 Pergamon Journals Lid

CALCULATION OF SUPERSONIC VISCOUS FLOWS USING THE PARABOLIZED NAVIER-STOKES EQUATIONS~" R. T. DAVIS and M. BARNETT University of Cincinnati, Cincinnati, OH 45221, U.S.A.

and J. V.

RAKICH~

N A S A Ames Research Center, Moffett Field,C A 94035, U.S.A.

(Received27 September 1984; in revisedform 21 May 1985)

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

Approximate forms of the full Navier-Stokes equations have been utilized for many years to analyze the flow of viscous fluids. In the early years of this century, Prandtl [ 1] developed boundary-layer theory, which has proven to be exceedingly robust, allowing the solution of many practical problems. Numerous other techniques have been developed to supplement the computationally costly full Navier-Stokes solution methods that have become relatively common since the advent of the high speed digital computer. One method which has received a great deal of attention over the last few years involves solving an approximate form of the l'qavier-Stokes equations which has become known as the "parabolic" or "parabolized" Navier-Stokes equations. This approximation falls somewhere between the full Navier-Stokes equations and the boundary-layer equations. The parabolic Navier-Stokes (PNS) equations can be viewed as a composite set of equations in the sense that they are valid in both the viscous and inviscid regions of the flowfield. All effects included in second-order boundary-layer theory are accounted for, in particular, the effect of boundary-layer displacement thickness growth upon the inviscid flowfield is included. This is crucial, as we shall see later, in permitting the PNS equations to be used to calculate flows where the viscous-inviscid interaction is strong. For such flows, ordinary boundary-layer theory or other procedures for weakly interacting flows are not appropriate. The boundary-layer equations can be recast into the so-called interacting boundary-layer equations and utilized to solve for flows with strong viscous-inviscid interaction. In this paper, we wish to restrict our discussion to the solution of the parabolic NavierStokes equations for supersonic flows. We will discuss the assumptions which lead to the PNS equations. As we shall see, the treatment of the streamwise pressure gradient term is of paramount importance in permitting the correct solution of both weakly interacting and strongly interacting flow problems. Characteristic and stability analyses are both undertaken to show how the streamwise pressure gradient term affects the nature of the governing equations. The results of these analyses coincide, both pointing out the need to somehow modify this term to permit the use of the parabolic Navier-Stokes equations to solve both weak and strong interaction problems. ~"This work was supported by NASA Grants NGT 36-004-800 and NCA2-OR 130-101. This paper is dedicated to the memory of the third author, John V. Rakich. C A F 14/3---B

197

198

R . T . DAWS et al.

If the parabolic Navier-Stokes equations are solved with the streamwise pressure gradient term fully included as a backward difference relative to the solution station, solutions exhibit what has become known as departure behavior, wherein the flow quantities branch exponentially from their expected values. The fundamental reason for this, as discussed more fully later, is that the equations are ill-posed for solution as an initial-value problem. The precise nature of departure solutions is determined in this study by considering the results of Lighthill's analysis of the sublayer equations [2] and the subsequently developed tripledeck theory. To obtain departure free solutions using the PNS equations for weakly interacting flows, one can drop a fraction of the streamwise pressure gradient where the flow is subsonic. For strongly interacting flows, a forward difference for a fraction of the streamwise pressure gradient is utilized in the present study. This satisfies the boundary-value nature of the problem and provides a mechanism for the upstream propagation of information, which is crucial to calculating such flows. The single sweep procedure for solving the parabolic Navier-Stokes equations for weakly interacting flows proposed by Vigneron, Rakich and Tannehill [3] is considered. The question of accuracy is discussed in detail and the conditions under which their approximation for the streamwise pressure gradient term is accurate are determined. The PNS equations are used in this study to solve the hypersonic viscous interaction problem and the results are shown to compare well with experiment throughout almost the entire range of the hypersonic interaction parameter, ~. In addition to the consideration of weakly interacting flows, the basis for a solution procedure utilizing the parabolic Navier-Stokes equations to solve for strongly interacting flows is developed. A new Alternating Direction Explicit (ADE) procedure is outlined which should permit efficient calculation of strongly interacting flows using a global iteration procedure for the PNS equations. This procedure has minimal computer storage requirements compared to full Navier-Stokes solvers since only the pressure is saved, and only in the subsonic portion of the boundary layer. 2. G O V E R N I N G

EQUATIONS

It will be convenient for the following discussion of the parabolic approximation to present the governing equations in their most general three-dimensional form. The NavierStokes equations for unsteady flow without body forces or external heat addition can be written in Cartesian conservation-law form as Ot+~x(E-Ev)+

(F-Fv)+

( G - Gv) = 0,

(2.1)

where Q r = {o, pu, or, pw, pet}.

(2.2)

Here O, u, v, w and et are respectively the nondimensional density, Cartesian velocity components and total internal energy (see [3]) as defined by U2 _[_ /)2 + W 2

et = e +

2

(2.3)

The subscript v is used to distinguish the viscous from the inviscid flux terms, all of which are given in [3-5]. The perfect gas equations of state for pressure p and temperature T are given by p = (3/ - l)pe

(2.4)

T = "},M2p/o,

(2.5)

and

Calculation of supersonic viscousflows

199

where 3, is the ratio of specific heats, and Moo is the free stream Math number. Finally, to close the system of equations, the viscosity coefficient u is obtained from the Sutherland law, I+C

t~ = ~ T+C

(2.6)

T 3/2.

A general coordinate transformation is introduced of the form = ~(x, y, z), n = n(x, y, z),

(2.7a-c)

~" = ~(x, y, z),

where ~, 7, ~"are boundary fitted coordinates which may vary for different body geometries. is the predominately streamwise coordinate, ~ the predominately normal coordinate and ~"the predominately crossflow coordinate. Following the analysis of Viviand [6], eqn (2.1) can be rewritten in terms of the new coordinates, while retaining conservation law form, as O

O

Ot ( J - ' Q ) = "~ { J - ' [ ~ x ( E - Ev) + ~r(F - F , ) + ~,(G - G0]} +

+ TIy(F- Fo) + ttz(G -

Go)]}+ ~O {J-Lff:,(E

O

~ {J-'[n,,(E

- E O + f y ( V - Fv) + ~z(G -

os

- E~)

Go)I},

(2.8)

where J is the Jacobian of the coordinate transformation defined by

r)

-'

J = a(x, y, z) = k0(~, ,1, g').l

"

(2.9)

In practice, the Jacobian and the derivatives of the coordinate transformation are evaluated numerically; see Ref. [3]. 2.1. P a r a b o l i c N a v i e r - S t o k e s e q u a t i o n s

For flows with slow variation in the streamwise, ~, direction, it is possible to neglect the derivatives of the viscous terms with tittle loss of accuracy. This approximation is consistent at high Reynolds number with the boundary-layer approximations for viscous terms as long as the coordinate system is aligned properly. If, in addition, the flow is assumed steady, one obtains what are commonly known as the parabolic Navier-Stokes equations. Actually, these equations are not truly parabolic since, as will be discussed in the following section; additional restriction must by made in order for the equations to be formally parabolic or of the initial-value type in the mathematical sense. With these assumptions, eqn (2.8) reduces to

a~ + aP

0"-~

aO

~ + " ~ - 0,

(2.10)

where ft. = J - ' [ ~ x E + ~ r F + ~zG], f f = J-~[~lx(E - Ev) + 7 1 y ( F - Fv) + ~lz(G - Gv)], d = J - t [ ~ x ( E - Ev) + ~ y ( F - Fv) + ~z(G - Gv)].

(2.1 1a-c)

Further, the derivatives appearing in the viscous terms Ev, Fv and Go are transformed with the expressions

200

R. T. DAVIS et

0

al.

0

0

r o- , 0

0

0

+ r,o0

a

0

(2.12a-c)

2.2. Thin-layer approximation

For boundary-layer type flows, where gradients in the ~/direction are dominant, it can be shown that the derivatives in the ~"direction of the viscous terms are negligible. In this case one obtains ft. = J-~[~xE + ~yF + ~zG], f f = J-~[n~(E - EO + n r ( F - F~) + nz(G - Go)], = J - ' [ ~ x E + ~yF + ~zal,

(2.13a-c)

and the following relationships for first derivatives can be used: 0

0

0

0

~ y ~ l"lyo'-~ , a c9 az = rlz &7 "

(2.14a-c)

The thin-layer approximation is used for the majority of applications of the PNS (Parabolic Navier-Stokes) technique because it simplifies the coding and, therefore, runs more efficiently. The more complete form of the PNS equations (2.10)-(2.12) contain crossderivatives, which are difficult to evaluate using ADI (Alternating Direction Implicit) methods and are therefore usually neglected. However, the second derivatives with respect to ~ can be included without difficulty. For this reason the ~'~"derivatives are often included in a problem, resulting in an approximation between thin-layer and full PNS. However, from an order of magnitude viewpoint, there is no advantage in retaining these extra terms without retaining the cross-derivatives. 3. ANALYSIS OF THE "PARABOLIC" EQUATIONS For an inviscid flow, it is clear that if the flow is supersonic the problem is well posed as an initial-value problem. However, in a viscous flow problem, even if the inviscid part of the flow is supersonic, the presence of a boundary layer on a solid surface will produce a subsonic region which provides a path for upstream influence. Attempts to solve these equations with a marching technique will, therefore, typically lead to what has become known as a departure solution, wherein the flow variables branch exponentially from the expected solution. On the other hand, if the pressure on the surface is specified, it is known from boundary-layer theory that the equations can be marched from an initial data surface as long as the flow does not separate. In boundary-layer theory, it is the pressure interaction between the outer inviscid flow and the viscous boundary-layer flow (displacement thickness) which produces the upstream influence corresponding to the subsonic layer near the wall in the full equations. A proper treatment of the interaction phenomenon will allow one to march boundary-layer equations into reversed flow regions [7, 8]. In previous attempts to march approximations to the Navier-Stokes equations containing the proper terms to model simultaneously the inviscid and viscous regions of flow, various approximations have been applied to the pressure terms in the G-momentum equation in

Calculationof supersonicviscousflows

201

order to eliminate the departure solutions. Davis [9], in early work on supersonic blunt body flows utilized the pressure behind the bow shock for viscous shock-layer flows and iterated globally to obtain converged solutions. This method was later improved by Srivastava, Werle and Davis [10]. Lin and Rubin [1 1], following a suggestion by Cheng [l 2], attempted to remove departure solutions by placing the pressure gradient term in the G-momentum equation on the right side of the equation and evaluating it with a backward difference. This approach was also followed by Lubard and Helliwell [13] and Schiffand Steger [14]. A slightly different approach was taken by Vigneron, Rakich and Tannehill [3], and that method is reviewed here. It is convenient for this analysis to consider the governing equations in two-dimensional Cartesian coordinates. The question considered by Vigneron et aL [3] is what fraction, o~, of the pressure gradient can be retained in the G-momentum equation without loss of a formally parabolic set of equations. Thus, eqn (2. l 0) is written as OE_____~*+ O(F - Fv) _ Ox Oy

OP Ox '

(3.1)

where g *T =

p/,/, OIJ2 ~- 60p, ,OUl),

flU, at/V, tOY2 "t- p,

F T =

"/'

"Y

p -I- ~ (1/2 nt_ 1)2) U ,

p d- ~ (Z/2 "t- 1)2) 1) ,

Ree O, uy, -~ vr, uu, + ~ vvy + (~, _ ~-M2pr

,

p T = { 0 , (1 - ~o)p, O, 0 } .

(3.2a-d)

Here, Pr is the Prandtl number and Re is the Reynolds number. We wish to determine under what conditions eqns (3. l) can be marched as a well-posed initial value problem, that is, when the characteristics are real and positive. The result of this analysis will determine what fraction 00of the pressure gradient can be included without ill-posed (elliptic) behavior. Note that the source term OP/Ox on the right side ofeqn (3.1) is considered a known function, and does not affect the characteristics, but this term controls the nature of the numerical solutions which can be obtained for these equations. For the present analysis, it is simpler to write eqn (3.1) in terms of primitive variables, U r -- {p, u, v, p}. This can be done since the characteristics are independent of the particular variables used. Thus, we write eqn (3. l) in the form ou B ou o2u A -~x + -~y = C--oy2 ,

(3.3)

where A=

I

0u 0 -a2u

pu p 0 0

0 o

001 0 pu 0

o~ 0 1 0 ' c~u

,

l0 o

~'P 0 p2p---~

0

o-~l

B =

[i o -a2v

and

pv 0 0

0 pv 0

0 1 v

'

¢, = ,y + w(l - ~).

(3.4a-d)

202

R.T. DAVISel al.

3.1. lnviscid flow In the limit Re ---. oo, C ---* 0 outside the boundary layer and these equations are wellposed as an initial-value problem in the x-direction if the eigenvalues o f (A-~B) are real. The eigenvalues ~, which satisfy the relation

IB

-

x,41 = 0

(3.5)

are X =

(1 + ¢)M~My +_-{(MxMv)2(¢ - l) 2 - 4[w(l - M 2) - ¢M2]} '/2 2(4~M 2 - 00)

(3.6)

The eigenvalues are real if w(1 - M 2 ) - ¢ M

2 <0

(3.7)

or

~o <

~m~ ( ' y - 1)M2x- M 2 +

(3.8)

1

Thus, the inviscid equations (Euler equations) have hyperbolic-parabolic characteristics if co satisfies eqn (3.8). Note, for later reference, that i f M v = 0, i.e. v = 0, then

~m~

o~ < - (3, - 1)M2x + 1

(3.9)

This inviscid limit was incorrectly given in Ref. [3] due to an algebraic oversight but the basic dependence o f o~ on M 2 was correctly stated. 3.2. Viscous flow W h e n Re is finite, the characteristics o f eqn (3.3) are governed by the highest-order derivatives. Thus, in this case, we require that the eigenvalues o f ( A - t C ) be real and positive. This leads to the characteristic polynomial; see Vigneron, et aL [3],

~ - (~b+Prr] "Y]Mx]+ 2] ")'M2x~ A(A-~){A2(¢M2x-¢O)+ A[Pr Pr )=0'

(3.10)

where A = Re

puh /z

(3.11)

The roots of(3.10) are A=0,

A = 4, A = [M~(1 + Pr)/Pr) - ~ ] ___ ifD 2('),M 2 - a~t) where D = 0~27: + k - - - - ~ r J [ L t = 1 +(5,1

7=-+(7P

,',~, Pr

lJ +

(3.12a-c)

}

(5'M2x)2 '

(3.13)

1)M 2,

(3.1 4)

l ) M 2.

(3.15)

Calculationof supersonicviscousflows

203

The roots of the quadratic (3.10) are difficult to evaluate, but the expression simplifies for a Prandtl number of unity. In this case one obtains

{ , A=

3"M~

(3.16)

3"M2x -- w t "

The last root gives a condition on w such that A > 0. The result is

3"M w < 1 + ( 3 ' - 1)M2

(3.17)

which is the same as the inviscid result for the case where My = 0. The significant point is that for both the inviscid and viscous analyses, w --~ M 2 results in real eigenvalues. This, plus the condition u > 0, insures that the equations are parabolic. The fact that the results (3.8) and (3.17) are almost identical is not surprising. Equation (3.17) reflects the fact that upstream influence occurs principally through the inviscid terms. The viscous terms have already been modified to remove their contribution to ellipticity. Therefore, in analyzing a set of parabolized equations it is only necessary to analyze the inviscid terms in order to detect the imaginary eigenvalues indicating upstream influence, a fact which makes the analysis much simpler. 4. FINITE DIFFERENCE METHOD In computations of steady boundary-layer type flows, it has long been a standard procedure to use implicit finite difference methods to eliminate severe restrictions placed on the marching step size by explicit methods. For the same reasons, we adopt the Euler implicit method in conjunction with the Beam-Warming [15] factored algorithm in delta form. Similar algorithms have been used by many other authors; see [16, 17], for example. The PNS equations can be written in delta form as

o

at +

A'e +

oo,

A'G =

-[ °e+ a°l' Lan

a .l'

(4.1)

where i is the grid index defined by = iA~

(4.2)

and the A operator is defined by Ai( 4.1. F a c t o r e d

)=(

)i+1_(

)i.

(4.3)

implicit algorithm

Application of the factored implicit method requires the calculation of Jacobian matrices of ff and G with respect to the primary dependent variables. For the time dependent form of the Navier-Stokes equations, the primary variables are obviously Q; see eqn (2.2). However, the PNS equations are complicated by the fact that the /~ vector appears in the G-derivative. The Jacobians of P and G with respect to E are more complicated than those with respect to (2. For this and other reasons, Q is chosen as the primary dependent variable for the PNS equations which results in

OEaQ

r_oo,ae

a ocq N_

[ ap aG-li

(4.4)

Finally, in order to make the equations formally parabolic, as described in the previous section, part of the pressure is removed from/~, giving

204

R. T. DAVISet al.

o--0 o-g +

e=-

(4.5)

g + orj

where (4.6)

E* =ff,-P

and (4.7)

f i = J - l ( ~ x P t + ~yP2 + ~zP3),

where

PT =

{o, (1 - w)p, o, o, o},

p r = {0, 0, (1 -- w)p, 0, 0}, P3r = {0, 0, 0, (1 - w)p, 0}.

(4.8a-c)

To first order accuracy]" in A~ the implicit Beam-Warming factoring procedure now gives

joe*

o OCl(Oe,) 'FOe, +

_

o oen ,_ A[OE*\-

_

OF -

_

OG _ i

(4.9)

We note that/~*, F, G depend on ( through the metrics (x, (y, (z, and this gives rise to the term on the right side of eqn (4.9); thus the subscript (~ indicates that the derivative is evaluated with (~ held constant. However, these terms vanish if the marching surface is chosen to be the planes of x = constant, which is the case in the present work. For more general marching surfaces these terms should be retained. Careful inspection reveals that the last right side term is higher order and hence is not needed from an accuracy standpoint for the first-order accurate scheme. Equation (4.9) is now in a form that allows solution by means of sequential onedimensional sweeps in ~ and ~"independently. This is done in four steps:

1:

[og* L OQ + zx~~}-~ J~QI = RHS(4.9),

2:

A02 = ~

A0,,

3:

ro£, LaQ

o oel

4:

(~i+1

o oC*n

0/~*

-

(4.11)

+ A( ~-~ ~--~]AQ = A(~2, =

(4.10)

(4.12)

O i + AiO.

(4.13)

Note that there is no need to compute the inverse ofd/~*/0(~. Second-order accurate threepoint central differences are used to approximate the n and ~"derivatives, resulting in block tridiagonal matrices. These may be solved with a block tridiagonal solver. The one used in the present study is due to Steger [ 18]. 4.2. B o u n d a r y conditions

The physical boundary conditions on the body are U=I)=W=0, T = Tw

or

= 0,

(4.14a-d)

W

]- The methodcan be made second-orderaccuratein A~but that would not significantlyaffectthe presentstudy.

Calculation of supersonic viscous flows

205

where w denotes the wall and n is in the outward normal direction from the wall. The manner in which these boundary conditions may be applied is discussed in [19]. The outer boundary is often taken as the shock wave emanating from the leading edge, and the Rankine-Hugoniot equations are satisfied there. In early applications [3, 5, 14] the shock boundary conditions were applied in an explicit manner. It has been shown by Barnett, Davis and Rakich [ 19] how the shock conditions can be imposed in a fully implicit manner. This implicit approach improves the efficiency of the PNS code. The shock boundary conditions are first-order accurate in the streamwise direction, which is compatible with the present difference scheme. Extension of the boundary conditions to second-order accuracy would be necessary if the present method is to be made fully second-order accurate. The details of the boundary condition calculation are beyond the scope of the present paper and are fully described in Ref. [19] for two-dimensional flows. In the present paper, attention is focused on the fundamental stability and/or branching behavior of the finite difference solution of the PNS equations and the solution of problems for which strong viscous-inviscid interaction is present.

4.3. Linear stability analysis In this section, the numerical stability of the two-dimensional PNS equations is examined. Consider the flow to be in a Cartesian (x, y) coordinate system. We consider, in particular, the effect of the Px term,t which was removed from the left side of the PNS equations in order to make them formally parabolic. The governing equations are

aE* aQ aP aQ _ aFo 02Q + OQ ax OQ ax aQy ay 2 '

(4.15)

where P is given by eqn (3.2d). Differencing with the Euler implicit method, and using backward difference for Px, gives

OE* + 2Ax OF~]t.)i+l OQ Ay ~ OQy]'~j

Ax OFv t,oi+l

,oi+tx

(OP

a

OE*\ i OP i l - ~ )eJ - - ~ Q~- = o,

Ay z OQv '~J+' + ~J-" + - ~

where

x = iAx,

y = jay.

(4.16a-b)

Considering solutions of the form Qj = riexp(~L-lkjAy) one obtains the amplification matrix for the system. Vigneron et al. [4] have obtained a polynomial of degree 8 for the eigenvalues of the amplification matrix. This polynomial is Xa(X -

1)B(X,

X)J(X, X, Mx, oo)

=

O,

(4.17)

where B=(1

j = {X2[~o(3,- -21)- 2

4 - ~ X ) X - 1,

2 - 1)]_ ( l - °°)(T2- l)} X]+X[ 3'+ 1-2o:(3, Pr(3,- l)M#J

2

+

X=

4gAx

sin2(KAy~

pu ReAy 2

\2j

1

(s- 1)M~

"r -

2]

2 ~' (4.18a-c)

l" Note that we are using the symbol P as in eqn (2.2d) and the (1 - o~)pfactor is hidden in the notation. The wp terms are in the terms shown in eqns (3.2a-c).

206

R.T. DAVISet al.

Equation (4.19) has five obvious roots which are = 0

(a

triple root),

X=I, X -- 1/(1 + ~X),

(4.19a-c)

which are all less than or equal to unity. The remaining eigenvalues are roots of J which are difficult to find analytically. These were numerically analyzed in Ref. [4] and the results are shown in Figs. 3 and 4 of that reference. The following conclusions follow from these results: (l) If ~0 = "rM~/[(3, - 1)M 2 + l ] for Mx < l, and ¢o = 1 for Mx >= 1, and the source term Px is included as a backward difference, then there is a domain for small X and small Mx, where the eigenvalues of the amplification matrix are greater than unity. This means that the finite difference method is unstable for sufficiently small step size, a result which coincides precisely with the analysis of Lubard and Helliwell [13] and the characteristic analysis presented previously. (2) If the source term Px is dropped and the ~0 factor is retained, the eigenvalues are all less than unity. In this case the finite difference method is stable for all Ax, no matter how small, again in agreement with the characteristic analysis presented previously. Thus, while the differential equations can be made formally parabolic by placing Px on the fight side, the finite difference equations may be unstable if this term is approximated by a backward difference. In practice, one can sometimes include Px using a backward difference without exciting the instability, but it is assured that the solution will fail if the step size is reduced. The safest procedure is to drop P~ entirely if one wishes to use a one pass method for weakly interacting flows. As will be shown in the section on Triple-Deck Analysis this is, in fact, a good approximation for large Mach numbers and for cold wall temperatures as long as the flow is of a weakly interacting type. For strongly interacting flows, one needs to retain the Px term in order to achieve accuracy. This term is in fact the term which allows for upstream propagation. It should once again be pointed out that in the present analysis P is given by (3.2d) and only includes the (1 - o~)p term appearing in the x m o m e n t u m equation. The stability and characteristic analysis are forcing us to the same conclusion--that Px should not be backward differenced and that this term is the one which accounts for upstream propagation. This means that a difference scheme for P~ should be of the forward type to allow for the accommodation of a downstream boundary condition. A global iteration scheme is, therefore, required to include this term. This conclusion is in exact agreement with triple deck theory; see Stewartson [20]. The inclusion of the P:, term will be discussed in Section 6, Globally Iterated Solutions. 5. DEPARTURE SOLUTIONS At this point we wish to investigate in more detail the effect of the streamwise pressure gradient term in the parabolized Navier-Stokes equations. We will deal only with the case of two dimensions and supersonic outer flow conditions. The phenomenon of upstream propagation in a boundary layer with a supersonic outer flow was observed by Ackeret et al. [21 ] and Chapman et al. [22] among others who showed experimentally that in a supersonic flow the boundary layer can, during a shock/boundarylayer interaction, separate upstream of the point of intersection of the shock with the boundary-layer. This occurs despite the fact that the supersonic outer flow is governed by hyperbolic equations and the boundary-layer flow is governed by parabolic equations. The implication is that a mechanism exists for upstream propagation of information within a boundarylayer which is supersonic except in a thin region near the wall. The fact that upstream propagation exists is not surprising, and it can be explained by analyzing the interaction between the boundary-layer displacement thickness and the outer inviscid flow. The conclusions reached from that type of analysis can be extended to the parabolized Navier-

Calculation of supersonic viscous flows

207

Stokes equations since, in a composite sense, the parabolized equations represent a boundary layer interacting with an outer inviscid flow. It will be shown in the present study that the growing exponential behavior displayed by the parabolized Navier-Stokes solutions (departure solutions) is identical to the behavior discovered by Lighthill [2], who analyzed the free interaction which exists between a boundary layer and a supersonic outer inviscid flow. This analysis of Lighthill led to the development of triple-deck theory (see Stewartson [20]) since he found that the interaction causes a sublayer to develop at the wall which is the primary contributor to the interaction with the outer inviscid flow. If one constructs a mesh which is capable of resolving the triple-deck structure, in particular the lower deck region, then the parabolized Navier-Stokes equations contain the mechanism which allows for the interaction phenomenon and, therefore, growing exponential solutions will arise. This will be demonstrated with numerical calculations wherein the pressure gradient term will be fully retained and departure solutions initiated at a chosen streamwise location. The results will then be compared to the predictions obtained from triple-deck theory for the initial departure. The analysis which follows is for a fiat plate immersed in a supersonic flow approaching along the x-axis of the plate. The y-axis is normal to the plate. We will not develop tripledeck theory here, one is referred to Stewartson [20] for the development in detail. We will only quote the lower deck (sublayer) equations along with the interaction condition and boundary conditions. These equations are sufficient for the analysis of the sublayer which will allow us to develop the expressions needed for comparison with parabolized NavierStokes solutions. The equations for the lower deck (sublayer) are as follows (see Stewartson [20] for details): Continuity: Ou

Ov

~x+~y=O.

(5.1)

x-Momentum: 0U

0U

de

u ~ x + v Oy

02//

d x + @2 •

(5.2)

Interaction (for supersonic flow): P -

dA dx"

(5.3)

Boundary and matching conditions: u= v= 0

at

y = 0,

(5.4a-b)

u "~ y

as

x ~ -~,

(5.5)

u "-, y + A ( x )

as

y ---, ~ .

(5.6)

In the above equations u, v and P have their normal definitions and A is the sublayer thickness. The Mach number, Reynolds number and wall temperature dependence have been scaled out through nondimensionalization of the equations. The relationships between dimensional and dimensionless variables are given by U•

- -

u* I)*

u* X* -- X~

= Re-I/8~k

l/4CI/S(M

-- Re-3/8X3/4c3/g(m2

=

2

_

1)-l/8(Tw/Too)l/2u,

-

l )l/S(Tw/T~)l/2v,

Re-3/sx-5/4C3/S(M2

-

1)-3/8(Tw/To)3/2X,

208

R.T. DAVISel aL __Y* = R e - 5 / s X - 3 / 4 C S / S ( M 2

_ 1)-l/8(Tw/T,~)3/2y,

x* P* -- P * - R e - i / 4 X ' / z C ' / 4 ( M 2 -- 1)-'/4P,

(5.7a-e)

p~U*~2

where the * superscripts denote dimensional variables, whereas variables without superscripts are dimensionless. The Reynolds number Re is defined by n--

1"%l:.

_

_

~Dco ~ ° o ' A P

u*

'

(5.8)

where x~ is the location of the initiation of the interaction measured from the leading edge of the plate. The quantity X ~ 0.332 is the Blasius constant, and Cis the Chapman-Rubesin constant given by # * T ~ / ( t z ~ T*). The most interesting property of the above equations is that even though the lower deck equations with prescribed pressure gradient are parabolic and the outer flow is supersonic, the coupled set of equations admit solutions which show growing exponential behavior in the streamwise direction. This fact was first discovered by Lighthill [2] and his results will be summarized here. One is referred to Stewartson [20] for details. Within the framework of triple-deck theory, there exists a solution to the lower deck equations for the flat plate of the form u = y-p=

alekxf'(y),

v = a,keUf(y),

A = --a,ekX/k,

a l e kx,

(5.9a-d)

where am is a small arbitrary constant. Substituting into the momentum equation (continuity and interaction have already been satisfied through definition of the variables) we obtain, after assuming at is small and linearizing, the following equation: f"

- kyf' + kf = -k.

(5.10)

Notice that the upstream boundary condition is satisfied by the assumed form of the solution. The wall and edge conditions are given by f(0) = f ' ( 0 ) = 0

and

f ' ( ~ ) = 1/k.

(5.11a-c)

The parameter al no longer appears in the equations, indicating an eigenvalue type problem, and it can be shown that a solution exists to the equations for only one value of k (in addition to the undisturbed case) which turns out to be k -- 0.8272. The solution f o r f i s given by an Airy function; see Stewartson [20]. The important result is that k is positive indicating growing exponential behavior of the solution. According to this theory the pressure and wall shear, for example, should behave like P = ale °'8272x

(5.12)

r = 1 -

(5.13)

and 1 . 2 0 9 a l e °8272x,

where r has been nondimensionalized by its starting point value. The constant at can be interpreted in a forward marching scheme as an error or small perturbation which initiates the departure. Its value could be positive or negative leading to the possibility of a compression or expansion type free interaction. Notice that the x variable in these expressions is the dimensionless one, whereas the dimensional x distance is defined by eqn (5.7c). Due to the Re -3/8 term in the expression

Calculation of supersonic viscousflows

209

we see that in physical terms the branching will be rapid and that mesh spacings must be scaled properly in order to capture the branch. The mathematical significance of the branching behavior is that the problem is improperly posed as an initial-value problem and that some downstream boundary condition is required. However, in many problems involving weak interactions we can suppress the departure solutions and produce accurate results. This is certainly not the case if one wishes to compute, for example, separated flows with a parabolized set of equations. The cause of departures in this case should be allowed to remain. It is possible to solve such problems by using a global type of iterative process for the parabolized equations similar to that used by Werle and Vatsa [7] for supersonic boundary-layer interactions. We will discuss such a procedure in the next section. However, for the present let us assume that we wish to filter the growing exponential behavior and, therefore, limit the class of solutions to flows with weak interaction in the streamwise direction. The purpose of this portion of our discussion is to show that the technique for eliminating branches developed by Vigneron, Rakich and Tannehill [3] is an appropriate method for doing this. In the present problem their approach would be equivalent to replacing eqn (5.3) by dA P = -~o dx '

(5.14)

where ~ois the filtering parameter based on the local Mach number. If we look for branches, eclns (5.9) remain the same except that eqn (5.9c) is replaced by (5.15)

P = ooale ~ ,

which results in a governing equation for the branches of the form f"

(5.16)

- k y f ' + k f = -~ok

in place of eqn (5.10). The boundary conditions (5.11) remain unchanged. A careful examination o f the above set of equations will show that departure solutions are eliminated as 00 ~ 0, i.e. f - - 0 when o0 -- 0. We know that the form for 00 used by Vigneron, Rakich and Tannehill [3] is proportional to the local Mach number squared [see eqn (3.17)] and, therefore, should eliminate the branch since ~o -- 0 at the wall. We still need to explore the Mach number variation in the sublayer o f triple-deck theory to see how o~ will vary across the lower deck in a typical case. Near the wall the normal component of velocity v is small and can thus be neglected in determining the local Mach number Mx. Therefore, from the definition of Mach number and eqns (5.7a-e) we find U* 2 M 2 - - - Re-I/2X'/2C1/4(M2 "yR T *

1)-l/4u2M2.

(5.17)

Since u is a lower deck variable of order one in the sublayer, we see that Mx is extremely small in the sublayer region for high Reynolds number. In fact eqn (5.17) shows that Mx = 1 will occur in the middle layer for the following reason. For an undisturbed lower deck (prior to initiation of an interaction) u = y and since y is related to the physical variable y* through eqn (4.7d) by [ T, x3/2

we can easily show that for the undisturbed lower deck (i.e. u = y) Y__*_*= R e - ' / z C '/2 ( T ~ 3 / Z M ~

(5.19)

210

R.T. DAVISet al.

and, therefore, the location where Mx = 1 is clearly in the middle deck since y * / x * = O(Re -I/2) there. We can also see from eqn (5.17) that since 00 pc M 2, the value of o~ will be extremely small in the lower deck region, approaching zero at the wall and remaining small across the whole lower deck. Therefore, the approximation of o~ = 0 in the sublayer as a filter for removing branches is consistent with what one would expect from triple-deck theory. We still need to show that it is not a severe approximation, if handled properly, from an accuracy viewpoint. • Near the wall, we can show that the tangential component of velocity behaves in physical variables like u* = z_.~y , + 1 d p y,2

~,

~ ox

w-2-

+

" " ""

(2.20)

Placing o~ in the m o m e n t u m equation to eliminate part of the pressure gradient term results in a modification to the series expansion near the wall and results in

7"w

u* = - - y * #

+coo

dp[ y,2 dxw

+ .''

"

(5.21)

However, since o~ pc M E pc y,2 the second term now belongs to the terms of order y,4 rather than y,2. This error is not severe as long as the value of o~has returned to one before we leave the linear portion of the velocity profile near the wall; i.e. as long as U* ~ -

Tw y , /.t

(5.22)

we can use the filter with reasonable accuracy. Therefore, the method can only be accurate if the Mach number becomes greater than one in the linear portion of the velocity profile near the wall. In this case, one can produce accurate solutions using the filter even for flows with pressure gradient since the influence of the pressure gradient is a higher order effect in the wall region. As long as the pressure gradient is fully incorporated into the equations outside the linear portion of the profile, accuracy will be retained. One obviously needs to modify the method if separated flows are to be computed. In that case the very mechanism which we are eliminating, the branching phenomenon, is of first order importance and one cannot calculate separated flows without dealing with it directly. It is the mechanism that allows for upstream propagation. A study of the accuracy of solutions obtained by using the oJ parameter will be presented in the section on Results and Discussion. The effect of ~ will be studied by examining selfsimilar solutions to the compressible boundary-layer equations with and without the incorporation of o~. In addition, we will study solutions to the parabolized Navier-Stokes equations and show that the observed departure solutions are of the Lighthill type. 6. GLOBALLY ITERATED SOLUTIONS In their original work Vigneron, Rakich and Tannehill [3] suggested that after one splits the pressure gradient term in the streamwise m o m e n t u m equation into two parts, i.e. writing the term in Cartesian coordinates and nonconservation form as

op

op

a~ = '~ ~

op

+ (i - o,) a-~'

(6.1)

that the first term can be incorporated implicitly into the streamwise marching process for the parabolized equations and the second term can be evaluated with a global iteration process using previously calculated values of pressure. However, it was found that a backward difference for the second term leads to divergence (departure solutions). Based on the analysis of the triple deck and the resulting upstream influence (ellipticity) due to the Op/Ox term it

Calculation of supersonic viscous flows

2! 1

is more logical to use a forward difference for the second term. This approach has been used with success by Rubin and co-workers; see [23, 24] and Khosla and Lai [29] for example, on several problems including subsonic flows. The major modification which needs to be made to a PNS code in order to incorporate the second term in (6.1) is as follows. For simplicity x is a Cartesian coordinate as in (6.1) which may be interpreted as any streamwise coordinate. If i is the index for the streamwise solution station, write

(~xx)

i

Pi-P~-~ -

= ~

Ax

+ (1 - -

60)P'+'-P~ Ax

(6.2)

'

where the subscript j for the normal direction is not written since the equation is the same at all j locations. Note that the second term is only active in the subsonic flowfield since when ~o -~ l this term vanishes. Therefore, if one uses the Pi+~ term which has been stored from the previous iteration, one can march the parabolized equations and develop an iterative process by performing repeated passes. In Section 7, Results and Discussion, we will 'give results obtained from this iterative process. This process, if applied as a Gauss-Seidel scheme, is slowly converging for an obvious reason--the correction can propagate upstream only one mesh point in the x direction per iteration. An alternative method is now developed based upon eqn (6.2). If one rewrites this equation by introducing a fictitious time dependent term as

(oa--xx, 1 = ~o,... - ~

+ (1 -

o.r.i+... (01] , L ~ ~,

(6.3)

an Alternating Direction Explicit (ADE) method can be developed as follows. Equation (6.3) is written.using a two step procedure, letting n denote the time level, as follows, step l:

o("-"-

+

-'7

-,. A, "n),]

(6.4a)

and step 2:

+(,_

,,,)re',,,+: :,,,:)"+' _ (,,'+' L~, ~

~t

L J"

(6.4b)

In both steps all terms other than Op/Ox are written at time level n + ½, so that eqns (6.4a) and (6.4b) can be combined to give

-

~,

Xt

h

At

h"

(6.4c)

If this equation is solved for p~'+~, one finds that

p,/+t = fn(p,/+ll, p~+m, p,,). The expression is written in this form because this new step 2 equation is solved by explicitly sweeping from the outflow location at i = IMAX upstream to i = I. An outflow boundary condition is applied during this forward sweep. Information is propagated upstream in a considerably more rapid manner than the one point per iteration of the Gauss-Seidel type scheme. Step 1 of the ADE scheme is implemented as in the weak interaction (single sweep) method. The parabolized equations are marched from the inflow to the outflow station writing the streamwise pressure gradient as given by eqn (6.4a). This generates pn+l/2 which is used in step 2, eqn (6.4c), along with p~ from the previous iteration, to generate pn+l.

212

R.T. DAVISel aL

It is important to note that the method being used here for the parabolized equations has implications as to how pressure terms should be differenced in full Navier-Stokes solvers. One should see improvements in convergence rate for Navier-Stokes solvers if similar differencing ideas are used. 7. RESULTS AND DISCUSSION As a demonstration of the applicability of the PNS method for flows of practical interest, a number of complex flowfields have been computed, some of which have been presented in previous papers. These include a delta wing with sharp subsonic leading edges, an ogive cylinder and a blunt slab delta wing, all at angle of attack; see [3-5]. O f interest in the present paper are simple flows that illustrate the essential features of the stability and accuracy of the method. Therefore, in the following examples, we consider demonstrations of the effectiveness and accuracy of the ~0 filter factor, and an example of hypersonic viscous interaction.

7. I. Suppression of departure solutions and flow over expansion and compression corners As mentioned earlier in Section 5, it is postulated that the branching behavior observed in solutions of the parabolic Navier-Stokes equations is the same as that obtained from the solution to the lower deck equations; see Eqns (5.10) and (5.11). In order to demonstrate this we wish to show that during a departure solution the pressure and wall shear behave in the manner described by eqns (5.12) and (5.13). Considering those expressions we can write, letting k = 0.8272 and kl = 1.209, dP __ = kale k~, dx hence logl~l

= loglkal[ + kx.

Similarly

logl~xl = log[kklall + kx. I f a departure solution computed using the parabolic Navier-Stokes equations is indeed the same as the solution to the lower deck equations then we should find that a plot of logldp/dx] or logldr/dxl versus kx should have unit slope. Since x is the triple-deck variable defined by eqn (5.7c) we must extract the computational x from that expression. In this analysis, a slightly different expression for the triple-deck x was obtained in order to account for the use o f a Sutherland viscosity law rather than a linear viscosity law as well as to account for an initial velocity profile with a wall shear which might not be exactly the same as the Blasius value. This, however, does not alter the conclusions as long as we plot the dependent variables against a consistently defined triple-deck x variable. It is important to note that the triple-deck analysis is a linearized analysis. Hence, any behavior which is observed with the parabolic equations should only be expected to agree in the initial portion of the departure where the behavior is still essentially linear. Before presenting the results it should be mentioned how the departure solutions were obtained using the parabolic Navier-Stokes equations. At the x-station at which it is desired to initiate a departure solution the streamwise pressure gradient term is fully retained through the entire boundary layer, that is, the o~ factor is set equal to unity. Small errors (truncation and roundoff) already present in the solution are sufficient to trigger the departure solution provided the streamwise pressure gradient term is retained downstream of the initiation of the branch. In order to resolve the departure solutions it is imperative that the finite difference mesh provide adequate resolution of the entire boundary layer and in particular the lower deck. Hence, streamwise stepsizes must be of order Re -3/8 and at the wall the mesh spacing must be of order Re -5/8 in the normal direction.

Calculation of supersonic viscous flows

213

o

pw

1.00

Flat

Plate

S®-3 Re

= 50 m i l l i o n

Tm -

100OK

Location

i. 000

of

u vs.

I. 002

y plots

i. 004

of

Figure

i. 006

I. 008

i. 010 x

Fig. 1. Departure solutions, wall pressure vs. x.

o

Flat ~w 2 p~u

Plate

M~

= 3

Re

= 50 m i l l i o n

T

= 100°K

T

w T

~ o

O .75 1.00

1.2s c~ 1.50

oo?i

1.000

1.002

1.004

i

1.006

i

1.008

i

i. OlO x

Fig. 2. Departure solutions, wall shear vs. x. CAF 14/3--C

214

R.T. DAviset al.

The basic behavior of a departure solution is illustrated in Figs. 1 and 2. They show the wall pressure and shear respectively plotted against the physical x variable for four different wall to freestream temperature ratios. In all cases the freestream conditions are identical. The branches are deafly exponential with their rate being a function of wall temperature. Figure 3 shows the velocity profiles for Tw/To~ = 1.0 at four successive x stations starting at the location of the initiation of the departure solution and proceeding downstream. It can be seen that the acceleration of the flow is confined to a region very close to the wall. This is in agreement with the expected behavior which suggests that most of the changes in the boundary layer during an interaction are confined to the sublayer. Stewartson [20] found similar behavior with his triple-deck calculations. This behavior can be interpreted to first order as an upward shifting of the initial upstream velocity profile by the developing lower deck. The streamwise locations of the plotted u versus y velocity profiles are given in Fig. l on the curve for TJTo~ -- I. Figures 4 and 5 are plots of logl0[dpw/dx[ and loglo[dTJdx[ versus the physical x variable and illustrate the different rates of departure in the unsealed streamwise dimension. Figure 6 is the plot of most interest in this discussion. It shows the variation of log~oidpJdx[ versus kx, where kx is defined using eqn (5.7c). The most important observations to be made are, first, that all curves collapse in the initially linear portion to essentially the same slope and, most crucial, that this slope is very close to the value predicted by the lower deck analysis. (Note that the slope is 0.4343 instead of unity because we plot log~0 instead of log, of the pressure and shear derivatives.) Figure 7 is the comparable wall shear plot. Again, we find that the curves collapse when plotted against k.x, however, they are not in as good agreement with the predicted behavior as the pressure results. In addition, oscillatory behavior is observed just after initiation of the branches. The wall shear was found throughout this study to be much more sensitive to changes in o~ than the pressure. This is reasonable since the o~ factor modifies the streamwise momentum equation which is the principal equation in determining u, hence the shear. It is believed that the incompatibility between the solution prior to initiation of the branch, where Op/Ox was not fully included near the wall, and that

o

o 6 Acceleration

of u / U

Velocity

Profile During Departure

Solution

J l

Case I x-station i

i

J

(i)

1.00,

(2)

1.0021

(3)

1.0042

(4)

1.0063

initiation

of

departure

7

Z

7o

o

I

i o

o o

0.00

0.50 u

1.00

0.00

0.20

0.40

0.60 u

0.80

0.00

Fig. 3. u/U®vs. y during departure solution.

0.10

0.20

0.30

u

0.40

Calculationof supersonicviscousflows

T

o m

X

'

2i 5

w

0.75 1.00 1.25

o~

1.50

Flat M

Plate = 3

Re = 50 m i l l i o n T

-'4 I

T

i. 000

i. 002

~

i. 004

I

i. 006

= 100°K

I

i. 008

I

i. 0 1 0 X

Fig. 4. Departure solutions, log,0~wx[ vs. x.

which follows inclusion of Op/Oxtriggered the observed oscillations in wall shear. This may also account for the less impressive agreement with triple-deck theory for the wall shear as compared to the wall pressure. Figure 8 shows the ability of the parabolic Navier-Stokes solution to recover from a departure solution if the streamwise pressure gradient term is again suppressed after branching has begun. Even if this is done late in the branch it is observed that the solution recovers after a small overshoot to the original weakly interacting solution. Figure 9 demonstrates a calculation of flow over a sharp one degree expansion ramp at a freestream Mach number of 3 and a Reynolds number of 50 million. Departure solutions were initiated at several stations upstream of the corner. It can be seen that a very small perturbation in the location of the initiation of the branch causes a large change in the downstream solution. The interesting result is that if the proper upstream location is found at which to start the branch, the solution will proceed more or less smoothly through the corner and asymptotically approach the proper downstream condition on pressure. This is further evidence of the physical significance of departure solutions. Use of this behavior was made in a 1968 study by Baum [25] where the branching was employed to develop a shooting method to solve strong interaction problems. Based upon the above calculations and comparisons between triple-deck theory and the results of parabolic Navier-Stokes calculations, we conclude that a departure solution is the same as the growing exponential solution to the lower deck equations discovered by Lighthill. This is intimately connected with separation phenomena. Thus, it has been demonstrated that the parabolic Navier-Stokes equations have the terms necessary to solve for flows with strong viscous/inviscid interaction, including separation. The mathematical significance of the branching behavior is that the problem is improperly posed as an initial-value problem and that a downstream boundary condition on pressure

o --

I

M T

_E

o ~4

T

ooo .75 1.00 1.25 1.50

o 7-

Flat M

7

Plate

= 3

Re = 50 m i l l i o n T

] . )00

i. 002

i

I

1.004

1.006

= 100OK

i

L

1.008

1.010 X

Fig. 5. Departure solutions, logmlT~.l vs. x. o

T

o

w

T

X

0.75

m

1.00

o r-4 o

1.25

~o

1.50

Flat

Plate

M~

= 3

Re

= 50 m i l l i o n

T

= 100°K

o o

,theoretically predicted slope

I

I

F

iO.O0

0.50

i

1.00

i

1.50

i

2.00

Fig. 6. Departure solutions, logmLD~l vs. k.x;.

r

kx~

2.50

T

_~w

__I )4

T .75

o

1.00 ,-.4

o

1.25

,--I I

1.50 Flat M

Plate

= 3

Re = 50 m i l l i o n

I

T

= 100°K

theoretically predicted slope

c~ I

I

I

o v

l

I

I

I

1.50

2.00

I

I

0.00

1.0O

0.50

kx£

2.50

Fig. 7. Departure solutions, loglo[r~[ vs. kxl.

1

2

3

Flat P~U~

M

Plate

= 3

Re = 50 m i l l i o n T w / T ~ = 1.0 T

/

= 100°K

Uninterrupted Departure Solution

Case

Xinterruotion

(i)

1.0028

(2)

1.0042

(3)

1.0056

o ~ o

~

1.000

1.005

1.010

i

1.015

m

m

1.020

1.025 x

Fig. 8. Effect of interrupted departure solution, pressure vs. x.

R. T. DAVISet

218

al.

Od

Pw o o 6

1° Expansion

Corner

M = 3 Re = 50 m i l l i o n

(2) (3)

.994185400 .994185110

L

(4)

.994184960

I

(5)

.994178B00

i i

d •9900

.9925

o .994203000

IT =

o

x

(i)

Tw/T

~

= 1.0

Case

= 100°K

t

J

x O .9950

.9975

,

5

t

1.0000

t

1.0025

1.0050

1.0075 x

Fig. 9. Flow over expansion comer, pre~ure vs. x. is required. All of these facts suggest the possibility of solving for flows with strong viscousinviscid interaction and separation using the parabolic Navier-Stokes equations in conjunction with a global type ofiterative process similar to that used by Werle and Vatsa [7] for supersonic boundary-layer interactions. The simplest type of global iteration process on pressure is that given by eqn (6.2). The same I o expansion comer case was recomputed using identical step sizes and globally iterated using this method. The procedure consists of saving the pressure array for grid points in the subsonic layer and inserting these old values into the Pi+t terms in subsequent iterations. All of the terms are computed implicitly. As mentioned previously, the method converges slowly and improvements are possible by using an ADE approach as in (6.4). The results are shown in Fig. 10 for pressure and Fig. 11 for shear. The results compare well with those obtained from the shooting method. The next case chosen to test the global iteration method is that of a compression comer. A triple-deck angle for the compression ramp was chosen such that mild separation occurs 7o

1°

Expansion

Corner

M = 3 Re

= 50x106 xc

L/T®

.

ao

.-,-...

. . . .

T

- 1.0

= 100°K

I ion

~

• 990

i

°

Iterated

PNS Solution /

1.0

r,.

.992

.994

.996

.998

1.80

1.002

1.004

1.006

1.008

Fig. 10. Expansioncomer, comparisonof shooting solution and iterated PNS solution,Pwvs. x.

1.01

Calculation of supersonic viscous flows

219

?

,-d

! 1° ExPMin~i~n Corner

Rexc = 50xi06

I ]

°o•' = I

~% ~'

/

/

lO

,o

Tw/T = 1.0 T = 100°K

"<,

/

.....

o,

"

,

I

•,

Solution

%

~990

.992

.994

.996

.998

1.00

1.002

1.004

1.006

1.008

1.01

x - Distance from Leading Edge Fig. l 1. E x p a n s i o n c o m e r ,

comparison of shooting solution and iterated PNS solution,

~w vs. x.

(see Burgraff, Wefle, Rizzetta and Vatsa [26]). Figures 12 and 13 show the results of these computations. The calculations are done on a triple-deck scale and the sensitivity to Ax step size is shown. The results are in agreement with expected results from triple-deck and interacting boundary-layer theory. We have shown how a departure solution is related to Lighthill's eigensolution to the sublayer equations and discussed the implications of its existence. In most calculations made in the past using the parabolic Navier-Stokes equations only weakly interacting flows were considered. For such flows we wish to suppress the departure solutions while still producing accurate results. Thus, we would now like to discuss limiting the class of solutions to flows with weak interaction in the streamwise direction.

O. Smooth Compression Corner •

M® = 3

8. | 88-

Rex=l " lxlO8

~

~/T® 1.4 -

-

T= = 560°R a - 2.5 (~* = 1.37 °)

Solution,

8.875-

i -

Ref. [54]

PNS Solution, ~x = .0002

I

PNS Solution,

|

]

, 8.8588 -

~x "

.0004

IBL Solu PNS Solution Ax = .0002

PNS Solutlon,~ Ax = .0004

8. 8 2 5 --

o*

1.0 0.080--~ G

.g7

9.98

0.9g

I .89

I .91

I .92

I .93

x - Distance from Leading Edge Fig. 12. S m o o t h

compression ramp, comparison

of PNS

solutions and IBL solution,

Pw vs. x.

R.T. DAVISet al.

220

7.2. Accuracy ofPNS equations with o~.6lter It is possible to analyze the effect upon the accuracy of the solution due to modification of the streamwise pressure gradient in the manner of Vigneron, Rakich and Tannehill [3]. To accomplish this the self-similar form of the compressible boundary layer equations are considered. The self-similar form of the m o m e n t u m equation is written as

d=I dd~ l d2f) ~ 2 + ,r.d~2 -/3(0- f,2)= O. The term/30 corresponds to the streamwise pressure gradient term where 0 = T/Ti and/3 is the usual pressure gradient parameter. The quantity df/d~ is the dimensionless scaled u velocity component and I is the ratio of the product of density and viscosity at a given location referred to that at the edge of the boundary layer. The term /30 is now multiplied by ~o where ~0 is defined in the same manner as by Vigneron et al. [3], that is

-rm x l + (7-

1)M 2

and Mx can be written in terms of F and 0 as F

Mx=Mi-~, where Mi is the inviscid outer edge Mach number. When ~o is greater than unity in the equation, it is set equal to one. The modification described above was made to a self-similar compressible boundarylayer code and the results which were obtained are compared with those of the unmodified version of the same program. For all results presented, we chose a Prandtl number of unity and a linear viscosity law. The value of/3 is 0.60. The results are presented in Fig. 14 in the form of velocity profiles. The inviscid Mach number is 8 and various wall to freestream temperature ratios are presented as well as the case of an insulated wall. As the wall

XIO -4 Smooth Compression Corner

8.S-.

M®- 3 Rex'l " ix108 Tw/T® - 1.4 T = 560°R ~ = 2.5 (~* ~ 1.37°)

~ ~

8.2~ ~~

IBL Solution, Ref. [54] PNS Solutions ~x

- .0004

g.l

k~

1.0

-g. I g.g7

g.g8

g.gg

I .gg

1 .gl

I .g2

I .g3

x - Distance from Leading Edge

Fig. 13. Smooth comprcssionn ramp, comparison of P N S solutions and IBL solution,rw vs. x.

Calculation of supersonic viscous flows

Boundary Layer Solutions

~__~

221

TW/Te

Me - 8

0.50

Pr - 1 8

1.00

= 0.6

1.50

Linear Viscosity Law

8

Adiabatic S t a n d a r d Px .....

Modified

Px

I c M ffi

o

.I0

0.0

.20

.30

.40

.50

.60

.70

.90

.80

1.0

Fig. 14. Effect of wall temperature on self-similarvelocity profile.

temperature decreases the accuracy of the modified result improves until at a value of O,~l = 0.50 the two velocity profiles are virtually identical. Figure 15 shows the effect of o~on self-similar solutions at Mach number 8.0 as a function of wall temperature and/3. Finally Fig. 16 shows the effect of Mach number for a given wall temperature as a function of/3. These results imply that treatment of the pressure gradient term by a method such as that of Vigneron et al. [3] is more accurate for cold walls than for hot walls. Care should therefore be exercised as the wall temperature is increased. It is also noted that the larger the magnitude of the streamwise pressure gradient, the poorer the accuracy for a given wall to outer edge inviscid temperature ratio. The results of Fig. 16 confirm the point discussed Boundary Layer Solutions Me - 8 Pr = 1

Linear Viscosity Law

%

o,

..... Cas_~e

n,

o

-

1

Standard

Px

M°dified

Px

Tw/T_____~e Adiabatic

~

~

Wall

~

2.50

3

1.50

/

.20

.30

i ,

/

2

I

/

-- . - - . . ~ ~ -_ .~.

-~" " " " --

_ .... -~" -- "" ~

i

,°

o °0. 0

• 10

.40

.50

.60

.70

8-

.SO

Pressure

.90

1.0

Gradient Parameter

Fig. 15. Effect of wall temperature on sel~similar skin friction parameter, F~..

222

R.T. DAVlS et al. Boundary Layer Solutions Tw/T e

=

i

Pr - 1 o.

Linear Viscoslty Law

Me = 4 , ~

M e = 6N~' M e " 8N%

o J

0.0

• 10

• 0

• 0

.40

.60

.50

.70

.80

.90

1.0

8 - Pressure Gradient P a r a m e t e r

Fig. 16. Effect of M a c h n u m b e r o n percent error in sel~similar skin fnction parameter, F,..

earlier with regard to the greater accuracy of the results obtained using the ~ofactor as the sonic line approaches the linear portion of the velocity profile.

7.3. Hypersonic viscous interaction The PNS equations, it is argued, should model the complete viscid-inviscid supersonic flow for high Reynolds numbers, in the absence of axial separation. As an example, the viscous high speed flow over a fiat plate is one which the growth of the viscous layer produces a displacement of the inviscid flow, with a corresponding pressure rise. This type of flow received much attention a number of years ago, resulting in weak and strong interaction theories; see Hayes and Probstein [27]. The two asymptotic theories apply for extreme values of the hypersonic interaction parameter, ~ = C~/2M3R -1/2, but fail at intermediate values. Here C --- # T ~ / # ~ T , is the Chapman-Rubesin constant, which is evaluated using wall properties. The results of the PNS code for the hypersonic viscous interaction problem are shown in Fig. 17; see Fig. 9.5 in [27]. Two cases are shown, M = 5.8 and 9.6, corresponding to EXPERIMENT, REF. 9

Y////2 M . - 5.e ~\\\\'~ M,.- 9.e 6

PNS METHOD M~ " 5.8 ~--0 M~° ,. 9.6

STRONG INTERACTION

5

,¢,k~ WEAK INTERACTION

°'3'

"

,

~

"

~"

2 /y/" 1~ ~

I 1

a 2

I 3

. 4

. 5

' 6

' 7

' 8

Fig. 17. Comparisonof hypersonicpressure interactionexperimentand theory.

Calculation of supersonic viscous flows

223

experimental conditions. We note first, that for small values of ~ (weak interaction) the PNS results, interaction theory, and experiment all agree. As ~ is increased, weak interaction theory fails, while the PNS result varies smoothly with the general trend of the experiment. For large values of ~ (rarefied flow) there is fair agreement with strong interaction theory, although not as good as for a weak interaction. The difference between theory and experiment at large ~ is not surprising, considering the difficulty of experimental measurements under rarefied conditions and the fact that the PNS model does not include some effects which may be important in this regime. The important point to note is that the PNS method captures the essential features of the viscid-inviscid interaction over the entire range of continuum viscous flow. 8. C O N C L U D I N G

REMARKS

In this present paper, we have attempted to place the parabolic Navier-Stokes (PNS) method on a more sound theoretical basis by examining the fundamental assumptions of the method, and by demonstrating its compatibility with the "triple-deck" theory of Stewartson [20]. The analysis and examples have shown that stable (departure free) marching of the PNS equations for supersonic compressible flow is related to the filtering of the eigensolutions to the sublayer (lower deck) found in triple-deck theory. If allowed to, the PNS equations correctly capture the known features of so-called departure solutions corresponding to a longitudinal compression or expansion wave. By use of an appropriate filter on the pressure term in the x-momentum equation, the departure can be suppressed, giving the usual desired weak interaction solution. A correct manipulation of the pressure term along with a global iteration process can allow one to correctly compute strongly interacting and separated flows as well. The fundamental reason for departure solutions is that the boundary-value nature of the equations is ignored in the PNS approximation. A downstream boundary condition on pressure is necessary to make the equations mathematically well-posed. The downstream boundary condition has little influence on the strong interaction solution if the downstream boundary is located far enough away from the strong interaction, but its inclusion is of paramount importance for computing such flows. Therefore, in order to compute strongly interacting flows, the boundary-value nature of the problem must be accounted for by some sort of iteration or relaxation procedure and the axial pressure gradient term in the streamwise momentum equation must be properly differenced. At the present time, for strongly interacting and separated flows, the compressible Navier-Stokes equations are usually solved in their time dependent form. A goal of future research is to solve such flows more efficiently using PNS codes with an appropriate iteration procedure to incorporate the downstream boundary condition. The paper of Barnett and Davis [28] considers the development of such a procedure utilizing the principles laid out here. These procedures should also be useful in improving convergence rates of full Navier-Stokes solvers. REFERENCES I. L. Prandtl, Uber Fliissigkeitsbewegung bei sehr kleiner Reibung. Proc. Third Intern. Math. Congress, Heidelberg (1904). 2. M. J. Lighthill, On boundary layers and upstream influence. II. Supersonic flows without separation. Proc. Roy. Soc. London A 217, 478 0953). 3. Y.C. Vigneron, J. V. Rakich and J. C. Tannehill, Calculation of supersonic viscous flow over delta wings with sharp subsonic leading edges. NASA TM 78500 (June 1978). 4. J. C. Tannehill, E. Venkatapathy and J. V. Rakich, Numerical solution of supersonic viscous flow over blunt delta wings. AIAA J. 20, No. 2 (Feb. 1982). 5. Y. C. Vigneron, J. V. Rakich and J. C. Tannehill, Calculation of supersonic viscous flow over delta wings with sharp subsonic leading edges, AIAA paper 78-1137. I lth Fluid and Plasma Dynamics Conference, Seattle, WA (July 1978). 6. H. Viviand, Conservation form of gas dynamic equations. Recherche Aerospatiale No. 1 (Jam-Feb. 1974). 7. M. J. Wefle and V. N. Vatsa, Numerical solution of interacting supersonic boundary layer flows including separation effects. Aerospace Research Laboratory Report ARL-73-0162 (Dec. 1973). 8. R. T. Davis and S. G. Rubin, Non-Navier-Stokes viscous flow computations. Comput. Fluids 8, 101-131 (1980). 9. R.T. Davis, Numerical solution of the hypersonic viscous shock layer equations. AIAA J. 8, No. 5 (May 1970). 10. B. N. Srivastava, M. J. Werle and R. T. Davis, Viscous shock-layer solutions for hypersonic sphere cones. AIAA J. 16, No. 2 (Feb. 1978).

224

R.T. DAvIs et al.

11. T.C. Lin and S. G. Rubin, Numerical methods for two- and three-dimensional viscous flow problems. PIBAL Rept. 71-8, Poly. Inst. of Brooklyn (Apr. 1971). 12. H. K. Cheng, S. Y. Chen, R. Mobley and C. R. Huber, The viscous hypersonic slender body problem. RM 6193-PR, The Rand Corp., Santa Monica, CA (May 1970). 13. S. C. Lubard and W. S. Helliwell, Calculation of the flow on a cone at high angle of attack. AIAA J. 12, No. 7 (July 1974). 14. L. B. Schiff and J. U Steger, Numerical simulation of steady supersonic viscous flow. AIA.A J. 18, No. 12 (Dec. 1980). 15. R. Beam and R. F. Warming, An implicit finite-difference algorithm for hyperbolic systems in conservationlaw form..I.. Comput. Phys. 22, (1976). 16. T. Lindemuth and I. Killeen, Alternating direction implicit techniques for two-dimensional magnetohydrodynamics calculations. J. Comput. Phys. 13 (1973). 17. H. McDonald and W. R. Briley, Three-dimensional supersonic flow of a viscous or inviscid gas. J. Comput. Phys. 19, No. 2, 150-178 (1975). 18. J. L. Steger, Implicit finite difference simulation of flow about arbitrary geometries with application to airfoils. AIAA paper 77-665, Albuquerque, NM (June 1977). 19. M. Barnett, R. T. Davis and J. V. Rakich, Implicit boundary conditions for solution of the parabolized NavierStokes equations for supersonic flows, J. Comput. Phys. 48, No. 2 (Nov. 1982). 20. K. Stewartson, Multistructured boundary layers on flat plates and related bodies. Advan. Appl. Mech. 14 (1974). 21. J. Ackeret, F. Feldmann and N. Rott, Investigations of compression shocks and boundary layers in gases moving at high speed. Translated in NACA TN1113 from E. T. H. Ztirich, No. 10 (1946). 22. D. R. Chapman, D. M. Kuehn and H. K. Larson, Investigation of separated flows in supersonic and subsonic streams with emphasis on the effect of transition. NACA Rept. 1356 (1958). 23. S. G. Rubin, A review of marching procedures for parabolized Navier-Stokes equations. In Numerical and PhysicalAspects of Aerodynamic Flows (Edited by T. Cebeci), pp. 171-185. Springer-Verlag, New York (1982). 24. S. G. Rubin, Analysis of global pressure relaxation for flows with strong interaction and separation. Comput. Fluids 11, No. 4, 281-306 (1983). 25. E. Baum and M. R. Denison, Interacting supersonic laminar wake calculations by a finite difference method. AIAA J. 5, No. 7 (July 1967). 26. O. R. Burggraf, D. Rizzetta, M. J. Werle and V. N. Vatsa, Effect of Reynolds number on laminar separation of a supersonic stream. AIAA J. 17, No. 4 (April 1979). 27. W. D. Hayes and R. F. Probstein, Hypersonic Flow Theory, pp. 333-374. Academic Press, New York (1959). 28. M. Barnett and R. T. Davis, A procedure for the calculation of supersonic flows with strong viscous-inviscid interaction, to be presented at the AIAA 23rd Aerospace Sciences Meeting Reno, NV (January 1985). 29. P. K. Khosla and H. T. Lai, Global PNS solutions for subsonic strong interaction flow over a cone-cylinderboattail configuration. Comput. Fluids 11, No. 4, 325-329 (1983).