Numerical studies in high Reynolds number aerodynamics

Comp~en and Fiaids Vol. 8, pp. 1-30 Pergamon Press Ltd., 1950. Printed in Great Britain

NUMERICAL STUDIES IN HIGH REYNOLDS NUMBER AERODYNAMICS J. C. LE BALLEUR, R. PEYRETt and H. VWIAND O~ce National d'Etudes et de Recherches A6rospatiales (ONERA) 92320 CMtillon, France

(Received 3 January 1979) Abstract--Two numerical approaches are presented for the computation of viscous compressible flows at high Reynolds' numbers, tn the first approach, named global approach, the whole flow field, which includes viscous and inviscid regions, is determined as the solution of a single set of equations, which may be the full Navier-Stokes equations, or some approximate form of these equations. The second approach, named coupling approach, consists in solving two different sets of equations in their respective domains simultaneously; one of the two sets governs an inviscid flow whose boundary conditions are provided by the viscous effects, determined by the other set. The discussion of the global approach is centred on two particular features of the finite-difference method used: a discretization technique, directly in the physical plane with arbitrary meshes: and a mesh adaptation technique, which combines a coordinate transformation to fit the mesh system to particular lines in the flow, and a technique of dichotomy for mesh refinement. Numerical results are presented for an axisymmetric compression comer and a shock-boundary layer interaction on a flat plate, both in supersonic regime, and for a transonic nozzle flow. For the coupling approach, emphasis is given firstly to the improvement resulting from an interacting analysis where the viscous and inviscid computations are matched, and not only patched. It is shown that the parabolic problems associated with simple viscous theories are always replaced by elliptic problems, even for supersonic flows, and that "supercritical interactions" or "critical points", as defined by Crocco-Lees, are removed. Secondly, a new coupling method, fully automatized and capable of solving directly a well-posed problem for supersonic flow, is illustrated by examples involving shock waveboundary layer interactions and reverse flow bubbles; they concern flows over symmetrical transonic airfoils and supersonic compression ramps.

1. INTRODUCTION

Progress made in computational methods in fluid mechanics allow the determination of flows of increasing complexity. Thus much attention is now given in computational aerodynamics to the viscous-inviscid interaction phenomena which are frequently encountered in real flows. The work presented here relates to the computation of high Reynolds' number compressible flows, such as flows of current interest in transonic or supersonic aerodynamics, in the case when the inviscid flow region and the viscous flow region are more or less strongly coupled. There are essentially two numerical approaches to determine such flows: a global approach and a coupling approach, and this paper deals with some aspects of these two approaches which are being studied at ONERAL In the global approach, the flow field of interest is determined as the solution of a single system of equations valid everywhere, that is to say independently of the importance of viscous effects. Generally the full Navier-Stokes (N.S.) equations (or the time-averaged N.S. equations together with a turbulence model in the case of turbulent flows) are used, but it is equally possible in many cases to use approximate forms of these equations involving simplified viscous terms. However, it is a characteristic feature of the global approach that all the non-dissipative terms of the N.S. equations, i.e. the full Euler equations, should be included in the equations used, so that a continuous transition is automatically insured between the inviscid flow and the viscous layers. The usual division of the flow field into an inviscid flow region and viscous layers remains however a valid approximation for viscous-inviseid interaction problems at high Reynolds' numbers, even if the classical boundary-layer theory breaks down. It is known that this failure is not always consequent to Prandtl's equations but actually to the simplified matching analysis tAIso, Institut de M6canique Th6orique et Appliqu6e, Universit~ Pierre et Marie Curie, Paris. ~Work performed with financial support from DRET and STA6 of the French Ministry of Defence. Section 2, relating to the global approach, is contributed by R. P. and H. V. and Section 3 on the coupling approach, is contributed by L C. L. B.

2

J, C. LE BALLEURet ai.

which relates the viscous boundary conditions to some given external inviscid flow, without possible allowances to strong viscid-inviscid interaction. If no approximation is involved in the matching, then high Reynolds' number flows can be determined partly from an inviscid flow solution, partly from a viscous additional computation, for which Prandtl's equations as well as less approximate ones may be used. Thus, one is led to the coupling approach which is based on the coupled solution of two different sets of equations, one of them governing an inviscid flow with boundary conditions provided by the viscous effects, determined by the other set. The first part of the paper (Section 2) deals with the global approach, the emphasis being laid on two main aspects of the numerical method used: (i) discretization of the N.S. equations in arbitrary meshes directly in the physical plane, using a simple, explicit finite-difference scheme, and (ii) adaptation of the mesh to the solution in a viscous layer. The coupling approach is considered in the second part of the paper (Section 3), restricted to flows for which thin shear layer approximations are relevant, so that computing costs can be greatly reduced by using integral methods for the viscous effects. This approach is subject to both analytical and numerical difficulties. The analytical ones relate to the proper approximation of the N.S. equations in order to obtain a set of viscous equations, as well as an adequate interacting formulation. The numerical problems lie in the inviscid flow solution, the viscous flow solution, and the coupling techniques; of these three points, we shall consider almost exclusively the crucial point of the coupling techniques. 2. G L O B A L A P P R O A C H

Research on the numerical solution of the compressible N.S. equations has developed rapidly in the last ten years, and much progress has been done, in particular concerning the reduction of computing times through the use of implicit or hybrid schemes. General reviews or discussions on this subject can be found in Refs. [I-5] and presentations of some of the methods currently used or under study, for example, in Refs. [6.-14]. The method which is pres~ted here still makes use of an explicit scheme, with the corresponding advantage of simplicity, but it involves two particular features which relate respectively to the problem of discretization in an arbitrary mesh and to the problem of adaptation of the mesh to the solution in a thin viscous layer. The discussion will be essentially restricted to these two features. 2.1 Discretization directly in the physical plane In general, for the problems we are interested in, two types of methods are used for discretization of the equations in arbitrary curvilinear meshes: either (i) in the case of finite-difference methods, a coordinate transformation is carried out so that the transformed equations are discretized in a rectangular mesh in terms of the new coordinates, or (ii) in the case of finite-volume methods, the equations are written in integral form for an arbitrary control volume, and the discretization is carried out directly in the physical space using the finite volumes (or cells) determined by the mesh system. The present method is of the finite-difference type, but, as with finite-volume methods, the discretization operates directly in the physical space without a transformation of the equations. This method is an extension to the viscous case of a method used by Veuillot and Viviand [15, 16], for inviscid flows. We start from the N.S. equations for plane (a = 0) or axisyrnmetric (a = 1) flows, written in conservative form in an absolute frame of reference, using cartesian or cylindrical coordinates (x, y), where y is the radial distance in the axisymmetric case: dU 0F

0G+

a-7+W+

F = F~ - F2

U=

, \.pE /

_

a--0

G = Gt - G2,

pu2+/~

\(pE + p)u/

=

(2.1)

H = H1 - 1t2

(,uvt

V_U+V 0 0

~(pE + p)vl

Y

Y10

Numericalstudiesin highReynolds'numberaerodynamics

3

where

E=e+~uZ+v2),

p=(7-1)pe,

/ ~ = p - a A S. Y

In these equations p is the density, u and v are the components of the velocity, p is the pressure and e the specific internal energy. The fluid is assumed to be a perfect gas with constant specific heats of ratio 3'. The dissipative terms Fz, G2 and He are written as:

F2 =

/°

)

"rxxp "rxy

U'rx~t, + tYrxy- qx

o ,G2 =

)

"rxy

o ,H2 =

/

%yo

2~(~d~y - v/y)]

UTxy "{" D"/'yyp -- qy

u¢~y+ v'ryyp- qy

This form displays the dissipative terms for the case of plane flows,, F: and G:, which are expressed in terms of the plane viscous stresses ~ and ~'.o:

Ou A Ov %yp = (A + 2 ~ ) ~ + AOu

7x

where A and ~ are the viscosity coefficients. For axisymmetdc flows, we have: 1"= = ~'~ + Ady,

~'. = ,'.p + Av/y

and the extra-term (Av/y) has been added to the pressure p in FI and G~, which explains the definition of p. The other terms in F: and G~ have the same expression for plane and axisymmetric flows:

%y=/Z ~y Ox/ k Oe q~= c~Ox' qY

l Oe c~Üy

where k is the thermal conductivity coefficient and c~ the specific heat at constant volume. The eqns (2.1) will be discretized on a mesh of classical type, i.e. made of two families of lines, so that a mesh point M can be identified by means of two indices (k, l). ARhough the mesh is generally constructed through a coordinate transformation (especially in view of fitting the boundaries of the computation domain on mesh lines), here we need not explicitly introduce such a coordinate transformation, and we consider that the mesh is known only through the coordinates xu and Yu of the mesh points (k, I). In some problems, in particular in problems involving an unknown boundary (e.g. a bow shock) or in the case of the flow past a body in arbitrary motion, part or all of the boundary moves in time, and the mesh points also move in the absolute frame of reference; the motion of the mesh must also be taken into account in an unsteady flow problem if a technique of self-adjustment of the mesh to the solution is used. Therefore we shall consider Xk.t and Yk.t as functions of time which are known or which are determined at each time level depending on the problem. The change in time of the value of any function f at a given mesh point with coordinates x(t), y(t) (omitting the indices k, ! when they are not necessary) will be denoted af/O~"with ~-= t, and it can be calculated as a material derivative for a medium with local velocity .~, ~, hence:

c~"

dt

ox

dy

Thus the Navier-Stokes eqns (2.1) can be written in terms of the change in time of U

J.C. LE BALLEURet al.

4

when following a given mesh point in its motion: 0U

.,gU . a U + a F + a G + a H = 0

Or

X~x-Y'~y

OX

3y

~,.~) .....

Approximations of the space derivatives are now obtained by using finite-difference operators defined on the instantaneous mesh. Approximations to the first-order derivatives Ocb/~x, O~/Ùy, at a mesh point M, of a quantity • defined at the mesh points, can be derived from the values of ~ at M and at two other neighbouring points P and Q (such that the points M, P and Q are not aligned). By expressing • (P) and ~ Q ) through Taylor expansions at M, to first order, we get a system of two linear equations for OO/0x and a~/Oy at M, the solution of which is:

7 x ~t (O~)

"~Y M

YM~MO- y~xM~, ~(6y~)M = XMp ~ MQ XMQ~ MP

(2.3)

xMPYMO --XMQyMp

where Duo =4)(Q)-O(M), yue = y(P)-y(M), etc... These first-order accurate finitedifference approximations (2.3) represent the constant values of Ocp/Ox, 0q)/0y over the triangle MPQ if • is assumed to vary linearly with x and y over this triangle. From the four neighbouring mesh points, i.e. points N, S, E, W as shown on Fig. l, we define the following finite-difference operators of the type (2.3): 8x, 8y using the points M, N, E, and 8~, 6y using the points M, S, W. The discretization of eqns (2.2) is carried out by a two-step explicit scheme using the operators 8~, 8y at the first step, and the operators ~, 6y at the second step (other variants are obviously possible). Let the superscript n refer to the time level t, = nat, and (h'~'i) refer to the predictor values; then the discretized equations at point M(x~t, y~t) read: ~k,lrT"+l= U~t.+ At{,~rIU + ySyU- SxF - 6yG - ctH}~l

u.÷,

z u"

""+"

+ ~,~,u- 8,F- &G- o,n}t~"

(2.4.a)

(2.4.b)

These equations do not give a complete picture of the scheme yet, since F, G and H include dissipative terms F2, G2 and H2 which are linear functions of the first-order derivatives of u, v and e with respect to x and y. These derivatives in F2 and G2 will be approximated by means of the operators 8x, gy at the first step, and by means of the operators 8x, 8. at the second step, but the derivatives in It2 will be discretized with the same operators as F and G. Clearly such a scheme is an extension of the well-known MacCormack scheme [17], In the case of a uniform rectangular mesh in the physical plane, the present scheme reduces exactly to the MacCormack scheme as far as the non-viscous terms are concerned; the treatment of the

W ~

l+1

Fig. 1. Definition of finite-difference operators.

Numerical studies in high Reynolds'numberaerodynamics

5

viscous terms is slightly different since the MacCormack scheme uses centered differences for some of the viscous terms. It can be shown that the truncation error of the scheme (2.4) in the case of a scalar equation analogous to (2.1) is 0(AE~, At AE2, At:, AE2e) where: AE, = Max {IME+ MWI, IMN + MSl} AE2 = Max {IMEt, IMWI, IMNI, IMSl}, so that the scheme is of second order accuracy if AEI = 0(AE::). In problems in which flow properties vary very rapidly in space, short wave-length numerical oscillations are likely to occur. Such instabilities, which are outside the scope of the usual linear stability theory for initial value problems, may lead to the divergence of the calculations, and it is necessary to damp them out by means of an artificial viscosity. Thus a third step is added after the corrector step (2.4.b), for example in the form of the fourth-order smoothing term used by Gnoffo [18]: Un+I

k,I Is ~

un+!

t. k.I -EllUk+2.1+ Uk-E,I-4(Uk+11+

- ~:{ Ukj+:+ Uk.l-2- 4( Ukj +I + Ukj-i) + 6 Uk.I}'+m.

k.l.tln+!

(2.5)

The positive coefficients el., e: are of the order of 10 -3 in the computations reported below. 2.2 The problem of adapting the mesh to the solution One advantage of the global approach is that it does not require to mate a distinction between viscous and inviscid flow regions, as far as the equations are concerned. Nevertheless, computing time and memory limitations make it advisable to keep the number of mesh points as small as possible, which means that the mesh size should be approximately adapted to the local nature of the flow. The problem of adapting the mesh to the solution is a complex one when considered in all its generality. Interesting approaches for constructing such meshes as functionals of the solution have been presented by Gough et al. [19] and by Yanenko et aL [20] for special cases and should lead to important developments. Up to now the methods used for practical applications in aerodynamics are less general in their principle and appeal to some a priori knowledge of general properties of the solution. We consider the restricted but important problem of mesh adaptation in the case of a flow made of an outer inviscid flow region and a thin viscous layer near a wall, so that the mesh size in the transversal direction must vary by several orders of magnitude. The technique to be described involves two aspects: (i) the basic mesh system should be fitted not only to the wall, but also approximately to the edge of the viscous layer, so that the mesh refinement introduces the minimum number of mesh points, (ii) the mesh refinement should be easily modified during a calculation for adaptation to the solution. To begin with, we assume that coordinates (~, 77) have been defined in the (x, y) plane, such that the wall can be represented by an equation of the form: ~ = ~off). The outer boundary of the computation domain, which is located well outside the viscous layer, has an equation of the type: rt = ,~(O- Assuming 7/o <~E, the computation domain is the curvilinear rectangle a ~<~ ~
st=~

a<~b

7/- r/o = g,(~ ~; t), r/E -- r/o

(2.6.a) 0~<~<1

(2.6.b)

6

J. C, LE BALLEURel at.

with the following conditions on the function g~: ~g~/ari > 0 g,(s¢, O; t) = O, g~(~, 1: t) = 1 g,(#,B:t)=g(~;t), 0<3
(2.7)

where x(se; t) = (n8 - no)/(ne - 70),

0< x < 1

(2.8)

A second transformation can be made independently on the ~ and ~ variables to define the final coordinates X, Y: =/(X) ¢1= g(Y),

g(0) = 0,

(2.9.a) g'(Y) > 0.

(2.9.b)

The basic mesh system is taken uniform in the X, Y plane, with mesh sizes AX, AY which are appropriate to the inviscid flow. Possibly, f = X and g = Y, but the transformation (2.9) can also be used to make a rough adjustment to some features of the inviscid flow as well as to make a slight preliminary refinement in the ~-direction near the wall. Contrary to current methods, the mesh refinement required for a proper description of the viscous layer is not to be obtained through coordinate transformations of the type (2.6.b) and (2.9.b). The essential purpose of the transformation (2.6.b), besides the fitting of the boundaries = ,lo and ,1 = ~e, is to take into account the variation of X with respect to ~; thus if X were constant, it would be sufficient to take g~ = ~ (hence 3 = )¢). In general, ~, is found to vary in some interval [X~, X2], and 3 should be chosen in this interval (e.g., 3 -- (X~ + ~¢z)/2, so that it is possible to make g~ = ~ when X =/3. The image of the approximate edge of the viscous layer in the (X, Y) plane is the coordinate line Y = Ys, with Ys defined by g~(~ g ( Y s ) ; t ) = X(~ t). Knowing Ys, we can decide upon the smallest mesh size in Y, say AY~>, to be used in the viscous layer, depending on the total number of mesh points used to describe this layer. The corresponding degree of mesh refinement is characterized by the ratio AY~U/AY, and it obviously depends mainly on Y d A Y . This refinement is achieved through a technique based on dichotomy of the mesh size in Y. which has been previously used by Viviand and Ghazzi [21] and which we describe now. The computation domain is divided into a number R of zones of the form: YI"~~< Y ~< Y["

for the "zone r", r = 1 to R,

the mesh size in Y being constant, equal to A Y~" in each zone. Moreover, from a zone to the next one, the mesh size varies by a factor of 2, so that, if the zones are numbered consecutively with r increasing with mesh size, we have: A y~,+l~= 2AYe" = 2,AYe I~ Thus in a given zone, the variation of the mesh size in the physical plane results only from the coordinate transformation; this variation can be kept small enough so that the accuracy of the scheme is practically the same as with a uniform mesh. Having thus determined the complete mesh in the X, Y plane, in order to apply the finite-difference scbeme (2.4) we only need to know the cartesian or cylindrical coordinates xk., y~ of the mesh points. These can be calculated explicitly through eqns (2.9), (2.6) and the deflm'tion of the eurvilinear system. A matching procedure must be set up to link the solutions in two adjacent zones. For this, it is convenient to make adjacent zones overlap as shown on Fig. 2 where horizontal lines in the same column represent the mesh lines Y = YI" = YI'~+(I - DAF t", 1 = I to L , in zone "r". The numerical domain of dependence of a two-step scheme such as the scheme (2.4), at a mesh

Numerical studies in high Reynolds' number aerodynamics

l : I-r

(

4

I

3

7

l-I Zone (r+ll ,

•

~

Zone

2

r

zone (r.1)

Fig. 2. Mesh refinement by dichotomy, and matching of adjacent zones.

point M, is centered about M and spreads over two meshes on each side. Therefore, considering each zone independently of the others, the two lower and the two upper mesh lines in zone "r" must be considered as boundary lines, and the solution can be advanced in time from t to t + At tr) only on the lines with index 1 such that 3 ~
U~"= U~),)-2, _

U ~ = U~"+')

(2.10a)

For the value of U1') on the line I = L, - I, an interpolation is necessary; we use the following 4-point relation:

U~_, =11--{6-UI'+')+9[U~'+I)+ U~"+')] - U(4'+') }

(2.10.b)

For steady flow calculations, the method need not be time-accurate (or consistent in time), and it is not necessary to match values at the same time level. This freedom allows many possibilities concerning the values of the time step and the respective numbers of iterations in the different zones. Thus, with an explicit scheme, the time step in each zone is restricted to some maximum value, and it is possible to advance the solution in time in each zone with the maximum time step allowed in this zone. However, to avoid possible convergence problems, especially when the calculation is started from arbitrary initial conditions, it is recommended that the number of iterations in a zone be approximately inversely proportional to the time step used, so that the time levels in the different zones remain approximately the same. For example, w e c a n use" At (,) =

2'At.) ]

N (r) = 2R-r

i

(2.11)

where N (') is the number of iterations carried out in zone " r " for one iteration carried out in the outer zone "R". Such a variation of time step corresponds approximately to the maximum value in each zone if the stability is governed by an inviscid C.F.L. type criterion. A systematic procedure is easily set up to sweep the regions in some definite way, including the matching of adjacent regions. Such a procedure is illustrated on Fig. 3 in the case when the relations (2.11) are used. On this figure, the dots represent the time levels (measured along vertical axis at which the solution in a zone is calculated, each zone corresponding to a vertical line. The set of all the iterations carried out in the different zones for one iteration in the outer zone (r = R) is called a cycle; Fig. 3 shows a complete cycle for R = 5. The order according to which the time levels in the various zones change, is obtained by following the line connecting all the dots, starting from zone 1. In general, the resulting method will not be consistent in time with the unsteady equa'tions, unless special care is taken in the application of the matching

8

J.C. Lz BALLEURet al. ttme level

At($ )

Ae" ~ ^ , " )

: !~"

'i

Fig. 3. Cycle of time iterations for a mesh with 5 zones.

relations (2.10). Thus the simple rule of applying these matching relations as soon as new values are obtained, is not time-accurate, but it is adequate if we are interested only in the steady-state solution. An interesting feature of the present method is that the mesh refinement in Y by the technique of dichotomy can easily be changed at any time during a calculation, the basic mesh system being kept fixed. Indeed the mesh structure is completely defined by the number of zones, R, and the number of lines, L,, in each zone (r = I to R), and the corresponding between the mesh lines of the new structure and those of the previous structure is easily established in a quite general way. Through this correspondence, the numerical solution can be transferred from the previous mesh structure to the new one without loss of accuracy, the solution on the new mesh lines (if there are some in the new mesh structure) being determined by interpolation techniques. This procedure has been made automatic and incorporated into the computer code as a subroutine. However, in the present stage of development of the mesh adaptation technique, the choice of the parameters R and L, (r = 1 to R) is not yet made automatically, and it requires the user's intervention. 2.3 Numerical applications We present three numerical applications obtained with the method described above. All three applications make use of the scheme (2.4) and of the technique of mesh refinement by dichotomy, but the technique of adaptation of the basic mesh system to the edge of the viscous layer has been implemented only in the third application concerning a nozzle flow. These applications are relative to laminar and two-dimensional (plane or axisymmetric) flows over adiabatic walls. The fluid is a perfect gas with 3' = 1.4, and the viscosity is calculated according to Sutberland's law. 2.3.1 Supersonic ]low in an axisymmetric compression corner. We consider the flow past an axisymmetric body made of a hemisphere cylinder extended by a conical flare. The flare angle is/3 = 8", and the cylindrical part is two nose-radius long. The free-stream conditions are: Mach number Mz= 2.94; Reynolds' number, based on nose radius R, Re~.R = 2.2 x 105, and the Prandtl number is Pr = 0.723. The flow field is divided into an upstream region which includes the spherical nose and a part of the cylinder, and the corner region. Only the flow in the corner region is considered here; the flow in the upstream region, which is not influenced by the corner, has been calculated previously, [21], and the results are used to determine all the flow conditions on the UPstream boundary of the corner region. The computation domain is shown on Fig. 4 with the basic mesh which is made of (30 x 21) points; the variable ~ is the distance along the wall (referred to the nose radius) measured from the stagnation point. The bow shock is fitted as the outer boundary # = ~E(~:; t), where # is simply the radial distance. The Rankine-Hugoniot relations, together with the value of the pressure behind the shock obtained from the discretized equations, are used to determine the shock velocity and the flow properties behind the shock, as described in Ref. [21].

Numerical studies in high Reynolds' number aerodynamics fl J:21

t~ i I

/ I 1 ,I

t t 11f t if I 11

,t I" I I I,

I ,"r~"lY i

~

~=2.9,~ ~J

/

. J , . ~ T J...,~~ "

.-"-;i=i

5,

~

j~:t9 2 3 0 I

/

2.83

-5~/=I

~,5; ~==-~-, 3.,;11 I

I

/

1

. . . .

~"~

5.80

4.uo

I

I

I

0

2

xlR

3

Fig. 4. Axisymmetric compression corner flow. Computationdomain and basic mesh.

The mesh structure obtained by the method of dichotomy involves 6 zones, and the number Y = const, in zone r is: LI = 29, L2 = 5, L3 = 7,/-,4 = 5, L~ = 7, L6 = 19. Zone 6 is a part of the basic mesh shown on Fig. 4. Profiles of the axial velocity component u and of the pressure in the corner region in the vicinity of the wall (zone 1) are shown on Figs. 5 and 6. A small separated flow region is seen to form. Isobar and iso-Mach contour maps are presented on Figs. 7 and 8. Because of the smallness of the separated flow region, only one shock is present. Experiments, [22], show that a large separated zone exists for/3 = 10°. Attempts to calculate this case have led to convergence difficulties which have not yet been solved, even when using

Lr of mesh lines

i=14

0

0.2

0

02

t i=15 , i=16 ,

0

Q2

0

0.2

0

:.17

i=t8

Q20'

Q2

i=19

Fig. 5. Axisymmetric compression comer flow. Axial velocity component profiles.

297~ l (zone 1)

I

1.20

(

.

f I

--.

-

1.25

Fig. 6. Axisymmetric compression comer flow. Pressure profiles. CAF Vol. 8, No. I--B

I0

J.C. LE BAI.LEURe t

aL

4-

/ "'" 2.2 , •

"

.

/

.

/

/

/ " / / / / ~ / /

/

• I/X~.."

"

" /

J /./

" /

,~

/

j.J

....

15

.

]

.

~

.f J09

/ / /

xlR

Fig. 7. ~,~xisyrnmelric compression corner ~ow. Isobar contour map.

2J- I

zy

/
3~.

2.2/'~ ~ / /

/

//

,/

,/

../

/'J ~

1/

/ " "

~

I

2

''

3

0.2

I

I

4

5

Fig. 8. Axisyrametric compression comer flow. Iso-Mach number contour map.

a twice smaller mesh size in the ~-direction. It may be necessary to use a much smaller mesh size in ~ in order for the calculation to describe correctly the strong pressure gradients at the separation and reattachment points and to produce the large displacements of these points from the corner region to their final positions at convergence. The artificial viscosity term also plays a crucial role in the convergence process, since its efficiency in damping out the oscillations should not interfere with the real viscous effects. 2.3.2 Shock-boundary layer interaction on a flat plate. This problem is treated for the flow conditions which correspond to Hakkinen's experimental data [23]: M= = 2, Re .... = 2.96 x I(P

where x, is the distance from the leading-edge of the point of impact of the incident shock on

Numerical studies in high Reynolds' number aerodynamics

11

the plate assuming inviscid flow. The Prandtl number is Pr = 0.72. The strength of the incident "shock is such that Pi/P= = 1.4, where Pi is the theoretical pressure downstream of the reflected shock in an inviscid flow. First the flow past the flat plate is calculated in a domain extending downstream to x/x, = 0.7, and outwardly in the free-stream. The incident shock is imposed on the upper boundary as a discontinuity in the free-stream conditions. The leading-edge shock as well as the incident shock are captured. The mesh was approximately adapted a priori to the leading-edge shock and to the boundary layer edge, using a transformation of the type (2.6) (where g ~ are simply the cartesian coordinates x, y), but where the upper boundary 71 = "OEis not imposed but is chosen in the free-stream region in such a way that not only the boundary-layer edge, but also the leading-edge shock ~---¢15h follow approximately mesh lines 7j = const. In this first calculation, the incident shock exits the computation domain on the downstream boundary without interacting with the boundary layer. The interaction of the incident shock with the boundary layer is studied in a domain, shown on Fig. 9, extending approximately from x/x$ = 0.5 to x/xs = 1.5. All the flow conditions on the upstream boundary (including the incident shock) are known from the previous calculation. The upper boundary of the computation domain is located in the free-stream beyond the leadingedge shock and downstream of the incident shock: in fact, the upper boundary could have been located below the leading-edge shock, but the gain in computing cost would not have been appreciable and the treatment of the boundary conditions would not have been as simple as with the known free-stream conditions. The basic mesh, made of 36 x 27 points, is shown on Fig. 9. As for the calculation in the upstream region, the adaptation of the mesh to the leading-edge shock and to the boundary layer'edge was made a priori. Furthermore some mesh refinement in y is introduced in the basic mesh by a stretching function of the type (2.9.b). As a consequence, it is sufficient to use three zones to refine the mesh further near the wall by the method of dichotomy. The number of mesh lines in each zone is: L~ = 17, L: = 13, L3 = 20, zone 3 being a part of the basic mesh. Figures 10 and I1 show isobar and iso-density contour maps, which exhibit the separation shock and the reattachment shock. The streamlines in the separated flow region are shown on Fig. 12, and the asymmetry of the separation bubble can be noticed. Figure 13 shows the distributions of pressure (referred to the free-stream stagnation pressure Pi~) and of tangential shear stress (referred to o=U~2) at the wall. The numerical results are compared to the experimental data of Hakkinen et al. [23]. The discrepancy can be reduced by a general shift of the numerical results in the x-direction: this can be explained by a shift of the computed incident shock due to the numerical errors in the region of intersection of the

y/x,

0.,5...,.

M=2"~, J_

0

I

'"

"

0.5

1

15

Fig. 9. Shock-boundary layer interaction on a fiat plate. Computation domain and basic mesh.

12

J.C. LE BALLEURet al.

025i

~

,1.20

5

12 R

Fig. 10. Shock-boundary layer interaction on aflat plate. Isobar contour map.

ylxs

0.50-

j l . . . . L - I " ~ / . ~

1.;'5

0.25.

O.5

t

5

,q

1.5

Fig. !I. Shock-boundary layer interaction on a fiat plate. [so-density contour map.

ylxj

002 i I I i

i aotJ

0

00

09

I

1.I

Fig. 12. Shock-boundary layer interaction on a flat piate. Streamlines in the separated flow region,

leading-edge shock with the incident shock which are both smeared by the method of shock-capturing. Also the pressure level at the upstream boundary, which is a result of the calculation in the upstream region, is slightly above the experimental value. 2.3.3 Transonic flow in an axisymmetric nozzle. We consider now the subsonic-supersonic flow in an axisymmetric nozzle made of a conical convergent of 45° ~ e and a conical divergent of 15° angle, joined together by a surface with a circular meridian section of radius R' such that R'/R = 0.625, where the reference length R is the nozzle throat radius. The Reynolds'

Numerical studies in high Reynolds' number aerodynamics

0.2

13

p/p .

•

•

.

.

•

• e x p e r i m e n t (Nakkinen)

2 ~ Ct.IO ~ o

°*~°

. present

.

calculation

I• +

°

• ..:.. _1

•

.

... x/X8

". I 1

I O.5

.,I 1.5

Fig. 13. Shock-boandary layer interaction on a flat plate. Distributions of pressure and shear-stress at the wall.

number, based on R and on reservoir conditions (sound speed co, density po, temperature To) is R e ~ = 104.

The computation domain and the basic mesh made of 41 x 21 points are shown on Fig. 14. The distortion of the basic mesh is the result of the adaptation of the mesh to the edge of the viscous layer according to the method of Section 2.2. Here ~, n can be taken to be simply the cylindrical coordinates x.y. The function g~ in eqn (2.7) is chosen parabolic in 4, i.e.: 1

s' = p(l - p) ~[x- 13~ + (f3 - x)~].

The parameter X is found to vary from 0.015 to 0.1; the constant ~ is taken equal to 0.075, that is to say that the approximate edge of the viscous layer is the line ~ = 0.075. The final coordinates X, Y are the distance.along the nozzle wall for X and Y = ~. The mesh refinement by dichotomy is obtained through 6 zones, and the mesh structure is defined by: L~ = 1I,/.,2 = 8, L3 = 10, L4 = 5, L5 = 5, ~ = 20. The artificial viscosity can be set equal to zero for this problem. Initial conditions are obtained from an inviscid flow calculation using a pseudo-unsteady potential method [16]. In this inviscid calculation, the flow direction is imposed on the upstream boundary as being that of a conical flow, but the modulus of the velocity results from the calculation. Then, in the viscous flow calculation, the no-slip condition is established gradually

•

y/R 45 °

IS"

17 ..~1

I~,

AI

I I I

l i ]/

d

1

I II

1

T

~=1

II

,

oT 21

3?

47

Fig. 14. Transonic flow in an axisymmetric nozzle. Computation domain and basic mesh.

14

J.C. LE BALLEURet

uL

by letting the wall velocity decrease slowly to zero, while the adiabatic wall condition is imposed from the beginning. Of course the mesh adaptation technique can be started only after that the viscous layer has sufficiently formed. Concerning the conditions in the viscous layer on the upstream boundary, calculations have shown that the influence of these conditions is felt only a short distance downstream, so that the viscous flow at i = 10 is practically independent of them. Thus it is possible to assume no upstream boundary-layer, i.e. more precisely to impose the inviscid flow solution at all mesh points on this boundary, except at the wall; however with these upstream conditions, the flow gradients in X are large and some oscillations occur near the upstream boundary. We found that better results could be obtained by imposing an upstream boundary layer, the thickness of which was determined after a preliminary calculation by extrapolating smoothly upstream the values of X determined from the numerical solution for i ~> 10, and by using the mesh adaptation technique up to i = 1. The total mass flow rate through the nozzle is not imposed, but the flow conditions on the upstream boundary are re-considered from time to time in order to let the mass flow rate adjust itself freely. This is achieved by using a technique similar to that of Lavai [24] for inviscid flows, in which the upstream flow conditions are modified so as to make the entering mass flow rate equal to the mass flow rate calculated in a section in the supersonic region. In this modification, the flow direction is maintained fixed. Figures 15-17 show profiles of Mach number M, of total enthalpy H (referred to co2) and of temperature T (referred to the reservoir temperature To), in the viscous layer at four stations. The mass flow rate (referred to pocoR2) is found equal to 1.686, and to 1.788 in the case of an inviscid flow; we see that the influence of the viscous effects on the flow rate is not negligible at this Reynolds' number. Finally it can be mentioned that practically the same results were obtained for this nozzle flow by using a thin layer approximation of the discretized N.S. equations. In this approximation, the calculation of the dissipative terms (and only them) makes use of finitedifference operators 8,v, 8~, derived from 8x, 8y (eqns (2.3)) by neglecting qb,o in the numerator if we assume that MP belongs to an X = const, line and MQ to a Y = const, line. In the present application, we have xMe ~ 0 and we get simply 8xv¢ = 0 and ,Syvqb= 8Jl, (~MP/YMP. =

i

1

e,i

~--~.;

-005

o39,

082

j, M 0

2

2

0

Fig. 15. Transonic flow in an axisymmetric nozzle. Mach number profiles.

,y. l

,.!

2. t

!\

\

_)_. 25

2.5

25

2.5

Fig. 16. Transonic flow in an axisymmetric nozzle. Total enthalpy profiles.

Numerical studies in high Reynolds' number aerodynamics

15

02.

i= fO] =-0,4[

01

I~ -oo~

24 Z3.I

34 082

I

,

r/), o --

o~

~

o9

7

o9

7

o:9

- i

Fig. 17. Transonic flow in an axisymmetric nozzle. Temperature profiles.

3. COUPLINGAPPROACH

3.1 Interacting thin shear layer approximations for finite Reynolds" numbers Many coupling approaches have been investigated, see for example [25]. Here, we develop an interacting layer method, including reverse flow bubbles and shock wave-boundary layer interactions, without requirement of consistency for infinite Reynolds' number R. However, the interacting Prandtl equations, and then the following extended approach, are known to be consistent, not only with boundary layer theories, but also with several strong interaction asymptotic analysis, as for example the laminar triple deck, see [26, 27].

3.1.1 Interacting thin shear layer: Patching approach. The idea, first suggested by CroccoLees [28], is to generalize the Prandtl assumptions of a viscous layer which is effectively relevant of equations without normal pressure gradient, parabolic in the streamwise direction. The non-interacting problem is unclosed without an "external" boundary condition, and for long time the closure has been the pressure, as for weak interaction asymptotic theories. One knows then the Goidstein square-root singular behaviour of the solutions near the point of zero skin friction, also reproduced by integral methods [27]. It has been now proved that such a singularity has not to be associated to a breakdown of the approach, or to a detachment of the layer from the wall, but only to the choice of pressure as closure. If the pressure is free, for example with one another boundary condition for closure, or with an inviscid interacting boundary condition, regular and realistic solutions with reverse flow are obtained [29, 30]. However, the interacting Prandtl layer problem is no more parabolic, since upstream influence is possible by way of an unknown closure boundary condition, the interacting upstream pressure evolution. A schematic representation of this ellipticity is given on Fig. 18, which summarizes the concept of viscid-inviscid patching: we have to define a free surface, which is a conventional edge y = 8(x) of a thin layer, on each side of which are solved equations that are locally representative of the flow. Normal pressure gradient is then effectively absent of the shear layer. Along the free surface, the aerodynamic quantities are continuous, their partial derivatives may be discontinuous, and mass transfer is present from the inviscid flow to the

IN VI3CID FL 0 W

Up3tream mf/uence by boundary conditions

Patch/ho

R LAYE

Jr

Fig. 18. Interacting thin shear layer: Patching. o,

..

16

J.C. LE BALLEURet al.

shear layer. The patching relation is the inviscid closure boundary condition along the free surface, issued from the continuity equation. In the following, we consider the usual curvilinear intrinsic x. y along the wall, with the density and velocity components denoted O, u. v for the inviscid flow, t~, a. ~ for the viscous flow. Then, if the wall curvature is negligible, we have: u ¢.,.~)= dx

(6 - 6") P

_1~.~)

(3.1)

with

{

8(*~[pultx.8) = p, u :inviscid

Y°

{[pu](~8)- [/~a]<~,} dy t~, a : viscous

Such a 6*(x) displacement formulation would not involve any approximation, if it was the same for 6a(x, y). A fundamental question of this patching approach is the following: is the ellipticity of the steady Navier--Stokes equations also recovered when the inviseid flow is supersonic, and that the upstream feed-back of influence may not be a purely inviscid one? From the initial work of Crocco-Lees, an important investigation has been done, see[25,27,31-33], the answer being positive for a so-called "suberiticar' layer, negative for a "supercritieal" one. In the subcritical ease, the coupling of an inviseid hyperbolic problem with a parabolic shear layer, leads to an initial values problem which is ill-posed, because a divergent growing of initial perturbations is possible through branching solutions, see for example [34]. This ill-posed nature is removed if one boundary condition is moved from upstream to downstream, which insures the ellipticity. Furthermore, the upstream influence of the downstream condition decreases so quickly, that it may generally be simply added as an extra downstream-condition. In the supercritical case, only a pseudo-ellipticity may be involved by way of weak discontinuous solutions. Pressure jumps are then present at the wall to insure sudden changes of the supercritical shear layers into subcritical ones, when upstream influence is necessary, see [35]. Such non-physical jumps have raised controversy about their interpretation, about the existence of supercritieal phenomena in non-integral solutions of the shear layers. We summarize here some new elements, detailed in Ref. [27], that lead us toward an improved interacting approach, presented in the next section. It has been early proved that supereritical questions are consequent to the simplified pressure field/~(x) = p(x, 8) of the shear layer. With an integral method, equation (3.1) gives a linear form between flow direction O = (v/u) and pressure gradient at y = 8(x): D, O = ~"~*dd~x + D3

(3.2)

D. is zero at separation or reattachment, where the pressure gradient is not trivial, but D2 is zero at critical points, where O = ( / ~ D i ) is needed. The linear form (3,2), or simply:

8*(dp/dx) = - (D3/Dz) to insure ~ a r i t y .

ctn

O = B-~r- + C (IX

D,

B = -'8" D1

(3.3)

is determinant for the branching solutions, stable with x if B < 0, unstable if B > 0. For an attached layer, the unstable case may be said subcritical, or elliptical. At the separation point however,/91 changes in sign, as consequently B, so that the suberitical case involves a stable branching in the reverse flow regions. There, ellipticity is then consequent to the regular reattachment condition #*(dp/dx) -- -(D31Dz), adjustable with a possible moving of the separation point, as allowed by unstable branching upstream of it. Finally, Fig. 19 gives a schematic summary of the total suberitical ellipticity.

Numerical studies in high Reynolds'number aerodynamics Pressure

"

Elliptieity of

~'~

//

$tmble

I

.

R

"

! 1 re~ttsch~nt- ax'--

Unstable br#nchmg Separ#tion

17

I

I

xj

Re#ttachment

Fig. 19. Supersonic separation over a compression ramp: Scheme of branching solutions for an initially subcritical layer.

With local Prandtl equations, a relation in analogy with (3.3) subsists [32], involving the local viscous Mach number M: B = lira ~ _----7-=~_ _ d .v ~-.*0 J,

ypM ~

(3.4)

For symmetrical wakes, the limit e ~ 0 is not singular, and the far flow, with ~r > 1 everywhere and B < 0, shows that supercritical patching does exist with local Prandtl equations. For boundary layers, a conclusion is not so easy because the limit is singular, and that divergent integrals B and C only cancel each other in relation (3.3). In the turbulent case however, available asymptotic theories [36, 37] involve a rapid growing of M inside of sublayers, where the flow deflexion from the wall is negligible. It is then possible that the limit e-* 0 has not to be associated to a vanishing Mach number at the wall, and that relation (3.4) would have an asymptotic meaning. If for laminar boundary layer, at least, the wall region is always determinant of the deflexion O, see [33], the branching question however requires only to look at a small perturbation from a given solution O0, dp0],

77j

and to write linearly:

(3.5) where B1 is the contribution from the subsonic sublayer and B: from the supersonic outer region. If an analytical expression for BI is not available, all numerical works show that B~ > 0, and that the supercritical branching B < 0 results from a balance with B2 < 0, given by relation (3.4) as: B: =

~l - ~ ¢ : , ~~ ay

fy

y,: sonic line

(3.6)

A determinant conclusion is that any supersonic boundary layer patching, which is subcritical because [B21
18

J. C, LE BALLEURet aL

turbulent supersonic separation over an adiabatic flat plate at increasing Mach numbers. subcritical configurations are obtained until approximately M = 1.2. If we look at the inviscid flow, and no more at the wall pressure, the continuous compression waves induced at the edge 6 coalesce at some distance into a separation shock-wave. As the Much number increases until M = 1.2, also the steepness of the branching compression becomes higher, and the coalescence point moves toward the 8 edge, to reach it at M = 1.2 (Fig. 20). From an inviscid point of view. the appearance of pressure jumps at the edge 6, for the supercritical range M > 1.2, may be interpreted as consequent to the penetration of the separation shock-wave inside of the shear layer of thickness & The unsatisfying aspect of patching is then purely a viscous one. because the lack of normal gradient does not allow the pressure to spread inside of the layer, and requires that global conservation laws with jumps will replace locally the Prandtl equations, tn our opinion however, caution is needed to implement conservation laws for the global kinetic energy, as usually done with integral methods, simply because shock-waves inside of the shear layer may involve infinite dissipation rates. A conservation of the 8 patching thickness will seem physically more convenient. At last, regardless of the jump relations, it has to be kept in mind that this trouble appears at low Much number for turbulent flows, M = 1.2, which excludes for example to build a general patching theory for transonic airfoils without a satisfactory answer to those questions. 3.1.2 Interacting thin shear layer: Matching approach. The patching approach is not, however, the only way to introduce the viscid-inviscid interaction [25-27, 38]. For example, the inviscid computation looks easier, either if there is no mass flow along its boundary, and one speaks then of displacement concept, either if the inviscid boundary is no more a free surface, but simply the wall. (a)

Approximate coupling relations

A first insight of such an extension of the inviscid flow inside of the shear layers is to consider an analytical report of inviscid boundary conditions from 8 toward the wall y = O: v

v

(3.7)

which, associated to relation (3.1), gives a slip condition if Y = 8', a wall mass flow if Y = 0, and corresponding displacements 8*:

~,8.~

dx + . . . .

(3.8)

[ pv l,~.o)= d [paa*l + ....

(3.9)

,.,.,Ivul L

J(x.Y)

J

C°maress;~/

6

6

,d u m p ~ Pl'aMtl eqe,ttiene / f /

/ / / f P ' f

l f

,

XL

J pr~atl ~qu~t~ns " f f f ~ f

lf

f

f

f

l

X

Fig. 20. Supersonic separation with the patching model: (a) Subcr/tical (b) Supercritical.

Numerical studies in high Reynolds' number aerodynamics

19

The mathematical equivalence of relations (3.1), (3.8)--(3.10) is however approximate, and for example caution is needed for an asymptotic meaning. Consistency as R ~ oo is insured with the second order theory, in which only first order terms (subscript 0) are then needed for 8": 8*~)[poUo]~.o~= R -"/~)

Jo

y = R-tll2)y

{[pou0]~x.0)-[Aao],~.~} dy

(3. l 1)

We have assessed that relations (3.8), (3.9) are also consistent with the asymptotic triple deck theory [27]. Despite this, a report of boundary conditions at the wall has large consequences on the patching approach, and represents in fact a first step toward a general analysis of viscous flows, where the thin shear layer approximation never induces supercritical behaviours, nor pressure jumps at the wall.

(b) Present non-asymptotic matching Here, the viscous and inviscid computation fields are fully overlapped, as in asymptotic theories, and we speak of a "matching approach", for finite Reynolds' numbers. The viscous solution, summarized on Fig. 21 by a(x, y), is split into three parts. The first one is an inviscid predictor, which is supposed to be representative of the far field (y > 8), by adjusting a non-zero normal velocity at the wall. The second part is a full viscous investigation. However, it is required to formulate it as an error analysis to the inviscid approximation u(x, y). Such an error may be defined for example in a multiplicative mode with ¢(x, y) as viscous variable instead of a(x, y): a(x, y) = u(x, y). ¢(x, y) (3.12) The error has only to vanish at large y, (with ~, = 1) and requires no more a conventional viscous edge 6. In the viscous layer, the error is important for some quantities, as the tangential velocity [u(x, y) - a(x, y)], rather small for others as the pressure [p(x, y) - p(x, y)]. We have thirdly to couple parts I and II, without weak interaction approximations, which looks like a reciprocal question: know the inviscid field u(x, y), to estimate the error g,(x, y) from a viscous computation, then adjust the normal inviscid velocity v(x, O) to vanish the error at large y. A schematic representation of the matched viscous and inviscid solutions is given on Fig. 22. The matching condition on v(x, O) may be expressed with a generalized displacement 8*(x): d

(3.13)

8*~[pu]c~,o)

=

f= {[pu]~,,) - [pt~]~,,)}

dy

Jo

Those relations are very much the same as relations (3.9), (3.10) but here, there is no approximation, only integration of the continuity equation both for viscous and inviscid fields. v i s c o u s SOL u.o sPurr Nc

/NV/SC/D FLOW ~ Over/appi'n--eo LVISCOUS CORRECTOR.~ PREDICTOR-u{x,y) | fields |5(x,y~u(x,y).~(x y) t I

"1"

y--. -

velocffy:v (xjo )

.(x,y)--1

J

-

at the wag

FULL COUPLING~.. Comput,tmn of v(x~l I Problems Z and ~ are strongly matched {no approxtmatmn of weak ~ter~ctwn Fig. 21. Interacting thin shear layer: Matching at finite Reynolds numbers.

20

J. C. LE BALLEURet al.

/ 1

Invi~ u

)

Invl$Clff L .lhock-wer¢~.

~Wall

~,

.~fpressore 9r~nt :

~

u L

Fig. 22. Schemeof viscous-inviscid matching at a given x-station.

Allowance to a non-approximate pressure field is present in 8*(x) by y-variations of viscous term ~ , and inviscid pu. Such a matching would be convenient to couple inviscid computations and Navier--Stokes solutions, with benefit from assuming the shock-wave capturing qu~,tion only into the inviscid part. Presently approximations are considered, for the viscous estimation to be as analytical as possible, and require both low computing times and simple turbulence modelling. Considering only viscous flows with a dominant direction, thin shear layer approximations are introduced. Then, if strong normal pressure gradients are still possible in the viscous layers, the pressure field becomes only slightly different from the inviscid-one. Defining now generalized defect thickness in the same way as for 8*(x), it is possible to extend formally the classical integral methods used in the patching approach, based for adibatic flows on integral momentum, kinetic energy or entrainment rate. F'mally, the viscous numerical computation reduces to global differential equations, involving two independent viscous defect thicknesses, and the inviscid quantities along the edges of the computation field, walls or center-lines of wakes, where x-coordinate is measured:

dS* + 8* dx

apu

pu ax

=v u

Continuity (Matching)

8*+2Oau +O~ _ = Cz

dO ¢ dx

I

d8

dS*+878*OpU=AE

dx

dx

dO*

u

dx+

ax pu

pax

x-Momentum eq.

2

Entrainment eq. if turbulent, with )t = )t(x)

ax

0* 0ou 3 _ 8 * - 8 * a u = ; t O + O p

X

ax

Kinetic energy eq. if laminar, with )t = 1, Op = 0

[8"~ + O~,J[puZ],.o) = J o {[PuZ]~"- [Pa2]"} dy with

+ O*~][pu3]~.o) =/~ {[pu~]~.., - [~a3]~.,} dy 8i(~)U(~.o)=

*

f:

{u(~.y)- a(~.,} dy

, ~ l r a,H a y -I A(x~E(x):entrainment rate at~ = p-'ff[_0(u- a)layJ,,.~)

2 C aa,

)t~x~O~: integral dissipation rate )t@ = ~ ~p(x>: extra term, negligible if

J0 ~"~yoy =0

(3.14)

Numericalstudiesin highReynolds'numberaerodynamics

21

Additional differential equations restitute a simple x-dependent lag-modelling for fluctuations and turbulent stresses, which influence the mean integral equations in A(x), with A = I for equilibrium flows. A detailed investigation of this viscous integral closure, including reverse flow, may not be given here. Two points will only be noticed. Firstly, when large normal pressure gradient occurs, as for example near a turbulent supersonic separation, the generalized equation for kinetic energy is no more formally the same as in Prandtl integral theories. An extra term Op has to be added to the global dissipation rate to deal with non-vanishing normal variations of (@/Ox). Such a complex implementing is not required to generalize the entrainment equation, which does not integrate informations across the whole layer, and is rather like a collocation with the local momentum equation near the outside edge 8. For turbulent flows, an entrainment equation is then used here rather than a kinetic energy one, not mainly for viscous closure improvement, but because of matching considerations. The highest interest of the matching approach is lastly assessed, for thin shear layer approximations, with the question of supersonic branching and ellipticity. Shortly, we found [27] that the investigation given for the patching approach may be approximately restarted if integral B' is substituted to the integral B of relation (3.4), with:

M(y): inviscid

B' = lim f~ M2 - " ~ " • -.O'E y p - - ~ - ~ t l y

M(y): viscous

(3.15)

A first important improvement is that the conventional 8 thickness has no more influence. Secondly, the term (1 - f F ) of relation (3.4), which generates negative values of B2 in relation (3.6), has been replaced by (M 2- ~2), essentially positive because of viscosity. Consequently, a total removal of supercritical branching behaviours may be waited for, and has been observed with the present integral solution. When a displacement concept is applied, which is an hybrid method between patching and matching, relation (3.15) becomes: B*=limLj ° ~oyj+j8 I

f=M

~dy

(3.16)

For largely supersonic inviscid Mach numbers, the local viscous Mach number M may still induce negative values of (1 - Jr/:) where y < 8", so that a supercritical branching is possible, as in the patching approach, especially for turbulent layers. Then, the matching approach has always to be selected, not only because the inviscid flows extended to the walls are easier to compute than with free boundaries, but firstly because the removal df supercritical branching leads always to a full ellipticity, as sketched on Fig. 19, and eliminate the jump question. Consequently, never a shock-wave has to reach the inviscid computing boundaries, walls or wakes center-lines, where the compressions are always spread over a length scaled by the local shear layer thickness. The shock-waves associated to a turbulent separation appears inside the computation field from a coalescence of compression waves, but the matching approach allows this coalescence may stand inside of the shear layer, see Fig. 23. The inviscid pressure

Shoe

I •

ii

i$~~""

_

]

R

"I

pressurep(x,o)

\Viscous wall presSUex~(X,O)

Fig. 23. Matching approach: Inviscid sketch of a supersonic turbulent separation. Pressure at the wMl.

22

J.C. LEBALLEURet al.

evolution at the wall, always continuous, issued from the present matching, is generally a good prediction of the viscous pressure. However, when a large inviscid streamline curvature is induced by the viscous interaction, as for supersonic turbulent separation, an approximate correction is finally added. We suppose for this that the inviscid pressure ptx, y) and the viscous one/~(x, y) satisfy inside of the shear layers: ap = K<~)pu2 cgy

__, ~ap = K~pu"

t3.17")

where K(x) is an averaged local streamline curvature. This leads to an integral correction: ptx.0~-/~x.0j = [8* + O]tx~K~x~[pue]~x.o~

(3.18)

3.2 lnvi$cid computation with boundary conditions of strong coupling to a thin shear layer The numerical solution of the previously outlined matching approach, using an integral viscous formulation, reduces then to an inviscid problem with known geometry, but unusual boundary conditions. The problem is not direct, with flow angle O prescribed, nor inverse with known pressure p. The boundary condition is a set of differential equations between p, O, and viscous quantities. With some elimination, the inviscid closure relation between p and O may be extracted from (3.14): AIO = A28*-~xx + A3

(3.19)

x is measured along the boundary. At, A2, A3, t$* depend on the remaining differential equations, on the initial viscous quantities, and on the inviscid evolutions Otx) and Ptx~. The equation (3.19) is highly nonlinear; At vanish at separations and reattachments; A2 is never zero because the matching eliminates the supercritical phenomena. The small scale 8" leads to branching solutions, whose importance has been discussed before, and which require to add extra-downstream conditions. Numerically, the matching is then very difficult, and may lead to solve unsteady or pseudo-unsteady equations [39], as for the global approach. Here, we consider only steady iterative methods where, either the inviscid solution is simple enough to repeat cycles based on exact viscous and inviscid solutions, either the inviscid solution uses a relaxation technique and provides yet an iteration loop where the iterative matching must be integrated. In this case, the matching looks like a relaxation for the boundary condition (3.19), and the increase of computing time required by viscosity is limited. For transonic or supersonic flows without separation, the first difficulty is to insure stability to the viscous-inviscid iteration, especially near shock.waves or incipient separation, and usually, questionable smoothings or arbitrary pre-selected underrelaxati0ns are added, see [25]. A second question is that poor treatments of the shock-waves boundary layer interactions are frequently involved, the compression at the wall being spread only by the inviscid artificial viscosity, without allowance to the non-linear interacting boundary layer evolution, as well to possible reverse flow bubbles. In these cases, non-conservative inviscid methods are often unduly preferred to more rigorous conservative solutions. A third question is to allow for separated bubbles. The thin shear layer approximation then requires "'inverse" techniques for viscid-inviscid iterations, in order that never, during the cycle, the pressure will be prescribed to the separating, separated or reattaching parts of the viscous layers. Without this, singular behaviours would be induced as for non-interacting computations. Those inverse techniques. where the viscous pressure is prescribed to the inviscid flow, are however difficult to apply when large regions with ordinary boundary layers are present, as it may be expected from the classical weak interaction theory. Then, mixed numerical techniques, direct or inverse according to regions, are to be elaborated. The new questions are those of an appropriate switching rule, if possible not predetermined, and of an inviscid problem that involves mixed boundary values, whose compatibility at juncture is not insured during the iterations.

Numerical studiesin high Reynolds'numberaerodynamics

23

For supersonic flows, the apparent simplification of marching techniques, involving a simultaneous solution of the boundary condition (3.19), has been used in the past, especially with simple-wave flows. We know however that the branching induces an ill-posed problem, which requires a trial and error iteration to satisfy a downstream condition. Practically, such a resolution is difficult, and almost impossible as soon as the beginning of the strong interaction region is not known before, or that successive interactions has to happen, as it is generally. Now, a new method [26, 40] makes it possible to prescribe directly a selected downstream condition, and the whole upstream influence, including eventually multiple separated bubbles, is found as a result, in the same way as in the global approach. The present method, detailed in Ref. [41], involves the same progress for supersonic range, and more, is relevant of general supersonic flows with two waves families. For transonic or subsonic flows, the present numerical technique avoids the previously cited limitations, and provides a fully stable and automatized procedure. 3.2.1 lterative numerical method (stead), flow) Looking at relation (3.19) as to a viscous solution, successive viscous and inviscid solutions are considered. They are said direct or inverse, according to the inviscid boundary condition prescribed at cycle (n + 1): Direct:

Inverse:

On+' = ~-~l~'[~x]" +,43 AI

(3.20)

['dp] n÷l Aj -n A3 L~J = 8"A2 (~ " 8'A2

(3.21)

For the direct technique, considering the inviscid flow is an operator associating a pressure pn(x) to a known angle On(x), and adding relation (3.20), a matching-loop operator F is obtained: O"+1 ix) = F[O~]

(3.22)

In the same way, when At may be zero, relation (3.20) has to be avoided, and with an integration of relation (3.21), an inverse matching loop may be involved for the pressure: p~'~J = G[p~]

(3.23)

(a) Direct iteration (small values of 8*/0) Unseparated viscous layers are considered here, but the viscous interaction may yet be strong, since no approximation on the coupling is included in relation (3.22). The stability of the iteration is however not insured, and depends on the influence matrix which is numerically associated to the linearized operator F. For that reason, a different iteration, based on a linearized Newton technique for eqn (3.22), has also been suggested [42]. However, it needs an effective evaluation of influence matrices. In the present method, only an approximate analytical estimation of the eigenvalue of F is used, in order to add a computed underrelaxation (or overrelaxation) to the iteration (3.22), and to insure an automatized convergence. We have found, see [41], that harmonic distributions for a perturbation [O~x~-O~x~] are approximately eigenfunctions for the linearized operator F, and that the corresponding eigenvalues are: a/Jg* ~n = x/i-Z_ Me

M: local inviscid Mach number a; harmonic frequency

(3.24)

/~ is a viscous layer constant which may be analytically related to B' in relation (3.15), or to (A,/AO in relation (3.19). The highest eigenvalue is used to compute an underelaxation coefficient ~o(x), see [41], strongly dependent of x, as it is for/zr~ A small mesh size increases the highest possible frequency a and requires a smaller value of

J.C. LE BALLEURet ai.

24

to(x), the reference length being the layer thickness 8*(x). The transonic range requires also smaller o~(x), as well as separating layers for which AI decreases, and consequently i/~i increases. At separation or reattachment, A) is zero and/~, as ~o, become infinite, so that the stabilizing underrelaxation to(x) has to become vanishingly small. (b) Semi-inverse method (large values of 6"/0, including separation) This requirement to vanish to at separation enhances that the direct iteration is not well-conditioned for separating or separated regions. Reciprocally, the eigenvalues #~r associated to a linearized operator G in the inverse technique (3.23) are found to be approximately #zt =(1/#zo). Consequently, the inverse iteration is well-conditioned for separating-separated regions, and requires no underrelaxation at separation (to = 1). For weak interaction regions however, small values of the viscous term/t lead to small/~O and large/zs. This indicates that the inverse technique is there poorly-conditioned, requires small values for oJ, and that a direct iteration is there better. The conclusion seems that the best solution is to switch each other the direct iteration for @(x), and the inverse iteration for p(x), according to the local properties of the viscous layer, summarized in/~ or simply in the usual shape parameter 8*/0. In each case, the stability may be insured by computing an underrelaxation to(x) from the knowledge of tip(X) and tit (x). Such a technique is however very difficult to handle when repeated inverse regions are present, especially to implement it into a transonic relaxation method. The reason is that we have to solve an inviscid problem with mixed boundary values, and that, until the convergence is not reached, large discontinuities between viscous and inviscid pressure may occur at the inverse-direct switching nodes. Consequently, in the present approach, we have elaborated a new "semi-inverse" technique for separating-separated regions. In this way, the inviscid flow is everywhere solved as a direct problem, with flow angle @'(x) prescribed, which gives an inviscid predictor of the pressure gradient dpn/dx, see Fig. 24. A switch toward an inverse computation at separating regions is operated only for the viscous calculation, which gives there a viscous predictor d/in+l/dx. The semi-inverse method is an automatized and stable technique to correct On(x) with [0,+1_ @n](X), knowing the error F~x)= [(d/~n+i/dx)- (dpn/dx)], in order that E(x) vanish for large n. For detail, see [41]. Shortly, the first idea is that such a small correction also vanishing at large n, is relevant of a linear approximate analysis, which may be easily an harmonic one, because harmonic perturbations for @(x) or p(x) have been found to approximate locally the eigenfunctions of the matching operator F and G, and also those of the inviscid direct or inverse operators. The second idea is to use this linear harmonic analysis to simulate an inverse matching iteration, that it will be possible to stabilize by computing an underrelaxation. This inverse simulation answers to the simple question: how to perturb the inviscid flow direction O n (x), in order that the inviscid pressure gradient dp"/dx becomes the viscous one d/~"+l/dx? Finally, very simple corrections have to be operated independently in each computing node of the inviscid boundary: M
o, ,~ , ['d/~"+l dp"] dx - dx J

[On+i-O n]=C,(M,o , n , a ) [

(3.25) M>I:

[On+l - O"] = C2(M, 8", B, a) [-'~f-x - dx ~ J

The dependence of the functions Ci and 6"2 with respect to the harmonic possible frequencies is at last dropped out, by selecting the smallest of the corrections [O n+l- O"](~). (c) Supersonic implementations (Downwind.differencing for matching dp/dx) According to relation (3.24), the approximate eigenvalues ~o or #zs are real for subsonic flows, imaginary for supersonic ones. As a consequence, during the convergence of the underrdaxed iteration (3.22), perturbations are not only vanished, but also convected, see [41]. The interesting point is that the convection depends on the sign of the viscous constant/i, so that for an attached layer, convection is towards upstream for a subcritical layer, towards downstream for a supercritical-one. Then, as only the subcritical case has to be considered in

Numericalstudiesin highReynolds'numberaerodynamics [_~ INVISCID DIREcTFLOW

25

CORRECTION ='~(×) SIMULATION OF AN INVERSE INVISCID FLOW CALCULATION

L--~ BOUN:A~;R s:AYER~

=

"'J

ATION IUNDER-RELAX n4,1 n W ~ n I

Fig.24. New"semi-inverse"iterativemethod. the present matching approach, it is equivalent to say that the upstream influence of an extra-downstream condition is automatically insured by the iterative technique, provided that a cut-off will not be involved by the numerical differencing, as it may be for example with a full viscous and inviscid upwind-differencing. If uncentered schemes are to be used for the pressure dp/dx in relation (3.19), downwind-schemes are needed for a subcritical matching. This conclusion, associated to the usual upwind-differencing for supersonic inviscid flows, gives a numerical insight into the viscid-inviscid supersonic ellipticity. Presently, the viscous pressure gradient d#/dx is matched to the inviscid supersonic pressure p(x) through a partly downwind differencing: (0 < X < 1): ["~]i =(1 - x )~p,+J - +p,-12~Lx+xP'*I~'~Lx p

(3.26)

3.2.2 Examples of application. The present examples have been selected to elaborate the numerical technique and the viscous approximations. Further improvements are now being implemented, firstly to reduce the inviscid approximations, and secondly to extend the method to wakes and lifting airfoils. However, those first results are yet interesting with respect both to experiments and to global Navier-Stokes simulations. (a) Supersonic compression ramp (adiabatic) Here, the inviscid computation is reduced to a simple-wave Prandtl-Meyer flow, which may be rather approximate because of the turbulent separation or reattachment shock waves. The computation is initiated from a flat plate flow, where the wedge angle/3 is zero. During a few first iterations, /3 is rapidly increased until its final value. This one is not prescribed at the beginning, only to avoid the infinite pressure gradient at the wedge, before the viscous spreading has begun. A different initiation with prescribed/3 would however be possible. As an extra-downstream condition of weak-interaction recovery, zero pressure gradient dp/dx is prescribed at the last node, far enough downstream. A well-posed problem is then solved, where separation bubble and upstream influcence are automatically raised by the iterations, either direct or semi-inverse, according to the instantaneous local shape parameter of the viscous layer. Successive separated bubbles would be solved without any new implementation, for example with successive wedges. The example of Figs. 25 and 26 has been investigated with an interacting approach by WerleBertke[43], with Navier-Stokes simulations including turbulence modelling by ShangHankey [44], and experimentally by Law [45]. The wall is adiabatic, the upstream Mach number is 2.96, the angle /3 is 25°, the Reynolds' number at the wedge is 107. For this high Mach number, the turbulent boundary layer would induce supercritical jumps for any usual integral technique, as well with a patching approach, than with a displacement concept. The present matching approach provides a fully continuous pressure evolution at the wall, given on Fig. 25 after that the curvature viscous correction of relations (3.18) has been applied. The present computation supposes the layer is turbulent from the leading edge. Agreement with experiment CAF Vol, 8, No. ]--C

26

J.C. LE BALLEURet aL

i P 5~ ~o-o v

1

i

=

t"

0.9

, _ . " "~'- ~ - ' - "

i Inviscdi ! i

~:IS ,:)t °

./

0.95

x

'

1

1.05

• E x p e r i m e n t (L:,w) . + + + Present i n t e r a c t i n g .. I n t e r a c t i n g

!

"

I

""

1.10

method

method

of Wefle-Bertke

Equililyium Frozen turbulence

I Navier-Stokes (Shang.Hankey)

Fig. 25. Pressure along a supersonic ramp (turbulent).

o.oo2t °'031~C~

8 ..I"-..

o.oo~+ t-~ o.oI~ ~

~

.-.

. . . . :-:-'-"-

=

o.95st'~..-"Rv ,b5

,.,o

Fig. 26. Skin-friction and viscous thickness along the ramp of Fig. 25.

is believed very encouraging, as well as the middle position between two Navier--Stokes simulations with different turbulence modelling. The questionable negative behaviour of the skin-friction predicted by the viscous method will be easily improved in the future, and have negligible influence on the pressure. The Ftg. 27 gives an application of the present numerical method to a laminar separation over an adiabatic compression ramp, obtained with a viscous integral closure issued from similar solutions, see About [46]. (b) Symmetrical supercritical airfoil The inviscid flow is now calculated by solving the transonic small disturbance equation for the potential. The method has been implemented into an inviscid relaxation code, involving the 2.6

•

° •` ~ ,e•e / o ee

N.=*:

2.~ R =s.8.Io" rw.r.w T,,=I3I*R

I

..'"

,

,,~J ~IP

" .-

A~,4aATIC WAll

" '-'° 7 L 1.4"

"Reottachmont

3

I o.~

o'e

L

i2

7:6

~

x ~.,

Hg. 27. Wall pressure evolution for supersonic laminar ran=,p-induced separation. O, Experiment (LewisKulx)ta-Laes); ...... , Navier-.Stokes solution (Carter); .... , Interacting method of Wede-Vatsa, - ......, Present interacting method.

Numerical studies in high Reynolds' number aerodynamics

2"/

Murman--Cole technique with conservative schemes. Without viscous effects, an approximately normal shock-wave is present near the wall. We hope from the present approach, not only to match the viscous layer upstream and downstream the shock-wave, but also to deal with the viscid-inviscid matching at the foot of the shock-wave. Then the lower part of the normal shock-wave has to curve towards upstream, and to end into compression waves that provide at the wall a continuous compression, scaled by the local thickness of the viscous layer. Generally however, the computational grids used for common inviscid calculations are considerably too coarse to resolve those phenomena, so that a choice for numerical cost has to be done. Either the inviscid grids are used, then the viscid-inviscid interactions are, for most of them, only weak interactions, and the computing cost is almost the same as for inviscid codes. However the viscous resolution is incomplete, and, for example, local control-volume approximations have to be added to limit the error. Either the mesh size has to be reduced locally for x as y, until to be a little smaller than the viscous layer thickness, then the strong viscid-inviscid matching may be involved everywhere, see Figs. 28 and 29, and the solution depends almost no more on the grid. Those computational grids are then intermediate ones between inviscid and Navier-Stokes grids. Consequently, they involve an increase of computation cost, even for an inviscid solution. It has to be noticed, however, that such a resolution in the x-direction is unduly not involved in many Navier-Stokes simulations of airfoils. Presently a local refinement of x-mesh size is moved during the computation with the sonic point at the wall. Then, a continuous compression at the wall is captured, as well as the expected shock-wave curvature, and as very important values for dS*/dx, even without separation. Downstream of the strong compression the Math number is rather near one, as it is known from experiments. However this result is presently not consequent to an arbitrary use of non-conservative schemes, but to the real viscous effects. Finally, Fig. 30 gives a comparison of

y/c

/l JlllllJIJ I

o.

/I I llilllllIrl, Jlrl r I I lllllI|

Jlllll

I

x/c 0.30

cf Pie

0.50

0.70

H- 6"/0

/" H I

GIc 0,004. 0.10-

10 / -

0.002. ~,05.

/

5

0.30

C+

+ 6

" - ~ ~ /

0.50

1,1~,,

~ A . - - ~

~

. ~,

~'~

F~. 28. Skock-inducedseparationover a symmetrical 18%circular arc airfoil at Mo= 0.?88. R = 4 x IP (transonic small disturbances approximation).

28

J.C. LE BALLEURet at.

Y/C

T

0.2!

Shock-wave

0.1

l,;!l[1~:lit! "

..-.-,.. . . . 6~c

', [ ]

0.5

0.6

~J'~--

0.7

0.8

Separation

,~6

/~

1

Fig. 29. Transonic shock-wave boundary layer interaction region.

,Kp - 1 M o = 0.7425z/4e~Y/5

1

-1

'

M o = 0.788 I

/

Inviscid " "

-2~--- t-- I ' ~ - - - -

o

0:4

ols

M = 1

£2

'-

o Experiment : McDevitt • Deiwert .

+ + ..

~"~7.~

Present computations Navier - S t o k u (according to turbulent

modelling) : Oeiwo't Fig. 30. Pressure over an 18% circular arc symmetrical airfoil.

wall pressure both with experiments[47], and with the envelop of Deiwert's Navier-Stokes simulations for different turbulence modelling[48]. The airfoil is an arc-circular one of 18% thickness, and the viscous layer is supposed to be turbulent. Although the present inviscid small disturbance approximation will be questionable because of large local streamline curvature at the shock-wave viscous interactions, agreement with experiment as with the global approach is noticed for small trailing edge separation. For large shock-induced separation however, only a qualitative prediction is obtained, as with the global Navier-Stokes simulations. Improvements are hoped from a full potential inviscid solution. Acknowledgements--The authors

wish to thank H. Hollanders and W. Ghazzi for their important contributions to the numerical applications relating to the global approach, and J. J. Chattot for providing the T.S.P. inviscid code for the coupling approach.

Numerical studies in high Reynolds' number aerodynamics

29

REFERENCES 1. B. S. Baldwin. R. W. MacCormack and G. S. Deiwert. Numerical techniques for the solution of the compressible Navier-Stokes equations and implementation of turbulence models, AGARD-LS-73 on Computational methods for Inviscid and Viscous Two-and Three-Dimensional Flow Fields, Von Karman Institute, Rhode-Saint-Genrse, Belgium, pp. 2.1-2.24 0975). 2. R. Peyret and H. Viviand, Computation of viscous compressible flows based on the Navier-Stokes equations, AGARDograph AG-212 (1975). 3. H. Viviand, Traitement des Probl/~mes d'interaction fluide parfait-fluide visqueux en 6coulement bidimensionnel compressible h partir des 6quations de Navier-Stokes, AGARD-LS-94on Three-Dimensional and Unsteady Separation at High Reynolds' Numbers, Von Karman Institute, Rhode-Saint-Genrse, Belgium, pp. 3.1-3.21 (1978). 4. R. W. MacCormack, Status and future prospects of using numerical methods to study complex flows at high Reynolds' numbers, AGARD-LS-94 on Three-Dimensional and Unsteady Separation at High Reynolds' Numbers, Von Karman Institute, Rhode-Saint-Gen/~se,Belgium, pp. t3.1-13.14 (1978). 5. S. I. Cheng, A critical review of numerical solution of Navier-Stokes equations, Lecture Notes in Physics, 41, Progress in Numerical Fluid Dynamics, (Ed. by H. J. Wirz, pp. 78--225), Springer-Verlag (1975). 6. Yu. A. Berezin, V. M. Kovenja and N. N. Yanenko, Implicit numerical method for blunt-body problem in supersonic flows, Computers and Fluids, 3, No. 2/3, pp. 271-281 (June 1975). 7. R. W. MacCormack and B. S. Baldwin. A numerical method for solving the Navier-Stokes equations with application to shock-boundary layer interactions, AIAA Paper 75-1 (1975). 8. G. S. Deiwert. Numerical simulation of high Reynolds' number transonic flows, AIAA J., 13, No 10, pp. 1354-1359 (October 1975). 9. C. M. Hung and R. W. MacCormack, Numerical solutions of supersonic and hypersonic laminar compression comer flows, AIAA J., 14 No 4, pp. 475-481 (April 1976). 10. R. W. MacCormack, An efficient numerical method for solving the time-dependent compressible Navier-Stokes equations at high Reynolds' number, NASA-TM X-73 129 (July 1976). I I. J. S. Shang, Implicit--explicitmethod for solving the Navier-Stokes equations, AIAA J., 16, No 5, pp. 496.-502 (1978). 12. W. R. Briley and H. McDonald, Solution of the multidimensional compressible Navier-Stokes equations by a generalized implicit method, J. of Comput. Phys., 24, No 4, pp. 372-397 (1977). 13. R. M. Beam and R. F. Warming, An implicit factored scheme for the compressible Navier-Stokes equations, AIAA J., 16. No 4, pp. 393--402 (1978). 14. J. L. Steger, Implicit finite-difference simulation of flow about arbitrary two-dimensional geometries, A/AA J., 16, No 7, pp. 679--686 (1978). 15. H. Viviand and J. P. Veuillot, Mrthodes pseudo-instationnaires pour le calcul d'rcoulements transsoniques, Publication ONERA No 1978--4(1978). 16. J. P. Veuillot and H. Viviand, A pseudo-unsteady method for the computation of transonic potential flows, AIAA Paper 78-1150 0978) and ONERA TP No 1978-47. 17. R. W. MacCormack, The effect of viscosity in hypervelocity impact cratering, AIAA Paper No 69--354 (1969). 18. P. A. Gnoffo, Solutions of the Navier-Stokes equations for supersonic flow over blunt bodies in a generalized orthogonal coordinate system, M. S. Thesis, Polytechnic Institute of Brooklyn, (1974). 19. D. O. Gough, E. A. Spiegel and J. Toomre, Highly stretched meshes as functionals of solutions, Lecture Notes in Physics, 35, pp. 191-1%, Springer-Verlag (1975). 20. N. N. Yanenko, V. D. Lisseikin, V. M. Kovenia, The method of the solution of gas dynamical problems in moving meshes, 36me Colloque International sur les Mrthodes de Calcul Scientifique et Technique, Versailles (1977). 21. H. Viviand and W. Ghazzi, Numerical solution of the compressible Navier-Stokes equations at high Reynolds' number with applications to the blunt body problem, Lecture Notes in Physics, 59, pp. 434--.439, Springer-Verlag (1976), also ONERA TP No 1977-17. 22. D. M. Kuehn, Laminar boundary-layer separation induced by flares on cylinders at zero angle of attack, NASA Technical Report R-146 (1962). 23. R. Hakkinen, I. Greber, L. Trilling and S. Abarbanel, The interaction of an oblique shock-wave with a laminar boundary layer, NASA Memo 2-18-59 W (1959). 24. P. Laval, Mrthode instationnaire de calcul de rrcoulement transsonique dans une tuyrre. Note Technique ONERA No 173 (1970). 25. J. C. Le Balleur, Calculs couplrs visqueux--nonvisqueux incluant drcollements et ondes de cboc en ~,oulement bidimensionnel. Agard Lecture Series No 94, on three-dimensional and unsteady separation at high Reynolds' numbers, V.K.I. (1978) (or TP ONERA 1978-5). 26. V. N. Vatsa and M. J. Werle, Quasi three-dimensional supersonic viscid-inviscid interactions including separation effects. Univers. Cincinnati, Ohio, Tech. Pep. AFFDL-TR-75-138(1975). 27. J. C. Le Balleur, Couplage visqueux-non visqueux: analyse du problrme incluant drcollements et ondes de choc. La Recherche A~rospatiale No 1977-6 (or TP ONERA 1977-162). 28. L. Crocco and L. Lees, A mixing theory for the interaction between dissipative flows and nearly isentropic streams. J. Aero. Sci., 19, No l0 (1952). 29. J. M. Klineberg and J. L. Steger, On laminar boundary layer separation, AIAA Paper No. 74-79, Washington (1974). 30. J. E. Carter, Inverse solutions for laminar boundary layer flows with separation and reattachment, NASA TR-R-447 (1975). 31. R. W. Garvine, Upstream influence in viscous interaction problems, The Physics of the fluids. II, pp. 1413-1423 (1968). 32. S. Weinbaum and R. W. Garvine, On the two-dimensional viscous counterpart of the one-dimensional sonic throat, J.F.M.. 39, p. 57--85 (1969). 33. King Mon Tu and G. Weinbaum, A non-asymptotic triple-deck model for supersonic boundary layer interaction, AIAA J.. 14, No 6, p. 767-775 (1976). 34. M. Holt, Numerical methods in fluid dynamics, Springer-Verlag, Series in Computational Physics (1977). 35. J. M. Klineberg, Theory of laminar viscous-inviscid interactions in supersonic flow, Cal. Inst. Tech., Ph.D. Thesis (1%8). 36. G. L. Mellor, The large Reynolds number asymptotic theory of turbulent boundary layers, Int. J. Engineering Sci., 10, p. 851-873 (1972).

30

J.C. L~ BJ~.LEU~eta/.

37. R. E. Melnik and R. Chow, Asymptotic theory of two-d/mensionaltra/ling-edse flows, NASA SP-347, Part 1. Langiey (1975). 38. L. Crocco, Laminar separation, AGARD CP-I~ (Supplement) on flow separation, G~ttingen (1975) (or TP ONERA 75-113). 39. W. R. Briley and H. McDonald, Numerical prediction of incompressible separation bubbles, J.F.M., 69, pt 4, pp. 631-656 (1975). 40. M. J. Werle and V. N. Vatsa, New method for supersonic boundary layer separations, A/AA J., 12, No I I (1974). 41. J. C. L¢ Balleur, Couplage visqueux-non visqueux: M~thode num~'iquc et application aux ~oulements transsoniques et supersoniques, La Recherche 3,£'rospatialeNo 1978--2(or TP ONERA 197g-55). 42. G. W. Brune, P. E. Rubl~'t and T. C. Nark, A new approach to inviscid flow boundary layer matching, AIAA Paper No 74-601, Paio-Alto (1974). 43. M. J. Werle and S. D. Bertke, Application of an interacting boundary layer model to the supersonic turbulent separation problem, Rep. AFL 76-4-21, Univers. Cincinnati, Ohio (1976). 44. J. S. Shang and W. L. Hankey, Numerical solution for supersonic turbulent flow over a compression ramp, AIAA J., 13, No 10 (1975). 45. C. H. Law, Supe,tsonic turbulent boundary layer separation, A/AA J., IZ (1974). 46. G. About, Doctor Thesis to be published, on: Supersonic laminar separation, new improved intc~al methods, University of Lille (1979). 47. B. McDevitt, L. Lionel, J. R. Levy and G. S. Deiwert, Transonic flow about a thick circular airfoil, AIAA Paper 75-878 (1975). 48. G. S. Deiwert, Computation of separated transonic turbulent flows, AIAA J., 14, No 6 (1976).