Numerical method of characteristics for three-dimensional supersonic flows

2 NUMERICAL METHOD OF CHARACTERISTICS FOR THREEDIMENSIONAL S U P E R S O N I C FLOWS* P. I. CHustmJNt Computing Centre of the U.S.S.R. Academy of Sciences, Moscow

Summary. The pt,esent paper lives a review of different schemes of the numerical method of ctmractefist~ for calculating three.41imemional steady supersonic I ~ flow about bodies moving at incidence. Some other methods for numerical solution of the supet'sonk flow problem are outlined. The schemes of method of ¢lmractedst~ pmpoNd for calculation of unsteady two-dimenaional flows are also immmted as far as they can he extended to the case under ~msideration. Four schemes of the method of characteristics are presented in detail, as developed recently in works by Soviet amhorJ for the calculation of three-dimensional flow about a body and checked in practical application,. Some results of numerical solutions are liven for illustration. This paper considers mahfly the flow of a perfect las, and the two last sect/ons deal with the llenerai/zation for the case of equilibrium and non-equ/libfium flow of real las.

1. N U M E R I C A L METHODS FOR SOLUTION OF THREE-DIMENSIONAL PROBLEMS IN GAS DYNAMICS

The study of supersonic three-dimensional gas flow about bodies at incidence is of great practical interest, and it is one of the actual problems in modem gas dynamics. The investigation involves different methods: experimental work, finding analytical solutions for equations governing this problem and numerical calculations on electronic digital computers. Experimental investigations of gas flows in supersonic wind tunnels imply considerable difficulties, especially at hypersonic velocities. Obtaining flows with high supersonic Mach numbers and with high stagnation t Enslilh trmmiation ed/ted by S. Neumark, R.A.E., Famboroulh.

t.Chta of Gas ~ m l ¢ .

~nment. 41

42

P. I. CHUSHKIN

temperatures requires expensive experimental installations. In some cases we are unable to simulate the flow in wind tunnels at all. Besides, due to complexity of measurements, we do not obtain good accuracy in these experiments and, as a rule, we get only the main parameters of the flow, rather than those of the flow fields which are of great interest. Mathematical calculation of three-dimensional steady supersonic flow about bodies consists in solving the boundary value problem for the system of non-linear partial differential equations depending on three variables. Due to complexity of this system, it does not seem possible to find accurate analytical solutions. All available analytical solutions are based on various simplifying assumptions, for instance, that of linearization or that of potential flow. These assumptions become nonsensical in the presence of large perturbations and at hypersonic velocities. In classical analytical methods it is very difficult to take into account physical and chemical reactions which take place in gas at high temperatures. In the general case the complete equations of gas dynamics can be solved only with the aid of numerical methods based on the use of high-speed electronic computers. Due to the development of such numerical methods, great progress has been made recently in the solution of many important problems in theoretical gas dynamics. It should be noted that the numerical methods enable us to determine the flow fields in the entire region under consideration. At present, there exist several effective and universal methods for numerical solution of gas dynamic problems depending on two variables, such as the method of finite differences, the method of integral relations, and the method of characteristics. By means of these numerical methods many calculations have been carried out for steady plane, axisymmetric, conical and unsteady one-dimensional flows of perfect and real gases, with such phenomena as viscosity, heat conductivity, radiation, and magnetic effects, taken into account. A review of applications of the above-mentioned numerical methods to gas dynamic problems with two independent variables was given in ref. 1. The change-over from the two-variable problems to those with a larger number of independent variables increases considerably the computational difficulties. First of all, the volume of computations increases sharply and the algorithms of calculation become more complicated; this leads to certain I/mitations in carrying out the numerical solutions by computers. The modern high-speed electronic computers enable us to calculate efficiently three-d/mensional steady and two-d/mensional unsteady gas flows. However, at present, the techrucal efficiency of computers appears

~-DIMENSIONAL

SUPERSONIC FLOWS

43

to be insut~cient for numerical solutions, with proper accuracy, of equations with four independent variables, such as describe three-dimensional unsteady gas flows. In numerical solution of multi-dimensional problems of gas dynamics there arise additional difficulties connected with the occurrence in such problems of more complicated physical features of the stream. For example, in a supersonic jet discharged from a three-dimensional Laval nozzle, or in supersonic flow about bodies at incidence, secondary shock waves may appear inside the flow region, and the calculation becomes more difficult because of their essentially three-dimensional character. In the three-dimensional hypersonic flow about blunt bodies the structure of the entropy layer (where the rapid change of functions takes place) becomes more complicated in comparison with the two-dimensional one, due to the presence of the cross flow. This may require the application of more exact numerical schemes or, in any case, a considerable refinement of the computational network. In this review of the three-dimensional method of characteristics for the calculation of supersonic steady flow about the bodies at incidence, we shall first consider briefly some other three-dimensional numerical methods. Only such methods will be discussed which have been applied to numerical computation of flow about bodies. However, attention will be given also to the numerical schemes proposed for calculating unsteady two-dimensional gas flows. We consider only essentially three-dimensional flows, and therefore do not intend to treat numerical solution of such problems as supersonic flow at incidence about plane and conical bodies or about triangular wings of zero thickness. Nor will this survey discuss works which investigate mathematical aspects of numerical methods and which are not applicable directly to gas dynamics. For numerical integration of partial differential equations depending on three variables, the method of finite differences is successfully used, in which the determination of functions at mesh points of a fixed network is reduced to the solution of a high-order set of algebraic equations. A method of finite differences of conventional type, for calculating three-dimensional supersonic gas flow about smooth bodies, has been developed by Babenko and Voskresenskii. <2~ The solution is carried out following layers perpendicular to the body axis. The non-linear set of finite difference equations is solved for each layer with the aid of iterations. The boundary conditions on the body surface and on the shock wave for each layer are satisfied by the factorization technique. The development of this method is given in the book <-~ which contains tables for supersonic flow

44

P . I . CHUSHKIN

past cones at incidence. By using this method, D'yakonov <4, s~ has calculated the supersonic flow region for blunt bodies at incidence, in the stream of perfect gas or of gas in thermodynamic equilibrium. In the general case of unsteady motions of gas, a number of shock waves appear in the flow field. As a consequence the direct use of finite difference schemes in this case becomes difficult. In calculating such flows some artificial dissipative terms are introduced into the original system of equations. This is achieved either by adding artificial viscosity (see, for instance, the paper by von Ncumann and Richtmyer ~6>) or by constructing a special finite difference scheme (see, for example, the papers by Lax (v) and Godunov(S)). In both these approaches spreading out of the shock waves into certain narrow regions takes place which makes possible continuous calculation without special treatment of discontinuities. Godunov, Zabrodin and Prokopo v<') have developed the fin/to difference scheme in which difference formulae are obtained by appi/cation of conservation laws to each cell of the network. They used this scheme for calculating unsteady axisymmetric supersonic flow about a sphere with a detached shock wave. Here the solution of the corresponding steady problem is found at the fimit for time increasing indefinitely. The Lax finite difference scheme, proposed initially for one=dimensional unsteady problems, has been extended by Lax and Wendroff (2s) and recently by Bohachevsky and Rubin (x°) to the case of unsteady streams with two space variables, and by Bohachevsky and Mates
THRI~=DIMEN~IONAL SUPERSONIC F L O W S

45

Another interesting approach, called particle-in-cell method, was developed by Harlow (~) and his colleagues for numerical calculation of unsteady gas flow. According to this method, the flow region is divided into cells by a certain network, and each cell is filled with a definite number of particles. The calculation is carried out by steps in time, in two stages. During the first stage the flow is calculated in the Eulerian presentation, and in the second stage the movement of particles is considered. Ewans and Harlow ~14~ have applied the particle-in-cell method to the calculation of unsteady supersonic flow past a circular cylinder. The disadvantage of this method is its limited accuracy which can be improved only by using high-speed computers with large storage capacity. A wide use has been made of numerical methods in which the gas dynamic equations are approximated with respect to one or two variables, and then the solution is reduced to integrating a system of ordinary differential equations. One method of such a type is the numerical method of integral relations proposed by Dorodnitsyn.c~ The idea of applying the method of integral relations to problems of three-dimensional flow has been put forward by Belotserkovskii and Chnshkin.~le) Minailos cv~ has used this method for calculating conditions in the three-dimensional region in front of the blunt-nose part of a body, in supersonic gas flow. Another efficientscheme of the method of integral relations for solving the same problem has been worked out in the monograph by Belotserkovskii et al.( ~ which presents some numerical results. The application of the method of integral relations for investigating threedimensional purely supersonic flow past bodies at incidence has been outlined by Katskova and Chushkin.c1.~ A numerical method which is in effect the modified method of straight lines, and which is also reduced to the integration of ordinary differential equations, has been worked out by Telenin and Tinyakov.¢~°.2v Using this method, they have calculated the three-dimensional flow in subsonic and transonic regions in front of smooth blunt bodies at incidence. This method was later extended by Telenin and Lipnitskiicm~ to the unsteady three-dimensional case, with unsteady perturbations assumed small. W e may note one more method proposed by Swigart ¢~ for numerical calculations of supersonic flow past a blunt-nosed body. This method, of inverse type (the body has to be determined corresponding to a given shock wave), is based on the construction of truncated double series inside the shock layer, the coeWcients of which series are determined successively by solving ordinary differentialequations.

P. 1. C H U S H K I N

The numerical methods of the above-discussed type are distinguished by their simplicity. They take into account the behaviour of the flow in the shock layer contiguous to the blunt nose and allow one to calculate this flow effectively. However, the extension of these methods to the purely supersonic three-dimensional domains of sufficiently great size may lead to difficulties due to the violation of the region of influence. Unfortunately, there exists as yet no experience in performing practical calculations in thl.q case. In investigating the purely supersonic three-d/mensional flows the particular subject of this survey) or unsteady multi-dimensional motions, the system of gas dynamic equations (the dissipative processes being neglected) will belong to the hyperbolic type. Thus, real characteristic surfaces will exist for this system. Consequently, it is purposeful to apply the numerical method of characteristics to solving such problems. The method of characteristics presents a number of advantages in comparison with other numerical methods. Firstly, the original equations are essentially simplified on the characteristic surfaces, and the number o f partial derivatives is reduced. Here the region of influence is determined rigorously, and thus the propagation of perturhations is adequately taken into account This circumstance gives a clear physical sense to the method of characteristics and enables one to consider rigorously some features of the flow, e.g. the rarefaction waves or the formation of shock in a wave of compression. For quasilinear hyperbolic equations with three independent variables, the existence and uniqueness of the solution of Cauchy problems by the method of characteristics have been proved. By using this method, it is easier to study the stability of numerical solution, and the stability criteria themselves become simpler. Finally, in the method of characteristics the particle paths or stream lines are always constructed, along which physico-chemical processes take place in a gas. This simplifies the integration of equations describing these processes. At the same time, the method of characteristics in practical application has some disadvantages. These are the complex configuration of characteristic surfaces, the large volume of computational work, the comparatively sophisticated and laborious computational algorithm, the difficulties in calculation of discontinuities inside the flow field. In this connection, it is reasonable to apply the method of characteristics to such hyperbolic problems of gas dynamics, in which one has to deal only with a small number of shock waves. The calculation of steady supersonic gas flow about bodies belongs to just such problems.

THREE-DIMENSIONAL SUPERSONIC FLOWS

47

The method of characteristics permits the construction of the characteristic networks in various ways. The characteristic surfaces can be issued forward (downstream) and the mesh points of the characteristic network can be found in the course of the numerical solution as the points of intersection of the characteristic surfaces. Using a network of such a type it is possible to account rigorously for the propagation of weak discontinuities. In this case, however, the mesh points of the characteristic network are placed in the flow field very irregularly, and this creates some inconveniences in practical calculations. Particularly serious di~culties of such kind arise in multi-dlmensional problems. The above considerations make it preferable to use a characteristic network of the so-called inverse type. Here the calculation proceeds layerwise. Each layer is represented as a coordinate surface on which mesh points are placed at locations fixed in advance. The characteristic surfaces are projected backward (upstream) from the mesh points of the layer to be calculated in the direction to the previous one, where the flow parameters at the corresponding points of intersection are found by interpolations. The computational network of inverse type was proposed first by Hartree ~ for unsteady one-dimensional flows. Such a computational network presents a number of advantages. It can be easily modified in the course of solution. The coordinates of mesh points are known in advance. This reduces the number of data to be stored in the memory of an electronic computer and yields the results at fixed points on the coordinate planes, which is convenient in practice. A disadvantage of the network of inverse type consists in the necessity of carrying out interpolations which require more computer time and may give rise to additional computational errors. It may be noted that, when using the network of inverse type, after the completion of calculations in each layer, it is possible to construct successively the full characteristic surfaces. By choosing this or that characteristic network, one can construct calculation schemes both of the first and second order of accuracy with respect to the step size. In the latter case, the iteration process should be introduced in determining the flow parameters at the point under consideration and all coefficients of finite difference equations should be averaged according to the data at points already known and at those now calculated. The use of schemes of the second order of accuracy saves computer time because larger steps may be chosen, but it is associated with a more complex computational algorithm. However, the schemes of the first-order accuracy are often preferred for the sake of simplicity.

48

P. I. CHUSHKIN

2. D E V E L O P M E N T METHOD

OF THE THREE-DIMENSIONAL OF CHARACTERISTICS

The numerical method of characteristics has been widely used for solving various gas dynamic problems dependent on two variables. Calculation schemes have been developed for supersonic gas flows with arbitrary thermodynamic properties (perfect gas, real gas with equilibrium or nonequilibrium processes). Such schemes have been constructed in the form specially adapted to electronic computers. The present state of the twodimensional method of characteristics is described in the surveycl~ which has been already mentioned. Though the mathematical theory of quasi-linear partial differential equations of hyperbolic type was developed long ago, the practical generalization of the method of characteristics for the case of more than two independent variables has encountered certain difficulties. This is why systematic computations with the aid of the three-dimensional method of characteristics have been really started only recently. The three-dimensional equations have some characteristic properties basically different in comparison with the case of two-dimensional equations. It is known that partial differential equations with two idependent variables, when written in the characteristic form, change into ordinary differential equations which can be integrated along the characteristics. In the case of three independent variables, however, the original equations written in the characteristic form remain partial differential equations and contain partial derivatives along two directions lying in the characteristic plane. Therefore, to perform numerical integration in this case, one should use essentially two-dimensional finite difference schemes. Furthermore, in the two-dimensional case, through each point of the hyperbolic region pass two characteristic lines and a stream line which also possesses the characteristic properties. In the three-dimensional case, however, there exists a one-parametric family of characteristic surfaces passing through each point, and their envelope is a conoid which determines the region of dependence and the region of influence of the point. The lines along which the characteristic surfaces touch this conoid are called bicharactcristics, and they also form a one-parametric family. As a result, when trying to solve numerically a system of threedimensional partial differential hyperbolic equations, it is possible to choose various types of computational schemes with various elementary cells, whereas in the two-dimensional case the calculation scheme is determined uniquely in a natural way. In fact, a number of authors have proposed

THItEEoDIMEI~IONAL S ~ N 1 C

FLOWS

49

many schemes of the method of characteristics both for the steady threedimensional and for unsteady two-dimensional flows. The number of computational schemes in the three-dimensional case is even greater due to the fact that there exists here a possibility to combine the two-dimensional method of characteristics with the finite difference method. Thus schemes of the so-called semi-characteristic type make an appearance. In these schemes, the original system of partial differential equations depending on three variables is reduced to a two-dimensional system of differential equations by means of approximations with respect to one of the variables. The resulting approximating two-dimensional system can be solved by the numerical method of characteristics developed for the equations with two independent variables. The semi-characteristic calculation schemes are distinguished by their greater simplicity in comparison with the full three-dimensional method of characteristics. However, the region of influence is not considered rigorously in semicharacteristic schemes, and therefore they can be naturally applied only for the calculating comparatively smooth flows. Let us now pass on to the description and brief analysis of actual numerical schemes of the three-dimensional method of characteristics, following mainly the chronological order. The first numerical schemes of the method of characteristics for threedimensional gas flows were proposed prior to the introduction of digital electronic computers in the late 1940's and early 1950's, by Ferrari, t26> Ferri, ts~> Coburn and Dolph, ts°> Sauer ts°~ and Sychev.t3°~ All these schemes belong to the semi-characteristic type and, in all of them, two families of characteristic surfaees are considered which intersect the third family of non-characteristic surfaces. All these authors limited themselves to the case of potential flow, and therefore the above papers are only of historical interesL A profound and wide investigation in the theory of characteristics of general gas dynamic equations was carried out by Rusanov ~3v in 1951 in his dissertation, and the results were published in the form of an article in a journal ~ considerably later. Rusanov studied the characteristic properties of systems of quasi-linear partial differential equations of general type. He derived the equation of the normal to the characteristic surfaces, and the compatibility relations. The question of inter-dependence between characteristic relations was also investigated. The results obtained were applied to two particular cases, viz. to the systems of equations governing the unsteady three-dimensional and the steady three-dimensional gas flows, respectively.

50

P . I . CHUSHKIN

For numerical calculations of the three-dimensional supersonic steady flows, Rusanov has developed in detail the characteristic network which is formed of elementary tetrahedrons and represents the natural generalization of the triangular characteristic network in two-dimensional case. The elementary tetrahedron is constructed here in such a way that the point n o w under consideration is its apex, while its three faces which are characteristic surfaces are tangential externally to the Mach cone issued in the Upstream direction from this point. The base of this tetrahedron is placed on the surface with the known data. A tetrahedral scheme of the above type was also proposed earlier by Thornhill (3a) and discussed later by Ferri. (a*) For the three-dimensional method of characteristics, Thornhill gave a further numerical scheme of the tetrahedral type, in which a new point is found as the intersection of three Mach cones issued in the downstream direction from the points with known data. In such a way, the bicharacteristics serve as the tetrahedron edges, and the tetrahedron itself will be inscribed inside the Mach cone issued in the upstream direction from the new point to be calculated. However, it became clear later on that this tetrahedral bicharacteristic scheme was unstable. Holt (m) has generalized the work by Coburn and Dolph (28) for the case of non-isentropic three-dimensional gas flows, and developed for them a prismatic numerical scheme belonging in effect to the semi-characteristic type. Here, the elementary cell is taken in the form of a prism whose two triangular bases are constructed on two pairs of bicharacteristics issued from computational mesh points, and one lateral rectangular face is placed on the surface with known data. In order to determine the parameters of the flow at a new computational point, Holt used not only relations along two bicharacteristics and the stream lines but also a supplementary finite difference equation along a non-characteristic line, viz. the prism edge. This numerical scheme, constructed in non-orthogonal coordinates, is very complicated for practical realization, and it is difficult to adapt it to the calculation o/" shock waves. In any case, it has not been used for actual computations of the three-dimensional flow about bodies. Bruhn and Haak (a6) considered unsteady two-dimensional and threedimensional gas flows. They presented the corresponding equations of motion at each point within the flow in the local system of coordinates formed by the trajectory of the particle, its principal normal and the binorreal. In such intrinsic coordinates the equation of the characteristic hypcrconoid and the relations along bicharacteristics were derived. The authors have developed a numerical algorithm for solution and used a computa-

THItJ~-DIMENn31ONAL SUPERSONIC FLOWS

51

tional network of inverse type with the successive time layers and with characteristic conoids issued in the upstream direction. Bruhn and Haak applied their numerical scheme of the method of characteristics to the calculation of non-stationary supersonic flow in a two-dimensional Laval nozzle. This scheme is quite adequate for the problems where shock waves do not appear, but it is not very suitable for flows with shock waves. For studying three-dimensional hyperbolic problems of gas dynamics, Butler c:7~ proposed a new numerical scheme of the method of characteristics. The equations, in the tetrahedral scheme, contain not only finite difference expressions along three bicharacteristics but include also the derivatives in non-characteristic directions which have to be calculated again by means of difference formulae. Buffer uses the network of inverse type and, in calculating an inner point in the flow field, issues from it four bicharacteristics and a particle line directed towards the preceding layer. Certain linear combinations of the finite difference analogues are then constructed for the compatibility relations along four bicharacteristics and for the transformed continuity equation along the particle line. This makes it possible to eliminate the derivatives along non-characteristic directions at the unknown point, which is undoubtedly advantageous. The scheme has some disadvantages, first of all the necessity of a great number of interpolations for determining functions at four intersection points of the bicharacteristics with the preceding layer; this increases the computer time and results in additional numerical errors. Further, the calculation of the boundary points becomes too cumbersome. Buffer has developed his scheme for unsteady two-dimensional problems. An example of hypersonic flow about a delta-shaped body which he presented was computed in effect as a non-stationary motion of an elliptical piston, and then the hypersonic analogy for slender bodies was used. Richardson (ss) has applied the given scheme to the numerical solution of the unsteady problems with axial symmetry. In 1961 Fowell'3'~ made a survey of the numerical method of characteristics for calculating three-dimensional steady supersonic gas flows about the bodies at incidence. He has analysed five different types of computational schemes as suggested by different authors (all of them discussed above), and estimated their advantages and disadvantages. Fowell gave preference to the Thornhill scheme in which the new point is found by the intersection of three Mach cones. For this scheme, considering as the basic unknown functions the pressure and two angles determining the direction of the velocity vector, Fowell derived all the necessary formulae for cal-

52

P . I . CHUSHKIN

culating mesh points of various types. However, as mentioned above, this numerical scheme proved to be unstable. Podladchikov ~ ~ used the Rusanov tetrahedral scheme, and was the first to work out this variant of three-dimensional method of characteristics on an electronic computer. He has brought in a number of improvements and developed a practical computational procedure for solving elementary problems in the method of characteristics for the cases where the new point to be computed is placed either on the shock wave, or inside the flow field, or on the body surface. The tetrahedrons were issued in the downstream direction forming a characteristic network analogous to the conventional network in the two-dimensional flows. However, for the multidimensional case, the manner of constructing a characteristic network as a whole is essentially significant. Podladchikov presented a convenient scheme which determined the sequence of including the mesh points into further computations. He considered various examples of supersonic flows past bodies at incidence. One more numerical scheme of tetrahedral type, for studying supersonic three-dimensional steady flow about sufficiently smooth bodies, has been given by Minostsev. ~ ' ¢> He applied the characteristic network of inverse type in which the computations are carried out following consecutive layers corresponding to constant values of one of the independent variables, and the Mach cones are issued in the upstream direction from the points fixed in advance on these layers. The unknown gas dynamic functions on each layer are represented by interpolational expressions in two other independent variables. Having in view calculations of smooth flows, the continuous approximations are used. The pertinent approximating formulae give all the information about the field of flow and enable one to determine both the gas dynamic functions themselves and their derivatives. The computational algorithm of this scheme is comparatively simple, and the practical calculations showed that, for very smooth bodies, it is sufficient to take the interpolation formulae of not too high order. As a rule, in computational schemes of the three-dimensional method of characteristics, the compatibility relations are used in such a form that they involve derivatives along two directions. For solving the three-dimensional steady problems of gas dynamics, Magomedov ~"'~> developed a numerical scheme in which only derivatives along bicharacteristic directions remain. This simplifies considerably the computational process and gives it a high logical orderliness. In the paper, ~"~ such a scheme was constructed for perfect gas flows, and later on, ~u~ this scheme was somewhat

THREE-DIMIENSIONAL SUPERSONIC FLOWS

53

modified and generalized for equilibrium gas flows with an arbitrary equation of state. These papers contain practical examples of the supersonic flows about different bodies at incidence. In the above-mentioned scheme, for calculating the point placed on a shock wave, the compatibility relation along one bicharacteristic and the geometrical condition written in the finite difference form are taken, in addition to the laws of conservation. Magomedov <~s~proved that, in the three-dimensional numerical method of characteristics, for determining the flow parameters at points located on a shock wave or on a free surface, it is necessary and sufficient to take only one combination of gas dynamic equations, e.g. one compatibility relation along a certain bicharacteristic. In a number of schemes of tetrahedral type aforementioned<3s,,0. ,,~ it was proposed to calculate points on the shock by using compatibility relations along two bicharacteristics, and thi.~ could lead to errors. Indeed, in the Minostsev ( ~ scheme for this case, an increase of errors in the values of functions on the shock wave was observed, especially in the values of the peripheral velocity component which exactly characterizes the threedimensional nature of the flow and is a very sensitive parameter. When the procedure was replaced by one utilizing only one bicharacteristic, the errors of such a type were eliminated. Sauerweinc4~) generalized the numerical method of characteristics for unsteady three-dimensional flows of chemically reacting multi-component media possessing perfect electrical conductivity. In order to make the system of partial differential equations governing the flow of such a type be hyperbolic, all terms connected with dissipative phenomena have to be rejected; in particular, in magneto-gas-dynamic problems, the effects of electrical resistance must be neglected. In the above paper, the characteristic equations and compatibility relations are derived for the case of unsteady three-dimensional magneto-gasdynamic flow with multi-component non-equilibrium thermodynamics. The characteristic hypercones in this case have a very complex appearance and represent the combination of slow, fast and Alfv~n waves. As a result, it is necessary, in numerical computations, to apply the inverse-type characteristic network and carry out layer-wise computations. In order to take into account all the perturbations in the region of dependence in the magneto-gas-dynamic problems, it is necessary to take not fewer than thirteen points on the base o f a hypercone in the hyperplano with the known data for calculating unsteady two-dimensional flows, whereas for the unsteady threo-dimensional flows at least sixteon points must be taken. It turns out, however, that when using such a number of points, the number of avail-

54

1'. I. CHUSHKIN

able compatibility relations will exceed the number of the unknown functions, and therefore these relations have to be solved by the method of least squares. The practical application of the method of characteristics to numerical solutions, even in the case of two-dimensional magneto-gasdynamic problems, produces serious difficulties, and Sauerwein confined himself to calculating a rather simpler example of a supersonic flow of perfect gas past a circular cylinder whose surface changes its form in time. As mentioned above, considerable computational simplifications may be obtained by using semi-characteristic schemes in which we encounter equations with only two independent variables. Within the last few years, two efficient numerical schemes of such type have been proposed for steady three-dimensional supersonic flows of gas in the state of thermodynamic equilibrium. From the set of three-dimensional gas dynamic equations Moretti <4s) excluded derivatives with respect to the third independent variable (which was chosen in accordance with the actual shape of the body), by introducing a number of interpolation planes and by approximating derivatives by means of the difference of the corresponding functions on two adjacent interpolation planes. Although the layer-wise computational network was used here, the straight-line segments of two-dimensional characteristics in each interpolation plane were constructed in the downstream direction towards a new layer, and the iterations (ensuring the second order of accuracy) were not carried out in this work. c4s~According to this scheme, the supersonic flow past blunt cones and delta wings was calculated; some results of these computations were given also in the paper by Moretti and Byrne. <4'~ Recently, in a short note, Moretti ~5°) gave information about some improvements of this scheme. Katskova and Chushkin ¢51~have worked out a semi-characteristic scheme which has some definite advantages. With the view of studying supersonic three-dimensional flow about bodies of revolution or similar at incidence they chose the cylindrical system of coordinates. With the aid of trigonometric approximations in the angular variable, the corresponding cross derivatives were excluded. The solution is constructed along layers perpendicular to the body axis, the characteristics are issued into upstream direction from points fixed in advance, and the iterations process is used in order to improve the values of the functions being sought. Such a computational scheme is of the second order of accuracy and enables one to apply continuous trigonometric approximations; this shortens computing time considerably, as a smaller number of layers and interpolation planes is

THREE-DIMENSIONAL SUPERSONIC FLOWS

55

required. In addition, a new variable is introduced which straightens the shock wave and the body, and this simplifies the determination of parameters at the boundary points. Katskova and Chushkin cx°' ~) applied their semi-characteristic scheme to investigating flows past both blunt-nosed bodies and bodies with inlets. In carrying out computations based on any numerical method, an essential role is played by the problem of stability of solution. In the three-dimensional method of characteristics the stability condition depends on the type of the computational scheme, which may be simplicial (when in computing a new mesh point a minimal number of points with the known data takes part) or non-simplicial. Inasmuch as the investigation of stability is made only for linear equations, the stability criteria are only "local" and, to a certain degree, of heuristic nature. The stability studies of three-dimensional schemes of the numerical method of characteristics were carried out by Sauerwein and Sussman~> and Heie and Leigh ~u~ who compared the theoretical analysis with the results of numerical calculations. Sauerwein and Sussman considered the stability of simplicial schemes. The necessary and sui~cient condition of convergence and consequently of stability for the finite difference schemes of such type is the CourantFriedrichs-Lewy criterion. As known, this criterion requires that the region of dependence of the finite difference scheme should include the region of dependence of the original partial differential equation. Sauerwein and Sussman pointed out that this criterion is not satisfied by the Thornhill-Fowell characteristic scheme which thus turns out to be unstable. Heie and Leigh investigated the stability of non-simplicial schemes and showed that, in this case, it is necessary to apply the more severe Nenmann's condition which takes into account not only the position of the initial points but also the method of inserting them into the finite difference scheme. The Neumann's criterion requires that the eigenvalues of the matrix of the set of finite difference equations should not be, in absolute magnitude, greater than unity. The non-simplicial schemes will be stable only in the case when they satisfy not only the Courant-Friedrichs-Lewy criterion but also the Neumann's criterion. It is known that, in some non-simplicial schemes of the numerical method of characteristics, the necessary Courant-Friedrichs-Lewy condition is not satisfied; nevertheless, in practical calculations, the numerical instability was not observed there. This is explained by the fact that using in these schemes the interpolations on the surface with the known data may enlarge the region of dependence of the finite difference scheme to such

56

p. I. CHUSHKIN

a degree that it becomes stable. In this case the rigorou anaslysis of stability depends on particular form of such interpolations. To conclude the present survey, we may note the papers by Burnat c55~ and Burnat, Kielhasinski, Wakulicz,~Se) where the three-dimensional method of characteristics is applied to the study of gas flows at supersonic velocities. They developed a numerical scheme of the tetrahedral type and constructed computational algorithms for the points inside the flow field, on the body and on the shock wave (strong or weak). However, the given simplest example (axisymmetrical flow past a cone) is of a preliminary nature and unfortunately does not permit to make a conclusion as to the effectiveness and accuracy of this scheme in the general three-dimensional case. In a recent note by Sauer, ~57~ a scheme of characteristics for threedimensional potential flows has been suggested, in which the solution is determined on the planes x = const. The author proposes to take two pairs of characteristic planes parallel to y- and z-axes respectively. However, this scheme has not been checked by any calculations of flow about bodies. Heie and Leigh~u~ have also analysed the stability of a non-simplicial scheme proposed by Brong. ~Ss~This scheme uses a pentahedral elementary cell whose faces are characteristic surfaces passing through the rectangular base where the initial data are known. Here, the computational algorithm involves interpolations. Since no numerical results are available, it is difficult to appraise the practical usefulness of this scheme; it is stable but seems too complicated. Sauerweinc¢) says that Strom cS°~ and Tsung ~°°~ carried out calculations with the aid of the three-dimensional method of characteristics in their dissertations. It is a pity that these two papers, as well as the paper by Brong, ~ss~ have not been available to the author of the present survey. In the following sections, we shall discuss in detail four numerical schemes of the method of characteristics which have been developed lately in the Soviet Union and which are being successfully used for systematic calculations of the three-dimensional steady supersonic flow past bodies at incidence. We shall consider in turn: the Rusanov-Podladchikov direct tetrahedral scheme, the Minostsev inverse tetrahedral scheme with continous interpolations, the Magomedov scheme using only bicharacteristic directions, and the Katskova-Chushkin semi-characteristic scheme. Special attention paid to these schemes in the given survey is justified by the fact that they are effective, have been quite recently realized in practice, and are obviously less known to the foreign reader. We shall consider these schemes for the case of flows of perfect gas, and then indicate a generaliza-

THREE-DIMENSIONAL S ~ N I C

FLOWS

57

lion for the case of equilibrium and non-equilibrium flows of real gas. Each of these schemes will be illustrated by certain numerical results. To begin with, we shall formulate the problem of three-dimensional supersonic gas flow and give the basic equations of the method of characteristics for this case. 3. F O R M U L A T I O N OF THE PROBLEM AND BASIC EQUATIONS OF THE THREE-DIMENSIONAL METHOD OF CHARACTERISTICS

Let us consider a supersonic steady iso energetic flow of an inviscid non-heat-conducting perfect gas, with constant specific heat ratio y, past a given body (Fig. 1). The three-dimensional nature of such a flow may be conditioned either by the direction of the flow whose velocity vector V . forms an angle of incidence 0c with the body axis, or by the body shape itself (if it is not a body of revolution or a two dimensional body). For calculating supersonic flow past a body, one should determine a bow shock wave and find the fields of the basic gas dynamic functions in the region between thi.q shock wave and the body surface. For simplicity, we assume that the body has a symmetry plane parallel to the velocity vector V . , in particular it may be a body of revolution. Furthermore, let us confine ourselves to such cases where other shock waves do not exist inside the flow region. The numerical solution of the flow problem under consideration is normal]y carried in two separate stages. First of all, the flow about the nose part of a body is calculated. As a result, we should obtain the necessary initial data in the flow region on a certain space-like surface, i.e. on such a surface which does not intersect at any point the corresponding characteristic cones. Further, for calculating the purely supersonic zone, it is necessary to solve the Cauchy problem with the initial data already known. The determination of the flow parameters near the nose of the body depends on its form. For blunted bodies, various numerical methods have been developed which make it possible to obtain the necessary initial data. If the blunt nose has a spherical shape, then such data can be determined (up to certain values of the incidence) by rotating the relevant numerical solution for the flow past the sphere. In the case of a conic nose, it is replaced by a certain small cone, and the corresponding flow parameters can also be calculated numerically or taken from the existing tables.~S, el, 62) Finally, for a body with a duct, in the case of an attached shock wave, the initial data on the leading edge of the body are obtained

58

P . I . CHUSHKIN

by considering the flow past a corresponding wedge at certain effective values of incidence and Mach number of the incident flow (see re/'. 78). For calculating steady supersonic flow in a certain region, one must integrate the system of equations of gas dynamics which, for an inviscid non-heat-conducting gas, has the form ~7~V = O, ~(V. v) V + v p = O, V . v S = O,

(3.1)

where V, p, O, S denote the velocity, pressure, density and entropy function, respectively. In the case of a perfect gas with constant specific heat ratio y, it is possible to express the density 0, enthalpy h and speed of sound a in terms of pressure and entropy by the following formulae:

O=Ptlvexp['l-Y 1, y

a2 =

S~ /'

h-

Y y-1

P o'

/ (3.2)

r , - -p. "

Q All quantities in the above equations will be considered as dimensionless, assuming as reference quantities a certain linear dimension (e.g. the radius of curvature of the blunt nose), the free stream velocity V. and density Q.. Then, in particular, the pressure will be referred to the quantity o, VS,. and the enthalpy to the quantity V_t. If, instead of V.., another reference velocity is sometimes used, this will be specially mentioned. The solution of the present problem can be carried out in various systems of coordinates. For instance, one may choose the Cartesian coordinates x, y, z connected with the free stream, directing the axis x along the velocity V . (Fig. 1). However, in such a system, it is difficult to represent the data on the body surface when changing the incidence. For bodies of

at

Flo. 1. Scheme of flow p u t a body in Cartesian coordinates connected with um/isturbed stream

THRI~-D~IONAL

°~'f~

~

59

SUPERSONIC FLOWS

'

-¢=cons,

FIo. 2. Scheane of flow past a body in cylindrical coordinates connected with the body

revolution, it is convenient to take the cylindrical system of coordinates x, r, 1Oconnected with the body (Fig. 2). The axis x is directed along the body axis, and it is assumed that I0 = 0 corresponds to the windward side and 1O-- n to the leeward side. Here, using inverse-type schemes, it is reasonable to straighten the region of flow to be calculated between the shock wave r = r,,(x, 1O) and the b o d y r = r b (x, ~), by introducing the normalized variable _ r - r b ( x , ~P) ,

(3.3)

e

where

= rw(x, lp)- rb(x, 1O).

(3.4)

Obviously ~ = 0 corresponds to the body, and ~ = 1 to the shock wave. Then, on each plane x = coast, it is possible to construct the network whose mesh points are formed by intersection of the surfaces ~ = const and meridional planes ~ = const. Let us note that, when using the variable ~ for obtaining the solution on a new plane x -- const, the calculation must be started by determining the coordinate r,, of a new point on the shock wave which is contained in the formula for ~. On the body surface, we introduce the condition that the normal velocity component vanishes,

v . n = o,

(3.5)

where m is the unit normal to the body surface. The following known relations expressing the conservation laws should be satisfied on the shock wave p+eV,,Zf p.+e.vL,

h+~

=h.+ V, = V , . .

V~,

(3.6)

60

P. I. CHUSHKIN

Here, the indices n and • refer to the normal and the tangent to the shock wave, and the index ~o to the undisturbed flow. The quantifies without the index ~ are the values of functions immediately behind the shock wave. If one considers, for instance, the Cartesian system of coordinates x, y, z, then it is possible to express, from (3.6), the velocity components V~, Vy, V: along the corresponding axes immediately behind the shock wave as follows cos (a, x) vx = v ~ . . - A p cos (n, V®) ' cos (n, y) cos (n, V . ) '

Vy = V , , - - A p

(3.7)

cos (u, z) V~ = V , . - A p cos (n, V . ) ' and also the density o, pressure p; and from one of the formulae (3.2) entropy S,

A.v~ 2

¢) = l - ( y -

2

l) AL(1 - V~..)/(y+ 1)'

p=p.+~p, Ap---2 (V~.. 1 ), y+l

_ P"

M~

1

S _ _ _ _ _ ~ l in p

~,MS.. '

y- 1

(3.8)

Ov

Here, M.. = V . . / a . is the Mach number of undisturbed flow and A.. = V../acr its velocity referred to the critical speed of sound aer, where AL = (y + 1) M~ 2 + ( y - 1) ML If the shock wave equation is taken in the form x = xw(y, z), then the direction cosines of a unit normal, i.e. its components along the coordinate axes, will be r

/~Xwx z

/~Xw\Z~-tlZ

n x = -- [ 1 + ( - - ~ - ) + ( - - ~ ) 0xw ny = - - n , c - ~ ,

J

0xw n : . = --nx ~z "

, (3.9)

It is obvious that (using the dimensionless variables under consideration) the velocity component of the incident flow, normal to the shock wave, will be Vn. = V_ .n = cos (n,V_). (3.10)

TIIREE-DIMIENSlONALSUPEItSONICFLOWS

61

Thus, under the given conditions of the flow, the gas dynamic functions immediately behind the shock wave are expressed by two derivatives, to be determined from the equation of its surface. If the shock wave equation is taken in the cylindrical coordinates • = r~,(x, ¥), then analogous formulae will apply. In particular v.= =

(3.11)

where tD _- ~O ~

@, rLl+

cos ~,+ sin g cos ~0+ ~ 1 mOr, sin ~ sin ¥, r, O~

ox!

/ "

The velocity components u, ~, w, along the corresponding directions x, r, V' connected with the body, will then be u =cos~

Ap 0rw t~) 0x '

v ---- --sin= c o s ¥ + ~ ,

(3.12)

A o 0¥w w = sin = sin ~0--~r"" t and the formulae (3.8) for ~, p and S will naturally remain unchanged. In connection with solving the set of gas dynamic equations (3.1) by the method of chara~,ristics, let us recall some concepts of the theory of s y s t e m s of partial quasilinear differential equations. The fundamental role in this theory is played by the concepts of the character~tic hypersurface and characteristic hypercone which, in the case of three independent variables, become simply the characteristic surface and cone. As known, the characteristic surfaces are defined by the property that the Cauchy problem has no unique solution on them. The characteristic conoid is an envelope of the family of characteristic surfaces passing through the given point. The lines along which the characteristic surfaces touch the characteristic conoid are called bicharacteristics. At its apex, coinciding with the point under consideration, the characteristic conoid is tangent to the corresponding characteristic cone. The normals to the characteristic surfaces at this point form the cone of normals. On the characte~tic surfaces, the original system of quasilinear equations can be replaced at each point by a certain number of independent characteristic relations which are also called the compatibility relations.

62

1,. I. CHUSHKIN

For the hyperbolic system (3.1) governing the steady supersonic gas flows, the equation of the cone of normals splits into two equations V-n=a,

n2 = 1

(3.13)

and V . n = 0.

(3.14)

The equations (3.13) determine the cone of normals to characteristic (wave) planes and show that thi.~ cone of normals (real only for V > a) will be a cone of revolution with the axis directed along the velocity vector V and with the semi-apex angle v = arc cos (a/V). The direction of bicharacter~tics in this case is defined by the vector I, which will be 1 = V-an. The characteristic cone here will be also a cone of revolution with the axis directed along the vector V and the semi-apex angle/~ = arc sin (a/V). The equation (3.14) determines the cone of normals to the stream planes passing through the velocity vector V and shows that this cone degenerates into a plane perpendicular to V. The corresponding characteristic cone then degenerates into a stream line which is tangent to the velocity vectorV. Rusanov c~') derived the compatibility relation for characteristic (wave) planes in the following form, convenient for numerical computations

eaA,. a , , v - s , a , a , =

a v-s a,,p,

(3.15)

1,2).

(3.16)

where Aj=V×sj,

Bj=n.Aj

(j=

Here sl and ~ are two independent vectors in the characteristic plane with the external normal n. The differentials of pressure p and velocity V along sx and ~ are expressed in a general form for a certain function ~ ( x , y, z) in the system of coordinates x, y, z, as follows: d*' (rJ: = ~--~x s ' x + - - ~ - s J" + - - ~ - s ' ~

(3.17)

where six, s~, sj, are projections of the vector sj on the corresponding coordinate axes. The pressure p in compatibility relation (3.15) may be expremed in t ~ m s of entropy S, enthalpy h and temperature T which are connected by the thermodynamic equality T d S = dh - d p .

Introducing the full stagnation enthalpy Ve h,t - T + h ,

THREE-DIMENSiONAL SUPERSONIC FLOWS

63

we obtain Gt . d.t V + Bt(Td., S + d.thtt) = Gx . d.,V + Bx(Td.,S + d.~htt),

where c j = aXj+BjV.

For isoenergetic flow h0t = const and then, in the case of a perfect gas, the compatibility relation on characteristic planes will assume the form (in our dimensionless variables) 7 ~ i " d~V + atB~l~ S = ~ 1 " d~ V + a2Btd.,S.

(3.18)

In the stream plane, the following relations will take place d v S = 0,

~V- dvV + dvp = 0.

(3.19) (3.20)

The last relation represents the well-known Bernoulli equation (which holds along the stream lines) V2

2

+ h = const.

(3.21)

The Bernoulli equation can be rewritten in the following form: +h =

v'.L=

(3.22)

2 '

where V=.. is the maximal adiabatic velocity which, for a perfect gas, is connected with the critical sound velocity act by the known formula

v:.,=

q

ac, •

Rnsanov
64

P . I . CHUSHKIN

basic elementary problems of the method of characteristics: calculation for a po/nt inside the flow field, for a point belonging to the body surface and for a point located on the shock wave. Apart from this, some more problems may be encountered here: the calculation of a rarefaction wave or shock wave arisen at a sharp corner on the body surface, and the problem of formation of a secondary shock in a compression wave. When investigating the jet flows, one has to deal also with the problem of calculation for a point on a free surface. The solution of elementary problems of the method of characteristics consists, generally speaking, of two interconnected procedures. Firstly, we must find the coordinates of a new computational point (in the case of direct-type characteristic network) or the coordinates of the initial points with known data (in case of inverse-type characteristic network). Then the gas dynamic functions at the new point are calculated with the aid of characteristic relations and corresponding boundary conditions (in the case of boundary points). To improve the accuracy of the values obtained, it is possible to introduce an iterative computational process. We shall now turn to the specific description of four separate schemes of the numerical method of characteristics, for calculating threedimensional steady supersonic gas flows.

4. D I R E C T

TETRAHEDRAL

SCHEME

Let us consider first the direct tetrahedral scheme of the method of characteristics developed by Podladchikov <4°' tl) on the basis of Rusanov's investigation. ~sl~ We shall describe this scheme taking into account those modifications which have been introduced recently by its author. The solution of the problem is carried out here/n the Cartesian coordinate system x, y, z connected with the free stream (Fig. 1). As an elementary computational cell for internal points of the flow region a tetrahedron is chosen which is circumscribed to the Mach cone issued in the upstream direction from the point under consideration 0 (Fig. 3(a)). Thus the characteristic planes, passing through three vectors tl, t2, ts linking the known points 1, 2, 3, serve as the tetrahedron faces. The coordinates of the tetrahedron apex 0 are determined following the normals ~ , ms, is to its corresponding faces. The normals to the characteristic planes are determined by the following formula n, =

t~ X II~ A~

(i = 1, 2, 3),

65

THlCJEE-DIMIEN~ONAJ[, S U P E l t ~ N I C F L O W 8

where, according to (3.16) At = V t x t t , and gi = alAt+ biVt. The value of bt is given by the expression

bl =

-I- .~/'(A~-a~t~).

0

t,

2

2

tz

'z ~-.o~, (o)

5

(b)

3

S

y / o6

(c)

Fro. 3. Calculation of typical points in direct tetrahedral scheme Since it is possible to draw two characteristic planes through the vector t I, so two values for b t are obtained, differing by their signs. If the surface on which the vectors t~ fie is the space-type surface and if, going round in direction of tt, corresponds to positive going round the triangular base of tetrahedron, then the value of b: with plus sign should be taken. As to the quantity B t in thl.q case, it is not difficult to get from (3.16) that B t -- - b : . The normals a t to the faces of tetrahedron having been determined, it is possible to find its edges It, and hence to calculate the coordinates of the new point 0. For the vectors It, we have the following formulae na'tx it ---- ns-(nsXux)(~X~),

Is = Ix-tx,

is = l l - t l .

66

P . I . CHUSHKIN

For calculating gas dynamic functions at the apex of the tetrahedron, we use the compatibility relation (3.18) on each face of the tetrahedron, in which the differentials are replaced by the corresponding differences of values of functions at the new point 0 and at the point already known. In the compatibility relation (3.18), the vector t, is taken as the vector s~, and the role of the vector sz is played by the vector directed along the bicharacteristic from the centre of the vector tt to the tetrahedron apex 0. Furthermore, the condition of constant entropy S (3.19) along the stream line must also be used. It is assumed that approximately the stream line is straight and coincides with the axis of the Mach cone inscribed into the tetrahedron. Let us denote by the digit 4 the point of intersection of a stream line with the tetrahedron base. The value of entropy at the point 4 can be found with the aid of linear interpolation between its known values at points 1, 2, 3. We shall then obtain $4 = 2zSz+2zS2+~3S3,

where v ,,

it

,

~3

v

= Ih, Xl, l+ll.,,Xl l+ll Xl ,l.

We have therefore S0 = S~,

Taking the last equafity in account, the three compatibility relations given above form a system of equations which enables one to determine three velocity components along coordinate axes at the point 0. Let us consider now the calculation for a point on the body surface, referring to Fig. 3(b). Let the point 3 be the point on the body surface already calculated and 0 the point to be determined. Let us draw the characteristic (wave) plane 1-2-0 through the vector tz and the stream plane 4-3-0 through the vector t4 = ~(ts-tz). The point where the intersection line 4-0 of these two planes meets the given body surface will be the new point. For determining flow parameters at the point 0 just found let us put the condition of constant entropy on the body surface, and take the boundary condition (3.5), the compatibility relation (3.18) on the characteristic plane 1-2-0, and the characteristic relation in the stream plane tangent to the body and passing through the stream line 3-0. The last relation on the stream plane can be presented, in accordance with the work ¢~ and with regard to the condition S = const, in the form s l . a . , v = s=.a

v,

SUPERSONIC FLOWS

~-DIMENSlONAL

67

where the vector sl is directed along the stream fine 3-0 and the vector s2 is a certain other vector located also in the plane tangent to the body. The four above equations will sui~ce for calculating gas dynamic functions at the new point on the body surface. For calculating a point placed on the shock wave, Podladchikov, with regard to the investigation by Magomedov <~>, proposes the following scheme (Fig. 3(c)). The position of a new point 0 on the shock wave is determined by the intersection of three planes. The first plane with the norreal m0 is passed through the vector t tangent to the shock wave at the point 3 already calculated, assuming in the first approximation no = as. The second plane is a characteristic plane given normal nx and passing through an internal point I with the known flow parameters. For convenience ofcomputation, the point I is chosen in the common meridional plane -- const with the point 3. The flow parameters at the point I can be obrained by interpolation between the known internal points 5 and 6, while the characteristic plane itself can be drawn parallel to the vector connecting these two points. The third plane is the above-mentioned meridional plane, the normal to which ~ is also known. The position of the point 0 can be determined with the aid of these three normals and of the known vector 6. It is not difficult to obtain the relationships aria il =

ns-(nl×nz) (nl×nl), So = it-J-Is.

Let us consider now a geometrical condition on the shock wave whose equation is x = x~,0', z). Let us take a local Cartesian coordinate system s, t at the point 3, directing the s-axis from the point 3 to the point 0. Then it is possible to write approximately

,at as/,

at )o

at

(4.1)

But 82x~ 8t Os

8~x,, 8s 8t "

The last derivative can be represented in terms of differences at two points 2 and 4 on the shock wave. We obtain

/ a,x.

= / ax.

where t24 = t 2 - - t 4 .

_/ax.

(4.2)

68

V. I. CHUSHKIN

Substituting further (4.2) into (4.1), returning to the usual Cartesian coordinates x, y, z, and taking into account the formulae (3.9), we shall have, finally, the following geometrical relation

ny

ny

nz

nz

where the corresponding vector components are denoted by the indices x, y, z. In calculating a point on the shock wave, we still need the compatibility relation (3.18) on the characteristic plane under consideration. Moreover, the relations on the shock wave (3.6) are used here. If, at the point to be determined, the normal no to the shock wave (actually its two components, since ~ = 1) is given, then all the flow functions can be found from the formulae (3.7)-(3.10). However, for two components of the normal arbitrarily assumed, the geometrical condition (4.3) and compatibility relation on the characteristic plane (3.18) will, generally speaking, not be satisfied. In order to satisfy these two equations within the desired accuracy, it is necessary to choose the normal no by using a procedure based on the application of the Newton's method with simultaneous averaging of parameters between the points 0 and 3, and between 0 and 1. In performing calculations with this scheme of the method of characteristics, the manner of introducing new points into a computing algorithm plays a significant role. Let the solution be started from a certain given characteristic surface of the "second family" which is called so by the analogy with two-dimensional case. The whole of this surface i~' divided into a certain number of rings containing a number of points. The first ring is formed by the points lying on the body surface, and the last ring by the points lying on the shock wave. The gas dynamic functions on each ring are written down in the form of a matrix whose rows consist of the coordinates of the given point, the velocity vector V and entropy S at this point. The number of rows in each matrix is equal to the number of points considered in a ring. Further, each matrix is put in correspondence with a certain point on the plane. Then the introduction of matrices into the computational algorithm will be analogous to that of the mesh points in the corresponding two-dimensional problem with the initial data on some characteristic of the second family. Let us illustrate the tetrahedral scheme of the numerical method of characteristics above described by several examples of supersonic flow about

TIIREE-DIMENSlONAL S U P E R S O N I C

FLOW'S

69

blunt bodies calculated by Podladchikov. All results will be presented in the system of cylindrical coordinates connected with the body and having the origin in the front point of the blunt nose (Fig. 2). The distances are referred to the radius of blunt nose, the adiabatic index being y -- 1-4. Below (in Figs. 4-6) some results are given for the cone with spherical nose and semiapex angle o~ = 9.5 °, at the free stream Mach number M , = 6 and the incidence = -- 5°. The sonic points are located in this case on the blunt nose, and therefore the characteristic surface of "the second family" with initial data was obtained by rotating the known numerical solution for supersonic flow about a sphere.

i i

p

~

~

"~. 3

I !

'

4.~"

~,,~,,,~

~90i ,8o.

2

3

4

x

L

5

6

F l o . 4. P r e s s u r e d i s t r i b u t i o n o n a b l u n t c o n e , oJ = 9.50, M . IO

7 = 6, = = 5°

!

O.g~

Cf

0.8 --

~

C¢

---tO.lO 0

0

9

0.7-~

0.08 ~

0.6--

~'~

o..-

7 o.o I

2 1,

Fro. 5. Variation of aerodymunic ~ t s M,. =6,===

3

4

alon8 a blunt cone oJ = 9.5°, 5°

70

r.z. cmJsugaN

Figure 4 shows the pressure distribution on the body along several generators in the planes ~0 = const. The pressure is referred here to the pressure in the undisturbed stream. As seen from thi,~ graph, the calculation on the leeward side (~0 = ~z)goes considerably further than on the windward side (F = 0). The reason is that the inclinations of the characteristic surfaces on the windward and leeward sides are rather different. This leads to considerable distortion of the rings with calculation points which were originally chosen in a uniform way. Such a distortion introduces great I

I

!

;

; 3.5

t

x:

i

/2~2.s ~2 .o

W 0.0

30

60

90 ~,,

Flo. 6. Variation

'20

tSO

80

de(J

of pefiphexal velocity on a blunt cone (o = 9.5 °, M . -- 6, ~,=

5°

errors and brings the computational process to a stop. This phenomenon becomes worse with the increase of the angle of incidence. For comparison the results by D'yakonov(') calculated with the help of finite-difference method are shown on this graph by circles. Figure 5 shows the variation of aerodynamic coefficients along the body under consideration. These coefficients are referred to the dynamic head and to the cross-sectional area of the body at the ~ven value of x. Curves are presented here for the coefficient of tangential force C,, the coefficient of normal force C,, the moment coefficient C m (with respect to the axis z), and for tho position of centre of pressure given by Cd = Cm/C,. The variation of peripheral component of the velocity w with the angular variable ~, for several values of x, is shown in Fig. 6. The value of the velocity w is referred here to the maximum adiabatic velocity ~ / m a x " The direct tetrahedral scheme permits the accurate calculation of the threedimensional supersonic flow past bodies with angular points on their

THltEE-DIMEN$IONAL SUPERSONIC FLOWS

71

contours. In Figs. 7-9 we shah present some computational results for the flow past a body whose configuration is formed by joining a spherical nose part with the semiangle oJ, = 80 ° and a conical aft part with the semiangle oJ = - 1 5 °. In this case the incident stream has the Mach number M . -- 4 and the incidence ~ = 5°. 1.5

I.O--

//

O.

/"

t

X

-0.5

-1.0

-I.5

F~3. 7. How-fieldpattern on the plane of symmetryfor the body with angular point on generator, at M . ffi 4, ¢c ---- 5° In the vicinity of the angular point of the body generator, a threedimensional rarefaction wave makes an appearance. Such a flow was investigated by Mikhailovcu~ who obtained a solution in the form of series. Using the results of thi~ work, a certain number of additional matrices is introduced on the fine of angular points. The computational scheme will then be similar to that for smooth bodies. Figure 7 shows the body investigated as well as the flow pattern in the plane of ~flmmetry. The crms-sections of the shock wave and of characs~trfaces of the "first and second families" in the meridional planes

72

P. I. CHUSHKIN

~o -- 0 and ~ -- :r are constructed here. Figure 8 presents the distribution of the relative pressure p/p.. (solid lines) and the relative density 0/o_ (dotted lines) on the body, for the same planes. Of some interest is Fig. 9 showing the pressure change in the plane ~ = 0 along a series of the "second family" characteristics designated by their numbers. As can be seen, when crossing the fan of the "first family" characteristics emanating from the line of angular points on the body, the pressure curves (as also those for other flow functions) show discontinuities of their derivative which decay slowly. In such a way, this physical 0.75 .

!

I

!

! o.IJ 02

i/

! i

I

iI

i t

!

i

i

j

i 2

FIG. 8. P r e ~ u r ¢ distribution (solid line) a n d density distribution (dotted line) along a body with angular point on generator, Moo = 4, ,, = 5 °

feature of the flow has been correctly rendered by the numerical scheme applied. The results given above and experience of calculations confirm that it is purposeful to use the direct tetrahedral scheme in those cases when one should obtain a solution in the zone stretched not too much downstream, hut it is necessary to take account of the region of influence accurately. To such problems belong, for instance, the three-dimensional flow in the vicinity of angular points on the wall (both for external flows about bodies and for internal flows in nozzles) and the formation of secondary shock waves in the flow field. Besides, the initial data found on a certain arbitrary surface in calculating the flow past the blunt nose at incidence can be re= computed with the aid of this direct scheme for a certain plane which is

73

THI~-DIMm~IONAL SUPERSONIC FLOWS

located lower downstream and serves as the plane of initial data for inverse-type numerical schemes. The defects of the scheme are its low effectiveness for large incidences and long bodies, as well as its considerable complexity. 12.~

IO.C

D Pw

,4 ,7

5.G

,9

2.5

0.5

0.75

i.o

1.25

I.~t

1.7*3

2.0

x

Fro. 9. Variation of pressure along charactericticz of the "second family" in the plane y =. 0, for the flow past a body with angular point, on generator M.

5. I N V E R S E

= 4, = = 5°

TETRAHEDRAL SCHEME INTERPOLATIONS

WITH

CONTINUOUS

This scheme of the three-dimensional method of characteristics has been proposed by Minostsev in his papers. ~a, ~) Here, in solving the problem of supersonic flow past a body of revolution at incidences the cylindrical coordinates x, r, ~ and the Cartesian coordinates x, y, z connected with the body (the axis x directed along the body axis) are considered simultaneonsly. The three velocity components V~, V~, V, along the Cartesian axes x, y, z and the pressure p are taken as the basic functions. The solution of the problem is carried out layer-wise, the layers being bounded by the planes x -- const. The normalized variable ~ defined by the formula (3.3) is introduced, and on each such plane the network (Fig. 2) is used whose mesh points are formed by intersection of the surface ~ = ~, (n = 0, 1 . . . . , ~ and the meridional planes ~0-- ~0k = k ~ / K ( k = O, 1. . . . . K). On each layer x -- cons1, gas dynamic functions accounting for their

74

P.I.

CHUSIIKIN

periodi~ty with respect to lp (the problem with a plane of symmetry is considered, 0 ~ ~0 ~ ~) axe represented by sections of the following double series (-~(x, ~, ~) =

7~ ~

x-x ~ a ~ (x)~" sin kv,

n~0

k--1

(5.1) ~ x , ~, v) = ~

~ b.~ (x)~" cos kv,

n--O

k--O

where odd and even functions in ~0 axe denoted by ('~ and (f~, respectively. For improving the approximations in ~, the zerm of Chebyshev's polynomials are taken here as mesh points for interpolation. For the shock wave to be determined (~ = 1), the following interpolational representation is applied: £

r,(x, ~o) = ~ ck (x) cos kv.

(5.2)

k--O

In computing flow parameters on a new layer xo = x , + Ax the M a c h cones are issued f r o m all mesh points on this layer in the direction towards the preceding k n o w n layer x = x , . The M a c h cones are constructed using the k n o w n data at x = x , . The equation o f these conical surfaces

(see, for instance, Ferri (u)) are of the form

( ~ + ~ - a , ) ( ~ ) , + ( ~ + ~ - a , ) (Ay), + ( ~ + ~ - a D (Az)a-2 r ~ v , Ax Ay - 2VxV, A x A z - 2VyV, Ay Az = O,

(5.3)

where

Ax = xo-x,

Ay = y o - y ,

AZ = Zo-Z.

Here, the sl)m, 0 denotes the values at the mesh point to be computed, the quantities Vz, Vy, V,, a are taken at the point 0', where x = x,, and the valu~ of ~ and ~ are the same as for the point 0. Substituting x = x . in the equation (5.3), we shall find the intersection line of the Math cone with thi~ plane x = x.. It is convenient to write the equation of this intersec~on line ¢~v, z) = 0 in local polar coordinates R, 0 on the plane x = x . whose pole is located at the point 4 where the stream line i u u a t from the mesh point 0 to be computed crosses this plane x = x , (in this case the measuring of 0 is made in the same way as that of,p). This point has the coordinates

Y' = Y°-(--~-Z)o, Ax,

z, = Z o - (

"V-~)o, Ax.

THR~-DIMBI~IONAL

75

SUPI~t$ONIC FLOWS

The equation of the intersection line then becomes Ax

R = v~, ( ~ - a * + ~

{v'(V'-,r') [ ~ + d ) -

v.g,+ v.~r']~

- ( V , 4 V ~ , - V ~ V . ) g - V~ V.(V, sin @+ V, cos@)

-(~-~ (v,, sin o +

(5.4)

v,, cos 0)},

where g=

V~, sin ~)- V. cos O.

Here the quantities with suffix 4 refer to the relevant point whereas the quantities without a numerical suffix to the point 0'.

(o)

(b)

~

0

I

(c) I~o. 10. Cakulation of typical points in the invene tetrahedral scheme with conttmmus intmpolatiom

In this tetrahedral-type numerical scheme the three bicharacteristics 0-1, 0-2, 0-3 and the stream line 0-4 (Fig. 10(a)) are taken for calculating a point inside the flow field. Having chosen any bicharacteristic defined by the angle 0~, it is possible to find the quantity R! from the equation (5.4) and then to determine the intersection points of this bicharacte~tic with the plane x = x,, viz. zt = z,-.R~ sin Ot.

The compatibility conditions (3.15) are valid on the characteristic planes, where the vectors sx and ss lie on the characteristic plane containing a selected b ' r . h a r a a m ~ . Let us take as the vector st the unit vector directed along the tangent to the curve R = R(O) defined by the equation (5.4). It is not dlmcult then

76

P. I. CHUSHKIN

to get the following expressions for components of the vector Sl in Cartesian coordinates: $1 x =

O,

,.~ly - -

~y

,

Slz

=

~Z

'

where ~ y , z) = 0 is the equation of the curve R = R(O) in Cartesian coordinates. As the vector s, we shall take the vector passing along the bicharacteristic in the direction towards the point 0 to be computed. Obtaining the components of the vector s~ in Cartesian coordinates also presents no difficulties. With the chosen directions of the vectors sl and ss, it is easy to find the vector ~ of the external normal to the characteristic plane, which is included in the compatibility relation (3.15). At the intersection points of bicharacteristics i -- 1, 2, 3 with the plane x = x., the vector at will be determined by the formula siXs2 n , = [slXs~[ "

(5.5)

Along the stream line 0-4 the condition of constant entropy (3.19) and the Bernoulli equation (3.20) are satisfied; they can be rewritten as follows o dp-~,p do = o, oVdV+dp = 0.

(5.6) (5.7)

Representing the compatibility relation (3.15) along three bicharacteristics and the equation (5.7) along the stream line in the finite difference form, we obtain the following system of four algebraic equations (i = 1, 2, 3, 4) for the values of the velocity components V~o, Vyo, V~e and the pressure po at a new mesh point 0 C,'Vo+~)@o = Et,

(5.8)

where, for i = I, 2, 3: Cl=~ta~A~,

~t=

--Bl~,

and for i = 4 C a = o,V4,

fD~= 1,

EL = p4+C*'Vt. Here, the components of a unit normal vector at are found from the equation (5.5), A~ and B# from (3.16), and the differentials of velocity components along the direction sl are calculated from the formula (3.17). The

77

THREE-DI]t~EN~ONAL ~LJ'PEIt~qONIC FLOWS

derivatives ~ / a y and ~--j~/Sz are expressed in terms of the appropriate derivatives in the cylindrical system of coordinates, viz. $y .-- e aT:

Oz

_

r ~ 1 (~ cos ~,+ sin v) - -

e

~o ' cos V ~ - f

a~

--, 8~

r

where

~'=-V~

~

~0

+

"

The derivatives with respect to the variables ~ and ~ on the layer already calculated, which are contained in the above expressions, are in accordance with the adopted approximations (5.1) and (5.2). The derivative with respect to x on this layer is computed by a finite difference formula, and the quantity e is expressed by the formula (3.4). Upon obtaining the quantifies V~o, V,e, V:e and/7o from the algebraic equations (5.8), it is possible to find the density 0e at a new inner point from the condition of constant entropy (5.6) along the stream line, written in the finite difference form; the speed of sound a0 is determined from the last formula (3.2). For calculating flow parameters at a point located on the body surface (~ = 0) two bicharacteristies 1-0, 2-0 and the stream line 4-0 on this surface are taken (Fig. 10(b)). Also, the condition of vanl.qhing normal velocity (3.5) is used here which, for a body specified by the equation x, y, z) - 0, will evidently be

V . v / = 0. The gas dynamic functions at a new point on the body are determined from the same system (5.8) as for a point inside the flow field. However, the coefficients of this system for i = 3 will be as follows C3=vf,

~:=Es=O,

and since the body equation is known, the derivatives in ~Tfcan be taken at a new point xo = x . + Ax. Let us turn now to calculating a new point located on the shock wave surface to be determined. In this case (Fig. 1G(c)) one should consider the compatibility condition (3.15) along one bicharacteristic 1-0 (01 -- ~0k +~z) and the relations on the shock wave (3.6), where the last vectorial relatior~

'78

P . I . CHUSHKIN

can be presented as two equations: Or,, . ¥ ,) - Or, V~-(V, cos ¥ + V: sm ~ - = V~..- V~. cos~o Ox ' V, sin~

r,

0~

-g,

(5.9)

(c o s ~ m,r, Or.1

~.) = V ~ . ( sin~o cos, r,,

(5.10)

Here, the shock wave equation is taken in cylindrical coordinates in the form r = r,,(x, lp). The coordinate r . , of a new point on the shock wave is obtained from equation (5.9) written in the finite difference form ( V ~ - Vx.) zix r,,e = r~,-t ( V ~ , - V~.) cos ~0+ V~, sin ¥ ' where all the quantities in the right-hand side of the equation are evaluated for the known preceding point 0' on the shock wave, located in the considered meridional plane. Having found in this way the radiuJ rw in all the meridional planes tp = tpt (k = 0, 1. . . . , K), one may approximate the shock wave by the expression (5.2) and determine from it the corresponding derivatives Orw/~, at every point. The value of the derivative 8r,,/ax, close to the value of this derivative at the point 0' is then assumed approximately at the point 0 of the shock wave to be calculated. The normal component of the velocity k',. can then be computed by the formula (3.11), and the three velocity components Vx, Vj,, V:, the pressure p and density 0 by the formulae (3.7)-(3.8). Further, it is necessary to check the correctness of the assumed value of 8r,,/Sx at the given point on the shock wave. Such a check consists in testing that the compatibility relation (3.15), taken for the bicharacteristic 1-0, is satisfied. Repeating the computationalcycie described, we find the final value of Or,,/Sx with the desired accuracy by Newton's method. By using this scheme, Minostsev (*a~calculated the supersonic flow past bodies of different shape for the case ofperfect air (7 = 1.4), and examined this scheme numerically. The practical computations showed that, for the flows past smooth bodies, it is sufficient to take the approximating polynomials (5.1) of the orders K = 3 and Ot = 7. However, for flows with more rapid change of gas dynamic functions, it is necessary to increase the order of polynomials with respect to ~ up to Ot = 11 whereas the effect of K is considerably less. It turns out that, for the bodies with discontinuities of curvature, it is impossible to perform the continuous calculation

THIt~-DIMBNSIONAL SUPIHt~NIC PLOWS

79

of the flow field by using this numerical scheme. Here, one can compute separately the flow in regions divided by the weak discontinuity line, or the body contour at the point of discont/nnity of curvature can be smoothed analytically with respect to several derivatives. However, it is most convenient (if possfole) to approximate the entire body by some body close to it and having an analytical contour. Let us reproduce now some numerical results of Minostsev.

iI

Fro. I I. Flow-field pattm'ns in the plane of symmetry for a n ellipsoid, at • ffi=5° a n d various ~ Mach numtmrs M .

Figure 11 shows, for an ellipsoid of revolution with ratio of semi-axes equal to 1-5, the pattern of supersonic flow field at incidence 0c - 5 ° and at two Mach numbers M . --- 3 and /14_ = ~. The shock waves and sonic lines in the plane of symmetry are shown on this graph. The dotted line shows the results of Tinyakov's numerical solution t'~) serving as initial data. The computation of the field of supersonic flow was carried out here by the scheme of characteristics discussed above but transformed for the system of spherical coordinates convenient in the given case. The distribution of pressure related to the stagnation pressure behind the normal shock is plotted in Fig. 12 for the leeward side of this ellipsoid. The results of flow computations ( M , -- 20, g -- 5°) past another smooth body, whose contour is composed by a circle and a cubic parabola (at the matching point at x -- -0-342 the first three derivatives are con-

80

P. [. CHUSHKIN

tinuous) are shown in Figs. 13 and 14. In this case, the linear dimensions are referred to the radius of the blunt nose, and the pressure to the quantity o.aSffi. The shape of shock wave and the pressure distribution over the body in three meridional planes are shown in these graphs. The corresponding data calculated by D'yakonov ~6> by the finite difference method are plotted as small circles.

I 0-'--

\i

0.2 P

I! +

!

0.I

-0.5

-0.4

-0.3

-0.2

I

-0. I

+

0

0

x

FIo. 12. Pressure d i s t r i b u t i o n o n the l e e w a r d aide o f a n ellipsoid, at • -- 5 ° a n d v a r i o u s free-stream M a c h n u m b ~ Moo

I

f

2

0

-I

0

i

2 x

3

I 4

5

FIo. 13. S h a p e o f a s h o c k w a v e for a b l u n t - n o s e d s m o o t h body, at M ~ a n d ~, ~ 5 °

-- 20

~=DIMENSIONAL

81

SUPERSONIC FLOWS

Minostscv ~'m also considered the supersonic flow past a body with analytical contour, whose nose part is close to a spheric segment with the sendangle m., and the rear part to the "inverse" cone with semiapex angle w. In the system o f spherical coordinates r, O, ~, with the pole displaced from the symmetry axis, the equation o f the family o f such bodies will be

[rS+2r(xo cos O+y o sin 8 cos ~)+~+yo2]tl ~ +

arctan

(r cos O+xo) sin c o - s i n (m+o~,)

0.8

I'

II

I.

iI

1t

o':t\

0

=

r

'

P 0.4

,

\\\ i

2

/

3

x

Fio. 14. Pressure distribution along a blunt-no~'d ~mooth body, at M.. = 20 and ¢ = 5°

Figure 15 shows three bodies o f such a type for I = m - 40, co. = 30 ° and a series o f values o f o~ (the shock wave is presented here for the body with m = 0°). The computations of flow past such bodies were performed for M _ = oo and0c = 15 °. The pressure distribution, related to o..V2m~, is shown on the surface o f these bodies in the meridional planes ~0 = 0 (solid line) and ¥ = :z (dotted line). These computations verify that the inverse tetrahedral scheme with continuous interpolations gives the desired accuracy for the flows past sufficiently smooth bodies.

~2

P I. CHUSH'KIN

/

2.0

1.5

1.0

/

J

/

t~=O*

Y O.f,

I i i

io

-0.~,

-I.C 0

0,5

g

1.0

1.5

Fro. 15. Shapes of several blunt bodies with analytical contour 6. S C H E M E

USING ONLY BICHARACTERISTIC DIRECTIONS

The scheme of three-dimensional method of characteristics, in which only bicharactcristic directions are used, was developed by Magomedov
W

6 = arctan ~/(u=+v=) "

The angle ~ represents the angle which the velocity vector V forms with the meridional plane ~0 = const. It is obvious that u=cosflcos6,

v=sinflcos~,

w=sin&

THRBE-DIMEI~IONAL $ ~ N I C

83

FLOWS

0.14

0.12

GIO

~

aoe

=0 •

il

p

II II 0.1~!

h

? 0,02

",.,, f -20"

. . Ii

Q. . . . .

. . . .

=

N 0

,.8

I~

x

Fro. 16. Premure distribution on windward side (solid line) and on leeward side (dotted line) of blunt bodies with analytical contour, at M . - =o and e~ - 15° Having in view the generalization of the present scheme for the case of gas dissociating ;,I equifibrium conditions, the functions C~) and H are considered here instead of the pressure p and enthalpy h: ~ = lnp,

H=

lnh.

We introduce also the local Cartesian coordinates connected with the velocity vector V. The unit vectors k: in this system are chosen as follows:

kx =

V

I

-V'

k , = cos-----~ ~/~ '

8ki

8kx k, =

~6 "

Assuming that/~, ~, ~D, H are the basic functions to be determined, the system of gas dynamic equations (3.1) can be written in the following form:

c o s d k s . V O + k s . V d + Q ( M ~ _ l ) k x . V : ~ ) ~ sin/~cos d = 0, r

cos ~$kx. v/~+Qk,. V ~

cos ~ sin2 ~ = O, r

kx. v 6 + Q k , , v~C)+ sin ~ sin 8 = O, r

kx. v H - ~

1

kx. ~:D = 0,

(6.1)

84

P. I. CHUSHKIN

where P Q - 0V 2,

Z = O--~-h p

(6.2)

and, as usual, 0 is the density and M is the M a t h number. The Bernoulli equation (3.22) is also added to this set. For the perfect gas, we have simply I

7~

Q=

Z-

7,-I

"

The family of wave characteristics for the system (6.1) satisfies the equat.ion •v = c o s / d t x - sin p ( - s i n ~k2+cos ~ks), (6.3) where the vector x is directed along the bicharaeteristic, the parameter varies within the interval 0 ~ ~ -~ 2n, and p is the M a t h angle

Having chosen from this family two bicharaeteristies for q~ = ~/2 and q~ = 3~/2, i.e. • t = c o s / d t x + ( - l ) t sin #kz, (i = 1,2), Magomedov ~u~ derived for them the compatibility relations in the form cos~_~+(_l)

t Q0 dr]) ds~

( - l)l sin/zk3. ~Tb +

_

cos # c o s ~ c o s ~ [~sin ~ - ( -

1)~70]

(i = 1, 2).

(6.4)

Here d('~/ds t represent derivatives along the bicharacteristics de7 =

I n addition, some new symbols are introduced: ~7= tanS,

~= tan6

(6.5)

and O = t a n t~.

As seen, the compatibility relations contain, apart from derivatives along the bicharacteristics, also the derivatives in the non-characteristic direction i~. To eliminate these last derivatives, two more bicharacteristics from the family (6.3) are introduced, with 9~ = 0 and ~ = =, i.e. x~ = c o s / ~ k 1 + ( - 1)~sin/~k3

(i = 3,4).

~*DIM][NBIONAL

$~NIC

85

FLOWS

Hence

2sm~' so it is easy to obtain the equality

ks.~t:

-- 2 s i n o

In such a way; this expression permits to express, in the compatibility relations (6.4), the derivatives along the non-characteristic direction ks in terms of the derivatives along two additional bicharacteristics. Along the stream line 'rs -- kl the two following relations will hold: d/~ + -k

"~s

Q s" ~ )

d~ &s

=

sin ~ sin 6 •

I d~ =0. Z &s

'

(6.7)

(6.S)

In the scheme using only bicharacterisfic directions, the inverse-type computational network is taken with the layers x = const and with mesh points determined by constant values of the variables x, ~, V' (Fig. 2). Let us consider the solution of elementary problems of the method of characteristics in this scheme. In defining the parameters at a new internal point 0 in the plane x - x0, four bicbaracteristics 1-0, 2-0, 3-0, 4-0 and the stream line 5-0 are issued from this point towards the preceding layer x , = x o - A x (Fig. 17(a)). The coordinates of points where these lines meet the plane x = x , are I

!

2

(a)

(b)

(c) Fro. 17. ~lculaflon of typical points in the scheme using only biclu~cteristi¢

86

P . I . CHUSHK1N

computed as follows: rt = ro--~r,,

V* = Vo--AW,

(i = 1, 2, 3, 4, 5),

(6.9)

and, by using the functions ~ and ~ (6.5), we have Ax

Arl = )1AX, Ar5 = ~ Ax,

Ax

¢+(-lyO ~/(1+~ t-(-tyo¢

~v,=

(i ---- 1, 2),

(~ = 3, 4),

Ax

~vs = -7- ~ V(l + ~f),

where

d = 1-(-1)%0v'(1+~). We further present the compatibility relations (6.4) along two bicharacteristics i = I, 2 in the finite difference form, after first introducing in them the functions )7 and ~, and after eliminating the derivatives along noncharacteristic ~ o n s with the aid of the formula (6.6). We have then ~7o+A~'Do = +T,+A~'Dt+Bt

(l = 1, 2),

(6.10)

where

--t ~ l +,f) A~ = ( - 1) ~ ~/(1.+~), . . l + , f [ Lg) .~__] B, = ( - 1 r - - - j - [v,H T~) ~ ( ~ - n 0 V ( t +P)) . The operator L((-f') involved here is expressed as

Similarly, let us write down, in the finite difference form, the relations (6.7) and (6.8) which are s a t i r e d along the stream line. They will become

Ho = H s + l ( ~ o - ~ 0 ,

(6.11)

where

Four equations (6.10) and (6.1 I) make up the system for determining the values of the functions r/o, ~o, Co and He at the point 0. The following

THRE~DIMENSlONAL $ ~ N I C

FLOWS

87

computational algorithm is used here. Initially, some approximate values of the functions at the point 0 are assumed (they can be taken, for instance, from the point with the same coordinates ~e, ~ on the preceding plane x = x.). Substituting t h e ~ values in the formulae (6.9), the coordiuates ofaU the points 1, 2, 3, 4 and 5 are determined, and then the valuu of all the necessary functions at these points are found by two-variable interpolation, using the known data on the surface x = x.. We turn then to solving the system of equations (6.10)-(6.11). In the first approximation, the coefficients At and Bt(i = 1, 2) are calculated at the points 1 and 2, the coemcient Cs at the point 5 and, in the operator L(('~), the above assumed approximate value of ~ o is taken. Having obtained the values of the functions at the point 0 it is possible, generally speaking, to repeat a waln the above-described computational procedure by now averaging parameters along all the lines O - i , although the first order of the numerical scheme does not require thi~ (see Magomedov~tO). In calculating the point 0 on the body surface (Fig. 17Co)) defined by the equation • = re(x , ,p), the condition (3.5) becomes 8r~ . ~ _ , . . . . . . 8rb ~/= ~ + " r Vtt ~'~') ~ -

(6.12)

The computational procedure in this case will be similar to that for an inner point of the flow field, but here the compatibility condition (6.10) is considered only along one bicharacteristic t = I and, instead of the compatibility relation for t = 2, the condition (6.12) is applied. Let us turn now to determining the flow parameters at the new point 0 on the shock wave (Fig. 17(c)). The computational procedure is based here on an investigation by Magomedov
1 $r,, t - r. 8~o

If 0' is the known point on the shock wave in the preceding layer x = x . then, with accuracy to the first order in the step Ax, we have r~e = r,~r + qo, Ax,

(°,)

to = re, + - ~ - o"

(6.13)

88

P. 1. CHUSHKIN

Using the condition ar,,t aq ~x = ~-~' it is possible to compute the value of (St/Sx)G, with the aid of difference formula, knowing the data in the plane x = x.. In such a way, the quantities r,, and t at the new point on the shock wave can be determined by the equalities (6.13). Since, under the known conditions of the incident flow, the values of gas dynamic functions immediately behind the shock wave depend only on q and t, so it is necessary to pre-assign in addition an approximate value of qo at the point 0. Then, from the formulae (3.10), (3.8) and (3.2), we find Vn., p, h and later ~D and H. The expressions (3.12) for the velocity components u, v, w permit to obtain the values of ~/and ~ at the point 0 -

= r " -

(~

~ sin ot cos V + At, ~cos,t-qAp ' ~ sin ct sin ~ - t Ap cos ~-q

Ap) V ( I + ~ )

(6.14) "

The compatibility condition (6.10) along the bicharacteristic 2 - 0 (i = 2) serves for checking the correctness of the pre-assigned value of q0 and its subsequent improvement until attaining the required degree of accuracy. By the scheme using only bicharacteristic directions, Magomedov calculated the three-dimensional supersonic flows past various bodies (sharpnosed and blunt-nosed bodies of revolution, delta wing with blunted edges) in the cases of perfect gas and of gas dissociating in equilibrium conditions. We present here some iesults of these computations, relating to the case of a perfect gas (7 = 1.4). Figure 18 illustrates the flow past a cone with a spherical blunt nose and semi-apex angle to = 9.5 °, at Mach number 3 4 . = oo and incidence = 5 °. The distribution of pressure p, depending on the coordinate x related to the nose radius and measured from the origin of the conical part of the body, is plotted on this graph, for three meridional planes. The solid lines show the pressure on the body surface, and the dotted lines that on the shock wave. The function ( is also plotted here, which is associated with the peripheral velocity (( = w/x/(u 2 +v2)) and characterizes the transverse flow about the body. The results of computations by the method of characteristics are in good agreement with the D'yakonov (4) data calculated by the finite difference method and shown in the graph by small circles. The flow past a slender cylindrical body with a sharp-nosed part is presented in F i g . 19 for the case 3,/. = 5, x = 5 °. Here, the pressure

89

TIffin-DIMENSIONAL SUPERSONIC FLOWS 0.11

& i

o., \

m

~,o

?,

7.

/

~o., I

0

J

4

8

1

12

I

i6

I

20

X

F ~ . 18. V~rtadon of W e m ~

on • blunt cone, m - 9 . ~ (solid liae) ~ . d on the

shock wave (dotted line), and variation ot im'Iphend veloc/ty function (duheddotted line) on the b o d y , at M .

-- ~ , a --

f.O

0.8

lOp 0.6 ,,0 0.4

0.2

0

FIo. 19. Flow l ~ t a sharp.nce~ body, at M . -- 5 a n d a - 5°. Premure distri. butkm on tim body and shape of the shock wave

90

P.I. CHUSHKIN

distribution on the body and the shape of the shock wave for W = 0 °, 900 and 180 ° are given. As seen, on the pointed elongated bodies (as well as, for example, on a blunted cylinder), the pressure on the leeward side (F = 180 °) varies non-monotonicvlly and. within a certain interval of lengths, may be higher than that for ~p < 180 °.

--~

0

.

5 0.4G.

r. r

--

~ ~ ~

•

I

I

Z

~0.3

i

4

1

6

/o.z

8

I0

X

l ~ o . 20. Shape of the shock wave for a delta winil (Z = 70 °) with blunt edges, atM. = 6 , = -- I~

0.20

O.IE

Y X x

Z

0,,2 P

| 0.08

0.04

m~.~ • ~

"-.,-

1

2

I

4

I

6

I A

10

F l o . 21. Premure dietriimtion on a delta win8 with blunt edsu, at M . ~==0 °

,= 6,

The numerical remits for a more complex three-dimensional problem of flow past a triangular wing o f finite thickne~, with an angle of sweep X = 700, are shown in Figs. 20 and 21, for M . = 6, = -- 0 °. The nose of

~-DIMENSIONAL

SUPEI~ONIC FLOWS

91

the wing represents a part of a sphere, whereas the lateral edges are blunted as parts of circular cylinders. The flow pattern past such a body is of compficated nature and is determined by interaction of the axisymmetrical flow past the spherical blunt nose and that past the edges on the flat part of the wing. Figure 20 presents the shape of the shock wave in two planes of the symmetry: in the plane z = 0 passing through the wing edges (solid line) and in the plane y = 0 (dotted line). These results show that the shock wave at considerable distances from the nose remains axisymmetrical and then, in the plane z = 0 behind the inflection point, it runs almost parallel to the lateral edge. Figure 21 presents the pressure distribution on a delta wing along the lateral edge (solid line), and along the central stream line (dotted line) which for oc = 0 ° is the singular line (the coordinate x is related to the nose radius). It is seen that the pressure along the lateral edge approaches, rather rapidly, to the asymptotic value equal to the pressure on an infinite cylinder with the corresponding sideslip angle. The experimental data obtained by Michel and Duong-Vinh-Hung (st) for the lateral edge at the free stream Mach number M . = 7 are also plotted on the same graph. Taking into account some difference in the values of Mach number Moo, the agreement of experimental and computed results is quite satisfactory. Various applications of the numerical scheme using only bicharacteristic directions confirm its effectiveness and convenience.

7. S E M I - C H A R A C T E R I S T I C

SCHEME

Let us consider now the sem/-characteristic scheme by Katskova and Chushkln(6X) intended for computing the three-dimensional supersonic flows past the smooth bodies. This scheme uses the cylindrical coordinates x, r, ~0 connected with the body (Fig. 2). The three velocity components u, v, w in this coordinate system, the pressure p and the entropy S are taken as the basic gas dynamic functions. The region of flow between the shock wave • = r,,(x, tp) a n d t h e b o d y surface r = rb(x, lp) is straightened by replacing the variable r by the variable ~ (3.3). The original system of partial differential equations of the problem (3.1), in which the continuity equation is previously reduced to the form O a s ~ . V + V . v p = O,

92

v.i. CHUSHKIN

will then be written in the variables x, ~, ~oas follows: Su

ax

1 (

8u

Sv

8w\

1 (

8p

+7 a~+~+~,-~-)+~- u ~ + ~ au

~p)

+~,=o,

au l {ap + a ap

"~+"~+~t~

7 ~ ) +~'' = o,

8v 1 8p u ~8v+ ,,-~ +-ff ~+¢,~ -- o,

(7.1)

aw . p 8p . ~ u -aw~ +. ~-~-+-ff ~ + ~ , = o, 8S

u~-~+

where

-~

I

Sr.

$rb~

arb

ax

ax]

ax'

~a

w(a~_w)

=

+ ,

=

O,

[lot. )Or ,l

1~

t'=-r

+

av

,

avJ

1 8

e = r,,(x, v)--rb(X, V),

i~, = 7

8S

v ~ - ( + Os

v = --(;tu+v+#w),

I

~ * " ' ~ ' - , av' ~, = ~ l ( W aw + , a , +VW ) r

r

-~

9~

'

w 8S

05 -

r 8~o" T i m transformed system contains the derivatives with respect to all three variables x, ~, ~,. In order to exclude the derivatives with respect to ~ from the system (7.1), the corresponding functions are approximated in this variable. For this purpose, let us take the system of the meridional planes ~o = ~k = = k~/K (k = O, 1. . . . , K) in the region 0 ~ ~0 -~ ~. Furtl2er, we represent the odd functions (-f (i.e. w) and the even functions ('jr (i.e. r, u, v, p, S) by the interpolational trigonometric polynomials in V' K-1

( ' ~ x , ~, ¥) =

~

ak(x, ~) sin k v,

k--1

(7.2) K

(-~(x, ~, V) = ~ bk(x, ~) COSk v. k--O

THREE-DIMBNBIONAL $ ~ N I C

FLOWS

93

The coefficients a t and b t are expressed in terms of the values of the approximated functions (-Fj and (~1 in the interpolation meridional planes, and are of the form • --1

•

= E

bt = E dtf ,.

J--1

k--O

From the approximations (7.2) the expressions for derivatives with respect to tp are as follows:

= E ekf'~l, k

J--1

---k

J

ftlt~l.

(7.3)

0

Let us note that determining the numerical coefficients c~, dk~, etj and

f~, depending only on K, presents no diflictdty. The paper Cry contains tables ofthese coefficients for K = 2, 3, 4. Hence, after using the approximations (7.2), the three-dimensional system (7.1) reduces to the approximating system of differential equations in two independent variables x and ~. However, this two, dlmensional system will include a greater number of dependent variables, namely the values of the basic functions in all the meridional planes ~ = ~0k, which should be determined. This approximating two-dimensional system, in the region where ua-t (~ +/~w)l >- a 2, 1 +pl will belong to a hyperbolic type and have two families of the real characteristics. Introducing two new functions ~ and ~ = l(v't'Pw)'--u

fl=N/[-~-+(l+P')(~-l)]'

(7.4)

the differential equation of characteristics of the first (i = 1) and second (l = 2) families will become d~

~-t

j

--- A,.

(7.5)

The compatibility relations for these characteristics may be written in the following form: d ~ - ~ t ~ + N , dx = 0 (i = 1, 2), (7.6)

94

P.I.

CHUSHKIN

where 3t = (- 1)* eu2 , N'=I[

Oauam ~-(-1)tflui-a' - O ' u'(-(-l)'a'fl~-a 2

O,+O,~-uw--dl:x]. One more family of lines which, by analogy with the axisymmetric case we shall conventionally call "stream lines", possesses the characteristic properties. The differential equation of "sla'eam lines" has the form -~

(7.7)

(~+~) -- As,

=

and the following relations will hold along them:

= O,

(7.8)

du+G dv+L d p + Q dx = 0,

(7.9)

dS + E dx = O,

(7.{o)

B dv-dw+Cdx

where B =#, L = I oU '

I

C=u(#Os-O,), Q = I(O2+¢Os),

G=¢, E =

O5 U

We may note that the equations (7.5)-(7.10) resemble by their form the corresponding equations for the axisymmetrical ease, and the latter can be obm_ined from the former if we put/~ = w = 0. This computional scheme uses the inverse-type network with mesh points formed by the lines ~ = const and ~p = eonst, and the solution is built up following layers of x = const. Unlike the t h r ~ schemes described above the semi-characteristic scheme now considered is of the second order of accuracy with respect to the step Ax, comprising the iteration process. Such iterations are performed simultaneously for all points having the same coordinate ~ and located in all meridional planes ~0 = ~s,. Let us consider now the practical computational algorithms for calculating partitular types of points. As usual, let us discuss first the computational procedure for an inner point in the flow field (Fig. 22(a)). We begin by taking any one meridional

95

THREE-DIMENSIONAL SUPERSONIC FLOWS

plane ~ = ~k. Let us issue characteristics of the first family 1 - 0 , of the second one 2-0, and a "stream line" 3-0, upstream from the point 0 on the plane xo = x . +Ax. The points of intersection of these lines with the plane x = x., in accordance with the equations (7.5) and (7.7), will have the coordinates ~1 = ~o-AiAx (i = 1, 2, 3). (7.11) ,O

4,

4

0 2 3

(a)

~0

0

(b)

(c)

Fro. 22. Calculation of typical points in semi-characterictic scheme

We find the values of required functions on the layer x = x., at the points 1, 2 and 3, by quadratic interpolation with respect to one variable ~. Then, from the finite difference analogs of the compatibility relations (7.6), for t = 1, 2, we express the quantifies po and to l p0 = px+~s_-~-~l [(x-(z-~l~px-p2)+(N2-Nx)Ax],

(7.12)

~o = ~ l + ~ x ( p o - p l ) - N x A x , and, from (7.10), we compute the entropy

So = Sa-EAx.

(7.13)

Using two relations (7.8) and (7.9), written in the finite difference form along the "stream line", and the formula (7.4) for ~, we define the values of three velocity components

B(¢o[us-L(po-ps)-Q6~x]-~s}+(1 + G¢o) (ws+ C6a¢) Wo =

1 + G¢o'4-B#o (7.14) uo =

I 1 +G~o [us-LO, o-ps)-QAx+G(v3+t~owo)],

t:o = uo~o-#owo. Such a computational algorithm is used successively in all planes tp = ~0k, for all points having the value ~ = to now under consideration.

96

i,. I, CHUSHKIN

Thus, the first iteration in solving the equations (7.11)--(7.14) is carried out. A similar cycle is repeated in further iterations which converge rapidly (correct digits are obtained after three iterations). In the first iteration the coefficients At and the derivatives with respect to ~0 appearing in the quantities O~ are taken from the relevant points 4 on the preceding plane x = x., and the remaining coefficients in the finite difference equations (7.11)-(7.14) are computed correspondingly at the points I, 2 and 3. In further iterations, the averaging of all coefficients is carried out, and the above-mentioned derivatives with respect to ~ are found from the formulae (7.3) from the data of the preceding iteration. The determination of flow parameters at a point on the body surface (~ = 0) is similar in many respects to the calculation for the inner point just discussed. Here, the boundary condition (3.5) at the new point 0 will become vo = AoUo+Vo+t~oWo = O.

Hence, allowing for the first formula (7.4), it follows to = --40.

(7.15)

Let us now issue from the point 0 (Fig. 22(b)) a characteristic of the second family 2--0. It will intersect the plane x = x. at the point ~ = = --AsAx, where all the values required are found by quadratic interpolation. From the compatibility relation (7.6) for i -- 2, represented in the finite difference form, we obtain 1 p0 = p2+~-~2 (Co-~2+NeAx).

(7.16)

Taking the relations (7.8)-(7.10) along the "stream line" 3-0, we come back to the same equations (7.13) and (7.14) for the values of So, uo, v0, w0 at the point on the body surface. These two last equations, together with the equations (7.15) and (7.16) for ~o and Po, form the algebraic set which is solved by iterations for all ~ = ~ok simultaneously. The computation for a point on the shock wave (~ = 1) is also carried out by iterations in all meridional planes simultaneously. Let us describe first the computational procedure for one plane V' = Tk. Here, at the new point 0 (Fig. 22(c)), we assume the values of derivatives Orw/Ox and Or,,[O~ (the latter remains unchanged during one iteration), close or equal to the corresponding values at the adjacent point 4 on the shock wave. It is possible then to find the radius of the shock wave r.o

--

[/at. I +l~, a x I],~J1

r.+ L ~, 0x ]o

2

Tlllt~-DIMI~N~ONAL $ ~ N I C

FLOW~

97

Therefore, knowing the values of the above-mentioned derivatives, the basic gas dynamic parameters at the point 0 are computed with the aid of formulae (3.8), (3.11) and (3.12). Further, we issue from the point 0 a characteristic of the first family 1--0 which, according to the equation (7.5), will intersect the layer x = x. at a point with the coordinate ~x =

1-AxAx,

where, as usual, all functions are determined by quadratic interpolation. It is then necessary to consider the compatibility relation (7.6) along this characteristic (i = 1) ~ o - ~ x - ~ x ( P 0 - p x ) + N x A x = O.

(7.17)

This relation is used for checking the validity of the assumed approximate value of ar,,/Ox, and thlg is changed as long as to have the equation (7.17) satisfied. The computational procedure described above is repeated successively for all planes ~0 = ~Pk,obtaining in this way the first iteration of the solution for the points on the shock wave. In subsequent iterations, the computational algorithm remains the same but then, in determining Or,/~o and other derivatives with respect to tp by the formulae (7.3), the values of functions should be taken from the preceding iteration and the coefllcients Ax, ~x and NI are averaged. Katskova and Chushkin appfied the semi-characteristic scheme described above for investigating the three-dimensional supersonic flow past blunt cones, and through inlets, in their papers Cxs~and ~s~~s~ respectively. Some results of these calculations will be presented below for illustrating the given numerical method. Figures 23-25 refer to the case of flow past sphere-nosed cones, for the Mach number M . = = , incidence ~, = 10° and ~, = 1.4. Here, the data are presented for either three (small circles in the graph) or five (solid lines) meridional planes. The origin of coordinates is placed at the front point of the blunt nose, and the radius of this blunt nose is taken as a linear dimension of reference. The pressure distribution along three generators of the blunt cone, with the semi-apex angle oJ = 20 °, is given in Fig. 23. The angular points of the pressure curves correspond to points of junction of the blunt nose with the conical part of the body, and are conditioned by the discontinuity of body curvature at these points. D'yakonov's results, Cs~ obtained by the finite difference method, are also plotted for comparison in this graph, as dotted

P. I. CHUSHKIN

98 0.281

0.24

0.20

0,~6 l

--

P

O. tZ

0.08

-

0.04

i 2

0

Flo. 23. P a ~ u r c

] 3

~ 4

' 5

x

distribution a l o n g three g e n e r a t o r s o f a b l u n t cone, to = 20 ° at M, = oo, = = 10°

_

,/,, ,-/2

/

4:

I ¢

0

2

4

6

2

0

¢AJ - ' 0 e X

0.10

_

/

x

4

6

~:20 ° x

=2 . 0

0.25

Y P 0.05

x=3.9

x:2.t

°x.-6. 2

,,l 0.=

[

o.,5

C

{

0.5

i

F I G . 24. F l o w past a b l u n t c y l i n d e r , oJ = G°, a n d p a s t a b l u n t c o n e , u, = 2 0 ° a t

M.

= o= and = = l 0 °. C o m p a r i s o n o f s h o c k wave s h a p e s a n d p r e s s u r e profiles, at tp = 0

T H R E E - D I M I ~ S l O N A L SUPERSONIC FLOWS

99

lines. The curve with small crosses refers to the case of axisymmetrical flow past the body at = = 0 °. Figure 24 presents the comparison of flows past the blunt cylinder co -- 0 ° and the blunt cone co = 20 °. It is seen that, for the blunt cone, there appears rather quickly an inflexion of the shock wave on the windward side (~0 = 0) The same graph gives also the pressure profiles across the shock layer on the windward side of the body, at several values of x. They have different character for the blunt cylinder and the blunt cone.

~,20 •

x '3.9 "~' 5.9 S

// /// l

-12.5

I

-I0

I / / i

-7.5

0

/ /

O.S

I.O

F r o . 25. E n t r o p y profiles o n a b l u n t c o n e , m ~ 20 °, f o r ¥ -- 0 ( d o t t e d line) a n d *f - ~ (solid line), at M . ffi 0% ec ffi 10 °

The development of the entropy layer along the blunt cone with co -- 20°, in the planes ~0 -- 0 (dotted line) and ~p = ~ (solid line), is shown in Fig. 25. It is clearly seen here that the thickness of the entropy layer decreases rapidly with increasing coordinate x on the windward side. Next, three figures illustrate the supersonic external flow past conic inlets at incidence, with the attached shock wave. The solid fines in the graphs correspond to the windward side, the dotted lines to the leeward one. For an inlet with semi-apex angle co -- 20 °, at incidence ~ -- 10°, and with various values of the free stream Mach number M . , the shock waves in the planes ~p --- 0 and tp = ~ are presented in Fig. 26. It is interesting that, at high values of M . , the shock wave on the leeward side is

I00

P . I . CHUSHKIN

located closer to the body than on the windward one. Such a phenomenon is observed also for sharp-nosed cones. Figure 27 shows the pressure distribution on the inlets with various angles ¢o, for M . = 4 and ~ = 5 °. Here, even within comparatively small M~-'3

/

/ / /

/ /

/~/

I

. --..-_--

.4 /6 .

I0

!

---.~.,o M~ : 3 Fxo. 26. S h o c k w a v e for an inlet, ~ -- 20 °, for Ip = 0 (solid line) a n d V = ~ (dotted line), at • -- 10 ° a n d various M a t h n u m b e r s Moo

°(

:

o.4

;~ ,

p 0

.

0

:

~'~

3

2

4

:so" i

~

i

6

8

i

[0

,2

x

FIG. 27. Pressure distribution on the w i n d w a r d side (solid line) a n d leeward side (dotted line), for three inlets, at M . . ----- 4 a n d ~ -- 5 °

T H I t ~ - D ~ I O N A L

$

~

N

I

C

101

FLOWS

distances from the nose, the pressure tends to its asymptotic value obtained for the sharp-nosed cone with the corresponding semi-apex angle oJ, this value being marked in the graph by small circles (these data are taken from the tables
f3

0 . 3 ' ~'

~3

tO" 5e

P

0.2

}.

.

.

.

O . . . .

--O~

0. I

I0" m

0

m

0.2

0.4

0.6

0.8

1.0

F~3. 28. Premme pmfilm on an inlet, oJ = 20°, for x m 14, y - 0 (solid line) and y :. :: (dotted line), at M , = 4 and various The pressure profiles across the shock layer are given in Fig. 28 for the inlet with oJ = 20 °, in the cross-section x = 14 (x is referred to the entrance radius of the inlet). The Mach number in this case was M . = 4, and the incidence ~, assumed the values -, = 0 °, 5°, 10° and 15°. In this graph, the small zircks correspond again to the data (3) for the sharp-nosed cone with the same w. The computational results given above show the effcctivness of the semicharacteristic scheme. A sul~,~ent accuracy is obtained here already with 5-7 mcridional planes and 20-25 mesh points on each radius. Because of the secondorder accuracy of this scheme, it is possible to use large steps with respect to x. In choosing thi~ step, one should also take into account the requirements of numerical stability. These considerations, together with great simplicity of the semi-characteristic scheme, makes it rather attractive. Let us note also that, for increasing efficiency in applying thk scheme, one should choose the coordinate system in accordance with the character of the flow to be computed. For instance, for a delta-shape wing with blunted edges, it is reasonable to apply elliptical coordinates when considering the appro][imations of the kind (7.2) with respect to the appro-

102

P. L CHUSHKIN

priate angular variable. In this case, it is possible also to take a combination of two coordinate systems, the Cartesian one on the fiat part of the wing, and the cylindrical one on the blunt part. 8. G E N E R A L I Z A T I O N FOR THE CASE OF E Q U I L I B R I U M GAS FLOWS At hypersonic Right speeds, intense heating of gas takes place behind the shock wave, and this causes various physico-chemioal processes in the gas, such as excitation of internal degrees of freedom of molecules, dissociation, ioniTation, radiation. Up to the temperatures of the order of 1000°K, the gas may be supposed to be perfect, and to have constant specific heats ratio. At higher temperatures, it becomes necessary to take into account the variation of specific heats. When the temperature exceeds 2000°K, the physico-cbemical changes in the gas begin to play a more and more substantial part, essentially affecting the flow field. The physico-chemical processes in a gas can occur at finite speeds, i.e. in non-equilibrium conditions. In this case, it is necessary to consider certain parameters character/z/rig the irreversible proccues, and to have differential equations defining the rate of change of these parameters in time. Under certain conditions, the physico-chemical reactions in a gas flow can occur at a very fast rate in comparison with the characteristic time, during which a gas particle remains in the given region of flow. One may then assume that, in every point of the flow field, the gas will be in the state of chemical and thermodynamical equilibrium although, strictly speaking, the complete equilibrium takes place only at infinitely high rates of reactions. In the case of the equilibrium flows, the thermodynamical state of the gas will be determined only by any two of thermodynamical parameters, for instance, by the pressure p and the temperature T, or by the pressure p and the enthalpy h. In the problems of equilibrium gas flows, it is also necessary (apart from the numerical integration of gas dynamic equations) to compute the thermodynamic functions for the mixture of reacting gases. It is possible here to apply several alternative techniques. One may solve the system of non-linear algebraic equations of chemical eqnifibrium directly. The equilibrium constants and the molar enthalpies must be represented preliminarily as functions of temperature. However, in this method, one must employ a laborious iterative procedure for solving the algebraic system.

~-DIMENSIONAL

SUPERSONIC FLOWS

103

Another technique, developed by Molodtsov~eS~ for calculating air flows, consists in reducing the original equations of chemical equilibrium to ordinary differential equations, by way of differentiation. In this case, it is necessary to know only the approximations of molar enthalpies. The system of equations obtained is integrated numerically together with the system of gas dynamic equations. This technique is more economical in comparison with the first one, but it still requires much time and has one disadvantage connected with the possible occurrence of trivial solutions. For multi-component gas media, such as air, the equations of chemical equilibrium are too compficated. Therefore, many authors have calculated, and presented in the tabular form, the thermodynamical properties of air in equilibrium state. Such tables contain various thermodynamical functions: pressure, temperature, density, enthalpy, entropy, internal energy, specific heats, speed of sound, concentrations of air components, etc. In these tables, there are two independent variables for which these or those thermodynamical quantities are chosen. In the U.S.A., tables of thermodynamical properties of the equilibrium air have been compiled, for instance, by Gilmore, c6° Hilsenrath e t al. ~6~ In the Soviet Union tables of thermodynamical functions of air, in the range of pressures from 0~01 to I000 arms, and for the temperatures up to 20,000°K, have been calculated by Predvoditelev e t al.
(8.1)

and the expression for enthalpy as h = ~,

T).

(8.2)

Hausen (~s) has constructed such approximations for the large range of pressures and temperatures. The errors of his approximate formulae are up to 4.5~ and, besides, they incorporate many exponential functions which increase the computational time appreciably. Similar approxima-

104

P. I. CHUSHKIN

tions were given also by Mikhailov, (~1' 7~) but their relative errors reach 10~o, as compared with the tables. (") Bade (vs) obtained approximate expressions for the density as functions of pressure and enthalpy, which are valid only for the pressure interval from 0.1 to 100 atms, and lead to errors up to 7-8~o. Two simple analytical equivalents of the tables/u) in the range of temperatures from 1000°K to 20,000°K, have been proposed by Kvashnina and Korobeinikov. (v° Their formulae express the enthalpy in terms of entropy and speed of sound, and are especially convenient for computating isentropic flows. These approximations have a comparatively small average error and their maximal errors at individual points do not exceed 5-7~o. The most accurate approximations of the type (8.1) and (8.2) for the thermodynamical functions of equilibrium air, on the basis of the tables, (ts) were developed by Naumova. (~s) These analytical equivalents in which an additional approximation of Mikhailov's results (v~)is used, are obtained by double series expansion in terms of Legendre polynomials. The approxi= mations by Naumova cover the pressure range from 0.001 to 1000 atms, and the range of temperatures from 2000° to 16,800°K is divided into three sub-ranges. These approximating expressions have the maximum error of only 0-084~ and have been often applied in the numerical calculations of the equilibrium gas flows. More simple but less accurate analytical relationships of the type (8.1)-(8.2) have been obtained by Kraiko. c)6) In the range of pressures from 0-001 to 1000 arms, and temperatures from 400 to 20,000°K, the enthalpy and density computed by these formulae differ from the tabular data by + 3 ~ and ± 1.5~, respectively. Recently, some authors have proposed approximations of a different type, based on the idea of choice of such thermodynamical functions which depend weakly on the variation of some special parameters. Such thermo= dynamical functions may then be represented by linear or parabolic interpolational formulae over a number of subintervals. This permits to save both computer time and memory substantially. In this case, using only 100-200 computer storage cells, it is possible to ensure maximum errors of approximation of about 1-2~o. The approximations of such type have been developed by D'yakonov. <~) He introduces a special thermodynamical function, viz. the effective isentropic exponent Ye = ~°a~/P • If the relation y, = y,(p, e) is available, then the equations for equilibrium gas flows with the basic functions V, p, o required will have the same form as for the perfect gas flows. However, the function y, calculated by the tables (") varies with the

THREE-DIMENSIONAL S U I P ~ I q l C

105

FLOWS

pressure and temperature in a very complicated and non-monotonic manner. D'yakonov considers the dependence of this quantity 7'o on the pressure p and a certain parameter q which is expressed in terms of pressure and density, and is taken for the air in the following form

where the dimensionality [P/0] = mS/se~ and p is the dimensionless pressure related to the pressure p , = 1.01325 bar. Figure 29 shows the varia1.4

P- I 0 0 0

I.$

/ ),o /

,i0. I °°'

V T 0

i

I

12

1~3. 29. Variation of the effectivespecific heat ratio y,. as function of p and q for the equilibrium air tion of the parameter 7,, as a function of the variables q and p. As seen, the depelidence on q is sat~'actory and on p rather weak and monotonic. These circumstances allow the function y, to be approximated with respect to q by quadratic parabolae over several separate intervals, and linearly with respect to log p. It turns out that, for the air, the enthalpy, the molecular weights, the electron concentration and some other quantities depend on the variables p and q also in a satisfactory way, and they are easily approximated within the whole range of the tables (st) with relative errors about 1 - 2 ~ . D'yakonov ~i) has used his approximations for calculating (by the finite difference method) the three-dimensional equilibrium flow past blunted bodies. Sinchenko (see lul) has also proposed analytical representations of thermodynamical functions, based on introducing special variables Z -~ P

and

O = ~h- .

(8.3)

106

P . I . CHUSHKIN

These quantities are considered as functions of the independent variables ~)=lnp

and

H=lnh,

which have already been employed in Section 6. The equations Z = Z(~), H),

=

(8.4)

H)

(8.5)

ate the equation of state and the expression for enthalpy, similar to the corresponding equations (8.1) and (8.2). For representing approximately the relations (8.4) and (8.5) in the case of equilibrium air, the following analytical formulae are suggested: Z = 90+f.Q1, • = ( Z - I) (~o+/W~),

(8.6)

where ~90, f21, ~o, ~1 are functions of the argument H - H., and 4

H . = 7. cJJ-lnrh, j-0

f = ~)+In I0 Kp,

Kp and Kh ate certain coefficients which depend on the parameters of undisturbed stream V. and ~.. The range of variation of the argument H--H., covmed by the tables, <6a~ is divided into five sub-ranges, in each of which the functions ~o, f2,, ~0, ~1 are represented by polynomials. The coefficients c~ in the approximations (8.6) are taken in accordance with the data from the tables above mentioned. When solving numerically the problem of equilibrium supersonic gas flow past a body at incidence, it is necessary to add to the basic system of gas d y m u ~ equations (3.1) the equation of state and the expression for e n t i ~ p y which can be taken, e.g. either in the form of (8.1)-(8.2) or in that of (8.4)--(8.5). In the case of equilibrium flows the boundary condition (3.5) on the body surface, and the relations (3.6) on the shock wave, remain the same as for a perfect gas. Here as before, the formulae (3.7), (3.9)-(3.12) for velocity components are vafid. However, instead of expressions (3.8) the following equations should be taken

t,

107

THREE-DIMENSIONAL SUPERSONIC FLOWS

where the equation of state and the expression for enthalpy must, of course, be taken into account. In deriving the equations for three-dimeusioual characteristics, the basic system of equations of gas dynamics (3.1) was considered without any assumptions about the gas being perfect. Therefore, for flows of a real gas in the state of thermodynamical equilibrium, the corresponding results of Section 3 will remain valid. In particular, for equilibrium flows, the compatibility relations (3.15) on characteristic planes, and the relations (3.19) and (3.20) along the stream lines, will hold. In the :light of the above, four numerical schemes of the three-dimensional method of characteristics, described in Sections 4-7 for perfect gas, are easily extended to the case of equilibrium gas flows. For the sake of brevity, we confine such a generalization only to the scheme using only bicharacteristic directions and to the semi-characteristic scheme. The scheme using only bicharacteristic directions was developed by Magomedov<'5~in such a way that it could be applied to equilibrium flows. Here, it is necessary to add only the equation of state (8.4) and the expression for the enthalpy (8.5"). Let us note that the function Q and the Mach number M, involved in the computational formulae, are easily expressed through the quantifies Z and • (8.3) and the Bernoulli equation (3.22). In calculating points on the shock wave, the functions ~7and ~ are found again from formulae (6.14). However, for computing Ap and h, we now must use the expressions (8.7) together with the approximations (8.6), and it is convenient to solve these equations by iterations with respect to density ~ = pZ/h. In the semi-characteristic scheme by Katskova and Chushkin c51~ for equilibrium gas flows, the equation of state and the expression for the enthalpy were taken in the form of (8.1)-(8.2). In this case, the temperature T will be added to the basic functions and, among the reference quantities which serve for introducing the dimensionless variables, the gas constant should be considered in addition. The temperature T is determined through other thermodynamical functions by the following differential relation: 1 ~h [ ~ h ' ~ -1 and we shall have the following formula for the speed of sound a: =

0

ap

"

(8.9)

Here the appropriate partial derivatives can be found from the equation&

108

1'. I. CHUSHKIN

(8.1) and (8.2). Hence, the expressions (8.1), (8.2), (8.8), (8.9) should replace the formulae (3.2) for the perfect gas. In this scheme, for the equilibrium case, the equation for temperature (8.8) represented in the finite difference form along the "stream line" is added to the system of equations (7.12)--(7.14), in calculating an inner point of the flow field. An analogous equation is taken, naturally, for computing a point on the body surface. In defining the flow parameters at a point on the shock wave, the formulae (8.7) are applied which, together with the equations (8.1) and (8.2) and for the given value of Vn., are solved by Newton's method and permit to determine the values of p, h, T and o. However, in order to calculate the entropy, the thermodynamical relation (8.8) is used again and is written in this case along the shock wave. It should be noted that analytical expressions of the type (8.4)-(8.5) are more convenient than the expressions (8. I)-(8.2) because they allow to disregard the temperature. A rearrangement of the semi-characteristic scheme for the equation of state in the form of p = ~p, h) presents no difficulties. Let us show now for illustration, some numerical results for the flow past a body at incidence in the case of the air dissociating in equilibrium conditions. Such a three-d/meusional supersonic flow was calculated by Magomedov ~ts~ by means of his scheme using only bicharacteristic directions. He analysed the flow past a cone with a spherical blunt nose and semi-apex angle oJ -- 9.5 ° moving at incidence ~ = 5 °, at altitude 60 km and at speed V.. -- 7500 m/sec (which corresponds to the Mach number M . ~ 23.5). The initial data in this case were obtained from the appropriate numerical solution for supersonic flow past a sphere which had been computed by Lunev, Parlor and Sinchenko. c~ Figure 30 shows the pressure distribution on the body surface as a function of the coordinate ~0, for several values of x. The pressure profiles across the shock layer for three meridional planes at x = 4 are plotted in Fig. 31. Here, the solid fines correspond to the equilibrium flow and the dotted lines to the perfect gas flow with the isentropic exponent ~, = 1.4 past the same blunt cone at M . = o~. Analogous enthalpy profiles are presented in Fig. 32 for x = 20. The above results make it possible to appreciate the effect of equilibrium dissociation on the flow past blunt cones at incidence at hypersonic flight speeds. The rapid change of enthalpy profiles near the body surface (the behaviour of the density and velocity components is analogous) confirms that the main effect of physico-chemical processes in the equilibrium flow takes place in vicinity of body surface.

SUPERSONIC F L O W S

THItEE-DIMENSIONAL

0

i

109

!

I

0.08

0.04

O0

35 0

J 90 ~, deq

45

155

180

F~o. 30. Variation of lpremure on blunt cone, 0) = 9.5 °, in a stream of air d i u o c i l t inig in equilibr/um conditions, at V _ = 7 ~ 0 m/see and ¢ = $o

I

oo~

'

/

I o.o~

i I~,.o

•

J ,

0.04

.

o.o~1

t

0.02

0

1,.,,/2

I

---'t

oo,

.' : .J

'

_:_

!

--

-~"

-

/

-

-:-

,

I

,,

0.2

0.4

0.6

~

t

~/

_ "

J

i

0.8

t.o

FIo. 31. Comparison of p r e ~ u ~ profiles for a blunt cone, oJ = 9-5% at ¢ -- ~,

x .= 4, in a flow of equilibrium air (solid line) and l~'r/'ect gas (dotted line)

II0

P. I. CHUSHKIN

I 0 30

~

.......

h 0.18

0.0~

--

0.2

94

0.6

08

,

Fro. 32. Comparison of enthalpy profiles for a blunt cone, to : 9.5 °, at • = 5°, x - = 20, in a flow of equilibrium air (solid line) a n d perfect gas (dotted line)

9. G E N E R A L I Z A T I O N OF NON-EQUILIBRIUM

F O R ~' T H E C A S E GAS FLOWS

In the non-equilibrium case, when physico-chemical processes in the gas flow occur at finite speeds, the calculation of the gas flow becon~es considerably more complicated. Here, unlike the equilibrium case, one cannot confine oneself to the application of thermodynamic tables or relevant analytical expressions depending on some two thermodynamic variables. In the non-equilibrium flow, it is necessary to consider in addition a system o f equations describing the behaviour of media relaxing in time. Let irreversible physico-chemical processes in a gas be characterized by m parameters ct (i = l, 2 . . . . . m). It can be assumed that the role of such parameters is played, for instance, by m a u concentrations of gas components, energies of internal degrees of freedom, etc. The thermodynamic state of the gas is determined by these m parameters, the pressure p and the temperature T of the translational degrees of freedom of a certain component. Thus, the equation of state of gas in non-equilibrium conditions will be o = oO~, T, c), (9.1)

THREE*DIMENSIONAL SUPERSONIC FLOWS

III

where for brevity the whole set of the parameters ci is designated by c. The expression for enthalpy can be represented in the similar form h = h(p, T, c).

(9.2)

Let us take the system of equations of gas dynamics (continuity, momentum and energy equations) in the present case in the form ~TQV = 0, 00r.xT)V + ~Tp = 0, oV.vh-V.~Tp

(9.3)

= O.

It is necessary to add to equations (9.1)-(9.3) some differential equations governing the variation of parameters c~ dc, _ F~p, T, c) dt

(i = 1, 2, . m), " "'

(9.4)

where F l is a known function. The actual form of these equations depends on the nature of the irreversible processes (chemical reactions, excitation of internal degrees of freedom, and so on). Let us deduce now the equations of characteristics and compatibility relations along them, for non-equilibrium three-dimensional supersonic gas flows. The equation of energy conservation from the system (9.3) can be transformed, by using equations (9.1), (9.2) and (9.4), in the following way: a2 de de -.62,

dt

dt

where f2=a 2 ~ /'

, hi,\

(9.5)

Here, the partial derivatives with respect to corresponding variables are designated by primes and subscripts. The quantity a is the "frozen" velocity of sound S, c I

which is expressed by the same formula (8.9), but now all derivatives involved are calculated from the relations (9.1) and (9.2). Further, using the transformed energy equation, it is not difficult to rewrite the continuity equation from the system (9.3) in the form 0a1~TV+V.vp = -f2.

112

p. L CHUSmr,IN

Hence, it is seen that the transformed continuity equation and energy equation in non-equilibrium conditions will differ from the corresponding equations for perfect gas or for the equilibrium case only by terms in the right-hand part. These equations, together with the remaining unchanged momentum equation, form a system which entirely determines velocity V, pressure p and density 0; therefore,the differentialequations of characteristics for non-equifibrium flows will be formally the same as for flows of a perfect gas or of a gas in equilibrium conditions. However, in the compatibility relation which holds along characteristics, some changes will take place in the non-equih'brium case, which are conditioned by the presence of the right-hand part in the transformed equations of continuity and energy. Turning to the work by Rusanov, ~ ) in which the characteristic properties of general systems of quasilinear partial differential equations had been investigated, we may obtain compatibility relation along characteristics for the present case of three-dimensional steady non-equilibrium gas flows. Here, after some transformations, we obtain an expression differing from (3.15) by an additional term, namely,

=

aAl.

,v- B1 ,p + Qo,

where sl and s~ represent, as before, two independent vectors lying in a characteristic plane with the external normal n, and quantities Aj and Bj are defined by the formula (3.16). In addition Q = n-(sl×s2),

and the quantity ~2 is found from the equality (9.5). The derivation of characteristic compatibility relation for unsteady three-dimensional flows of multi-component media reacting in nonequilibrium way, is available in another form in the Sauerwein's paper. ~4" The streamline in non-equilibrium flows with m parameters which determine irreversible processes is a characteristic of (In + 1) th order. As previously, the Bernoulli equation (3.21)will hold along a streamllrle but the entropy will now not remain constant there. For calculating temperature, the energy equation from the system (9.3) is applied, in which the enthalpy h is eliminated with the aid of expression (9.2), i.e.

dT+ ~.~.(h;_l) dp +--~--r, I .~"x h;, dc, -- 0.

(9.6)

Besides,m equations (9.4)serve as characteristicrelationsalong the streamfines, and these express the change of non-equilibrium parameters ct in the directionof these lines.

THR~oDIMIgNSIONAL ~ N I C

FLOWS

113

All the above-mentioned relations, which hold along characteristics and streamlines, together with the equations (9.1)-(9.2) and formula (8.9), are sufficient for cakulating unknown functions V, p, ~, T, h, c~ at the mesh points of the characteristic network inside the field of non-equilibrium gas flow. For defining the flow parameters at the boundary points, we must resort to the conventional conditiom, viz., the boundary condition on the body surface and the relations on the shock wave. It should be emphasized that, as long as the shock wave has zero thickness, the parameters ct characterizing irreversible processes remain constant when crossing it. Thus, the relations on the shock wave in a non-equilibrium flow will be exactly the same as for the equilibrium case. Keeping all this in mind, it is not difficult to generalize the various numerical schemes of the three-dimensional method of characteristics for calculating supersonic gas flows in the prmence of non-equlh~rium physico-chemical phenomena. However, in computing the supersonic non-equilibrium flows by the method of characteristics, some specific difficulties may arise due to the behaviour of equations (9.4) as regards the parameters c i. The point is that the fight-hand parts of these equations have the following structure:

F~p, T, c) = ~ ,

T, c)f~p, T, c).

In the general case, the right-hand parts of the equations (9.4) may represent the sum of t e ~ s of similar type. Here, the function cpl(p, T, c) is determined by the speed of the tth irreversible process. If the characteristic relaxation time is very large, i. e. the physico-chemical processes are frozen in gas, the functions ~l -" O. If, however, the characteristic relaxation time is very short, i.e. the proc esses in gas occur at very large speeds, then the functions ~l -" ~ . In the latter (equilibrium) case, it is natural to assume that fKp, T, ¢) = O. The characteristic relaxation time for various irreversible processes, taking place in the region of flow under consideration, may vary within rather large limits because, generally speaking, some of the parameters ct may approach the equilibrium state, while others, on the contrary, may tend to be "frozen". Therefore, when integrating numericaUy ordinary differential equations (9.4) along the streamlines, the maximum admissible step of integration should be determined from requirements of computational accuracy of the most rapid process. In such cases, even using numerical schemes of high order, it may be required to apply a very small step of integration along streamlines. However, still more serious complications, connected with the stability

114

P . I . CHUSHKIN

of a numerical solution, make appearance when approaching equilibrium conditions. These difficulties are due to a special behaviour of relaxation equations (9.4) near equilibrium state because they are, in fact, equations with small parameters accompanying the highest derivative. Here, the functions f~p, T, c) tend to zero and begin to depend strongly on variation of temperature T and parameters ct. As a consequence, introducing small errors in values of T and ct leads to a growth of error in integrating equations for ca, and this in turn incwases the errors in T. Thus, a numerical instability arises. This question is of such importance in calculating nonequilibrium flows that it will be discussed in more detail. For suppressing instability in using standard numerical integration schemes of explicit type for ordinary differential equations, an exceptionally small step is required. Close to equilibrium, the size of step, determined by stability requirements, can be lower by several orders of magnitude than that which could be chosen in the basis of sufficient accuracy in approximating the functions. Even in the case of applying the Rungv--Kutta numerical scheme, which provides the fourth-order accuracy, the necessary size of step may turn out so small that the numerical integration may become practically unrealistic. First attempts to overcome these computational difficulties near equilibrium state cannot be considered as satisfactory. Cleaver ~9) suggested to apply, for nun~rieal integration in this case, the Euler first-order scheme which, though ensuring stability, decreases greatly the accuracy of approximation. Pavlova (e°) and Eschenroeder, Boyer and Hall (81) lineariz~d the relaxation equations with respect to equilibrium solution, but this approach considerably complicates practical computations. A simple and efficient finite-difference scheme, for calculating gas flows close to equilibrium conditions, was developed by Katskova and Kraiko.(~, sa) Here, the differential equations (9.4) are integrated following a scbeme of second-order accuracy, the functionsf~p, T, c) being represented by two terms of a series in terms of that parameter ct which tends to the equilibrium value. Then, for calculating this parameter c t, it is possible to obtain an expression ensuring numerical stability near equilibrium. As a result, the step of numerical integration can be increased many times as compared with the standard scheme. This scheme was proposed initially for computating flows with one nonequilibrium parameter c~. Later, Galyun and Kraiko ~s~) investigated this scheme and extended it for the general case. They suggested, for all relaxation equations, to expand the functionsf~p, T, c) in series in terms of all non-equilibrium parameters c~ and also of temperature T. In this case,

THREE-DIMENSIONAL SUPEIL~)NIC FLOWS

115

the system of linear algebraic equations must be solved to find the unknown values of T and ct. Numerical stability in computating supersonic flows near equilibrium is then obtained, even if the step of integration is increased by several orders of magnitude. Another method of solving relaxation equations, for the conditions close to equilibrium, was proposed by Moretti. ~ss~This method is based on replacing these equations, at each step of integration, by the corresponding system of linear non-homogeneous differential equations. The solution of this system is written in a closed analytical form so that, in this case, no problem of providing numerical stability arises. Moretti checked his scheme in a simple aerodynamic one-dimensional problem, viz. calculating combustion of hydrogen in the air at constant pressure. In linearizing equations of chemical kinetics, the author used essentially the specific properties of the flow under consideration, assuming for simplicity that the speeds of forward and backward reactions, the logarithmic derivative of density, and the sum of concentrations of all components, are constant during each step of integration. In this problem, an admissible step size was 200 times greater than that in the Runge-Kutta method. Let us note, however, that, when finding an explicit solution of the system of linear differential equations in the Moretti scheme, it is unavoidable to compute complex eigenvalues and eigenvectors of a high-order matrix. With large number of kinetic equations, these operations are very laborious. Therefore, DeGroat and Abbett ~ss~suggested to approximate the exponential solutions of the system of the linear differential equations in the Moretti scheme by polynomials of second order with respect to independent variable, i.e. time. The coefficients of these polynomials are found from a system of linear algebraic equations. The computational experience has shown that the step size in such a scheme can be increased, say by one order as compared with the Runge--Kutta method. Treanor ~s?~devised a numerical scheme for integrating ordinary differential equations, which considerably improves the Runge-Kutta scheme. He accounts especially for the behaviour of solution under the conditions when the derivative is strongly affected by dependent variable. Indeed, such a case takes place for the relaxation equations in approaching to equilibrium. In this new scheme, the solution involves exponential functions at each step and, for small steps, the new formula of numerical integration reduces to the Runge--Kutta formula. In an example given by Treanor for calculating equations of chemical kinetics behind the strong shock wave in the air, it was possible to increase the step 25 times as compared with the Runge-Kutta scheme, without causing numerical instability.

116

P.I. Cl/USHKII~

One more method of stable numerical calculation of gas flows with physico-chemical processes has been developed by Kamzolov and Pirum o v : ss) These authors apply an implicit finite-difference scheme. A solution of the system of fin/te-difference equations obtained is found by Newton's method, the initial values for the unknown functions being determined by extrapolation using two points previously calculated. The computation of one-dimensional flow of a chem/cally reacting multicomponent medium in a nozzle shows a good effectiveness and stability of the method. U~g this or that numerical scheme, one can compute supersonic non-equilibrium gas flows. In the two-dimeusional case, a number of authors used the method of characteristics and obtained numerical solutions of various problems with non-equilibrium physico-chemical phenomena. As regards the case of three-dimensional flows of media reacting in a non-equilibrium way, the practical calculations are only starting. Each of the four numerical schemes of the three-dlmensional method of characteristics discussed previously may, in principle, be extended to the case of real gas in non-equilibrium state. Here, we shall confine ourselves to such a generalization for the semi-characteristic scheme described in Section 7.

In thi~ scbeme, for non-equilibrium flows, the differential equations of characteristics (7.5) and the compatibility relations (7.6) along them preserve their form, but the expression for O1 will involve an additional term conditioned by the presence of irreversible processes. The quantity O1 will n o w be

=--

T

-~

pa~ 8~

~-v +

.

.

The differential equations of "streamlines" (7.7) and the relations (7.8) and (7.9) along them will not change in general. However, instead of the equality (7.10) for the entropy, which is not considered in non-equilibrium case, the equation of type (9.6) should be taken, serving for computating the temperature T. In the semi-characteristic scheme the equation of such a type holds along the "streamlines" and can be represanted in the form

where

:

+ h-7,.Z h:,.

~DIMBNSIONAL

SUPERSONIC FLOWS

117

The relaxation equations (9.4), for parameters c~characterizing irreversible processes, are written along the "streamlines" as follows: dot+ Oi+l dx -- 0

(i = 1, 2 . . . . . m),

14

(9.8)

where Oi+l-- W,-F,,

W, = w 8c,

New derivatives with respect to ~p, which are contained in the quantifies Os and Os+t, are eliminated in the usual way by means of approximations (7.2). In computing an inner point in the region of non-equilibrium flow, the pressure p and the velocity components u, v, w are determined as previonsly by the formulae (7.12) and (7.14). The temperature T and the parameters ct are found from the equations (9.7) and (9.8), represented along the "streamlines" by a finite-difference scheme of the second order of accuracy. The last equations are applied also for computating a point on the body surface. The procedure of determining the functions at a point on a shock wave remains fully the same as in the equilibrium case. If an arbitrary parameter c t approaches the equilibrium conditions then, for ensuring convergence of iterations carried out in the given semi-characteristic scheme, the procedure proposed by Katskova and Kraiko Csz' sa~ can be applied. For computating the parameter ct from the corresponding equation (9.8) considered along the "streamline" 0-3, we then obtain the expression

c,o = c = - T

\

Us

where and the subscript * me~n~ that the function is taken for the values of the arguments p~ To, cM (i ~ / ) and c=. The formula is of the second order of accuracy and makes it possible to carry out calculations both near to and far from equilibrium conditions, whereas at qp, -- ~ this formula reduces to a corresponding expression for equih~rium flow. When required, it is po.~dble to use for calculations the general formulae obtained as a result of expanding in series not only with respect to c t but also with respect to other arguments as it was proposed by Galyun and Kraiko! st~

I 18

P.i. CHUSHKIN CONCLUSION

Coming to the end of this survey of applying the numerical method of characteristics to the calculation of the three-dimensional supersonic flows past bodies flying at incidence, one may make the following conclusions. At the present time, several efficient numerical schemes of the method have been developed and realized on electronic computers, which permit to calculate the supersonic steady flow past bodies of revolution or other bodies of a simple shape (for instance, a delta wing with blunted leading edges) with sufficient accuracy. Considerable experience has been accumulated for the use of these schemes, both in the case of a perfect gas and for a real gas in thermodynamical equilibrium state. The schemes of the three-dimensional method of characteristics, developed for calculating the external flows about bodies at incidence, are easily extended to the internal flows in nozzles and to the jet flows with supersonic velocities. The generalization of the three-dimensional method of characteristics for supersonic flows with non-equilibrium physicochemical processes may be considered as quite realistic, though practical computations in this case have only been started. It seems, however, that the practical application of the method to calculating three-dimensional supersonic flows past bodies of complicated form, when interference takes place and a series of shock waves arise inside the flow field, appears to he not probable at present. The technical possibilities of present-day electronic computers put also a limit to obtaining sufficiently accurate numerical solutions for unsteady three-dimensional problems of gas dynamics by the method of characteristics. AKNOWLEDGEMENTS

The author expresses his deep indebtness to Pro£ W. Fiszdon, D. Kiichemann and S. Neumark who have helped in publication of this survey. The author is very grateful to K. M. Magomedov, V. B. Minostsev and Yu. N. Podladchikov for presenting results of their investigations. REFERENCES I. O. M. l~t,oTsz~covs~ and P. I. Cxusmcn~, The numerical solution of problems in gas dynamics. In: Basic Developments in Fluid Dynamtcs (ed. M. HOLT), vo|. 1, pp. 11--126. Academic Prca, New York-London (1965). 2. K. I. B ~ [ o and G. P. VoeKl~mm~gR, A numct'ical method of csicu/atinj threedimensional flow imst bodies in supersonic gas stream. Zll. Vych. Mat. i Mat. Fiz. 1, 1051--60(1961).

THIUEE-DIM]ENBIONAL SUPERSONIC F L O W S

119

3. K. I. B ~ o , O. P. V o s r ~ , A. N. Lx'uu~ov and V. V. RUSANOV, Threedlmen~onal Flow about Smooth Bodies in Ideal Gas. Nauka, Moscow (1964). 4. Yu. N. D'Y~oNov, Three-dimensionalflow past blunt bodies in supersonic stream of perfect gas. Izv. Akad. Nauk $$SR, Mekhan. l Mashinostr. 150-3 (1964). 5. Yu. N. D'YAKONOV, Three-dimensional flow past blunt bodies allowing for equiiibrinm physico-chemical reactions. DokL Akad. Nauk S$$R 157, 822-5 (1964). 6. J. YON NrOMAZ~ and R. D. R.ICHTMYn, A method for the numerical calculations of hydrodynamic shocks. J. AppL Phya. 21,232-7 (1950). 7. P.D. LAx, Weak solution of non-linear hyperbolic equations and their numerical computation. Communs. Pure and AppL Math. 7, 159-93 (1954). 8. S.K. GoDtmov, Finite-difl'eaer~ce method for nmnerical computation of discontinuons solutions of equations of fluid dynamics. Mathem. Sb. 47(89), 271-306 (1959). 9. S. K. GoDtmov, A. V. ZAmLOOnqand G. P. l~oKot,ov, Difference scheme for twodimensional non-stationary problems of gas dynamics and calculation of the flow with a ~ shock wave. Zh. Vych. Mat. i Mat. Flz. 1, 1020--50 (1961). 10. I.O. BOnACtmVSKYand E. L. Rtmnq, A direct method for computation of nonequilibrium flows with detached shock waves. AIAA Journal 4, 600--7 (1966). 11. I. O. BOHACtmVSrYand R. E. MATlS, A direct method for calculation of the flow about an axisymrnetri¢ blunt body at angle of attack. AIAA Journal 4, 776-82 (1966). 12. F . H . H~u.ow, The particle-in-cell computing method for fluid dynamics. In: Fundameutai Methods in Hydrodynam/ca, (eds. B. ALD~X, S. FEn~ACH and M. ROT£mUmO). Vol. 3 of Methods in Computational Physica, pp. 319--43. Academic Preus, New York-London (1964). 13. M. W. EWANSand F. H. I'IZZLOW,Calculations of supersonic flow past an axially symmetric cylinder. J. Aeronaut. Sct. 25, 269-70 (1958). 14. M.W. EwAm and F. H. HAllOW, Calculation of unsteady supersonic flow past a circular cylinder. AILS Journal 29, 46--48 (1959). 15. A . A . DOtODm'rSYN, On a method of numerical solution of certain non-linear problems in aero-hydrodynamics. Proc. $rd All-Union Math. Confr., 1956, vol. 3, pp. 447-53. Akad. Nauk SSSR, Moscow (1958). 16. O. M. Bstozasio~ovsr, n and P. I. C s u s m t ~ , Numerical method of integral relations. Zh. Vych. Mat. t Mat. F/z. 2, 731-59 (1962). 17. A. N. Mmxn.ce, Calculation of flow over blunt bodies of revolution at incidence in supersonic gas flow. Zh. Vych. Mat. i Mat. Fiz. 4, 171-7 (1964). 18. O. M. Bstorsm~ovsr,n, A. B t n . m ~ v , M. M. ~ v , V. G. GRuvNrrsrai, V. K. Dusm~, V. F. IVANOV, YU. P. L t m ' ~ , F . D . PoPov, G. M. RYABINKOV, T. YA. TIMOFZZVA,A.I. Tot.vrYr,H, V. I'4. Fom~ and F. V. SHuo^zv, Flow past BiJmted J ~ i o s in Supersonic Gas Stream. Theoretical and Experimental Investltatiom. Vych. Tsentr. Aired. Nauk SSSR, Mmcow (1966). 19. O.N. I~TSKOVAand P.I. CmJsm¢.~, Three-dimensional supersonic flow about bodies. Paper presented to the VII-th Symposium on Advanced Problems and Methods in Fluid Dynamics. Jurata, Poland (1965) (to he published). 20. G. F. ~ and G. P. TmY~J~ov, Method of calculating three-dimensioual flow about bodies with detached shock wave. DokL Akad. Nauk $S$R 154, 1056-8 (19~). 21. G.P. TnqYAKOV,Investigation of three-dimensioual flow about blunted bodies. l~v. Akad. Nauk S$$R, Mckhan., 10--19 (1965). 22. G. F. TmJmm and Yu. M. ~ l m , Non-stationary supersonic flow about blunted bodies with detached shock wave. lzv. Akad. Nauk $$SR, Mekhan. Zhid. i~ 19-29 (1966). 23. R. J. SWXO~UtT,Hypersonic blunt-body flow fields at angle of attack. AIAA Journal 2, 115-17 (1964).

120

P . I . CHUSHK1N

24. D. R. I - I A J t ~ Some Practical Methods o f uaint Characteriatica in the Calculation o f Non-steady CompreaMble Flow. U.S. Atom/c Energy Comm. Rep. No. AECU2713 (1953). 25. P. D. LAx and B. Wm~'D~OFV,Difference schemes for hyperbolic equations with high order of a~na'acy. Communa. Pure andAppL Math. 17, 381-98 (1964). 26. C. ~ I n t e r f ~ e between wing and body at supersonic speeds-analysis by the mothod of c,~_~_~istics. J. Aero. Sci. 16, 411-34 (1949). 27. A. F~ltI, Elements o f Aerodynamics o f Supersonic Flows. Macmillan, New York (1949). 28. N. ContntN and C. L. DOLPH, The method of characteristics in the three-dimensional statiomu7 supersonic flow of a comprcssible gas. Proc. Symp. Appi. Math. 1, 55 (1949). 29. JL SAtnnt, Dre/d/mensinnale Probleme der Characteristikentheor/e partieller Differential-Gleichun~n. Z. aneew. Math. Mech. 30, 347-56 (1950). 30. V. V. S v c s ~ , Cakulation of pressure distribution on bodies of revolution at angle of attack in supersonic stream of gas. In: Collection o f Theoretical Works on Aerodynamics, pp. 127-39. Oborongiz, Moscow (1957). 31. V. V. RUSANOV, Calculation o f Supersonic Gas Flows by Method of Characteristics. Disse~. Cand. Tech. Sci., Joukowsky VVIA, Moscow (1951). 32. V.V. RUSXNOV, Characteristics of general equations of 8ns dynamics. Zh. Vych. Mat. l Mat. Fiz. 3, 508--27 (1963). 33. C. K. THOI~'Un-L, The Numerical Method o f Characteriatics for Hyperbolic Problems in Three Indepe~alcnt Variables. ARCR and M No. 2615 (1948). 34. A. F~Jtt, The method of characteristics. In: General Theory o f HlSh Speed Aerod y ~ m ~ a (ed. W. R. S l ~ u ) . Vol. IV of Hith Speed Aerodynamics and Jet Propulsion, Princeton University Press (1954). 35. M. HOLT, The method of characteristics for steady supersonic rotational flow in three dimensions. J. Fluld Mech. 1, 409-23 (1956). 36. G. BIUm~ and W. ~ Ein Chataktetmtikanverfahren fUr dreidimensionale instationiL~ GasstrOmungen. Z. Angew. Math. Phys. 9, 173-90 (1958). 37. D. S. Btrr~it, The numerical solution of hyperbolic systems of partial differential equations in three independent variables. Proc. Roy. Soc. A, 2.25, 232-52 (1960). 38. D . J . RICHXaX~ON, The solution of two-dimensional hydrodynamic equations by the method of characteristics. In: Fundamental Methods in Hydrodynamics (eds. B. AI.D~, S. Fmt~mACHand M. R ~ , ~ R O ) , Voi. 3 of Methods in Computational PhyMcs, pp. 295--318. Academic Press, New York-London (1964). 39. L. F. Fowm.l., Flow r"YeidAnalysia for l.£fting Re-entry Confi~urationa by the Method o f Charaeterlatics. IAS Paper, No. 61-208-1902 (1961). 40. Yu. H. Ih~LADCHn~OV, A contribution to the calculation of supersonic gas flows in space by method of c h a r a c t e r . DokL Akad. Nauk $SSR 163, 1092-5 (1965). 41. Yu. N. PODL~m::H~OV, Method of characteristi~ for calculating ~ i m c n s i o n a l supersonic ~ flows, lzv. Akad. Nauk SSSR, Mekhan. 3-12 (1965). 42. Z. D. Z~RYANOV and V.B. Mn,~TnV, Computational method for thr¢¢~ i ~ ' s o n i c ~ flow about bodies. Izv. Akad. Nauk SSSR, Mekhan. i MasMnostr., 20-24 (1964). 43. V.B. Mn,~s'rmv, Computational method for supersonic three-d/menaional flow past smooth bodies. Izv. Akad. Nauk $$8R, Mekimn. ZMd. I Gaza, 126-33 (1967). 44. K. M. MAOotm~v, Method of characterist~ for numerical calculation of threedimemional ~ flows. Zh. Vych. Mat. i Mat. Fiz. 6, 313-25 (1966). 45. K . M . MAOOMZDOV, Calculation of three-dimensional flow about blunt cones allowing for equilibrium physko-chemical processes. Izv. Akad. Nauk $SSR, Mekhan. Zhid. i Gaza, 130-7 (1967).

THREE-DIMENSIONAL SUPERSONIC FLOWS

121

46. K . M . MAootvmov, On calculating surfaces in three-dimensional methods of characteristics. Dokl. Akad. Nauk SSSR 171, 1297-1300 (1966). 47. H. SAtnmwlm~, The method of characteristics for the three-dimensional unsteady magnetofluid dynamics of a multi-component medium. J. Fluid Mech. 25, 17--41 (1966). 48. G. Moltwrrl, Three-dimensional supersonic flow computations. AIAA Journal 1, 2192-3 (1963). 49. G. M o a r r n and R.W. B ~ , Hypersonic flow field calculations of threedimensional and real gas flows. Ch. 8. AGARD. In: The High Temperature Aspects of Hypersonic Flow, pp. 149--62. Pergamon Press, London (1964). 50. G. MoRrrn, Some comments about three-dimensional supersonic flow computation, AIAA Journal 4, 1695 (1966). 51. O. N. KArSKOVAand P. I. Cm.rstw,~, Three-dimensional supersonic equilibrium flow of a gas about bodies at incidence. Zh. Vych. Mat. i Mat. Fiz. 5, 503-18 (1965). 52. O.N. KATSKOVAand P.I. CHvsmcm, Supersonic ducts at the angle of attack. In: RAeS Centenary and ICAS V Congress Proceedings, pp. 309-20, Pergamon Press, London (1967). 53. H. SAUeaW~N"and M. SU~t~L~N, Numerical stability of the three-dimensional method of cha~cteristi~. AIAA Journal 2, 387-9 (1964). 54. H. H m and D . C . I . ~ o t t , Numerical stability of hyperbolic equations in three independent variables. AIAA Journal 3, 1099-1103 (1965). 55. M. BUgNAT, Method of characteristics for quasi-linear equations of gasdynamic type. Bull. Acad. Polon. Scl., ScJrie Sci. Math., Astr., Phys. 10, 91-98 (1962). 56. M. BUaNAT,A. KtELBA,S~Srdand A. W^Kuucz, The method of characteristics for a multidimensional gas flow. Arch. Mech. Stosowanej 16, 179-87 (1964). 57. R. S^uee, Numerische Ermittlung dreidimensionaler Oberschalh~r0mungen ohne Symmetricannahmen. Wehrtechn. Monatsh. 63, 67-71 (1966). 58. E. A. BltON¢3,A Method o f Characteristics for Solution of Three-dimensional Supersonic Shock Layer Flows. General Electric TIS 64 SD271. 59. C.R. STaOM, The method of characteristics for three-dimensional steady and unsteady reacting gas flow. Ph.D. thesis, Univ. of Illinois (1965). 60. C. C. TSUNO,Study of three-dimensional supersonic flow problems by a numerical method based on the method of characteristics. Ph.D. thesis, Univ. of Illinois

(1960). 61. Z. KOI,XI., Tables of Supersonic Flow around Yawing Cones. Tech. Rep. No 3, M.LT. (1947). 62. Z. KoPAJ., Tables o f Supersonic Flow around Cones o f Large Yaw. Tech. Rep. No. 5, M.I.T. (1949). 63. V. N. M~ttan.ov, Calculation of three-dimensional rotational supersonic gas flow in vicinity of a curve, along which stream lines have angular points. Prikl. Mat. i Mekh. 27, 1083-9 (1963). 64. R. MICHEL and Duo~3-Vn,m-Huso. Flux de chaleur au bord d'attaque d'une Rile it forte fl/:che en hypersonique. ICAS I V Congress Proceeding, pp. 1015-39. Washington-London (1965). 65. V. K. MoLom'sov, A contribution to calculating equilibrium gas flows. In: Numerical Methods o f Solution o f Problems o f Mathematical Physics, pp. 207-13. Nauka,

Moscow (1966). 66. F. R. GILMOPJS,Equilibrium Composition and Thermodynamic Properties of Air to 24,000 °K. Rand Rep. RM-1543 (1955). 67. J. I-In_qnq~nu~l"tlet aL Tables o f Thermal Properties o f Gases. Bur. Standards, Cir. 564 (1955). 68. A.S. P l u ~ V O V ~ V , E.V. S ~ N X O , E.V. SAMUn.OV, I.P. STAKHANOV,

122

69. 70. 71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84. 85. 86. 87. 88.

1,. I. CHUSHKIN

A. S. PLam~NOV and I. B. R O ~ N S K I 1 , Tables o f Thermodynamic Functions o f Air. Akad. Nauk SSSR, Moscow, vol. 1 (1957), voi. 2 (1959), vol. 3 (1961). N. M. Kuzmia,sov, Thermodymsm/e Funetiom and Shock Adiabates for Air at High Teml~watues. Muhinostroenie, Moscow (1965). C. F. H A I q ~ , Approximations for the Thermodynamic and Transport Properties o f High-temperatwe Air. NASA TR R-50 (1959). V.V. MiKttan.ov, On analytical representation of thermodynamic functions of dimociati~ air. l&th. Sb. 28, 36--43 (1960). V. V. Mimiau~v, On analytical representation of thermodynamic functions of air. Inzh. Sb. 31,206-16 (1960). W. L. B ~ s , Simple analytical approximation to the equation of state of dissociated air. A R S Journal 29, 298-9 (1959). S.S. K v ~ A and V.P. Koaommmmov, Solution of some problems of air motion in pretence of dissociation and ionization. Izv. Akad. Nauk S$SR, Mekhon. i MaaMnostr. 2~ ~1 (1960). See also ARS Jourmd 31,997-1001 (1961). I. N. NAmmvA, Approximation to thermodynamic functions for air. Zh. Vych. Mat. 1 Mat. Fiz. 1, 295--300 (1961). A.N. ~o, Analytical representation of thermodynamic functions for air. l~Th. Zh. 4, 548-50 (1964). V.V. L t m ~ , V. G. PAVLOVand S. G. SINCI~NKO,Hypersonic flow past a sphere in air ~ t i n g in equilibrium conditions. Zh. Mat. i Mat. F/~. 6, 121-9 (1966). O. N. KA'rsKovA and P. I. CHustlm~, Supertonk flow about open-noted bodies at the ~ of attack. Izv. Akad. Nauk SSSR, Mek/um. Zhtd. t Gaza, 124--30 (1967). J. W. Cutavtm, The Two-d/mensiona/Flow o f an Ideal Diasoclatlnf Gas. The College of Aemnautim, Cmnfield, R. No. 126 (1959). L. M. PAVLOVA,Non-equilibrium gas flows close to equilibrium flows, lzv. Akad. Nattk S ~ R , Meklmn. 1 McuMamtr. 13-17 (1961). A.Q. ~ o m ~ l t , D. W. Bo~am and J. G. HALL, Non-equilibrium expansion of air with coupled chemical reactions. Phya. FlaMa $, 615-24 (1963). O . N . Ka'mcow, and A . N . K t a m o , Calculation of two-dinmmional and axisymmetrical tmpmmonic flows in the premnce of ha'eversible processes. Priki. Mekhan. i Tekhn. Fiz. 116-18 (1963). O. N. KATaltOVAand A. N. IO, UKO, Calculation o f Two-diraemional and Axisym. metrical Sup~rso~c Flows in the Presence o f Irreversible Processes. Vych. Tsentr Akad. Nauk SSSR, Moscow (1964). N.S. Gm.~rN and A. N. KRAmO, Calculation of non-equilibrium flows. Izv. AN $S$R, Mekhan. t Mashinostr. 41-47 (1964). G. M o u ' r n , A new technique for numerical analysis of non-equilibrium flow. AIAA Journal 3, 223-9 (1965). J. J. DmGitoAT and M. J. Amm'rr, A computation of one-dimensional combustion of methane. AIAA Journal 3, 381-3 (1965). C. E. Tmml~R, A method for the numerical integration of coupled first-order dtffet,tmdal eqnatiom with Meatly different time constants. Matha. Comput. 20, 39-45 (1960. V. N. ~ v and U . G . Pmut~v, Calculation of non-equilibrium flows in nozzke. Izv. A N SSSR, Mekhan. Zhid. i Gaza, 25--33 (1960.