Implicit numerical method for blunt-body problem in supersonic flows

Computersand Fluids,Voi. 3, pp. 271-281. PergamonPress, 1975. Printedin Great Britain

IMPLICIT NUMERICAL METHOD FOR BLUNT-BODY PROBLEM IN SUPERSONIC FLOWS Yu. A. BEREZIN,V. M. KOVENJAand N. N. YANENKO Computer Center, Novosibirsk, U.S.S.R.

(Received 25 September 1974) INTRODUCTION

Blunt-body problem in supersonic flows is a classical gas-dynamic problem. In the last decade many scientists have carried out calculations including viscosity and thermal conduction. Such phenomena are described by a very complex system of the partial differential equations, and numerical methods are necessary to employ. Many calculations of the two-dimensional problem were made by means the explicit schemes[I-6]. All the schemes are very easy in sense of realization and conditionally stable: h z~
h 2Rep}

lul+lvl+cV '

3"

The condition on time step z (or iterative parameter for steady problem) is a very serious one, particularly at moderate and low Reynolds numbers. In paper [7] a scheme has been suggested, which excluded the second condition connected with viscosity because the convective and viscous terms were taken on upper time level. However, the first condition remained (usual Courant condition in gasdynamics), because pressure gradients were taken on the lower time level. In papers [8, 9] the schemes have been treated, which excluded the second condition too, but the scheme[9] with Dufort-Frankel approximation to viscous terms was conditionally approximating one, and it is necessary to have z ~ h 2 in calculations. For solving of steady problems by relaxing the solution in time it is important to have absolutely stable scheme (without any limitations for iterative parameter). In this paper, we construct an implicit scheme with absolute stability and full approximation for two-dimensional (plane and axisymmetric) blunt-body problem in supersonic viscous and thermal conductive gas flows. 1. F O R M U L A T I O N

OF THE PROBLEM

To solve the blunt-body problem we will use curvilinear 'natural' coordinates connected with a body surface and depicted in Fig. 1. Let gas be an ideal one, with constant Prandtl number Pr and constant cp [co = y, viscosity/z

A

0

C

Fig. 1. 271

B

272

Yu. A. BEREZIN,V. M. KOVENJAand N. N. YANENK0

is equal to

g = g=(T/T=)%

(1)

Index oo denotes freestream functions 0.5 ~< oJ ~< 1.0, volume viscosity go is equal to go = - ~ g (Stokes' formula). Free stream is supposed to be uniform. System of equations for compressible viscous thermal-conductive flows can be written in the following non-dimensional form (v = 0 is a plane case, v = 1 is an axisymmetric case, external forces are absent):

t + at Os t Uon ~- al"a2

~ au 4 a au ] Ou+vOu+ Ou K v 2 + I Op_ 1 [1a2osa a2 ot a2 7s U~n a2 p O-'n- pRe Os ~--~-~-ffg~--~+F. Ov+ v Or+ Ov + Kuv + 1 Op_ 1 [ 4 O g Or+ 0 Ov -1 at a,O-s u-~ a2 pa20s pR-----e3a= Os a20s -~ng~n +F°j dT

v aT

aT

p

--+----+u--+"-';-~ at a2 0s On am a2p

am"v+o-~a,"a2u =p-~e

F.=a~-~gk~-ff---~]- ~

g(d2+ vd4)+ 2g K

1 0 g OT ~~__~nlZ_~n + dp Prr (12Os a:Os z,a3

-

,ou F: =

gdo-~-~: ~ g t ~ n

(Ou]2+vd2 ~b = 2g[d22+ \On/

a2

L \

a2 /

1.2 +~a, --~

' '

r r La2

Os

do= 1 0 u Kv d, = do + OV d3 = - - 1 -Oa, ala2 0s ' a2 Os a2 ' On ' d2 = a : l ( ~ + Ku), d, = ua3 + vd3, a m = R + n .cos ~, a2 = l + K n ,

a3=cos¢/am.

(2)

Here R = R (s) is the distance from body axis to body surface, ¢ = ¢ (s)--angle between the axis of symmetry and a tangent to a body contour, u, v are n, s - - c o m p o n e n t s of gas velocity, p is the density, Re is Reynolds number, K is a longitudinal curvature of a body surface. Density and velocity are scaled by freestream values p~, U~, and temperature is scaled by U®2/cp. Characteristic length L is a curvature radius of a body in a front critical point. State equation has the following form P = Y-lpT. 3'

(3)

Solution of the problem (1-3) we will find over domain D, limited by a body surface (n = 0, 0<~ s ~< sl), axis of symmetry AB and upper boundary Rl(s), which was chosen so that the disturbances from a body were dissipated up to the boundary. On the body surface we have no slip condition u = v=0

(4)

or

(5)

and one of the thermal regimes

OT --=0 0n

T = T ( 0 , s).

On the axis of symmetry

Op OU OT v = O. Os Os Os

. . . . . .

(6)

273

Blunt-bodyproblemin supersonicflows On the part 0 ~< s -< so of the boundary Rt(s) an undisturbed flow is given:

1

p = 1, u = - sin q~, v = cos q~, T = (y _ 1)M®2

(7)

o_~: o

(8)

at s >So Oe

where f = (p, u, v, T), e is directed along the axis of symmetry, M~ = U®/X/(7 - 1)M® is Mach number in a free stream. If we want to solve the problem in a part of the domain 0 ~< s ~< So, 0 <~n <~R~(s) it is necessary to have an additional boundary condition, for example,

Of-o

or

Oe -

0~[-

~

-

0

(9)

(see [10]). Effect of the boundary conditions (8, 9) on a flow in the domain must be studied by methodical calculations. For the sake of simplicity we transform the domain D into a unit square by means the following way:

(lo)

= S/S*, ~l = n/Rl(S),

where s* = s~, when the solution is found in a whole domain, and s* = So in a case of part of domain. Excluding pressure p by means state equation (3), we get the equations (2) in the following form (11)

s:t;/ --t:'7 O~, O2, B - - t h e matrices of fourth order with non-zero elements on the main diagonal:

fl,=~zovO

qzo~ o g o}

t a--~O¢ pRea2 0¢ a2 ~ ' {( zv) 0 1 [ 0 0 qz 0 gz 0 ] } 1~2= z ~ u - ~ O~ oRe q~z~ -~ lz ~-~-~ a2 On a2 an

~ qzo [zO ix a

q =

1

alxz~]}

ql =

4/3 0

e~=y-1;

e

y-1

------,

y

_ . 0 ¢. _. 1 .

ZO .

Os

s*'

,

PU Zl =

Or/ =

-

1

On R,(s)

,

z = ~ s - "oz~R, (s), R, (s)= OR, as" 03, O,--the matrices accounting for pressure in the motion equations and the terms div ti in the

274

Yu. A. BEREZIN, V. M. KOVENJAand N. N. YANENKO

remaining equations:

a~a2

I i lh=

00 zopKv al~ ic91 a2

zoeT a Kv

0

zoe

'

~Pao"@ 0a"zoe~T ~ a:o~/ at~a2 al ~

l

0

pz.......~,!

zp

a

al~a2071a]~a2 a~a2Or/a~

ez,____T~ p Or/

0

0

ezT O 0 a2p Or/ 0 elz~T O 7 a , ~a2 al"a---~ r~

0

ar/~ 03 ez, _ ez O [ a2

ze,T O a~a2 Or/ a, ~

0

7

Boundary conditions (4-7), (9) can be written as before with change n ~ 77, ~ ~ s because on the body surface z ---0, and condition (8) has the form o f _ Of Os¢+ of Or/=o" Oe @ 0 e Or/ 0e

(12)

Let us suppose that the solution existe and unique. For solving the problem we will use a relaxing method with (7) as initial conditions. 2. FINITE-DIFFERENCE SCHEME In the domain {0~
di--l.e(fj

1,e+l--~--,,e--l)].

(13)

We define also matrix difference operators lqj", B approximating differential operators l~j, B with the second order of accuracy: convective terms in the operators IL h, fl2 h are approximated by formula

aAf= a

a-

a [(33~ - 4~_, + ~_=) + - ~ h a

a

(4~+, - 3fi - ~+2)

(14)

the second derivatives are approximated by formula

AdAf = ~hz[(dj + dj+,)(~+,-.~) - (dj + d~_,)(fj - fj_,)]

(15)

and mixed derivatives--by formula (13). Approximation of the first derivatives is made according to the following rule: if a > 0 (see (14)), the terms in the equations for density and temperature have

Blunt-bodyproblemin supersonicflows

275

A_f - 3J~ -4)~_, +3~-2 2h

(16)

A+f - 43~+'2 ~ -.~+2

(17)

the form

in the equations of motion--

and vice versa in a case a < 0. We introduce finite difference operators ~jh with the first order approximation according to -

a+lal

aAf=

2h

([j -)~-')+

2h

0~+1 -fJ) (18)

h

h

For solving the system (11) let us consider the following scheme with splitting according to the physical processes and coordinates:

flih{h"+ B h ]'hn - F ~.n ,

(E +rfi,")~"+"'=-r (E

+ " T f i 2 h ) f n+l/2 = +

(n+,]4,

(19)

=

+ ,a2)y

÷1

=

/:+1 = f : + Here ~ (a = n + j/4) is the vector of the same structure, which has ]h~. After excluding of the fractional steps the scheme (19) becomes a scheme of universal algorithm:

efh "+1- [h"_ _ ~ l-i h[h. + Ph"-- Bh[~ ~ T

j=I

(20)

where C = 14I (E + zl)jh). Our finite difference scheme (19) or (20) approximates the differential 1=1

equations (11) with the second order of accuracy in a case of a steady problem. In a case of non-steady problems the scheme (20) approximates the equations (11) with the first order of accuracy 00" + h). Therefore in a case of non-steady problems we propose the operator C = I41 (E + rlajh) instead of t~. The scheme

C ]h~÷l - Jth"- - ~ f~,hfh" + Fh" - Bhfh" T

(21)

j=l

approximates the equations (11) with 0(r + h2). Operator is named by inconsistent operator. Choosing the operator as a stabilizator, we can realize the scheme (20) by means three-point sweepings, in contrast with the scheme (21) which is realized by five-point sweepings. For solving the finite-difference equations (20) or (21) it is necessary to have the points inside a body contour (see (13-16)) or to use another formula. For example, in the points near the boundary the first derivatives can be approximated by employing the symmetric finite differences. Condition of monotonicity of the symmetric schemes h < 2/~/a [11] is stronger than for schemes (20), (21), that is 4Aj+I+ Aj I> Aj+2, Aj =j~÷1-3~ (see (14)).

276

Yu. A. BEREZIN, V. M. KOVENJAand N. N. YANENKO 3. L I N E A R A N A L Y S I S OF S T A B I L I T Y

We will make this analysis for model linear equations with constant coefficients.

#+

0~+ d,-~-~ dz

le

Zo~ Ov

oV+bs O.oV oV O-----~ - c 0~:On

+ (1~3+ l~4)f = b, 0~2

\

/

-dr (T Op+OT'~ +du l a~ \~ o¢ o¢1 zoe~7" av

P[zlor] a2 r t o. OT1 z [TOp

OT-I

(22)

-e~LF an+~] el Lz,g as

/iZo2 Kf Zo~ z~ d~=--, ds= z~a--- b~=qR--~pa2, d=--,a2 az as'

[©+ q,z, 2], c = 21~ee~ _q a2 2

bs= ~

The system (22) was received by linearizating equations (11) with v =/zo = 0. We will give a spectral analysis of (22) in the following limit cases: (1) inertial terms are equal to z e r o (dl = dz = 0) and ~ 3 f = ~ 4 f = 0; (2) viscous terms are equal to zero (b~ = b2 = c = 0). In the first case the finite-difference schemes (20), (21) for the equation (22) has the form (E - rb,A,,)(E - bsAss)fh"+l- fh"= (b,A,, + bsA=s- c A,s)L" T

(23)

where Aij is operator approximating the second derivatives by formula (13), (15). Characteristic equation is (1 + B0(1 + Bs)(A - I) + B~ + B2 - C = 0 where B~ =(4~'/hjS)bj equation are equal to

sinS(kjhj/2), j= 1, 2; C=(,c/h,hs)c sin klhl .sinksh2. Roots of the 1+ C +

BIB2

Aj = (1 + B~)(1 + Bs) and IA[ = max IA, I ~ 1. i

We consider now a stability of the difference schemes in the second case when bl = b2 = c =0. Let transformation (10) is orthogonal one, free terms connected with curvilinearity of coordinates are small, and dj/> 0. Then we have the following system of linear equations C'(A - 1) + ~ = 0

1 ~N: nN2 1 = to 1 0

~N1 0) c2NIN2 N: 1 I '

(24)

277

Blunt-body problem in supersonic flows

~=

nN2

too

n

0 exTN2

0 too e~TN~

N2 ~Oo

tOo = A, + A2, to = (1 + A,)(1 + A2), c = ~/(7 - 1)T, 2r

Ao,=~sin 2

k,h,

,

. r

.

Bo,=Aoj+t~slnkjhj,

n

7" ~,

A j = dj[Boj - Ao,(COS kjhj - i sin k, hj)], j = 1,2.

Nj = zffAoj(COS kjhj + i sin k j h j ) - (Ao, - i ( r / h ) sin kjhj)]. Characteristic equation for the system (24) has the form

[to(X - 1) + too]{[to(X - 1) + too]3 +

c'N~2N22(A -

1)[to0t - 1) + 1]2 - c2(N~ 2 + N22) x [to(A - 1) + 1]2[to(A - 1) + tOo]}= O.

In the case of hypersonic flows c2.~ lu [2 the roots of the characteristic equation are equal to Ai=l-~o

(Do

and

[ X , [ ~ < l , ] = l . . . . . 4.

In another limit case c2-> lul 2 we have AI = 1 _ to_.-? to '

1

A2=A3= 1-~o, A,=I, Ix, l~
Stability of the finite-difference scheme with inconsistent operator C we study on the following linear equation

Ou+ Ou a-'t" a ~

=0' a>0

for which the scheme can be written as n+l

(1 + ~.a~,) u

n

-u

T

=-aAu"

where Au - u ~ - u j _ l h '

Au = 3 u j --4Uj-l+Ul-2 2h

If u "+1 = Au", t h e n

[A [2 = 1 + 2Co cos k h + Co2 sin 2 kh (1 + Co)2 (1 + Co)2+ c, 2 ~<(1 + co)2+ ci 2~<1' 2"r • 2 k h "r C o = - ~ a sin ~-, c~ = ~ - a sinkh

and the finite difference scheme is stable. N o w let us study an effect of curvilinear coordinates and nonorthogonality of transformation (10) on stability of schemes (20), (21). If in the equation (22) we conserve only terms, connected with mentioned reasons, the scheme (20) can be written as (b, = b2 = c = 0)

Yu. A. BEREZIN,V. M. KOVENJAand N. N. YANENKO

278

C fh"+'-/h°_o/:

(25)

'T

D =

l:

oa2

zen A2+ do

0

0

C1 =

0

zetl.a2+ [ r 2aao--~z r zen A2+

\-o

°: 7

etZTA2 a2

1

|

ze A:

"CPZA2

0

1

a2 -'cdo

"c:dozeA

rdo

1

r l Aa2 2+

0

ze

~2+

l

a2

do =KfJ,- n : = .7~ a2

p

Finite-difference operators A2: are defined by (16), and (17). For the scheme (25) we have the following characteristic equation (A - 1)2[(A - 1)2+ b2A 2] = 0 where

b = r4do2 + 4Z22c2(doZ + l)sinz K2h2 a2 n: 2 The roots of the equation are equal to 1

A,=A2=I,

A3.4=l___ib

and

[xl=maxlajl~
Hence, in the limit cases considered the finite-difference schemes (20), (21) for the system of linear equations (21) are unconditionally stable. 4. R E A L I Z A T I O N OF THE BOUNDARY CONDITIONS

For realization of the boundary conditions we construct a finite difference scheme similar to (20). For example, the scheme for density on a body surface n = 0 can be written as

pn+l/2 = ~ n+l/4

n

~pn+314

=

~pfl+l]2

~TZop A ~ n+3/4 Xtl~v

a2

~o~+' = ~o" ÷ ~ -

n

~z,o A~.

ph"÷l = ph~ + ~'o"÷~.

n+l

,

Blunt-bodyproblemin supersonicflows

279

By excluding fractional steps we receive the scheme

.+,.

ph r- p" + ~'ph"

u.+,-.)

A1

r- v," ~_z~A= h -r" Uh. = --ph"

AlVh" + ziAzuh"

)

or equivalent scheme with a full approximation

z,a,,,.o+,) = o, 1"

\a2

The scheme (20) for the boundary condition (12) can be written in the following form

~"+'"+raAl~ "+''=

(aA,+bA=)]h",

~"+"=

~.+,,4_

,__~..+,,=

..

~--+,,, = ~-+, = ~°+,,~, f;+, = f; + ~-°+,. It is easy to see that the schemes (26), (27) are absolutely stable. 5. N U M E R I C A L

CALCULATIONS

With the scheme (20) proposed we have calculations of steady supersonic viscous flows around cylinder and sphere. Two cases of domain boundaries were compared. In the first case (Fig. 1) and domain D was limited by a body surface, axis of symmetry and upper boundary Rl(s), where

R~(s)=a+b b-2 2

a cos (s ~).

(28)

Here a is a distance from a front critical point at s = 0 to the external boundary RI(0), b is a distance from a body at s = ~r to Rl(Tr). The curve R~(s) was chosen so that the coordinates in a vicinity of symmetry axis were orthogonal coordinates. In this case the flow near the symmetry axis was almost one-dimensional and depended on n only. With changing the parameters a, b, a we can stretch the domain to direction of a wake. In our calculations a = 3/4. The parameters a and b depend on Mach and Reynolds numbers. For example, at M = 6, 10 ~ 5 influence of back boundary of domain on the flow was only ~ 0.1 per cent near a body. In the second case we made calculations in a part of domain D (0 -< s ~< so). On the line s = so,

IC-

3-0~

"\ Ff" ~\

--U

--

//

0.5

F/

1", I 05

Fig. 2.

I.O

280

Yu. A. BEREZIN, V. M. KOVENJA and N. N. YANENKO

0 -< r~ -< 1 the condition (Of/Os)ls=,o= 0 was used. Comparison of solutions in a whole domain D and in a part domain D showed that the condition (Of/Os) = 0 influenced on a flow only near the back boundary. In the calculations the following parameters were used: I = 1 6 , K = 2 0 , w = 0 . 5 (or 0.75), P r = 0 . 7 2 , 3/=1.4. To check an influence of the scheme order on a solution we have made the numerical calculations of a flow around sphere (M~ = 2, Re= = 90, ~o = 0.75) by schemes of the first (see (18)) and the second order of accuracy (see (14, 16, 17)). In Fig. 2. the profiles of velocity and density on symmetry axis in front of the body are represented. Solid lines are for the second order scheme, dotted lines are for the first order scheme, circles are the results of paper [5]. It can be seen that using the first order scheme we get 'spreaded' shock wave due to numerical viscosity. In Figs. 3, 4 distribution of velocity and density on the symmetry axis in front of cylinder and in a section s = 0 and s = 7r[2 are represented (M = 2, Re = 45, oJ = 0.75). The condition of thermal isolation on the body is taken. The calculations were made by the scheme of the second order of accuracy and in a whole domain D. Solid lines are for our results and dotted lines are for results of paper[5]. The shock wave is spreaded on 5-6 grid points.

3.C

--

I0

/ <

-u

i\'/ #

e<

I0--

u 05

I 0

I-0

20

Fig. 3.

O~

_u

,

IO n

Fig. 4.

I

I

20

30

Blunt-bodyproblemin supersonicflows

281

We have also made calculations of hypersonis viscous flow around a cooled sphere at M = 40, Re = 200, ~o = 0.5, To = 16.05. In Fig. 5 velocity distribution in a section s = 0.576 (33 °) is represented. Our results were compared with results of paper [10] obtained by method of 'large particles' (dotted lines). In our calculations we used a part of domain D (So = 0.5777r) and upper boundary was taken by formula (28) with a = 0.5, b = 2. In hypersonic flows there is merging of a shock zone and boundary layer, and all the profiles are smooth (Fig. 5). There are strong increasing of density near the body, maximum of temperature is near the shock front in an ideal gas and almost equal to value in ideal gas.

IOI

C'

-u

O2

I

]

O4

n .?.5. Steady problem was solved by relaxing in time up to max I(l/o)(oo/ot)l ~ 10-3. In all the cases considered it was necessary 150-200 iterations for steady solution. Initial conditions were freestream conditions (7), and on a body surface-conditions (4), (5). Iterative parameter r was at first equal to r = h ~ , ] 2 (h~n = rain ( h , h2)) and after 50-70 steps was increased up to ~"~ 3hm~; then the calculation was made with a constant time step ~'. The change of iterative parameter gave us a possiblity to damp the strong disturbances near the body which could increase due to inconsistence of initial and boundary conditions at initial time interval. Control of accuracy was made by change of grid steps (that is by change of a boundary R , ( S ) location), by conservation laws, by comparison of parameters in the critical point and with another calculations. If we changed steps of finite difference grid twofold in the direction n, the solution on body surface and the axis of symmetry was changed in 3-4 per cent. Error in conservation laws was maximum for the density (~1.5-2 per cent) and it can be explained by small number of grid points over the domain ( I = 16, K = 20). REFERENCES 1. RichtmyerR. D. and Morton K. W., Difference Methods for Initial-Value Problems, 2nd Ed., Interscience,New York (1967). 2. McCormac R. W., The effect of viscosity in hypervelocityschemes, AIAA Paper, 69 (1969). 3. BrailovsnayaI. Yu., Differencescheme for numericalsolutionof the two-dimensionalunsteadyNavier-Stokesequations for compressiblegas, Dokl. Akad. Nauk SSSR, 160, 1042 (1965). 4. Tolstyx A. I., Method of numericalsolutionof the Navier-Stokes equations for compressiblegas in a wide range of the Reynolds numbers, Dokl. Akad. Nauk SSSR 210, 48 (1973). 5. PavlovB. M., Numericalstudy on blunt-bodyproblemin supersonicviscousflows,in Nekotoryeprimeneniya metoda setok v gazovoi dinamike, IY, 181, MSU, Moscow (1971). 6. HawlowF. H. and Hirt C. W., Recentextensionsto Eulerianmethodsfor numericalfluiddynamics,Zh. vychisl. Mat. mat. Fiz. 12, 656 (1972). 7. Podezhaev V. I., Numerical study on natural convection in fluidsand gases, in Nekotorye primeneniya metoda setok v gazovoi dinamike, IY, 86, MSU, Moscow (1971). 8. Scala M. and Gordon P., Solutionof time-dependentNavier-Stokes equations for the flow around a circular cylinder, AIAA 3" 6, 815 (1968). 9. VictoriaK. J. and WidhopfG. F., Numericalsolutionof the unsteady Navier-Stokesequationsin curvilinearcoordinates, Lecture Notes in Physics 19, 254, Springer-Verlag(1973). 10. BelotserkovskiiO. M. and SeverinovL. I., Conservative"flow"methodand numericalsolutionof blunt-bodyproblemin a viscous heat-conductivegas, Zh. vychisl. Mat. nat. Fiz. 13, 385 (1973).