A numerical method for solving the equations of a viscous shock layer

64 (20)-(23) have the form of partial differential equations with a small Eqs.(8)-(171, parameter accompanying the heighest derivative. The use of these equations enables one to construct efficient numerical schemes, since there are terms in the equations which are analogous to artificial viscosity, and to obtain a better description in the case of the rapidly occurring processes whih are accompanied by a rapid change in the characteristic physical parameters. We note that equations of similar structure have been postulated in /4/ with the aim of proving the unique solvability as a whole of the initial-boundary value problems which are described by the Navier-Stokes' equations in the case of an incompressible However, the equations obtained here differ both from those postulated in /4/ and liquid. from those used in /5/ and /6/.

Remark added in proof. It is known that the problem of the "sharpening" of the velocity distribution function arises in the numerical integration of the Vlasov equation. As a result it is necessary to resort to a filtration procedure which leads to the appearance of cross second-order derivatives inthevlasov equation (see /7/, for example). It is readily seen that the above-mentioned equations are a special case of Eqs.(6). REFERENCES 1.

2. 3. 4. 5. 6. 7.

HIRSCHFELDER J., CURTISS CH. and BIRD R., Molecular Theory of Liquids and Gases, /Russian translation/, IIL, Moscow, 1961. ALEKSEYEV B.V., Mathematical Kinetics of Reacting Gases, Nauka, Moscow, 1982. BRAGINSKII S-I., Transport phenomena in a plasma, in: Problems of Plasma Theory, 1, Gosatomizd., Moscow, 183-272, 1963. LADYZHENSKAYA O.A., Mathematical Problems of the Dynamics of a Viscous Incompressible liquid, Nauka, Moscow, 1970. ELIZAROVA T.G. and CHETVERUSHKIN B.N., On a computational algorithm for calculating gas dynamic flows, Dokl. Akad. Nauk SSSR, 279, 80-83, 1984. ELIZAROVA Q.G. and CHETVERUSHKIN B.N., Kinetic algorithms for calculating gas dynamic flows, Zh. vychisl. Mat. mat. Fiz., 25, 10, 1526-1533, 1985. KLIMAS A.J., A method for overcoming-the velocity space filamentation problem in collisionless plasma model solutions, J. Comput. Phys., 68, 202-226, 1987.

Translated

Vo1.27,No.3,pp.64-71,1987 U.S.S.R. Comput.Maths.Math.Phys., Printed

in Great Britain

by E.L.S.

OO41-5553/87 $10.00+0.00 01988 Pergamon Press plc

A NUMERICAL METHOD FOR SOLVING THE EQUATIONS OF A VISCOUS SHOCK LAYER* S.A. VASIL'YEVSKII,

G.A. TIRSKII

and S.V. UTYUZHNIKOV

A method is proposed for calculating supersonic flows around blunt bodies. This is based on the numerical solution of the complete two-dimensional Allowance for the upward transfer equations of a viscous shock layer. of disturbances through the stream into the subsonic regions is ensured by carrying out global iterations along thewhole of the segment of the Each global iteration is calculated by the shock layer considered. Questions regarding the stability of iterative process marching method. are considered as well as the correctness of the formulation of the The results of the calculations mixed problem on each global iteration. are in good agreement with experimental data. 1. Introduction. In this paper we consider the stationary supersonic flow around blunt bodies by a viscous In this case heat-conducting gas under conditions when the continuous medium model is valid. However, the the flow of the gas is described by a system of Navier-Stokes' equations. numerical solution of these equations for the above-mentioned class of problems, especially at fairly high Reynold's numbers, is an exceedingly difficult problem which makes considerable We note that the demands on computer time and core memory and is not always necessary. investigations of flow problems using the Navier-Stokes' equations with the help of the simplest explicit scheme of the method of establishment required many thousands and even tens of thousands of iterations (time steps) /l, 2/ while calculations using the more efficient implicit schemes in /3/ or /4/ required of the order of 300-800 iterations for just the windward surface of a sphere. In order to carry out practical calculations it is therefore necessary to use simplified equations which are less then universal but sufficiently accurate (in an asymptotic sense), and suitable for solving the important practical class of problems *Zh.vychisl.Mat.mat.Fiz.,27,5,741-750,1987

65 on the supersonic flow around bodies by simple and efficient numerical methods. One such simplified system of equations is formed by the complete equations for a viscous shock layer (v.s.1.) (see /5-g/, for example) which describe the supersonic flow around the windward surface of smooth blunt bodies at moderately small and large Reynold's numbers (Re,biOo) and contain all the terms of the Euler equations and of the Prandtl equations. An economical computational method is proposed for solving the v.s.1. equations in the present paper. This method is based on the idea of carrying out global iterations for the whole of the computed region which was proposed for the first time in /6/ for calculations in a thin shock layer during flow around a hyperboloid. This idea has been developed further It has been found thattheimplementation of this method in /0, lo-15/ and other papers. enables one to reduce significantly the demands made on computer time (by a factor of approximately 10 compared with implicit methods of establishment). The method is also highly economical in the demands it makes on computer memory since, apart from the position of the shock way, it only needs to remember the pressure field for the next global iteration. 2.

Formul,ation of the problem. 1. Let us consider the flow of an ideal gas around a smooth axially symmetric (v=i) or planar (v=O) blunt body around which the gas flows under zero angle of attack. The v.s.1. equations and the corresponding boundary conditions for this case can be obtained on the basis of /6, 8/ and /16/. We shall write the v.s.1. equations in a curvilinear orthogonal coordinate system: X is the length of an arc of the contour of the body, measured from the critical point and y is the distance along the normal from the surface of the body. This choice of coordinate system ensures that the v.s.1. equations in the region of the preboundary layer pass asymptotically into the Prandtl equations as Re,-m. The v.s.1. equations can then be written in the form

(2.1) (2.2)

-_

1

Fle,HtZrv

a dy

a”_” )f( (ff,?+ (ay H,R (2.3)

Fig.1

(2.4)

Here u and v are the physical components of the velocity along the 5 and y directions, H is the radius of curvatureof is the total specific enthalpy, o is the Prandtl number, R(x) the body, H, is a metric Lam& coefficient, r=r,+ycosa (see Fig.1) and Re,=p,U,R(O)/p.... In Eqs.(2.1)-(2.4) the density p is referred to the density of the approach stream p_, the velocity projections are referred to the velocity cf the approach stream u,, the pressure P is referred to p-u,*, the viscosity p to the value of the viscosity in the approach stream p_,H to the value of the total enthalpy in the approach stream and X, y, and also r, are referred to R(0). The generalized Rankine-Hugoniot conditions /16/ are the boundary conditions on the internal boundary of the shock wave for the v.s.1. equations: sin!3 u.=u,tgp.--k.COSP.'

H,=i..-

‘* Re,o sinp

P,=(l-k.) sin'~+(TM_')-'. is the Mach number of the approach stream and 7 is Here k,=p,-', fi.=B-a, M,-[p,U,'I(~P~)I"' the ratio of the heat capacities. Slipping and temperature jump conditions /6/ are used as the boundary conditions on the It is also assumed that the surface of the body is impermeable or an injection of body. The wall temperature, referred to specified intensity occurs from the surface of the body. the impact temperature To, was assumed to be equal to a certain constant value T=. The following relationships /16/ were used in order to determine the detachment and inclination of the shock wave:

USSR 27:3-E

66

“0

r*‘+‘/(v+l)=

spur’ dy+Q,,

(2.5)

0

(2.6)

Here the injection

parameter

Q.--j

as.

p,v,r,’ 0

Relationship (2.5) expresses the mass balance while relationship (2.6) is a geometrical relationship. 2. The v.s.1. equations, although they are substantially simpler than the Navier-Stokes' equations, are complex for numerical solution since they retain the elliptic nature of the In order to demonstrate this, we follow the increase problem in the subsonic flow regions. Let in the perturbations during the solution of the Cauchy problem for the v.s.1. equations. us represent the required functions in the form of the sum of the exact solution and a small short-wavelength perturbation W. We shall assume that the perturbation varies so rapidly At the throughout space that the terms with the highest derivatives are the decisive terms. same time the exact solution only changes slightly over the period of the harmonics of the perturbation and it may be set equal to a constant in the neighbourhood of each point being considered /17/. As a result, by neglecting the effect of small source terms, the following local system of equations, which we shall consider as a model system, can be obtained from the v.s.1. equations: Aaw+Bdw= ax

(2.7)

Pw

Ddyz.

ay

Here

6P WE

,“I

A=&

(

B=

u 0 l/P OurI’ 0 0

-a%

v

”

:,

;“”

- a%

P

n 1

8P

0

P

0

I‘

I/p

0

I’

’

0 ”

11

()

D=+

u 0

n

II

l/P 0

0 0 on

0

0

o

m -

at/of)

0

0 *

Y/UP.

where

a=(yPlp)“* is the local velocity of sound referred to U,. Let us now map the flow domain between the body and the shock wave on a half-plane. (S,.n),s=Z, n=!//y,. Then, from In order to do this, we transformtotheindependentvariables (2.71, we have

dW

(2.8)

A-+Cz=D,s,

as

where C=gA+$B,

D, = -l_D. Y**

For the Cauchy problem to be well-posed with respect to s, it is necessary in the case should be non-negative of system (2.8) that the real part of the spectrum of the matrix A-ID, In This is not satisfied close to the wall where subsonic local zones always exist. /17/. It is these subsonic regions one of the eigenvalues of the matrix ;l-'Dt is negative /15/. known that the Cauchy problem is also ill-posed in a subsonic non-viscous domain since the This implies computational eigenvalues of the matrix A-‘C have a non-zero imaginary part. instability when the v.s.1. equations are integrated numerically using the marching method. The determination of the inclination of the shock wave, as was shown in /ll/, also requires that account is taken of the upward propagation of perturbations through the flow. We note that the different regularization methods (see /lS, 19/, for example), which are directed towards the construction of marching computational algorithms, are only orientated In these and in many other papers, towards supersonic flows with narrow subsonic regions. other methods and, as a rule, the method of establishment, have been used to calculate the flow in the region of bluntness.

3. The method of global iterations. The current global The iterative process is constructed in the following manner. awas in Eq.(2.2) and at a fixed angle iteration is calculated with a specified derivative preceding of inclination of the shock wave which are obtained from the calculation of the The pressure field P'"'-" and the angle of inclination of the shock wave global iteration.

67 fj(‘“fl)

for the subsequent J relationships:

(m+l)-th

w+u

8.

global

iteration

are determined

=-7arctg( &[ %I,}

+(I-?)

from

the

P!“‘.

relaxation

(3.2)

is the pressure obtained as a result of the calculation of the m-th global Here Pcm) iteration, [dy,lti],, is the difference approximation of the derivative and t and T* are the relaxation parameters. An approximate analytical formula /16/ was employed in calculating the initial global iteration for an angle p. and, in doing this, the pressure was taken in the Newton approximation. 1. In order to prove that the global iteration method is correct it is necessary to show that the Cauchy problem and the boundary value problem are correct on each global iteration In the general non-linear case it is practiand, also, that the global iterations converge. callyimpossibleto carry out such an investigation. It may, however, be possible to carry out a corresponding analysis under certain simplifying assumptions in the case of the local system of Eqs.(2.8), for example, which arises in the treatment of small short-wavelength perturbations. Although such a treatment does not finally answer the question as to whether the method is correct or not, it is of considerable interest in establishing the characteristic features of problem /20/. Since, on each global iteration, i)Pldx is the "source" term in the momentum equation, the matrix A may be replaced by the matrix At: u

P u 0 0

0 o

A=& 1

-u%_h

” 0 u 0

0 0 0 u

when investigating whether the Cauchy problem is well-posed with respect iteration in the case of system (2.8). The eigenvalues of the matrix A,-‘C are real and have the form

to x

on each global

It follows from this that the system of equations for a non-viscous gas on each global iteration is of hyperbolic type and the Cauchy problem with respect to s is well-posed for it over the whole of the computed region. In order that the mixed problem (the initial problem with respect to s and the boundary value problem with respect to n) should be well-posed, it is necessary that the boundary conditions should be solvable with respect to the required functions on the boundary jointly with the characteristic relationships for the characteristics approaching this boundary (as s increases) withrespect to the required functions /20/. It can be shown that this requirement is satisfied if the following inequality holds: (3.3)

cos/L~]v,]la,

Condition where V, is the component of the velocity vector on the normal to the shock wave. (3.3) is always satisfied in the case of the windward surface of smooth bodies when there is hypersonic flow around them. In a number of cases such as, for example, the case of bodies with a sharp change in the inclination of the generatrix, relationship (3.3) cannot be satisfied. This corresponds to the fact that the family of characteristics in this region which was earlier approaching the boundary n=l becomes an emergent family. The mixed problem for the equationsofan ideal gas then becomes ill-posed. The terms with the highest derivatives play a decisive role in the treatment of the Cauchy problem for the complete system of Eqs.(2.8). The eigenvalues of the matrix A,-‘D, are:

h*,z=O,

PHl

&3=._

Re,puy.2

’

h, =

WHl Re,opuy.’

In the absence of return flows all the eigenvalues are real and non-negative. This means that system (2.8) is not completely parabolic on each global iteration. When this is so, the Cauchy problem is well-posed /15/. 2. Let us investigate the convergence of the global iterations when there is a fixed external boundary in the form of a rectangle R=[O, sJX[O, 11. The iterative process may be written in the form, whichis canonical for a two-layer iterative scheme, proposed by Samarskii: w(n+l)_w(n)

+ LW’“‘=O.

.L

(3.4)

7P

Here L is

the operator

of the v.s.1.

equations

and L,.is the operator

of the system of equations

68 on each global iteration. Inversion of the operator L, is significantly simpler than inversion of the operator L since it is associated with the solution of a well-posed mixed (initial-boundary value) problem. We approximate the differential operatorsLand L, by means of difference equations and write out the solution of the corresponding difference equations when there are homogeneous boundary conditions in the form of finite series (we assume that such a solution exists): w;~‘=So-8

2:

r,

X'"'exp[-i(k,ph,S-k,lh,)l. llEI

Here k is a two-dimensional wave vector which, on summation , traverses the lattice k,=2nz,. where z1 and z: are integers so that 0<.zlh,
k,=4nz,isu,

The convergence of the iterative process is intimately associated with the stability of the difference scheme in which tF plays the role of a step size in time. Stability of the difference scheme means that the matrices G(T,, k)” must be uniformly bounded with respect to k for all 04~,nGT (see /21/l. If all the gi in (3.5) are uniformly bounded with respect to k, it follows from the Bachenan criterion that, for stability, it is necessary and sufficient that Confirmation of condition (3.6) has been obtained in practical numerical calculations. In the case of a simplified (model) version of the difference approximation of the derivatives, the uniform boundedness of gi for any fixed step sizes h, hasbeen demonstrated in /15/. In order to investigate the convergence of the iterative process (3.2) we will represent ,(nl) c, 9 as was done previously, in the form of the exact solution and a small rapidly varying perturbation o. In the blunt region, o satisfies the relationship /15/

Here m-y,(M'-I). In the case of a hypersonic shock layer, m-p,/p,-(r-l)l(yfl) on the bluntness. It follows that account should be taken, in constructing the iterative process, of the change in sign of m in passing across the sonic line in the neighbourhood of the shock wave, the position of which is previous known. For instance, it is found from the CourantFriedrichs-Lewy condition that a homogeneous difference scheme of the "explicit angle" type for (3.2) is unstable (here T plays the role of a step size in time) and this has been confirmed by numerical calculations. In the numerical calculations a scheme of the second order of accuracy with central differences was used for [dy.ldxl h . When this is done the following constraint on r and h holds: Imrlh1<1. (3.7) The constraints (3.6) and (3.7) on 5 and 7, are substantially weaker than the constraints which arise from the requirement of stability in the method of establishment. The stability of the solution in the case of the local scheme provides one with grounds for believing that the solution will also be stable for the initial system of v.s.1. equations. This has been confirmed by numerous calculations.

4.

Numerical implementation of the method and examples of calculations.

1. Independent variables of the type ofDorodnitsyn employed in carrying out the numerical calculations:

+I

+

variables

(5, q)

(see /16/) were

-.ycosa r’,

Solution of the v.s.1. equations on the current global iteration is based on the decomposition of the initial system of equations into two subsystems and, in fact, into two parabolic Eqs.(2.2) and (2.4) and a hyperbolic subsystem containing the second-order Eqs.(2.1) and (2.3). Scalar pivotal condensation using Petukhov's method with a fourth order of approximation with respect to n. but with an approximation of the derivatives with respect to 5 different from that employed in /22/. was used in solving the boundary value problem on each step in E for the parabolic equations (2.2) and (2.4) for known P and v. The boundary value problem for the hyperbolic subsystem (2.1) and (2.3) was solved for P and v by a two-point matrix pivotal condensation /14, 15/. An analytical investigation of equations with "frozen" coefficients showed that condition (3.3) is necessary and sufficient for the stability of

69 the pivotal condensation method which was used /15/. The derivatives with respect to g, which occur in all of the v.s.1. equations, were approximated with the aid of an implicit difference scheme t2-lr /23/ which combines a second order of approximation with good stabilizing properties. The so-called "soft" boundary conditionsforthe pressure: d'P/dgl=Oand ?P/a~"=O were imposed on the outside section. Numerical calculations showed that the effect of the abovementioned boundary conditionsonthe solution is only propagated upwards through the stream for l-2 computational sections. A variable step size in n: An(f, n), which was chosen at each point depending on the change in the functions in its neighbourhood from the relationship

was used to calculate flows with large stream gradients (large Reynold's numbers). This enables one to carry out calculations over a wide range of Reynold's numbers (1004 Re,GlO") which encompasses flow regimes from a completely viscous shock layer to a non-viscous shock layer with a thin boundary layer around the body /24/. 2. Calculations of the flow of an ideal gas around smooth blunt bodies over a wide range of flow regimes: T,/T,GT,GI, M,B->l,10
---_ 0

z

I

z

J

Fig.2

ofthepressure

SO

90

60

w,degree Fig.3

Fig.4 The convergence

0

on the body

Fig.5 Pa)(z) and the detachment

of the shock wave

The Y!"'(2) for a loo-hyperboloid when M-=20, Re,=lOS and T,=0.1 are shown in Fig.2 for i=$ and and Yc' stationary solution is indicated by the solid lines, the value of P$

by the small crosses. 2 by the dashed lines and the results of the calculation for i=4 In Fig.3 the thermal flux distribution along the windward surface of a sphere, calculated when M-=6, T,-0.4,Re,=10*. 10' and iOn, is compared with the experimental data /25, 26/ and with the known solution of the boundary layer equations /27/. increase as Re, decreases, which is in good The valueslof ~m(o)=qur/qowhen 30"~;o
70

equations presented-in /3/. We note that our calculation for Re,=lO' is in good agreement with experiment /26/, while our calculation for Re,=108 is in good agreement with the solution presented in /27/. A calculation when Re,=lO' without allowing for slipping on the body is represented by the broken line in Fig.3. The values of Pm (the solid curves) and C,Rez (the dotted and dashed curves) along the surface of a loo-hyperboloid , calculated for M-=20, T,=O.1, Re,=lO*, 103,10' and 10' taking account of slipping on the body, are shown in Fig.4. The small crosses are data from nonviscous flow tables /20/. The results of our calculations for P, and y, when Re,=lO* are identical with the tabulated values from /20/ to within 1%. Fig.5 shows the emergence, at large Re, on the asymptotic form of the boundary layer, of the values StoRez (St=q,/(l--H,) is the Stanton number) and (C,/COSa)oRe? for the critical point of the sphere, calculated both allowing for slipping (the solid curves) and ignoring slipping (the broken curves), and the temperature jump when M,=6 and 2',=0.4. We note that previously presented calculations on the cirtical line of a sphere /8/ have shown that the resulting asymptotic form for large Reynold's numbers Re,
71

14. VASIL'YEVSKII S.A. and TIRSKII G.A., A numerical method for the solving the complete equations of a viscous shock layer, Preprint 3119, Inst. Mekhan. Moskovsk. GOS. Univ. (MGU) , Moscow,

1985.

15. VASIL'YEVSKIIS.A.,TIRSKII G.A. and UTYUZHNIKOV S.V., The method of global iterations for solving the complete equations of a viscous shock layer, Preprint 3138, Inst. Mekhan. Moskovsk. Gos, Univ. (MGU), Moscow, 1985. 16. TIRSKII G.A., On the theory of hypersonic flow over planar and axially symmetric bodies by a viscous chemically reacting flow of gas with injection, Nauchn. Trudy Inst. Mekhan. Moskovsk. Gas. Univ., Izd. Moskovsk- Gos. Univ. (MGU), Moscow, 39, 5-38, 1975. 17. KOVENYA V.M. and YANENKO N.N., The Method of Decomposition in Problems of Gas Dynamics (Meted reasshchepleniyav zadachakh gazovoi dinamiki), Nauka, Novosibirsk, 1981. 18. VORONKIN V.G., Calculation of the viscous shock layer on blunted cones, Izv. Akad. Nauk SSSR, Mekhan. zhidkosti i gaza, 6, 99-105, 1974. 19. KOVENYA V.M., CHERNYI S.G. and YANENKO N.N., Simplified equations for describing flows of a viscous gas, Dokl. Akad. Nauk SSSR, 245, 6, 1322-1324, 1979. 20. LYUBIMOV A.N. and RUSANOV V.V., Flows of a Gas around Blunt Bodies (Techeniyagasa okolo tupykh tell, 1, 2, Nauka, Moscow, 1970. 21. RICHTMEYER R. and MORTON K., Difference Methods for Solving Boundary Value Problems, /Russian translation/, Mir, Moscow, 1972. 22. PETUKHOV I.V., Numerical calculation of two-dimensional flows in a boundary layer, in: Numerical Methods for Solving Differential and Integral Equations and Quadrature Formulae (Chisl. metody resheniya differentsial'nykhi integral'nykhuravnenii i kvadraturnye formuly), Nauka, Moscow, 304-325, 1964. 23. CHUDOV L.A., A difference method for calculating flows in a boundary layer possessing the property of strong stabilization of high-frequency perturbations, in: Some Applications of the Method of Nets in Gas Dynamics (Nekotoryeprimeneniya metoda setok v gazovoi dinamike), Issue 1, Izd. MOskovsk. Gos. Univ. (NGU), Moscow, 196-210, 1971. 24. VASIL'YEVSKII S.A.andTIRSKII G.A., The effect of multicomponent diffusion and high approximation for the transport coefficients on the heat transfer during the hypersonic flow round a blunt body, in: Applied Pmroblemsinthe Aerodynamics of Aircraft (Prikl. voprosy aerodinamiki letatel'nykh apparatov), Naukova Dumka, Kiev, 100-103, 1984. 25. ZAVARZINA I.F., Experimental investigationsof local heat flows on a sphere and on the spherical blunting of an axially symmetric body, Izv. Akad. Nauk SSSR, Mekhan. zhidkosti i gasa, 4, 158-161, 1970. 26. HICKMAN R.S. and GIEDT W.H., Heat transfer to a hemisphere-cylinderat low Reynolds number, AIAA Journal, 1, 3, 665-672, 1963. 27. LEES L., Laminar heat exchange on blunt bodies at large supersonic velocities, in: Scientific Problems of Artificial Satellites /Russian translation/, Izd. Inostr. Lit., MOSCOW, 243-279, 1959. 28. GERSHBEIN E.A. and TIRSKII G.A., Flow of a viscous heat-conducting gas in a shock-layer in the neighbourhood of bluntness under intense injections, Nauchn. Trudy Inst. Mekhan. Moskovsk. Gos. Univ. (MGU), Izd. MGU, Moscow, 1, 1970.

Translated by E.L.S.