- Email: [email protected]
Numerical solution of some problems in gas dynamics
NUMERICAL SOLUTION OF SOME PROBLEMS IN GAS DYNAMICS* 0. M. BELOTSERKOVSKH, F. D. POPOV, A.I. TOLSTYKH V. N. FOMIN, A. S. KHOLODOV Moscow (Received
20 April
1969)
Notation s is the length of arc measured along the contour of the body; n is the normal to the body (the contour of the body corresponds to n = 0); t is the time; R is the radius of curvature of the body; r is the distance from the axis of symmetry of the body; t is the distance along the normal from the surface of the body to the shock wave, [=&c(s), & 1 -[; is the angle between the tangent to the contour of the body and the direction of the incident flow; u and u are the velocity components in the direction n and s respectively; w is the magnitude of the total velocity; p is the pressure; p is the density; h is the enthalpy; ‘I’ is the temperature; p is the molecular weight; p” is the coefficient of viscosity; Moois the Mach number of the incident flow; Pr is the Prandtl number; Re is Reynolds number, calculated for the parameters of the incident flow, the characteristic diameter of the body and the viscosity for h = w,‘; is the coefficient
*Zh.
u’jchisl.
Mat. mat. Fiz.,
10,2, 401-416,
1970.
158
of friction;
NumericuZ
solution
of some problems
in gas
dynamics
159
is the heat flow;
v is the frequency; hv’ is the absorption coefficient Bv is the Planck function.
taking induced emission into account;
Recently the study of problems in gas dynamics by means of numerical methods using high-speed computers has been intensively continued. Further improvement has been made by means of serial calculations by well-known methods, which have also led to the development of new numerical algorithms defined by the consideration of a number of new problems in gas dynamics of scientific and technical interest. The latter include problems connected, for example, with the consideration of the streamlining of blunt bodies by the flow of a gas at low supersonic speeds (1 < M, < 2); the determination of the characteristics of blunt bodies in the presence of discontinuities of the generator; the study of stream regimes in the presence of viscosity, physico-chemical transformations, radiation etc. Below we describe the results of some investigations of this type of problem using new numerical schemes of the method of integral relations, finite differences, combinations of them, and also build-up methods. The following four types of problems are considered: (1) the streamlining of an axially symmetric cylindrical the supersonic flow of a perfect gas; (2) the streamlining
of blunt bodies at low supersonic
butt and sphere by
speeds;
(3) the streamlining of blunt bodies by the hypersonic flow of a viscous perfect gas for low values of the Re numbers; (4) the streamlining of blunt bodies by a supersonic flow of an equilibrium dissociated and ionised gas taking selective radiation and absorption into account. Here the main attention is given to the results of the investigations. In the description of the numerical methods only the fundamental ideas of the construction of the corresponding schemes are presented and the necessary references are given. This paper forms the main content of a report by the authors at the Third Aii-Union conference on Theoretical and Applied mechanics (MOSCOW, January 1967).
160
0. M. Belotserkouskii
et al.
1. As is well known,. at the present time different numerical schemes for solving the non-stationary equations of gas dynamics have been developed which are also used successfully to obtain the stationary characteristics of the stream by the build-up method. Approaches of this type are now possible in practice due to the introduction of computers of comparatively large capacity. Moreover, it has to be borne in mind that in some cases the use of “classical” stationary methods, well developed earlier, encounters certain difficulties in the solution of a number of important problems of gas dynamics (for example, the investigation of transonic streams, obtaining the stream field for bodies with discontinuities etc). We consider one of these %on-stationary” approaches in the classical problem of the streamlining of an axially symmetric cylindrical butt by the steady supersonic flow of a gas. This problem has been solved by using the modified non-stationary Euler-Lagrange method (close to the discrete model of Harlow), which can provisionally be called the “large particle” method [ll. We briefly describe the computational process. The initial system of equations of gas dynamics is described for the equations of motion and energy, in Lagrange’s form d
z
s
pYdz
= -..
d
s
ai-+ (0 pEdz
pds,
t (0
=
s PV d% s (0
(1.1)
and for the continuity equation, in Eulerian form
s +zdr dp
=-
s
pVdS.
(1.2)
s
Here p, p, E are respectively, the pressure, density and total energy; V = {u, ul is the velocity vector; the integration in (1.1) is carried out with respect to the fluid (Lagrangian), and in (1.2) with respect to a fixed element. The equation of state must also be added to the system (1.1)~(1.2). perfect gas it has the form p=p(p,E--VZ)
=
For a
(x-1),p(E49/2),
where x is the adiabatic exponent. In the design region ABCDEO (Fig. 1) we introduce the cylindrical
system
Numerical
FIG. 1.
solution
of some problems
in gas dynamics
161
FIG. 2.
of coordinates I, z, and consider a fixed (Eulerian) computing net with rectangular cells Ar, AZ, using the simplest finite difference app~ximation of equations (1.1) and (1.2). The flow parameters at the instant n + 1 are calculated from their state at the instant n in three stages. At the first stage we compu~ the in~~~ia~ pa~rne~~ U, u, E of the fluid Lagrangian elements coinciding at the instant n with the fixed Eulerian cells. For example, for the projection of the equations of motion on the z-axis we have
and so on: At the second stage of this approach we use the continuity equation, from which we determine the values of the density of the Eulerian cells at the instantn +I:
Here the quantities A&f?+, j are the mass flows through the boundaries of the corresponding cells. In calculating the values of hMn it is assumed that the flow of mass across the boundary carries with it intermediate values of velocity and specific energy. Tlren, for example, for the value of AM:_,,*, j
162
0. M. Belotserkovskii
we obtain an expression
et al.
of the form
And finally, at the third stage, from the conservation laws, we calculate the parameters a, u, E of the flow on the time layer n + 1 for the cells of the fixed net. Here the final expression
can be written in the form
Here X = U, U, E and C = 0 or 1, respectively, the cell boundary considered.
if the flow flows out or in across
The region in which the calculations are carried out, and also the computing net are shown in Fig. 1. As initial data we here use the parameters of the unperturbed flow. The boundary conditions were taken as follows: on the left AB and upper BC boundaries the conditions in the incident flow of the gas were taken; on the axis of symmetry A0 and on the body OED the conditions of symmetry of the stream and the usual boundary conditions at a wall (no flow conditions) were used; at the right boundary of the region CD the stream parameters were extrapolated outside the domain considered. The state of the flow at succeeding instants up to the onset of the steady state is calculated cyclically from the initial conditions and the difference scheme taking into account the boundary conditions. The number of time steps necessary to obtain a stationary solution depends on the chosen space step and the stream regime. As a rule, the build-up begins last in the neighbourhood of a stagnation point. Figure 2 shows as an example the nature of the build-up of p in the neighbourhood of a stagnation point for M, = 4.1. Questions connected with obtaining time stability of the solution are central in this method. Here they are decided by the introduction of terms of the “artificial viscosity” type into the difference equations. At the present time schemes of a higher order of spatial accuracy based on representations of the method of integral relations are also being considered.
Numerical solution of some problems in gas dynamics
163
Calculations of the s~eam~ining by a su~rs~ic flow of a perfect gas = 1.4) at a zero angle of attack of a cylinder with a flat butt for M, = 1.2 + 03 were carried out by this method on the BESM-6 computer. Figure 1 shows the position of the shock waves and of the sound lines for different values of M,, The shape of the shock waves was determined as the curves on which the derivative of the density with respect to one of the spatial directions has a maximum. The results of an experiment for the case M, = 2 and 4.1 are shown there by dots. Figure 3 shows the density distribution along the axis of symmetry of the body. The values of the density at a direct compression shock and at a stagnation point are shown by dots. It is seen that the discontinuities caused by the jumps are “blurred” spatially over several cells, and for smaller values of the numbers M, the band of blurring of the wave is wider. We notice that the method explained above is more economical (in time and memory) than Harlow’s method of particles. Recently another explicit scheme of the build-up method using an approach and foundation of the problems similar to that explained above, has also been developed and tested on a number of problems in gas dynamics. The scheme is based on the use of compatibility conditions along the curves of intersection of the coordinate and characteristic surfaces and a fixed spatial net. It is explained and studied in sufficient detail in 121 for a first-order system of quasilinear equations of hyperbolic type, hence we here present only some computing formulae for the twodimensionaf non-stationary equations of gas dynamics and illustrate it by a number of examples.
164
0. M. Belotserkovskii
et al,
The initial system of equations comprises, firstly, the differential equations of motion and continuity in Euler’s form, written in matrix form, and secondly, the energy equation written along the trajectory of the particles (that is, in characteristic form), and also the equation of state. The equations of continuity and motion contain the derivatives of the unknown functions along the particle trajectories and two spatial directions. In consequence of their hyperbolicity in time, these equations can also be reduced to a normal kharacteristi c) form, with the result that in the one case they contain derivatives along some “characteristics” and one of the spatial directions, for example the x-axis, and in the other case, along “characteristics” and the other spatial direction. The initial equations are approximated by finite differences only along the particle trajectories, and the compatibility conditions only along the characteristic directions. From these 10 difference equations (only 4 of which are linearly independent), we eliminate the 6 spatial derivatives and determine the 4 unknown parameters at some point Q to be computed on the time layer n + l(n = 1, 2, 3,. . .) in terms of the parameters on the preceding layer n. For the simplest approximation the expressions point Q have the form P=m+F~+PsfP6-2Po 2
-!- P@o c + -PI--2 uQ
1
Pip-:
v,- 04 ux - & 2 -I- - 2 +%&&>
+
0 (At21;
‘1-k 2
2pQ aQ
vQ =
for the parameters at the
+t+
U2+f
At-
2o
’
+ ftQAt;
QQ hQ
Here h and a are respectively
=
ho+
(PQ-PO)
/pa
the enthalpy and speed of sound.
The form of the functions f,, f,, f, depends on the specific choice of the In particular, for a coordinate system and the nature of the non-stationarity. Cartesian coordinate system and a stationary stream f, = f, = f, = 0. The parameters with the subscript i (i = 0, 1, 2, 3, 4) are determined by interpolation at the points of intersection of the characteristics issuing from the point Q, with the plane t = nAt. The formulation of problems for performing calculations
by this scheme is
Numerical
solution
of some problems
similar to that explained above. The difference wave by means of known Hugoniot relations.
in gas dynamics
consists
165
in isolating the shock
This method was used to compute the streamlining of a sphere by the supersonic stationary flow of a perfect gas for different M, numbers, and also that of a number of other blunt bodies at a zero angle of attack. As an illustration the crosses in Fig. 3 show the results of calculations using this scheme. 2. At the present time two schemes: scheme I and scheme II, of the Dorodnitsyn-Delotserkovskii method of integral relations for calculating the supersonic streamlining of blunt bodies, have been developed in detail and are widely used. In scheme I the gas dynamic equations are approximated transverse to the boundary layer, and in case II along it. The resulting approximating system of ordinary differential equations is integrated in scheme I downstream, and in scheme II from the wave to the body. Certain difficulties arise in the use of these schemes in a number of cases, and it is sometimes practically impossible to find a solution. In scheme I these difficulties are connected with the presence of singular points at which it is necessary to satisfy conditions of regularity of the solution and numerically to (Lby-passn these points. Moreover, if the variation of the parameters transverse to the shock layer is of a complex nature, high approximations are necessary which lead to a large expenditure of computer time. In scheme II when the velocity of the incident flow decreases, a singularity begins to develop similar to that of case I, only now for the velocity component normal to the body. This leads to a position where the dependence of the solution for the initial data on the shock wave is increased as the values of M, decrease, and consequently the determination of the initial parameters for the solution of the boundary problem is hampered. It must be mentioned that beginning with some value of M, the continuous dependence of the solution on the data at the shock wave practically disappears altogether. Hence, to overcome the difficulties enumerated above a new numerical scheme based on the method of integral relations, which can provisionally be called scheme III [3l, has recently been developed. In this scheme the functions are approximated simultaneously in two directions and the initial system of partial differential equations is replaced by an approximating system of algebraic equations. Unlike the usual finite difference schemes, in this case it is possible to use a straight through approximation over all the nodes of the computing net for the divergent form of the initial system of equations. What are approximated are not the gas dynamical variables themselves, but some combinations of them in the i&grand. This
166
0. M. Belotserkouskii
d
et al.
Mao=2
1
5
-2.
-7
-
i
1
d
FIG. 5.
FIG. 4.
makes
it possible
comparatively
to obtain
small
a good representation
number of nodes
functions
for a
net, and consequently,
of unknown parameters.
for a small number
The use of the method of integral leads
of the initial
of the computing
to the approximating
system
relations
being
to construct
completely
closed,
the given scheme at the same time
the region of influence is exactly accounted for and all the boundary conditions of the original problem are satisfied to any approximation. The approximating system in scheme III is a system of non-linear algebraic equations of comparatively low order. The program for computing the streamlining scheme III has a simple logic and occupies a comparatively
of blunt bodies by small amount of
computer storage. by some iteration
equations is sought requires small
amounts
The solution of the system of non-linear method, and as is shown by calculations,
of computer
Calculatuons in the incident
time.
by scheme
III were carried
flow, but special
attention
out for a wide range of conditions
was paid to streams
at low supersonic
speeds (M,d 2). The results of calculations in this range of speeds was the best check of the advantages of scheme III (apart from the undoubted practical interest of these results). Calculations were performed down to M, = 1.05.
Numerical
solution
of some problems
in gas dynamics
167
FIG. 6.
The to 1, the different the same computer
calculations showed that scheme III remains stable even for M, close iterations converge fairly quickly, and the computing technique of the cases, for example, for M, = 1.05 and M, = 10, is quite similar and for size of approximation requires approximately identical expenditures of time.
Some of the results obtained are shown in Figs. 4-6. Calculations were carried out for the 4-th approximation both transverse to the shock layer, and also along it W x M = 4 x 4). Figure 4 shows the position and shape of the shock waves, the sound lines and the boundary characteristics of the first and second families for different values of M,. As is obvious, as the Mach number is decreased the picture of the flow changes noticeably. The sound lines deviate sharply from the axis of symmetry downstream, turning about the sound point, the position of which depends weakly on Moo. The region of the influence of the supersonic part of the flow on the subsonic part is bounded by two characteristics of different families, touching a sound line, the point of contact being shifted
168
0. M. Belotserkovskii
et al.
towards the shock wave as the velocity of the incident flow decreases, and the principal changes of the gas dynamic quantities occur only close to the surface of the body. Noticeable bending of the stream lines also occurs only in the immediate proximity of the body. It is interesting to notice that this type of minimal region of influence (bounded by two characteristics) holds for the streamlining of a sphere down to M, = 1.1. Figures 5-6 show a comparison between the computed data of F. D. Popov and the experimental data of V. G. Maslennikov, A. P. Bedin, G. I. Mishin and others [41 on the deviation of the shock wave along the axis of symmetry (c,,), and also the shape and position of the waves. It is obvious that the agreement between the computed data (full line) and experimental data (circles and triangles) is very good. Tables of the fields of flow of these streamlining modes are given in [31. 3. We now consider the problem of calculating the flows of a viscous gas for small values of the Reynolds numbers. In a number of cases the influence of viscosity can be so important that it is impossible to distinguish the stream regions corresponding to the Lexternal” non-viscous flow and to the boundary layer. It is here necessary to regard the flow as completely viscous and study it with the aid of the NavierStokes equations. The method of integral relations (all three types of schemes) is also used to investigate the streamlining of blunt bodies by the flow of a viscous gas. Here we first consider the system of “shortened” Navier-Stokes equations, in which there are no terms with second derivatives along the surface [3l. As calculations have shown, these terms play a negligible part right down to small values of the Re number, while their absence enables a unique solution of the problem to be obtained without artificially introduced closing boundary conditions downstream. The boundary value problem is here solved in a finite region, the boundary conditions on the external boundary of the region (“shock wave”) being obtained from a consideration of the equations of the structure of the shock wave to a first approximation in an expansion with respect to the Re parameter. For the efficient application of the method of integral relations the coordinates of the region of the viscous shock layer are stretched by means of a special transformation, the degree of deformation not being specified beforehand, but obtained during the solution as a function of the coordinate s. This approach enables a definite “smoothing” of the functions represented to be carried out. The approximating system for scheme I has the form S(z, s)dz/ds
= b(z, s),
(3.1)
Numerical
solution
of some problems
0.6
FIG.
7.
in gas dynamics
as 5
Re = 25;
Scheme III: * - sj = 0.25, 0.5, 0.77 0 _ Sj = 0.3. 0.6, 0.8 X _ Sj = 0.2, 0.4, 0.6, 0.8;
Scheme
m = 4;
I
a FIG. 8.
I - Re = 50; 2 - Re = 25; 3 - Re = 10; 4 - Re = 10; 5 - Re = 5.
169
170
0. hf. Belotserkouskii
et al.
where S is a matrix, and b is the vector of the right sides, the components of which are the unknown functions, that is, the values of u, u, p, h at the edges of the bands and the geometrical parameters. As in the case of a non-viscous gas, the conditions on the axis of symmetry and the conditions of regularity at the singular points, serve as the boundary conditions of this quasilinear system of equations. The app~ximat~g integral identities
system of equations for scheme II is obtained from the
in which the subscript j denotes the values of the vectors z on the rays s = sj Cj = 0, 1,. . . , rd. The i&grands are replaced by interpolation polynomials for even and odd powers of s. Calculations by scheme I were performed mainly for plane and axially symmetric bodies with a head insulated surface, for various values of the numbers M, and Re. Calculations by scheme III were performed both to compare the efficiency of the two schemes, and also to investigate the aerodynamic characteristics of spherical blunt bodies with a cooled surface. In the last case the calculations were performed for values of the exponent of degree o for a variation of viscosity p” ‘u ho, given by o = 0.5, 0.67, and 1. The system of algebraic equations into which the relations (3.2) are transformed, is solved by the method of local variations (for example, i31). Some results of the calculations of viscous flows are shown below. Figure 7 shows the values of the flow parameter at the edges of the strips in the case = 1.4, of scheme I and scheme III (a sphere, M, = 10, Re = 25, Pr = 0.72, Wz/&)n = * = Ok Here m is the number of rays and s, is the coordinate of a boundary ray. As is seen from the graphs, the results for three and four rays are practically the same and do not differ from the results for scheme I in the region up to the singular points. It turns out that the expenditure of computer time when scheme III is used is considerably less and the computing process is here better adapted to automation. Figures 8-9 show the values of the principal parameters of the flow on the surface of the body, obtained for different Reynolds numbers. Measurements along the contour of the body of the pressure p (sf (relative to pressure at the critical point p 1011,the enthalpy h, and the coefficient of friction Cf are shown
Numerical
FIG. 9.
solution
2 -hW =h,
of some problems
= 0.025;
in gas dynamics
171
2 - hW = 0.2; 3 - hW = h,.
in Fig. 8 for the case of a heat insulated sphere. Also shown is the variation of the relative heat flow Q(.s)/q (0) (q (0) is the heat flow at the critical point) along the surface of the spherical blunting, the temperature of which equals the temperature of the unperturbed flow. The dashed and dash-dot curves in Fig. 8 correspond to the limiting cases Re = 0 (the free molecular stream, diffuse reflection, the coefficient of accommodation equals 1, Moo>>1) and Re = 00 (ideal gas, M, = 10). It follows from the results obtained that the distributions of the relative values p (s)/p (0) and q (s)/q (0) are almost independent of the Reynolds number and are well approximated by the corresponding limiting dependencies. Hence the variation of the quantities p (0) and q (0) with the Reynolds number is of interest. In Fig. 9 the pressure at the critical point of a heat-insulated and cold sphere is represented by the ratio p (0)/p;, where pi is the stagnation pressure after a direct jump in an ideal gas; the dashed lines indicate the limiting values of this ratio for free molecular streamlining. It is obvious from these graphs that the quantity p (0)/p; (close to unity for moderate Reynolds numbers) increases sharply in the band Re < 10 as the rarefaction of the medium increases, and as Re + 0 approaches its own limit, depending on the cooling of the surface.
172
0. M. ~elotse~kou~k~i
FIG. 10.
et al.
FIG. 11.
Figure 9 also shows the variation with the Reynolds number of the heattransfer coefficient at the critical point Ch = Q(O)/(ho - hw) (h, is the stagnation enthalpy), obtained for two values of the exponent in the viscosity variation relation (o = 0.5 and 1). The value of o has considerable effect on the heat flow. However, independently of the viscosity variation, the heat transfer coefficient (as also the pressure) approximates to its limit in a free molecular stream when the Reynolds number is reduced. Therefore, the results obtained focally for small Reynolds numbers (Re = 11 are found to be in satisfac~ry qualitative agreement with the limiting picture of streamlining. Hence, it may be supposed that for the problem considered the region of applicability of the Navier-stokes equations extends considerably further than its formal limits. 4. In conclusion we consider the streamlining of a blunt body by the hypersonic flow of a selectively radiating and absorbing gas. For sufficiently great Mach numbers (M, > 25) the radiation can have considerable effect on the stream field in the shock layer. We will assume that the gas behind the shock wave is in ~uilibrium dissociation and ionisation; we
Numerical
neglect the viscosity
solution
of some problems
173
in gas dynamics
and heat conduction.
We use scheme II of the method of integral relations for calculating the streamlining of blunting. The initial system of equations consists of the equations of gas dynamics (where in the energy equation there occurs a term allowing for the radiation and absorption of energy) and the equation of the transfer of radiactive energy. The system here is integrodifferential. To construct the approximate (differential) representation of the equations of the transfer of radiative energy it is possible to use the P,-th approximation of the method of spherical harmonics, which is widely used in the theory of neutron transfer. By means of this approximation the transfer equation is replaced by two equations for the vector flow of radiative energy H and an integral in the direction of the intensity of radiation J,. The system of equations of the problem considered is written in dimensionless form just as in [61. The equations of gas dynamics:
P =
PWP
(~9
T).
The transfer equations:
&Jo, v= - 3k’(A&
-&Jo, y = - W%,
8 -I- Ess’C, E),
E = 2 [k,‘(J,,
v-
E,
4nB,)],Av+
Here
a=
E = (AZ - &‘g)r, f =
tr,
t =
pu,
(Ah - ,&‘t)r,
X=-Z/R+APsinO, pv,
Z=
zr,
A=i+CeIR,
F=P+P~~,
I=
m=
(AH - Edz)r,
H=
P+PU",
c=n/&(s), z = pUllr
8 = gr,
P=dc0se+~e,
Jo,v=
(4.2)
174
0. M. Belotserkovskii
et al.
0 FIG. 12. EfO,R,=lm,M,=33;-.-.-.EfO,R,=O.lm.M,=33;
0.5
l-E=E
FIG. 13. E&O; ---E=O. ---E=
0, Mm = 33.
Hv s
andHv5 are the components of the vector flow of radiative energy H, along the axes of s’ and E respectively; the primes denote differentiation with respect to s. It is known that air at high temperatures radiates and absorbs energy selectively, the absorption coefficients kv’ (p, T) depending in a complex manner on V, p and T. The coefficients obtained in [7l were used in the calculations. The range of variation of k,,’ was divided into n intervals (in the calculations n = 20),and the values of k,,’ were recorded at the nodal points for v, p and T. The values of kv’ required in the computing process for intermediate values of p and T, were found by linear interpolation with respect to p and T on a logarithmic scale. The system of equations (4.1M4.2) was solved as follows (the second approximation was considered). The equations of motion and continuity are the same as in scheme II ignoring radiation; the energy equations are either integrated twice, or approximately represented along the curves s, (the limiting
175
Numericalsolutionof some problems in gas dynamics
FIG.14. O-E
CO;
A-E =0,x/R = 3.
~h~ac~ristic} and s, = s/Z in terms of the corresponding values on the streamlines. The diffe~nti~ transfer equations are inbred twice from s = 0 to s = s, and s = s,, where to obtain a closed approximating system it is necessary to have 8 equations for each interval (8 unknown functions Hgf, Hat, Joi, i= 0, = 0). The remaining 4 relations are obtained by writing the total 1,2;&I derivatives of the unknown functions with respect to $ in terms of the known form of the partial derivatives of these fictions with respect to 5 and s. The boundary conditions for the gas dynamic functions are the same as those ignoring radiation. The boundary conditions for the quantities characterizing the radiation are assigned approximately on the assumption that the radiation in the shock layer does not enter from the incident flow; the surface of the body absorbs all the radiative energy of the gas, and radiates a quantity of energy negligibly small in comparison with the radiationof the shock layer. These conditions can be written in the form
Jvoi- 2&&i = 0
(shock wave)
3voj+
ZH,&i = 0
(body)
The numerical integration of the approximating system was carried out as in [3, 61.
’
As a typical We will discuss osme of the results of the calculations. example we consider the supersonic streamlining of a sphere of radius R, = 1 m by a flow of radiating-absorbing air for the following values of the parameters of
176
0. M. ~e~otserkovsk~i
yt al.
the incident flow: Moo= 33, p, = 0.0029 atm, T, = 257°K. Figure 10 shows the variation of the temperature profile from the wave (4 = 1) to the body (t = 0) along the axis of symmetry s = 0, and also along the lines s = s, (the limiting characteristic) and s = s/2. The dashed curve represents the temperature profile obtained for the case of a non-radiating gas (E = 0); the continuous line shows the temperature profile for calculations allowing for selective absorption and radiation (E f 01; the dash-dot curves show the results of calculations of selective de-excitation using the data given in El, and the dashdot curves with two dots show the results for de~x~itation using the data given in [sl (E f 0). Figure 11 shows similar results for the density profile. As we see, in all the cases considered the radiation leads to a decrease in the temperature and an increase in the density in the shock layer, Calculations have shown that the radiation has practically no effect on the pressure p and leads to a negligible decrease in the velocity w close to the surface. Therefore, radiation can have a considerable effect on the radiative heat flow at the surface of the body and has no noticeable effect on its aerodynamic characteristics. The departure of the shock wave from the body is somewhat decreased in consequence of the lowering of temperate and the increase in density in the boundary layer (Fig. 12). Figure 13 shows the variation of temperature in the shock layer (F= 0 is the wave, F = 1 is the body) in a gas mixture of 9% Co,+ 91% N,. Calculations show that the qualitative picture of the variation of p and 7’ remains the same as in the case of air. On the basis of the calculations carried out for the different methods of approximating to the energy equation, it can be concluded that the methods are very accurate. Estimates show that in a significant part of the region considered the error of the ~alc~ations amo~~d to a qu~tity of about 1%. The solution obtained for the region of blunting was continued into the supersonic zone by the method of characteristics especially adapted for calculating the flow of a radiating gas [31. The results were obtained on the assumption of volume de-excitation. Figure 14 shows the variation of temperature, density and pressure in the shock layer on the lateral surface of a cone with spherical blunting with an aperture angle w = 30° in the section x/R, = 3 (M, = 33, p, = 0.003 atm, T = 257*K, R, = 1 m). As we see, accosting for the effect of radiation leads to a fall in tem~ratu~ and an increase in density in a narrow zone situated close to
Numerical
solution
of some problems
in gas dynamics
177
the surface of the cone. In the remaining part of the shock layer the parameters of the gas remain unchanged. Translated
by J. Berry
REFERENCES
1.
RICH, M. A method for Eulerian No. LAMS-2826, 1963.
fluid dynamics.
Los
AIamos
2.
MAGOMEDOV, K. M. and KHOLODOV, A. S. Construction hyperbolic type on the basis of characteristic relations. Fiz., 9, 2, 373-386, 1969.
Scient.
Lab.,
of difference Zh. uychisl.
Rep.
schemes
of
Mat. mat.
3.
BELOTSERKOVSKII, 0. M. (editor) Streamlining of Blunt Bodies by a Supersonic Flow of Gas. Theoretical and Experimental Investigations. Trudy VTs Akad. Nauk SSSR. 1967.
4.
DUNAEV, YU. A. (editor). Aerophysical Studies of Supersonic Flows (Aerofizicheskie issledovaniya sverkhzvukovykh techenii). “Nauka”, Moscow-Leningrad, 1967.
5.
CHERNOUS’KO, F. L. variational problems.
Method of local
variations
for the numerical 5, 4, 749-754,
Zh. vy’chisl. Mat. mat. Fiz..
solution 1965.
of
6.
FOMIN, V. N. Streamlining of blunt bodies by a hypersonic flow of gas taking radiation into account. Zh. uychisl. Mat. mat. Fiz., 6, 4. 714-726, 1966.
7.
BIBERMAN, Radiation 1964.
8.
KIVEL,
B.
Aero Sci.,
L. M.. VOROB’EV, V. S., NORMAN, heating in hypersonic streamlining.
Radiation from hot air and its effect 28, 2, 96-102, 1961.
G. E., YAKUNOV,
Kosmich.
issl.,
on stagnation
I, T.
2, 3. 441-454,
point heating.
J.
20 April
1969)
Notation s is the length of arc measured along the contour of the body; n is the normal to the body (the contour of the body corresponds to n = 0); t is the time; R is the radius of curvature of the body; r is the distance from the axis of symmetry of the body; t is the distance along the normal from the surface of the body to the shock wave, [=&c(s), & 1 -[; is the angle between the tangent to the contour of the body and the direction of the incident flow; u and u are the velocity components in the direction n and s respectively; w is the magnitude of the total velocity; p is the pressure; p is the density; h is the enthalpy; ‘I’ is the temperature; p is the molecular weight; p” is the coefficient of viscosity; Moois the Mach number of the incident flow; Pr is the Prandtl number; Re is Reynolds number, calculated for the parameters of the incident flow, the characteristic diameter of the body and the viscosity for h = w,‘; is the coefficient
*Zh.
u’jchisl.
Mat. mat. Fiz.,
10,2, 401-416,
1970.
158
of friction;
NumericuZ
solution
of some problems
in gas
dynamics
159
is the heat flow;
v is the frequency; hv’ is the absorption coefficient Bv is the Planck function.
taking induced emission into account;
Recently the study of problems in gas dynamics by means of numerical methods using high-speed computers has been intensively continued. Further improvement has been made by means of serial calculations by well-known methods, which have also led to the development of new numerical algorithms defined by the consideration of a number of new problems in gas dynamics of scientific and technical interest. The latter include problems connected, for example, with the consideration of the streamlining of blunt bodies by the flow of a gas at low supersonic speeds (1 < M, < 2); the determination of the characteristics of blunt bodies in the presence of discontinuities of the generator; the study of stream regimes in the presence of viscosity, physico-chemical transformations, radiation etc. Below we describe the results of some investigations of this type of problem using new numerical schemes of the method of integral relations, finite differences, combinations of them, and also build-up methods. The following four types of problems are considered: (1) the streamlining of an axially symmetric cylindrical the supersonic flow of a perfect gas; (2) the streamlining
of blunt bodies at low supersonic
butt and sphere by
speeds;
(3) the streamlining of blunt bodies by the hypersonic flow of a viscous perfect gas for low values of the Re numbers; (4) the streamlining of blunt bodies by a supersonic flow of an equilibrium dissociated and ionised gas taking selective radiation and absorption into account. Here the main attention is given to the results of the investigations. In the description of the numerical methods only the fundamental ideas of the construction of the corresponding schemes are presented and the necessary references are given. This paper forms the main content of a report by the authors at the Third Aii-Union conference on Theoretical and Applied mechanics (MOSCOW, January 1967).
160
0. M. Belotserkouskii
et al.
1. As is well known,. at the present time different numerical schemes for solving the non-stationary equations of gas dynamics have been developed which are also used successfully to obtain the stationary characteristics of the stream by the build-up method. Approaches of this type are now possible in practice due to the introduction of computers of comparatively large capacity. Moreover, it has to be borne in mind that in some cases the use of “classical” stationary methods, well developed earlier, encounters certain difficulties in the solution of a number of important problems of gas dynamics (for example, the investigation of transonic streams, obtaining the stream field for bodies with discontinuities etc). We consider one of these %on-stationary” approaches in the classical problem of the streamlining of an axially symmetric cylindrical butt by the steady supersonic flow of a gas. This problem has been solved by using the modified non-stationary Euler-Lagrange method (close to the discrete model of Harlow), which can provisionally be called the “large particle” method [ll. We briefly describe the computational process. The initial system of equations of gas dynamics is described for the equations of motion and energy, in Lagrange’s form d
z
s
pYdz
= -..
d
s
ai-+ (0 pEdz
pds,
t (0
=
s PV d% s (0
(1.1)
and for the continuity equation, in Eulerian form
s +zdr dp
=-
s
pVdS.
(1.2)
s
Here p, p, E are respectively, the pressure, density and total energy; V = {u, ul is the velocity vector; the integration in (1.1) is carried out with respect to the fluid (Lagrangian), and in (1.2) with respect to a fixed element. The equation of state must also be added to the system (1.1)~(1.2). perfect gas it has the form p=p(p,E--VZ)
=
For a
(x-1),p(E49/2),
where x is the adiabatic exponent. In the design region ABCDEO (Fig. 1) we introduce the cylindrical
system
Numerical
FIG. 1.
solution
of some problems
in gas dynamics
161
FIG. 2.
of coordinates I, z, and consider a fixed (Eulerian) computing net with rectangular cells Ar, AZ, using the simplest finite difference app~ximation of equations (1.1) and (1.2). The flow parameters at the instant n + 1 are calculated from their state at the instant n in three stages. At the first stage we compu~ the in~~~ia~ pa~rne~~ U, u, E of the fluid Lagrangian elements coinciding at the instant n with the fixed Eulerian cells. For example, for the projection of the equations of motion on the z-axis we have
and so on: At the second stage of this approach we use the continuity equation, from which we determine the values of the density of the Eulerian cells at the instantn +I:
Here the quantities A&f?+, j are the mass flows through the boundaries of the corresponding cells. In calculating the values of hMn it is assumed that the flow of mass across the boundary carries with it intermediate values of velocity and specific energy. Tlren, for example, for the value of AM:_,,*, j
162
0. M. Belotserkovskii
we obtain an expression
et al.
of the form
And finally, at the third stage, from the conservation laws, we calculate the parameters a, u, E of the flow on the time layer n + 1 for the cells of the fixed net. Here the final expression
can be written in the form
Here X = U, U, E and C = 0 or 1, respectively, the cell boundary considered.
if the flow flows out or in across
The region in which the calculations are carried out, and also the computing net are shown in Fig. 1. As initial data we here use the parameters of the unperturbed flow. The boundary conditions were taken as follows: on the left AB and upper BC boundaries the conditions in the incident flow of the gas were taken; on the axis of symmetry A0 and on the body OED the conditions of symmetry of the stream and the usual boundary conditions at a wall (no flow conditions) were used; at the right boundary of the region CD the stream parameters were extrapolated outside the domain considered. The state of the flow at succeeding instants up to the onset of the steady state is calculated cyclically from the initial conditions and the difference scheme taking into account the boundary conditions. The number of time steps necessary to obtain a stationary solution depends on the chosen space step and the stream regime. As a rule, the build-up begins last in the neighbourhood of a stagnation point. Figure 2 shows as an example the nature of the build-up of p in the neighbourhood of a stagnation point for M, = 4.1. Questions connected with obtaining time stability of the solution are central in this method. Here they are decided by the introduction of terms of the “artificial viscosity” type into the difference equations. At the present time schemes of a higher order of spatial accuracy based on representations of the method of integral relations are also being considered.
Numerical solution of some problems in gas dynamics
163
Calculations of the s~eam~ining by a su~rs~ic flow of a perfect gas = 1.4) at a zero angle of attack of a cylinder with a flat butt for M, = 1.2 + 03 were carried out by this method on the BESM-6 computer. Figure 1 shows the position of the shock waves and of the sound lines for different values of M,, The shape of the shock waves was determined as the curves on which the derivative of the density with respect to one of the spatial directions has a maximum. The results of an experiment for the case M, = 2 and 4.1 are shown there by dots. Figure 3 shows the density distribution along the axis of symmetry of the body. The values of the density at a direct compression shock and at a stagnation point are shown by dots. It is seen that the discontinuities caused by the jumps are “blurred” spatially over several cells, and for smaller values of the numbers M, the band of blurring of the wave is wider. We notice that the method explained above is more economical (in time and memory) than Harlow’s method of particles. Recently another explicit scheme of the build-up method using an approach and foundation of the problems similar to that explained above, has also been developed and tested on a number of problems in gas dynamics. The scheme is based on the use of compatibility conditions along the curves of intersection of the coordinate and characteristic surfaces and a fixed spatial net. It is explained and studied in sufficient detail in 121 for a first-order system of quasilinear equations of hyperbolic type, hence we here present only some computing formulae for the twodimensionaf non-stationary equations of gas dynamics and illustrate it by a number of examples.
164
0. M. Belotserkovskii
et al,
The initial system of equations comprises, firstly, the differential equations of motion and continuity in Euler’s form, written in matrix form, and secondly, the energy equation written along the trajectory of the particles (that is, in characteristic form), and also the equation of state. The equations of continuity and motion contain the derivatives of the unknown functions along the particle trajectories and two spatial directions. In consequence of their hyperbolicity in time, these equations can also be reduced to a normal kharacteristi c) form, with the result that in the one case they contain derivatives along some “characteristics” and one of the spatial directions, for example the x-axis, and in the other case, along “characteristics” and the other spatial direction. The initial equations are approximated by finite differences only along the particle trajectories, and the compatibility conditions only along the characteristic directions. From these 10 difference equations (only 4 of which are linearly independent), we eliminate the 6 spatial derivatives and determine the 4 unknown parameters at some point Q to be computed on the time layer n + l(n = 1, 2, 3,. . .) in terms of the parameters on the preceding layer n. For the simplest approximation the expressions point Q have the form P=m+F~+PsfP6-2Po 2
-!- P@o c + -PI--2 uQ
1
Pip-:
v,- 04 ux - & 2 -I- - 2 +%&&>
+
0 (At21;
‘1-k 2
2pQ aQ
vQ =
for the parameters at the
+t+
U2+f
At-
2o
’
+ ftQAt;
QQ hQ
Here h and a are respectively
=
ho+
(PQ-PO)
/pa
the enthalpy and speed of sound.
The form of the functions f,, f,, f, depends on the specific choice of the In particular, for a coordinate system and the nature of the non-stationarity. Cartesian coordinate system and a stationary stream f, = f, = f, = 0. The parameters with the subscript i (i = 0, 1, 2, 3, 4) are determined by interpolation at the points of intersection of the characteristics issuing from the point Q, with the plane t = nAt. The formulation of problems for performing calculations
by this scheme is
Numerical
solution
of some problems
similar to that explained above. The difference wave by means of known Hugoniot relations.
in gas dynamics
consists
165
in isolating the shock
This method was used to compute the streamlining of a sphere by the supersonic stationary flow of a perfect gas for different M, numbers, and also that of a number of other blunt bodies at a zero angle of attack. As an illustration the crosses in Fig. 3 show the results of calculations using this scheme. 2. At the present time two schemes: scheme I and scheme II, of the Dorodnitsyn-Delotserkovskii method of integral relations for calculating the supersonic streamlining of blunt bodies, have been developed in detail and are widely used. In scheme I the gas dynamic equations are approximated transverse to the boundary layer, and in case II along it. The resulting approximating system of ordinary differential equations is integrated in scheme I downstream, and in scheme II from the wave to the body. Certain difficulties arise in the use of these schemes in a number of cases, and it is sometimes practically impossible to find a solution. In scheme I these difficulties are connected with the presence of singular points at which it is necessary to satisfy conditions of regularity of the solution and numerically to (Lby-passn these points. Moreover, if the variation of the parameters transverse to the shock layer is of a complex nature, high approximations are necessary which lead to a large expenditure of computer time. In scheme II when the velocity of the incident flow decreases, a singularity begins to develop similar to that of case I, only now for the velocity component normal to the body. This leads to a position where the dependence of the solution for the initial data on the shock wave is increased as the values of M, decrease, and consequently the determination of the initial parameters for the solution of the boundary problem is hampered. It must be mentioned that beginning with some value of M, the continuous dependence of the solution on the data at the shock wave practically disappears altogether. Hence, to overcome the difficulties enumerated above a new numerical scheme based on the method of integral relations, which can provisionally be called scheme III [3l, has recently been developed. In this scheme the functions are approximated simultaneously in two directions and the initial system of partial differential equations is replaced by an approximating system of algebraic equations. Unlike the usual finite difference schemes, in this case it is possible to use a straight through approximation over all the nodes of the computing net for the divergent form of the initial system of equations. What are approximated are not the gas dynamical variables themselves, but some combinations of them in the i&grand. This
166
0. M. Belotserkouskii
d
et al.
Mao=2
1
5
-2.
-7
-
i
1
d
FIG. 5.
FIG. 4.
makes
it possible
comparatively
to obtain
small
a good representation
number of nodes
functions
for a
net, and consequently,
of unknown parameters.
for a small number
The use of the method of integral leads
of the initial
of the computing
to the approximating
system
relations
being
to construct
completely
closed,
the given scheme at the same time
the region of influence is exactly accounted for and all the boundary conditions of the original problem are satisfied to any approximation. The approximating system in scheme III is a system of non-linear algebraic equations of comparatively low order. The program for computing the streamlining scheme III has a simple logic and occupies a comparatively
of blunt bodies by small amount of
computer storage. by some iteration
equations is sought requires small
amounts
The solution of the system of non-linear method, and as is shown by calculations,
of computer
Calculatuons in the incident
time.
by scheme
III were carried
flow, but special
attention
out for a wide range of conditions
was paid to streams
at low supersonic
speeds (M,d 2). The results of calculations in this range of speeds was the best check of the advantages of scheme III (apart from the undoubted practical interest of these results). Calculations were performed down to M, = 1.05.
Numerical
solution
of some problems
in gas dynamics
167
FIG. 6.
The to 1, the different the same computer
calculations showed that scheme III remains stable even for M, close iterations converge fairly quickly, and the computing technique of the cases, for example, for M, = 1.05 and M, = 10, is quite similar and for size of approximation requires approximately identical expenditures of time.
Some of the results obtained are shown in Figs. 4-6. Calculations were carried out for the 4-th approximation both transverse to the shock layer, and also along it W x M = 4 x 4). Figure 4 shows the position and shape of the shock waves, the sound lines and the boundary characteristics of the first and second families for different values of M,. As is obvious, as the Mach number is decreased the picture of the flow changes noticeably. The sound lines deviate sharply from the axis of symmetry downstream, turning about the sound point, the position of which depends weakly on Moo. The region of the influence of the supersonic part of the flow on the subsonic part is bounded by two characteristics of different families, touching a sound line, the point of contact being shifted
168
0. M. Belotserkovskii
et al.
towards the shock wave as the velocity of the incident flow decreases, and the principal changes of the gas dynamic quantities occur only close to the surface of the body. Noticeable bending of the stream lines also occurs only in the immediate proximity of the body. It is interesting to notice that this type of minimal region of influence (bounded by two characteristics) holds for the streamlining of a sphere down to M, = 1.1. Figures 5-6 show a comparison between the computed data of F. D. Popov and the experimental data of V. G. Maslennikov, A. P. Bedin, G. I. Mishin and others [41 on the deviation of the shock wave along the axis of symmetry (c,,), and also the shape and position of the waves. It is obvious that the agreement between the computed data (full line) and experimental data (circles and triangles) is very good. Tables of the fields of flow of these streamlining modes are given in [31. 3. We now consider the problem of calculating the flows of a viscous gas for small values of the Reynolds numbers. In a number of cases the influence of viscosity can be so important that it is impossible to distinguish the stream regions corresponding to the Lexternal” non-viscous flow and to the boundary layer. It is here necessary to regard the flow as completely viscous and study it with the aid of the NavierStokes equations. The method of integral relations (all three types of schemes) is also used to investigate the streamlining of blunt bodies by the flow of a viscous gas. Here we first consider the system of “shortened” Navier-Stokes equations, in which there are no terms with second derivatives along the surface [3l. As calculations have shown, these terms play a negligible part right down to small values of the Re number, while their absence enables a unique solution of the problem to be obtained without artificially introduced closing boundary conditions downstream. The boundary value problem is here solved in a finite region, the boundary conditions on the external boundary of the region (“shock wave”) being obtained from a consideration of the equations of the structure of the shock wave to a first approximation in an expansion with respect to the Re parameter. For the efficient application of the method of integral relations the coordinates of the region of the viscous shock layer are stretched by means of a special transformation, the degree of deformation not being specified beforehand, but obtained during the solution as a function of the coordinate s. This approach enables a definite “smoothing” of the functions represented to be carried out. The approximating system for scheme I has the form S(z, s)dz/ds
= b(z, s),
(3.1)
Numerical
solution
of some problems
0.6
FIG.
7.
in gas dynamics
as 5
Re = 25;
Scheme III: * - sj = 0.25, 0.5, 0.77 0 _ Sj = 0.3. 0.6, 0.8 X _ Sj = 0.2, 0.4, 0.6, 0.8;
Scheme
m = 4;
I
a FIG. 8.
I - Re = 50; 2 - Re = 25; 3 - Re = 10; 4 - Re = 10; 5 - Re = 5.
169
170
0. hf. Belotserkouskii
et al.
where S is a matrix, and b is the vector of the right sides, the components of which are the unknown functions, that is, the values of u, u, p, h at the edges of the bands and the geometrical parameters. As in the case of a non-viscous gas, the conditions on the axis of symmetry and the conditions of regularity at the singular points, serve as the boundary conditions of this quasilinear system of equations. The app~ximat~g integral identities
system of equations for scheme II is obtained from the
in which the subscript j denotes the values of the vectors z on the rays s = sj Cj = 0, 1,. . . , rd. The i&grands are replaced by interpolation polynomials for even and odd powers of s. Calculations by scheme I were performed mainly for plane and axially symmetric bodies with a head insulated surface, for various values of the numbers M, and Re. Calculations by scheme III were performed both to compare the efficiency of the two schemes, and also to investigate the aerodynamic characteristics of spherical blunt bodies with a cooled surface. In the last case the calculations were performed for values of the exponent of degree o for a variation of viscosity p” ‘u ho, given by o = 0.5, 0.67, and 1. The system of algebraic equations into which the relations (3.2) are transformed, is solved by the method of local variations (for example, i31). Some results of the calculations of viscous flows are shown below. Figure 7 shows the values of the flow parameter at the edges of the strips in the case = 1.4, of scheme I and scheme III (a sphere, M, = 10, Re = 25, Pr = 0.72, Wz/&)n = * = Ok Here m is the number of rays and s, is the coordinate of a boundary ray. As is seen from the graphs, the results for three and four rays are practically the same and do not differ from the results for scheme I in the region up to the singular points. It turns out that the expenditure of computer time when scheme III is used is considerably less and the computing process is here better adapted to automation. Figures 8-9 show the values of the principal parameters of the flow on the surface of the body, obtained for different Reynolds numbers. Measurements along the contour of the body of the pressure p (sf (relative to pressure at the critical point p 1011,the enthalpy h, and the coefficient of friction Cf are shown
Numerical
FIG. 9.
solution
2 -hW =h,
of some problems
= 0.025;
in gas dynamics
171
2 - hW = 0.2; 3 - hW = h,.
in Fig. 8 for the case of a heat insulated sphere. Also shown is the variation of the relative heat flow Q(.s)/q (0) (q (0) is the heat flow at the critical point) along the surface of the spherical blunting, the temperature of which equals the temperature of the unperturbed flow. The dashed and dash-dot curves in Fig. 8 correspond to the limiting cases Re = 0 (the free molecular stream, diffuse reflection, the coefficient of accommodation equals 1, Moo>>1) and Re = 00 (ideal gas, M, = 10). It follows from the results obtained that the distributions of the relative values p (s)/p (0) and q (s)/q (0) are almost independent of the Reynolds number and are well approximated by the corresponding limiting dependencies. Hence the variation of the quantities p (0) and q (0) with the Reynolds number is of interest. In Fig. 9 the pressure at the critical point of a heat-insulated and cold sphere is represented by the ratio p (0)/p;, where pi is the stagnation pressure after a direct jump in an ideal gas; the dashed lines indicate the limiting values of this ratio for free molecular streamlining. It is obvious from these graphs that the quantity p (0)/p; (close to unity for moderate Reynolds numbers) increases sharply in the band Re < 10 as the rarefaction of the medium increases, and as Re + 0 approaches its own limit, depending on the cooling of the surface.
172
0. M. ~elotse~kou~k~i
FIG. 10.
et al.
FIG. 11.
Figure 9 also shows the variation with the Reynolds number of the heattransfer coefficient at the critical point Ch = Q(O)/(ho - hw) (h, is the stagnation enthalpy), obtained for two values of the exponent in the viscosity variation relation (o = 0.5 and 1). The value of o has considerable effect on the heat flow. However, independently of the viscosity variation, the heat transfer coefficient (as also the pressure) approximates to its limit in a free molecular stream when the Reynolds number is reduced. Therefore, the results obtained focally for small Reynolds numbers (Re = 11 are found to be in satisfac~ry qualitative agreement with the limiting picture of streamlining. Hence, it may be supposed that for the problem considered the region of applicability of the Navier-stokes equations extends considerably further than its formal limits. 4. In conclusion we consider the streamlining of a blunt body by the hypersonic flow of a selectively radiating and absorbing gas. For sufficiently great Mach numbers (M, > 25) the radiation can have considerable effect on the stream field in the shock layer. We will assume that the gas behind the shock wave is in ~uilibrium dissociation and ionisation; we
Numerical
neglect the viscosity
solution
of some problems
173
in gas dynamics
and heat conduction.
We use scheme II of the method of integral relations for calculating the streamlining of blunting. The initial system of equations consists of the equations of gas dynamics (where in the energy equation there occurs a term allowing for the radiation and absorption of energy) and the equation of the transfer of radiactive energy. The system here is integrodifferential. To construct the approximate (differential) representation of the equations of the transfer of radiative energy it is possible to use the P,-th approximation of the method of spherical harmonics, which is widely used in the theory of neutron transfer. By means of this approximation the transfer equation is replaced by two equations for the vector flow of radiative energy H and an integral in the direction of the intensity of radiation J,. The system of equations of the problem considered is written in dimensionless form just as in [61. The equations of gas dynamics:
P =
PWP
(~9
T).
The transfer equations:
&Jo, v= - 3k’(A&
-&Jo, y = - W%,
8 -I- Ess’C, E),
E = 2 [k,‘(J,,
v-
E,
4nB,)],Av+
Here
a=
E = (AZ - &‘g)r, f =
tr,
t =
pu,
(Ah - ,&‘t)r,
X=-Z/R+APsinO, pv,
Z=
zr,
A=i+CeIR,
F=P+P~~,
I=
m=
(AH - Edz)r,
H=
P+PU",
c=n/&(s), z = pUllr
8 = gr,
P=dc0se+~e,
Jo,v=
(4.2)
174
0. M. Belotserkovskii
et al.
0 FIG. 12. EfO,R,=lm,M,=33;-.-.-.EfO,R,=O.lm.M,=33;
0.5
l-E=E
FIG. 13. E&O; ---E=O. ---E=
0, Mm = 33.
Hv s
andHv5 are the components of the vector flow of radiative energy H, along the axes of s’ and E respectively; the primes denote differentiation with respect to s. It is known that air at high temperatures radiates and absorbs energy selectively, the absorption coefficients kv’ (p, T) depending in a complex manner on V, p and T. The coefficients obtained in [7l were used in the calculations. The range of variation of k,,’ was divided into n intervals (in the calculations n = 20),and the values of k,,’ were recorded at the nodal points for v, p and T. The values of kv’ required in the computing process for intermediate values of p and T, were found by linear interpolation with respect to p and T on a logarithmic scale. The system of equations (4.1M4.2) was solved as follows (the second approximation was considered). The equations of motion and continuity are the same as in scheme II ignoring radiation; the energy equations are either integrated twice, or approximately represented along the curves s, (the limiting
175
Numericalsolutionof some problems in gas dynamics
FIG.14. O-E
CO;
A-E =0,x/R = 3.
~h~ac~ristic} and s, = s/Z in terms of the corresponding values on the streamlines. The diffe~nti~ transfer equations are inbred twice from s = 0 to s = s, and s = s,, where to obtain a closed approximating system it is necessary to have 8 equations for each interval (8 unknown functions Hgf, Hat, Joi, i= 0, = 0). The remaining 4 relations are obtained by writing the total 1,2;&I derivatives of the unknown functions with respect to $ in terms of the known form of the partial derivatives of these fictions with respect to 5 and s. The boundary conditions for the gas dynamic functions are the same as those ignoring radiation. The boundary conditions for the quantities characterizing the radiation are assigned approximately on the assumption that the radiation in the shock layer does not enter from the incident flow; the surface of the body absorbs all the radiative energy of the gas, and radiates a quantity of energy negligibly small in comparison with the radiationof the shock layer. These conditions can be written in the form
Jvoi- 2&&i = 0
(shock wave)
3voj+
ZH,&i = 0
(body)
The numerical integration of the approximating system was carried out as in [3, 61.
’
As a typical We will discuss osme of the results of the calculations. example we consider the supersonic streamlining of a sphere of radius R, = 1 m by a flow of radiating-absorbing air for the following values of the parameters of
176
0. M. ~e~otserkovsk~i
yt al.
the incident flow: Moo= 33, p, = 0.0029 atm, T, = 257°K. Figure 10 shows the variation of the temperature profile from the wave (4 = 1) to the body (t = 0) along the axis of symmetry s = 0, and also along the lines s = s, (the limiting characteristic) and s = s/2. The dashed curve represents the temperature profile obtained for the case of a non-radiating gas (E = 0); the continuous line shows the temperature profile for calculations allowing for selective absorption and radiation (E f 01; the dash-dot curves show the results of calculations of selective de-excitation using the data given in El, and the dashdot curves with two dots show the results for de~x~itation using the data given in [sl (E f 0). Figure 11 shows similar results for the density profile. As we see, in all the cases considered the radiation leads to a decrease in the temperature and an increase in the density in the shock layer, Calculations have shown that the radiation has practically no effect on the pressure p and leads to a negligible decrease in the velocity w close to the surface. Therefore, radiation can have a considerable effect on the radiative heat flow at the surface of the body and has no noticeable effect on its aerodynamic characteristics. The departure of the shock wave from the body is somewhat decreased in consequence of the lowering of temperate and the increase in density in the boundary layer (Fig. 12). Figure 13 shows the variation of temperature in the shock layer (F= 0 is the wave, F = 1 is the body) in a gas mixture of 9% Co,+ 91% N,. Calculations show that the qualitative picture of the variation of p and 7’ remains the same as in the case of air. On the basis of the calculations carried out for the different methods of approximating to the energy equation, it can be concluded that the methods are very accurate. Estimates show that in a significant part of the region considered the error of the ~alc~ations amo~~d to a qu~tity of about 1%. The solution obtained for the region of blunting was continued into the supersonic zone by the method of characteristics especially adapted for calculating the flow of a radiating gas [31. The results were obtained on the assumption of volume de-excitation. Figure 14 shows the variation of temperature, density and pressure in the shock layer on the lateral surface of a cone with spherical blunting with an aperture angle w = 30° in the section x/R, = 3 (M, = 33, p, = 0.003 atm, T = 257*K, R, = 1 m). As we see, accosting for the effect of radiation leads to a fall in tem~ratu~ and an increase in density in a narrow zone situated close to
Numerical
solution
of some problems
in gas dynamics
177
the surface of the cone. In the remaining part of the shock layer the parameters of the gas remain unchanged. Translated
by J. Berry
REFERENCES
1.
RICH, M. A method for Eulerian No. LAMS-2826, 1963.
fluid dynamics.
Los
AIamos
2.
MAGOMEDOV, K. M. and KHOLODOV, A. S. Construction hyperbolic type on the basis of characteristic relations. Fiz., 9, 2, 373-386, 1969.
Scient.
Lab.,
of difference Zh. uychisl.
Rep.
schemes
of
Mat. mat.
3.
BELOTSERKOVSKII, 0. M. (editor) Streamlining of Blunt Bodies by a Supersonic Flow of Gas. Theoretical and Experimental Investigations. Trudy VTs Akad. Nauk SSSR. 1967.
4.
DUNAEV, YU. A. (editor). Aerophysical Studies of Supersonic Flows (Aerofizicheskie issledovaniya sverkhzvukovykh techenii). “Nauka”, Moscow-Leningrad, 1967.
5.
CHERNOUS’KO, F. L. variational problems.
Method of local
variations
for the numerical 5, 4, 749-754,
Zh. vy’chisl. Mat. mat. Fiz..
solution 1965.
of
6.
FOMIN, V. N. Streamlining of blunt bodies by a hypersonic flow of gas taking radiation into account. Zh. uychisl. Mat. mat. Fiz., 6, 4. 714-726, 1966.
7.
BIBERMAN, Radiation 1964.
8.
KIVEL,
B.
Aero Sci.,
L. M.. VOROB’EV, V. S., NORMAN, heating in hypersonic streamlining.
Radiation from hot air and its effect 28, 2, 96-102, 1961.
G. E., YAKUNOV,
Kosmich.
issl.,
on stagnation
I, T.
2, 3. 441-454,
point heating.
J.