- Email: [email protected]
On the convergence on a shock wave of continuous calculation difference schemes which are invariant with respect to similarity transformations
U.S.S.R.
Comput. Maths. Math.
Phys,
Vol. 26,
OOdl-5553186 $10.00+0.00 Q 1988 Pergamon Journals Ltd.
No. 6, pp. 39-50, 1986
Printed in Great Britain
ON THE CONVERGENCE ON A SHOCK WAVE OF CONTINUOUS CALCULATION DIFFERENCE SCHEMES WHICH ARE INVARIANT WITH RESPECT TO SIMILARITY TRANSFORMATIONS* V. V. OSTAPENKO The convergence of continuous calculation uniform difference schemes on discontinuous solutions which are invariant with respect to similarity transformations is investigated. A definition of the conservativeness of a difference scheme as a purely internal property of the scheme is given and it is shown that, by virtue of its definition, conservativeness is a sufficient and, subject to a certain general assumption, necessary condition for the convergence of an invariant scheme on a shock wave. A method is developed for the theoretical estimation of the unbalances of non-conservative schemes on a shock wave on the basis of which estimates are obtained of the unbalances which arise during the calculation of the problem of the decay of a discontinuity using the nonconservative schemes of gas dynamics in Eulerian coordinates. The property of conservativeness is one of the most important properties of finite-difference schemes which approximate conservative systems of non-linear differential equations [l-4]. When this property is not present, it would be expected that the limiting discontinuous solutions of a difference scheme would not be generalized solutions of the approximated system. For instance, an example was constructed in [S] and [6] where a non-conservative difference scheme for the heat conduction equation, which ensures second-order accuracy in a class of sufftciently smooth coefficients, diverges in a class of discontinuous coefficients and it is shown in [7] and [8] that the limiting discontinuous solutions of non-conservative invariant schemes of gas dynamics in Lagrange coordinates do not satisfy the law of conservation of the total energy. Conservative difference schemes are constructed using integro-interpolation, variation and projection methods on the basis of the corresponding integral conservation laws [l-4] and, also, by means of the direct standard approximation of the divergent form of representing the system [4]. When this is done, the scheme is considered to be conservative when certain difference conservation laws, which are analogues of the corresponding integral conservation laws [l-4], are satisfied. With such an approach the conservativeness of a scheme is a property which is defined in terms of the conservativeness of the system of differential equations which is approximated. At the same time, since a difference scheme may be considered as an independent mathematical object, there is considerable interest in another approach which treats the conservativeness as a purely internal property of this scheme which is not directly related to which system is approximated or how it is approximated. In the development of this approach in this paper, a formal definition of the divergence of a difference operator and the concept of the conservativeness of a difference scheme associated with it is given as a purely internal property of such schemes. It is shown that, by virtue of the given definition, the conservativeness of a difference scheme, approximating a conservative differential equation, which is invariant with respect to similarity transformations is a sufftcient and, subject to a certain general assumption, necessary condition for its convergence on a shock wave. A simple criterion is obtained for checking actual difference operators for the property of divergence. A canonical form of representing a divergent difference operator is obtained by drawing an analogy between linear difference operators and algebraic polynomials. The method which has been developed enables one to estimate the unbalances which arise in the calculation of a shock wave using non-conservative schemes. Using this method we have obtained estimates of the unbalances arising in the calculation of a shock wave using a non-conservative scheme for the non-linear transport equation and in the calculation of the simplest discontinuous flows using the two-layer non-conservative schemes of gas dynamics in Eulerian coordinates.
1. Divergences of finite-difference operators 1. By T,hir we denote a shift operator, the action of which on each functionf(v(x,r)) T&f(v(z,
t))=f(v(s+ih,
is defined by the identity
l+jr)).
where v = (u ,, . . .,uN)is a vector-function which depends on the continuously varying arguments x and t,fis a specified scalar function and h and r are the spatial and time steps. Subsequently, we shall also use the abbreviated notation I): everywhere in the case when the steps h and z for the operator T,P are fixed. We will denote by AFT%a non-linear operator which is specified on a uniform template located within the (n, and n2 are integers), the action of which on each function v&t) is defined by rectangle [-R~, nJ’x 1-h ntl the _ ~~‘“‘[vl=~~‘~~‘~~of.,(v) m-l *-I (1.1) = &“,fi
“l-1 lZh.vychisl.
Mat.
mat. Fiz., 26, 11, 1661-1678,
1986.
Imr(v (z+i&
k-1
39
t+jn,r) ) ,
where --n,G&&nr, -n,Gj&32,, m--l, 2,. . . , M, k--l, 2,. . . , Km; f,,,&are known functions of v, a, are constant coefficients which are dependent on It and r and .f*d. Definifion I. We shall say that the difference operator (1.1) is divergent and that the difference scheme corresponding to it
AT [VI=0
(1.2)
is conservative if: (1) for all integers N,, Nr, M, and M2 such that .N,-N,Xn,, x0, z,, the sum
M,-ill,~?&,,
and for any real numbers
Nr =lr s = c c ht:,n’ Iv (softh, t&r) i-I*‘*f-.X,
1
(1.3)
after the reduction of similar terms is independent of the values of the function v(x,?) from the rectangle n-{30+(N,+nl)h~z~~~~(Nz-n,)h, t,+(M,+~)rat~t,+(M,-na)tj; 2)
Ah:‘“‘[Cl-0 vc * coni3t. It follows from this definition that, by analogy with the case of a divergent differential operator, the integral of
which over a region is written in the form of a curvilinear integral along its boundary, only terms remain the values of which are taken from a closed band of finite width along the boundary of the domain, as a result of the summation of a divergent difference operator over a sufftciently large domain. Remurk 1. Summation of a difference equation over a net region, in order to check it for the property of conservativeness, is applied quite frequentiy [1,2,4]. Therefore, in the present paper, the property of “summabiiity” lies at the foundation of the definition of the divergence of a difference operator and the conservativeness of a difference scheme. We also note that definition 1 may be considered as a fo~~ization of the requirement mentioned above in the introduction that a conservative scheme should satisfy a certain difference conservation law. 2. Let us now obtain a simple criterion for checking actual difference operators (1. I) for the property of divergence. In order to do this we give the following definition. Definition2. We shall cali the operator K_
the difference monomial. We shall say that two difference monomials %[v] and @,,,,[v] are analogous if numbers a and b exist such that the identity @,[v]-T~~~U&[V] holds. This relationship, which lies in the fact that two monomials are analogous, is an equivalence relationship and, by virtue of this, it subdivides all monomials occurring in the operator (1.1) into non-intersecting groups of anaiogous monomials. Let there be 2, such groups and let M, monomials occur in each I-th group. Then, the operator (1.1) can be written in the form (1.4a)
-nlGih+ir,Gz,, Lemma I. The condition
-nrGjrm+jrrGnl.
(1.4b)
Y.
E arm=0
Vl=1,2,...,L
(1.5)
is necessary and suf%ient for the divergence of the difference operator (f .4). Proof. Let us write the sum (1.3) in the form
Only monomials which contain values of the function v from the rectangle II remain in this sum. After reduction and taking account of (1.4b), we obtain
It follows from this formula that the sum S, vanishing identically is equivalent to condition (1.5). The lemma is proved.
41
As an example of the use of this lemma, let us consider the difference operator A[p, u]=IL(““(pu)I-‘/*u’““u(~~Jp,, where f”‘hof(tfr)+(l-o)f(t),
f,-r-‘[f(t+r)--f(t)]
, which approximates the divergent differential operator
d(pu) 1 ~--_-_*~_~-
dt
(1.6)
2
dp
d
dt
dt
This difference operator arises in passing from the difference equation for the internal energy to the difference equation for the total energy in the case of gas-dynamic difference schemes in Eulerian coordinates. Let us rewrite the operator (1.6) in an expanded form: A[& u]=-r-‘[ (or--hr)p~~-(l--cit-ar)pu~-((o*--X~)~uri f(l-o;-l,)p;ia+~.*p8’-2hlp;tP]; X(1-o,), ht=‘:,[~l(l-a*)+u*(l-a~)] has been where the notation f-f(t), f=f(th), ~FYWJ& L= %(i-u,) introduced. This operator only contains a single pair of analogous monomials: ii’ and pu*, Hence, by equating the sum of their coefftcients to zero and also the coefficients accompanying all the remaining monomials, we obtain a system of equations for determining the parameters u,, i = 1,2,3 for which the difference operator (1.6) is divergent:
20,~lfL-h,-Q, ~*==0*0s==0,
.o,--hs=O, 2aP=(l-o,)
l-a,-h,=O,
(i-%)~O.
Solving this equation, we find that u, = 0.5, us = 0, CT)= 1 (or u2 = I, u3 = 0). At the same time, the operator (1.6) is divergent and can be represented in the divergent form:
_ &zf-pu’ . 25 3. We shall refer to the different operator A~(TI-E)OAI+(T’-E)oA~,
E-T,“,
(1.7)
where A, and A, are certain operators of the form of (1.1) as a canonical operator. It is obvious that a canonical difference operator is divergent by virtue of definition 1. It is found that the inverse assertion also holds. ZXeorem I. Each divergent difference operator (1.1) can be represented in the canonical form of (1.7). Proof. Applying Lemma 1, let us write the divergent difference operator (1.1) in the form
and ‘I+are difference monomials. Let us transform (1.8) in the following manner:
(1.9) where
are linear operators It follows from differential operator, now set an algebraic
which are divergent, respectively, with respect to the variables I and x. formula (1.9) that each two-dimensional divergent difference operator (l.l), just like a divergent may be written in the form of a sum of one-dimensional divergent difference operators. Let us polynomial L,-
&IT, (-0
with the same coefficients a, to correspond with each one-dimensional
linear difference operator
42
This correspondence is an isomorphism of a ring of linear difference operators (Lm) with the operations of addition and composition on a ring of polynomials [Pm) with the operations of addition and multiplication. Hence, the wellknown results on the expansion of polynomials are automatically transferred to linear difference operators. Let us formulate the results we require in the form of a lemma. Lemma 2. A divergent operator L, is uniquely representable in the form L,,,- (T,-E) a&,-,. The proof follows from the fact that the polynomial P,,, (the sum of its coefficients is equal to zero) which corresponds to a divergent operator L,,, has a root x0 = 1 and, hence, according to Bezout’s theorem, is uniquely representable in the form P, = (x-l)p,_,. By transforming the operators A, and A@occurring in (1.9) with the help of this lemma and carrying out some simple calculations, we arrive at formula (1.7). The theorem is proved. It follows from this theorem that formula (1.7) may also be taken as a definition of a divergent difference operator and is equivalent to definition 1 (such a definition is used, for example, in [9]). At the same time, the difference analogues of integration by parts hold in the case of a divergent operator written in the canonical form (1.7)
[11. 2. Conservation law< for invariant difference schemes for continuous calculation on a shock wave 1. Let us assume that the difference scheme (1.2) approximates conservative differential equation
(in the sense of the definition presented below) a
where cp,@mC’. Definition 3. The difference operator (1.1) approximates the differential operator F and the difference scheme (1.2) approximates the differential equation (2.1) if, for all sufficiently smooth functions u(x,t) limhr:‘“‘[U]~F[U]. &T-r0
(2.2)
We note that theapproximation of an equation is usually [2] understood as a less strong requirement that Eq. (2.2) should be satisfied only on the smooth solutions of this equation. Let us next assume that, in accordance with the principle of invariance [lo], the difference scheme (1.2) (taking account of the fact that the time step r in it is chosen in accordance with Courant’s stability condition as being proportional to the spatial step h) allows of the similarity transformation Z’=%.z,
t’v&t,
h/-ah.
(2.3)
Allowing for this, we shall introduce the abbreviated notation A&(2,
t) I-O.
(2.4)
for it. We note that practically all difference schemes for one-dimensional planar-symmetric problems in gas dynamics which are used in practice allow of such a similarity transformation [ 1I]. The following lemma is central in the study of the convergence of scheme (2.4) on a shock wave. Lemma 3. Let h, and v&) (where .$= x - 0 t, @ = const) exist such that &[v&)]-0 for all EEB, Then, when h > 0,AJvr(t)] -0 is satisfied in the case of the functions
for all &&I. The proof transformation It follows 0 corresponds
follows directly from the fact that scheme (2.4) allows of the similarity E’-=a’& h’-ah in the self-modelling variable {. from this lemma that a whole family (2.5) of solutions v*(t) of this scheme with different step sizes !I > to each self-modelling solution v&) of scheme (2.4). At the same time. if
lim ME) - vo, t-.+-D
lim vha(t)I--m
VI,
vo+vi,
(2.6)
then, by passing to the limit as h - 0 with respect to this family of solutions (2.5), we obtain the limiting discontinuous solution, that is, the “rack” v’- (vO,E>O; v,, E
h(~k-b~(u)ll.:~-o
(2.7)
of Eq. (2.1) in the case of this solution does not follow by any means from the approximation of this equation by means of scheme (2.4) on smooth functions in the sense of definition 3, since the limiting values of the functions (2.5) when E-t= are independent of the scheme step size h.
43
2. Let us now investigate under what conditions the limiting discontinuous solution v* satisfies the conservation law (2.7). In order to do this, we shall write the operator A* (which, in the general case, is non-divergent) in the form of the sum of a divergent operator b* which approximates the divergent differential operator F and a nondivergent operator A, which approximates zero in such a manner that, just like scheme (2.4), the scheme &‘[v] -0 allows of the similarity transformation (2.3). The following theorem provides an answer to the question which has been posed. llrebrcm 2. The discontinuous solution v*, which is obtained by passing to the limit as h - 0 with respect to the family (2.5), (2.6) of solutions v&) of scheme (2.4), satisfies the conservation law +(2.8) [~(u)-~~(u)]l;;~Z, Z-r afn(E)lf% _I
when Zis a finite integral which is invariant with respect to A. Proof. Using Theorem 1, let us rewrite the divergent operator A** in the form ~-h-‘(T~-E)ohn,+7-‘(T’-E)onN,
(2.9)
where the operators AI, and lzll approximate the functions cp and Q respectively. Let us now substitute the solution v#) into scheme (2.4) which has been transformed with the aid of (2.9) and integrate the resulting identity with respect to f from +oo to -00. As a result, we have -_ -*
Introducing the notation AfE]=tlsJvr(@
1, we find
&I‘[-N-h,
-NJ,
& =[N-h,N].
Let us show that (2.12)
and obtain the first of these formulae. From (2.6), it follows that
limA[EJ-Ih~[~_~_v~fEll=I\n[v,~. 6-w” Since &, approximates cp and v, = const, then, taking account of the fact that the scheme &[v]-0 is invariant with respect to the similarity transformation (2.3), we obtain &Jn]-cp(vr) for all h. The second formula of (2.1 I) is proved in an analogous manner. Using (2.11) and (2.12), we transform the first term on the left hand side of (2.10): -. -I
-
lim AI&l- llm AteI-cp(v,)-I. t-r-%*+*
Starting from the fact that a shift of - 0 7 with respect to the self-m~elling in the variable f, we obtain l? (T’-E)Bb,[v,(E)1~=-6(9h)-~(v*)). 7 -*
variable t corresponds
to a shift of r
(2.14)
in an analogous manner to that used previously. Subtituting (2.13) and (2.14) into (2. lo), we arrive at formula (2.8). Its left hand side is finite and independent of Ir. Hence, the integral Zon the right-hand side of this formula is also finite and invariant with respect to h. The theorem is proved. It follows from this theorem that, if the difference scheme (2.4) is conservative, Z = 0, and the conservation law (2.7) is satistied in the case of its limiting solution v*. Let us now assume that, if A, * 0, scheme (2.4) allows of some solutions of (2.5) and (2.6) for which Z # 0 {in particular, when I is finite). This assumption is quite general since it follows from the fact that the operator A, is nondivergent by virtue of definition 1 that the integral I (see (2.8)) depends on the values of the function v&) at all points on the numerical axis, that is, on its values for all &Q;(-0, +=) . Under this assumption the conse~ativeness of the difference scheme (2.4) is also a necessary condition for its limiting discontinuous solution v* to satisfy the conservation law (2.7) of the approximated equation (2.1). In particular, it has been shown in f7] and [S] that the
44
conservation law for the total energy is violated in the case of the limiting solutions v* of non-conservative gas dynamic schemes in Lagrangian coordinates. It is interesting to note that the non-divergent operator A, can have as high an order of smallness with respect to /I as may be desired on smooth functions and, in spite of this, as follows from Theorem 2, it will lead, in the general case, to the violation of the conservation laws for the limiting solutions v* as h - 0. Remark 2. The result obtained may be considered as a direct extension of the well-known results [S, 61 on the convergence of difference schemes in a class of discontinuous coefftcients to the case of the convergence of continuous calculation schemes, which are invariant with respect to similarity transformations, on discontinuous solutions of conservative hyperbolic systems. 3. Let us now consider in greater detail the question of the convergence, as h - 0, of the solutions (2.5) and (2.6) to the discontinuous solution v*. In order to do this, we introduce the following notation: @I) is a monotonically decreasing positive function such that Iims(h)==O,
lii[k/e(k)
b-.*
h-O
R--R\[-~81,
3-0, (2.15)
llfll~-=M~) L.1
1,
where R is the set of real numbers. The functions (2.5) and (2.6) converge, as h - 0, to the discontinuous function v* on the whole of the R axis in a certain integral norm. However, by far the greatest interest lies in the convergence of these functions in the norm C. It can be readily shown that, for all t > 0, these functions converge in the space C[RJ. However, they do not converge in the whole of the space C[R\O] in which the function v* is defined. The following lemmas show in what sense it is nevertheless feasible to speak of the convergence of these functions in the norm C on the set R\O. Lemmu 4. Let there exist he > 0 and k 2 0 such that Then:
Proof. Assertion (1) is obvious. Let us now prove assertion (2). From the condition of the lemma, it follows that there exists h, such that, for every ei > 0, there exists N(h,,e,) such that, for every e>N(h,, E,) , we have lu,,,(E)-u’l<.s,/]~l’. Using (2.5) and replacing the variable f = h f/h,, we obtain Iv~(%‘)-v’I<(k/k,)‘el/l%‘l*.
(2.16)
This inequality is satisfied subject to the condition that 1r / >hN( h.,, eJ/h,. It fohows from (2.15) that a S > 0 Whence, aRowing for (2.16), we obtain that, for each exists such that, for every h 4 6, we have e(h) >hN(h,, et)&. e2 = q’ht, a Sexists such that, for every R < S and every ffr/>s(h) , }~~(%‘)-v*l<[kle(k)]*e~. The lemma is proved. It follows from this iemma that the functions (2.5) and (2.6), which satisfy these conditions, converge as h - 0 to the function v* in the following cases: (1) for all e > 0, t = const, in the space C(RJ and the order of this convergence is greater than !I’. (2) in the single parameter family of spaces C[RJ which pass into the space C[R\O] as h - 0. In this sense one may speak of the convergence of these functions in a norm C on the set R\O. The order of this convergence is greater than [h/e(h) I”. The following lemma is proved in a similar manner. Lemmo 5. Let there exist ho > 0, i L I and k 2 0 such that the following are satisfied: 1) for all
PER there exists v~(E)=d’v&Q’;
2) lim %%l’ (%)-0. I-r Then: (1) for a11k > 0 and
2) 2
&ER , there exists vhio (g) such that Iii ply;‘) (g) -0; r-r-
~8k~b)iV-‘JUYA(f~c%,ila,,,,~o.
It follows from this lemma that, when k > i and h - 0, the functions vb”‘(n converge to a function which is identically equal to zero in the following cases: (1) for all t > 0 in the space C[R,] and the order of this convergence is greater than hX-‘. (2) in the single-parameter family of spaces C[Rdr,] subject to the condition that the quantity /P/e’(h) is bounded as h - 0. The order of this convergence is greater than k”-‘le*(k). As an example of the use of Lemmas 4 and 5 we shall cite the solutions obtained in [4, p.631 of the differential equations of gas dynamics in Lagrangian coordinates with a linear artificial viscosity which, after the formal replacement of the coefftcient of viscosity Y by h, assume the form of (2.5) and (2.6). It can be shown that these solutions satisfy the conditions of Lemmas 4 and 5 for all k L 0 and i 1: 1. These solutions as we11as their derivatives
45
of all orders therefore converge, as h - 0, in the family of spaces C[R,,,] and the order of this convergence is equal to h’ (greater than hk for all k). In concluding this section we note the following important fact. If a certain solution v,(C) from the family of solutions (2.5) and (2.6) has oscillations on the shock wave front (in the neighbourhood of the point 5 = 0) then all the solutions from this family have oscillations of the same amplitude on the wave front. This is associated with the fact that the solutions ~~(6)when h = h, are obtained from the solution u&) as a result of compression of the latter by a factor of h,,/h along the f axis towards the origin of coordinates 6 = 0. As h - 0, all oscillation of the solutions (2.5) and (2.6) are thereby contracted towards the origin of coordinates while their amplitude remains constant. This means that, in the calculation of a shock wave using the corresponding difference scheme, it is impossible to attain any substantial reduction in the amplitude of the oscillations of the numerical solution solely by means of a reduction in the step size of the scheme h.
3. A method for the theoretical estimation of the unbalances of a non-conservative scheme on a shock wave 1. The violation of the conservation law (2.7) of the approximated equation (2.1) on the shock wave in the case of the non-conservative difference scheme (2.4) leads to the development of unbalance 6v,-vl-; in the determination of the variables behind the wave front and 6LB--6-a in the determination of the velocity of propagation of the wave front (; and @ are quantities associated with the exact relationship (2.7)). Let us assume that Z, Se!, and 6v, are quantities of the same order of magnitude which are small compared with v, and @. By varying Eq. (2.8), allowing for this, we obtain
This relationship, which is written in the form 6v &iv,(M, I), is conveniently considered as an unusual Hugoniot condition for the unbalances on a shock wave. For the approximate determination of the integral Z, let us suppose that Eq. (2.1) occurs in the system of equations
of the form of (2.1) which is closed, possibly by taking account of certain algebraic relationships G(u)--0
(3.3)
of the type of equations of state in gas dynamics, while the difference equation (2.4) occurs in the difference scheme which approximates this system and allows of the similarity transformation (2.3). With the help of the first differential approximation of this scheme [lo], we find the viscosity w which is characteristic for it which, in the general case, is written in the self-modelling variable in the form
where 6 is a certain sufficiently smooth function of two arguments such that e(u, 0) -0. By introducing this viscosity as an artificial viscosity into system (3.2) written in terms of the self-modelling variable f and integrating it after this with respect to 6 from +w to the current value off with the boundary conditions
lh
E-+0
u 6) -vo,
lhl
lim - = 0, re+r dt
(3.4)
we obtain
f&h;) Then, since a k L 1 exists such that the operator Z will be approximately defined by the formula
-Idu)-mmll.:. An/Z? approximates
(3.5) a certain differential operatorf,
the integral
where ;I(0 is the solution of (3.5),(3.3) with the boundary conditions (3.4) which correctly transmits the profile of a smeared out running shock wave [4,p.62]. Further simplifications are possible if the system (3.5),(3.3) taking (3.4) into account, like the analogous system of equations in gas dynamics [4], is solvable in the form
u,-u,(un,93),
t=4,2,.
(3.7)
. . ,N-i,
where, for all @ > 0, the function ~(uN, S?) is continuous and does not vanish in the interval (USA,%;I).The integral (3.6) can then be rewritten approximately in the form (3.8)
where the functionf, is obtained fromfby replacing all the fun~ions u,(z) and their derivatives in the latter by functions of uHand @ with the help of (3.7). Formula (3.8) will subsequently be the basis for determining the unbaiances. Remark 2. The method described above is based on a purely heuristic assumption concerning the “closeness” of the solution $0 of system (3.Q (3.3) to the corresponding solution of the difference scheme. Nevertheless, the unbalances in actual differences schemes of gas dynamics are determined using this method with a sufficiently high accuracy [7,8,12]. 2. As an example let us determine the unbalances which appear during the calculation of a shock wave using a difference scheme which is non-conservative when g # 0.5 Vt+V,,,t%=vVS, where Y-v(r,
t), YC==I-*(T~-E)*V, u,-~-~_‘(T,--E)~u; q=~-,*u., r--Ah,
which approximates
v-&h,
(3.9) u~~,-@E+(~-~)Z’_,~V,
A, c-const,
(3.10)
the non-linear transport equation
~+u$-0. Allowing for (3.10), this scheme admits of the similarity transformation Rewriting scheme (3.9) in the form -vvs+v,+(u’/2)r=-h(p-0.5)
and applying Theorem 2 to it, we obtain that its limiting ~~continuous relationship
(3.11)
(2.3).
(3.12)
(I&)*
solution U* (when u, = 0) satisfies the
At the same time, since here there is a fixed sign function, which is identically different from zero, under the integral sign, 1st 0 when j3 # 0.5. This means that the non-conservative scheme (3.9) on the shock wave does not satisfy the conservation law of equation (3.11). With the help of the first differential approximation we find that the characteristic viscosity in the case of scheme (3.12), which is the sum of an artificial and an approximate viscosity, has the form
By writing the viscosity in the self-modelling variable and transforming it after this taking account of (3.10) and also of the fact that, by virtue of the differential approximation of scheme (3.9), the equality
holds, we obtain
By introducing this viscosity as an artificial viscosity into the right hand side of Eq. (3.11) written in terms of the variable t and integrating it after this with the boundary conditions (3.4), WCobtain
au
-P
4
1
--u(rn-u). 2&h
Now, by applying formula (3.8), we find I=_
[email protected], 2Cih
5
&se
U (2s -z&u)&=
2 fP- 0.5)P 3Cl
Whence, taking account of (3. l), we have 2(~-0.5)SY ?r
8v,--2&B-
(3.13)
*
The quantity @ on the right-hand side of this formula is replaced by the quantity 3 which is close to it. 3. As an example of the use of this formula let us consider the numerical solution _ofthe problem of the propagation of shock waves which are specified by two different boundary conditions: ; =I 2 * @ = 0.5 (Fig. I, lines l-3) and G = 1 * 0 = 0.5 (lines 4-6) with respect to a uniform background v,, = 0. In this case, since u, = i and, consequently, au, = 0, we obtain from (3.13) the formula for determining 6 @ : (3.14)
32.6
3625
w
w.t5
,5p
3x1
.z
FIG. 1 b-0.1,S-4& linesZ-3- for A- b.2,C-0.5, T-40. 4-S- for A= 0.4,T-80; lines I-when@dt.~, (conscrvativc scheme), j 2 - when @i, 3 - when B-0, 4 -when b-s/,, C-V,, 5 -when B-S/s, C-.x0, 6 -when fi--I/,, C==‘/,. The unbalance (3.14) in the velocity of propagation
of the shock wave leads to the emergence of an unbalance 85~S-S--693T
(3.15)
in the distance traversed by the shock wave after a time ?‘(where S-62’ is the distance traversed by the wave front during the numerical calculation and S-&’ is the distance traversed by the shock wave when there is no unbalance). TABLE 1
The results of a comparison of the theoretical unbalances Ss,, mkuiated using formulae (3.14) and (3.15). and the unbalances &Sc which arise during the numerical caiculation are shown in Fig. 1 and Table I. Retnurk 4. In Fig. 1 the configuration of the shock wave front is shown by means of a vertical line. The broken lines are the positions of the fronts computed taking account of the theoretically calculated unbalance SS, (formulae (3.14) and (3.15)). The dotted and dashed lines are the shock waves which are obtained in the numerical calculation using scheme (3.9) with the corresponding parameters. The coordinate of the front of the wave in the numerical calculation is determined by the point of intersection of the corresponding graph with the line p = 0. The quantity A’S is defined bv the formula
4. The unbalances of non-conservative gas-dynamic schemes in Eulerian coordinates 1. Let us consider the system of differential equations of gas dynamics in Eulerian coordinates: (4.la)
9=-P&
e-Pe,
T-_(r-lke
(4.lb)
48
with a linear artificial viscosity w--vpi3u/as.
g=pfw,
(4.2)
By integrating this system in the self-modelling variabie t with the boundary conditions (3.2), we obtain, after transformations and allowing for the fact that r+,= 0 [4],
where
We approximate system (4.1),(4.2) by the differential difference scheme (4.4a)
(4.4b) Taking account of example (1.6), we obtain from this scheme the difference equation for the total energy:
where the unbalance of the Cst order with respect to T, which vanishes with the weights has the form 6E,=-r
[ (o,-a,)gl
G
5,=5,,
ape;,
5,-5,~0.5,
+ (0&.5)g@l’ $
1, while the unbalance of the order of rr has the form &E, - _$F
(&)‘.
By applying Lemma 1, it can be shown that the sum 6E, + 6E, is non-divergent for any choice of the time weights u,. The following theorem follows from this when account is taken of the fact that the divergence of its differentialdifference analogues is a necessary condition for the divergence of a difference scheme. Rkeorem 3. Among the two-layer difference schemes for which the schemes (4.4) are the differential~iffe~n~ analogues, there is not one scheme which is conservative. It follows from Theorem 3 that, among the above-mentioned schemes, there is also not one which is completely conservative. This result can be considered as a generalization of the analogous result in [9] to the case of an arbitrary spatial approximation. 2. Let us approximate system (4.1) by a difference scheme for which scheme (4.4) is the differential-difference analogue:
et--
(4Sb)
(~(51310,1)t_g(.,)~~~,
where f-Z’_& &=Tsa~a, g,--(2h)-‘( T,-T-J og. In this difference scheme not only the time approximation is non-conservative but also the spatial approximation (that is, the differential with respect to time and the spatial difference analogue of thii scheme is also nonconservative). A scheme with such a spatial approximation arises, for example, in the Euler direction of the combined Euler-Lagrange method proposed in [ 131. Under the conditions (3.10), the difference scheme (4.5) is invariant with respect to the similarity transformation (2.3) and, by virtue of this, Theorem 1 is true for this scheme. By transforming formulas (3.1), which correspond to scheme (4.5), we obtain
47P* 6% + (r+,)a’m
- -
By applying formula (3.6), we find, after some transformations,
2(7-f) D’(7+-i) (V-9) that
(4.6b) z
49
(4.7)
I=I,-kls,
where
(4.8a)
h--o,-o~-o,+0.5, +a
xP=u‘--o*,
h,-u‘-o,-u,+0.5,
(4.8b)
We shall refer to the quantity L, as the unbalance in the time approximation of the difference scheme (4.5) since this quantity is completely defined by the differential-difference scheme (4.4) which corresponds to this scheme. Similarly, we shall refer to the quantity Zs as the unbalance in the spatial approximation since it is completely determined by the differential with respect to time and spatial difference analogue of the difference scheme (4.5). It thereby follows from (4.7) that, in the approximation of the method being considered, the total unbalance of the difference scheme is equal to the sum of the unbalances of its time and spatial approximation. By supposing now that its artificial viscosity (4.2) is the characteristic viscosity in the case of scheme (4.5), we obtain from (4.8), making use of (4.3) and (3.8), that n zT_
(r+i)Tn% j 2v
- - (Yy;vrTnO
a(n) (Q-‘l)
(?-Pt)drl
;h,l,+LEi.+h,l,),
where
(4.10) 3. As an example of the use of the formulae which have been obtained, we shall give a calculation, using an explicit scheme (4.5) with the weights a ,-c~~=a~==cr~=u~=O,~,==-a,-I, u, -0.5*~,--lr-0, h-,-O.5 , of the problem of the decay of a discontinuity which arises when a diaphragm which separates cold gas (X) from a warmer gas (H) is instantaneously removed, where nOx-nOn-l, pox-O, p,,%- 2, r-2. When this is done, a shock wave will propagate through the cold gas with a velocity d = 1.13755 while a rarefaction wave will propagate through the warmer gas. The unbalances Sf,‘, 6D on the shock wave are related by Eqs. (4.6) where, when account is taken of (4.7), (4.9), and (4. lO), z=
-
(r+OD 24v
(h-r@)
[ (2no’+Snof)-4’) (no-q) -6~9
inn+] .
The unbalances on the rarefaction wave are related by the equations [ 121
where
c,=(yp,qo)”
is the initial velocity of sound. They are obtained as a result of a variation of the equations
which the gas-dynamic parameters on the rarefaction wave satisfy [3, p.1881. From this, by equating the unbalances of the pressures and velocities of the gases on the boundary of contact, we obtain a closed system for determining the unbalances 6~~-6u,~=6u~~, 6ps=6p+8p,‘, 6D, 6qix, 6qt’. After these have been found, the unbalances in the internal energy are obtained from the equation aa,-( ~-1)-1(@6~,+~6p,),which follows from the equation of state (4. lb). USER26:6-D
TABLE 2
1.1
0.0147 -0.0158
:i
0:s 9.4 0.8
-Fzz -0:0380 -0.Oli8
The results of a comparison of the theoretically computed unbalances with the unbalances which arise in the numerical calculation (h-0.1, z-0.02, v-0.05) are shown in Table 2, where
The author thanks S. M. Shugrin for useful remarks. Tronsiotedby E. L. S.
REFERENCES 1. SAMARSKII A. A., Theory ofDifference Schemes (Teoriya raznostnykh skhcm), Nauka. Moscow, 1983. 2. GODUNOV S. K. and RYABEN’KII V. S., Dijference Schemes (Raznosmye Skhemy), Nauka, Moscow, 1973. 3. ROWDESTVENSKII B. L. and YANENKO N. N., Systems of Quasi-Lineor Equations (Sistcmy kvasilineinykh uravnenii). Nauka, Moscow, 1978. 4. SAMARSKII A. A. and POPOV YU. P., Difference Methodsfor Solving Problems in Gas Dynamics (Raznwtnyc metody resheniya aadach gazovoi dinamiki), Nauka, Moscow, 1980. 5. TIKHONOV A. N. and SAMARSKII A. A., On the convergence of difference schemes in a class of discontinuous coefficients, Dokl. Akod. Nauk SSSR, 124.3, X9-532.1959. 6. TIKHONOV A. N. and SAMARSKII A. A., On uniform difference schemes, Zh. vychisl. Mot. mut. Fiz., 1, I, 5-63, 1961. 7. OSTAPENKO V. V. and SAPOZHNIKOV G. A., On non-conservative difference schemes of gas dynamics, in: Numerical Mefhods of the Mechunics of u Conrinvous Medium (Chisl. metody mekhan. sploshnoi sredy), Vol. 14, No. 2, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1983, pp. 103-l 13. 8. OSTAPENKO V. V., Estimates of the unbalances in the gas-dynamic parameters which arise in the calculation of a running shock wave using non-conservative difference schemes, in: Numerical Methods of the Mechanics of o Conrinuws Me&m (Chisl. metody mckhan. sploshnoi sredy), 14,6, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1983, 130-W. 9. KUZ’MIN A. V. andMAKAROV V. L., On an algorithm for constructing completely conservative dilfcrcncc schemes. Zh. vychisl. Mar. mot. Fir, 22, 1. 123-132, 1982. 10. SHOKIN YU. I., The Method of Diff erential Approximation (Meted differcntsial’nogo priblizheniya), Nauka, Novosibirsk, 1979. 11. YANENKO N. N. ET AL., Classification of the difference schemes of one-dimensional gas dynamics by the method of differential approximations, in Numerical Methods of the Mechanics of u Continuous Medium (Chisl. mctody mekhaniki sploshnoi sredy). Vol. 12, No. 2, Inst. Twr. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1980. pp. 123-159. 12. OSTAPENKO V. V., Estimates of the unbalances in the gas-dynamic parameters which arise in the calculation of the simplest discontinuous solutions using non-conservative difference schemes, in: Numerical Me:hodr of the Mechanb of o Continuotu Medium (Chisl. mctcdy mckhaniki sploshnoi srcdy), Vol. 15, No. 6, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1984, 110-124. 13. FRANK R. M. and LAZARUS R. B., A mixed method using the variable of Euler and Lagrange, in: Computationof Merhods in Hydrodynamics (Vychisl. melody v gidrodinamikc), Mir, Moscow, 1967.55-75.
U.S.S.R.
Cornput. Maths
Mafh.
Phys.. Vol. 26, No. 6,
Printedin Great Britain
pp. 50-59, 1986
0041-5553/86 E10.00+0.00 0 1988 Pergamon Journals Ltd.
A METHOD OF COMPUTING IDEAL GAS FLOWS IN PLANE AND AXISYMMETRIC NOZZLES WITH CONTOUR BREAKS* A. N. KRAIKO, N. I. TILLYAEVA and S. A. SHCHERBAKOV A method for computing inviscid and thermally non-conducting gas flow in a smooth or unsmooth nozzle, including one with a sudden contraction, is described. It combines a difference scheme that preserves the approximation on any mesh (in the case of an unsmooth generator its lines condense into a “sonic” or “supersonic” break), the method of characteristics, and analytic solutions, valid in the neighbourhood of the break and on the axis of symmetry. When using the method of characteristics to compute the supersonic flow
lZh.
vychisl. Mat.
mat Fiz.. 26,
11, 1679-1694, 1986.
Comput. Maths. Math.
Phys,
Vol. 26,
OOdl-5553186 $10.00+0.00 Q 1988 Pergamon Journals Ltd.
No. 6, pp. 39-50, 1986
Printed in Great Britain
ON THE CONVERGENCE ON A SHOCK WAVE OF CONTINUOUS CALCULATION DIFFERENCE SCHEMES WHICH ARE INVARIANT WITH RESPECT TO SIMILARITY TRANSFORMATIONS* V. V. OSTAPENKO The convergence of continuous calculation uniform difference schemes on discontinuous solutions which are invariant with respect to similarity transformations is investigated. A definition of the conservativeness of a difference scheme as a purely internal property of the scheme is given and it is shown that, by virtue of its definition, conservativeness is a sufficient and, subject to a certain general assumption, necessary condition for the convergence of an invariant scheme on a shock wave. A method is developed for the theoretical estimation of the unbalances of non-conservative schemes on a shock wave on the basis of which estimates are obtained of the unbalances which arise during the calculation of the problem of the decay of a discontinuity using the nonconservative schemes of gas dynamics in Eulerian coordinates. The property of conservativeness is one of the most important properties of finite-difference schemes which approximate conservative systems of non-linear differential equations [l-4]. When this property is not present, it would be expected that the limiting discontinuous solutions of a difference scheme would not be generalized solutions of the approximated system. For instance, an example was constructed in [S] and [6] where a non-conservative difference scheme for the heat conduction equation, which ensures second-order accuracy in a class of sufftciently smooth coefficients, diverges in a class of discontinuous coefficients and it is shown in [7] and [8] that the limiting discontinuous solutions of non-conservative invariant schemes of gas dynamics in Lagrange coordinates do not satisfy the law of conservation of the total energy. Conservative difference schemes are constructed using integro-interpolation, variation and projection methods on the basis of the corresponding integral conservation laws [l-4] and, also, by means of the direct standard approximation of the divergent form of representing the system [4]. When this is done, the scheme is considered to be conservative when certain difference conservation laws, which are analogues of the corresponding integral conservation laws [l-4], are satisfied. With such an approach the conservativeness of a scheme is a property which is defined in terms of the conservativeness of the system of differential equations which is approximated. At the same time, since a difference scheme may be considered as an independent mathematical object, there is considerable interest in another approach which treats the conservativeness as a purely internal property of this scheme which is not directly related to which system is approximated or how it is approximated. In the development of this approach in this paper, a formal definition of the divergence of a difference operator and the concept of the conservativeness of a difference scheme associated with it is given as a purely internal property of such schemes. It is shown that, by virtue of the given definition, the conservativeness of a difference scheme, approximating a conservative differential equation, which is invariant with respect to similarity transformations is a sufftcient and, subject to a certain general assumption, necessary condition for its convergence on a shock wave. A simple criterion is obtained for checking actual difference operators for the property of divergence. A canonical form of representing a divergent difference operator is obtained by drawing an analogy between linear difference operators and algebraic polynomials. The method which has been developed enables one to estimate the unbalances which arise in the calculation of a shock wave using non-conservative schemes. Using this method we have obtained estimates of the unbalances arising in the calculation of a shock wave using a non-conservative scheme for the non-linear transport equation and in the calculation of the simplest discontinuous flows using the two-layer non-conservative schemes of gas dynamics in Eulerian coordinates.
1. Divergences of finite-difference operators 1. By T,hir we denote a shift operator, the action of which on each functionf(v(x,r)) T&f(v(z,
t))=f(v(s+ih,
is defined by the identity
l+jr)).
where v = (u ,, . . .,uN)is a vector-function which depends on the continuously varying arguments x and t,fis a specified scalar function and h and r are the spatial and time steps. Subsequently, we shall also use the abbreviated notation I): everywhere in the case when the steps h and z for the operator T,P are fixed. We will denote by AFT%a non-linear operator which is specified on a uniform template located within the (n, and n2 are integers), the action of which on each function v&t) is defined by rectangle [-R~, nJ’x 1-h ntl the _ ~~‘“‘[vl=~~‘~~‘~~of.,(v) m-l *-I (1.1) = &“,fi
“l-1 lZh.vychisl.
Mat.
mat. Fiz., 26, 11, 1661-1678,
1986.
Imr(v (z+i&
k-1
39
t+jn,r) ) ,
where --n,G&&nr, -n,Gj&32,, m--l, 2,. . . , M, k--l, 2,. . . , Km; f,,,&are known functions of v, a, are constant coefficients which are dependent on It and r and .f*d. Definifion I. We shall say that the difference operator (1.1) is divergent and that the difference scheme corresponding to it
AT [VI=0
(1.2)
is conservative if: (1) for all integers N,, Nr, M, and M2 such that .N,-N,Xn,, x0, z,, the sum
M,-ill,~?&,,
and for any real numbers
Nr =lr s = c c ht:,n’ Iv (softh, t&r) i-I*‘*f-.X,
1
(1.3)
after the reduction of similar terms is independent of the values of the function v(x,?) from the rectangle n-{30+(N,+nl)h~z~~~~(Nz-n,)h, t,+(M,+~)rat~t,+(M,-na)tj; 2)
Ah:‘“‘[Cl-0 vc * coni3t. It follows from this definition that, by analogy with the case of a divergent differential operator, the integral of
which over a region is written in the form of a curvilinear integral along its boundary, only terms remain the values of which are taken from a closed band of finite width along the boundary of the domain, as a result of the summation of a divergent difference operator over a sufftciently large domain. Remurk 1. Summation of a difference equation over a net region, in order to check it for the property of conservativeness, is applied quite frequentiy [1,2,4]. Therefore, in the present paper, the property of “summabiiity” lies at the foundation of the definition of the divergence of a difference operator and the conservativeness of a difference scheme. We also note that definition 1 may be considered as a fo~~ization of the requirement mentioned above in the introduction that a conservative scheme should satisfy a certain difference conservation law. 2. Let us now obtain a simple criterion for checking actual difference operators (1. I) for the property of divergence. In order to do this we give the following definition. Definition2. We shall cali the operator K_
the difference monomial. We shall say that two difference monomials %[v] and @,,,,[v] are analogous if numbers a and b exist such that the identity @,[v]-T~~~U&[V] holds. This relationship, which lies in the fact that two monomials are analogous, is an equivalence relationship and, by virtue of this, it subdivides all monomials occurring in the operator (1.1) into non-intersecting groups of anaiogous monomials. Let there be 2, such groups and let M, monomials occur in each I-th group. Then, the operator (1.1) can be written in the form (1.4a)
-nlGih+ir,Gz,, Lemma I. The condition
-nrGjrm+jrrGnl.
(1.4b)
Y.
E arm=0
Vl=1,2,...,L
(1.5)
is necessary and suf%ient for the divergence of the difference operator (f .4). Proof. Let us write the sum (1.3) in the form
Only monomials which contain values of the function v from the rectangle II remain in this sum. After reduction and taking account of (1.4b), we obtain
It follows from this formula that the sum S, vanishing identically is equivalent to condition (1.5). The lemma is proved.
41
As an example of the use of this lemma, let us consider the difference operator A[p, u]=IL(““(pu)I-‘/*u’““u(~~Jp,, where f”‘hof(tfr)+(l-o)f(t),
f,-r-‘[f(t+r)--f(t)]
, which approximates the divergent differential operator
d(pu) 1 ~--_-_*~_~-
dt
(1.6)
2
dp
d
dt
dt
This difference operator arises in passing from the difference equation for the internal energy to the difference equation for the total energy in the case of gas-dynamic difference schemes in Eulerian coordinates. Let us rewrite the operator (1.6) in an expanded form: A[& u]=-r-‘[ (or--hr)p~~-(l--cit-ar)pu~-((o*--X~)~uri f(l-o;-l,)p;ia+~.*p8’-2hlp;tP]; X(1-o,), ht=‘:,[~l(l-a*)+u*(l-a~)] has been where the notation f-f(t), f=f(th), ~FYWJ& L= %(i-u,) introduced. This operator only contains a single pair of analogous monomials: ii’ and pu*, Hence, by equating the sum of their coefftcients to zero and also the coefficients accompanying all the remaining monomials, we obtain a system of equations for determining the parameters u,, i = 1,2,3 for which the difference operator (1.6) is divergent:
20,~lfL-h,-Q, ~*==0*0s==0,
.o,--hs=O, 2aP=(l-o,)
l-a,-h,=O,
(i-%)~O.
Solving this equation, we find that u, = 0.5, us = 0, CT)= 1 (or u2 = I, u3 = 0). At the same time, the operator (1.6) is divergent and can be represented in the divergent form:
_ &zf-pu’ . 25 3. We shall refer to the different operator A~(TI-E)OAI+(T’-E)oA~,
E-T,“,
(1.7)
where A, and A, are certain operators of the form of (1.1) as a canonical operator. It is obvious that a canonical difference operator is divergent by virtue of definition 1. It is found that the inverse assertion also holds. ZXeorem I. Each divergent difference operator (1.1) can be represented in the canonical form of (1.7). Proof. Applying Lemma 1, let us write the divergent difference operator (1.1) in the form
and ‘I+are difference monomials. Let us transform (1.8) in the following manner:
(1.9) where
are linear operators It follows from differential operator, now set an algebraic
which are divergent, respectively, with respect to the variables I and x. formula (1.9) that each two-dimensional divergent difference operator (l.l), just like a divergent may be written in the form of a sum of one-dimensional divergent difference operators. Let us polynomial L,-
&IT, (-0
with the same coefficients a, to correspond with each one-dimensional
linear difference operator
42
This correspondence is an isomorphism of a ring of linear difference operators (Lm) with the operations of addition and composition on a ring of polynomials [Pm) with the operations of addition and multiplication. Hence, the wellknown results on the expansion of polynomials are automatically transferred to linear difference operators. Let us formulate the results we require in the form of a lemma. Lemma 2. A divergent operator L, is uniquely representable in the form L,,,- (T,-E) a&,-,. The proof follows from the fact that the polynomial P,,, (the sum of its coefficients is equal to zero) which corresponds to a divergent operator L,,, has a root x0 = 1 and, hence, according to Bezout’s theorem, is uniquely representable in the form P, = (x-l)p,_,. By transforming the operators A, and A@occurring in (1.9) with the help of this lemma and carrying out some simple calculations, we arrive at formula (1.7). The theorem is proved. It follows from this theorem that formula (1.7) may also be taken as a definition of a divergent difference operator and is equivalent to definition 1 (such a definition is used, for example, in [9]). At the same time, the difference analogues of integration by parts hold in the case of a divergent operator written in the canonical form (1.7)
[11. 2. Conservation law< for invariant difference schemes for continuous calculation on a shock wave 1. Let us assume that the difference scheme (1.2) approximates conservative differential equation
(in the sense of the definition presented below) a
where cp,@mC’. Definition 3. The difference operator (1.1) approximates the differential operator F and the difference scheme (1.2) approximates the differential equation (2.1) if, for all sufficiently smooth functions u(x,t) limhr:‘“‘[U]~F[U]. &T-r0
(2.2)
We note that theapproximation of an equation is usually [2] understood as a less strong requirement that Eq. (2.2) should be satisfied only on the smooth solutions of this equation. Let us next assume that, in accordance with the principle of invariance [lo], the difference scheme (1.2) (taking account of the fact that the time step r in it is chosen in accordance with Courant’s stability condition as being proportional to the spatial step h) allows of the similarity transformation Z’=%.z,
t’v&t,
h/-ah.
(2.3)
Allowing for this, we shall introduce the abbreviated notation A&(2,
t) I-O.
(2.4)
for it. We note that practically all difference schemes for one-dimensional planar-symmetric problems in gas dynamics which are used in practice allow of such a similarity transformation [ 1I]. The following lemma is central in the study of the convergence of scheme (2.4) on a shock wave. Lemma 3. Let h, and v&) (where .$= x - 0 t, @ = const) exist such that &[v&)]-0 for all EEB, Then, when h > 0,AJvr(t)] -0 is satisfied in the case of the functions
for all &&I. The proof transformation It follows 0 corresponds
follows directly from the fact that scheme (2.4) allows of the similarity E’-=a’& h’-ah in the self-modelling variable {. from this lemma that a whole family (2.5) of solutions v*(t) of this scheme with different step sizes !I > to each self-modelling solution v&) of scheme (2.4). At the same time. if
lim ME) - vo, t-.+-D
lim vha(t)I--m
VI,
vo+vi,
(2.6)
then, by passing to the limit as h - 0 with respect to this family of solutions (2.5), we obtain the limiting discontinuous solution, that is, the “rack” v’- (vO,E>O; v,, E
h(~k-b~(u)ll.:~-o
(2.7)
of Eq. (2.1) in the case of this solution does not follow by any means from the approximation of this equation by means of scheme (2.4) on smooth functions in the sense of definition 3, since the limiting values of the functions (2.5) when E-t= are independent of the scheme step size h.
43
2. Let us now investigate under what conditions the limiting discontinuous solution v* satisfies the conservation law (2.7). In order to do this, we shall write the operator A* (which, in the general case, is non-divergent) in the form of the sum of a divergent operator b* which approximates the divergent differential operator F and a nondivergent operator A, which approximates zero in such a manner that, just like scheme (2.4), the scheme &‘[v] -0 allows of the similarity transformation (2.3). The following theorem provides an answer to the question which has been posed. llrebrcm 2. The discontinuous solution v*, which is obtained by passing to the limit as h - 0 with respect to the family (2.5), (2.6) of solutions v&) of scheme (2.4), satisfies the conservation law +(2.8) [~(u)-~~(u)]l;;~Z, Z-r afn(E)lf% _I
when Zis a finite integral which is invariant with respect to A. Proof. Using Theorem 1, let us rewrite the divergent operator A** in the form ~-h-‘(T~-E)ohn,+7-‘(T’-E)onN,
(2.9)
where the operators AI, and lzll approximate the functions cp and Q respectively. Let us now substitute the solution v#) into scheme (2.4) which has been transformed with the aid of (2.9) and integrate the resulting identity with respect to f from +oo to -00. As a result, we have -_ -*
Introducing the notation AfE]=tlsJvr(@
1, we find
&I‘[-N-h,
-NJ,
& =[N-h,N].
Let us show that (2.12)
and obtain the first of these formulae. From (2.6), it follows that
limA[EJ-Ih~[~_~_v~fEll=I\n[v,~. 6-w” Since &, approximates cp and v, = const, then, taking account of the fact that the scheme &[v]-0 is invariant with respect to the similarity transformation (2.3), we obtain &Jn]-cp(vr) for all h. The second formula of (2.1 I) is proved in an analogous manner. Using (2.11) and (2.12), we transform the first term on the left hand side of (2.10): -. -I
-
lim AI&l- llm AteI-cp(v,)-I. t-r-%*+*
Starting from the fact that a shift of - 0 7 with respect to the self-m~elling in the variable f, we obtain l? (T’-E)Bb,[v,(E)1~=-6(9h)-~(v*)). 7 -*
variable t corresponds
to a shift of r
(2.14)
in an analogous manner to that used previously. Subtituting (2.13) and (2.14) into (2. lo), we arrive at formula (2.8). Its left hand side is finite and independent of Ir. Hence, the integral Zon the right-hand side of this formula is also finite and invariant with respect to h. The theorem is proved. It follows from this theorem that, if the difference scheme (2.4) is conservative, Z = 0, and the conservation law (2.7) is satistied in the case of its limiting solution v*. Let us now assume that, if A, * 0, scheme (2.4) allows of some solutions of (2.5) and (2.6) for which Z # 0 {in particular, when I is finite). This assumption is quite general since it follows from the fact that the operator A, is nondivergent by virtue of definition 1 that the integral I (see (2.8)) depends on the values of the function v&) at all points on the numerical axis, that is, on its values for all &Q;(-0, +=) . Under this assumption the conse~ativeness of the difference scheme (2.4) is also a necessary condition for its limiting discontinuous solution v* to satisfy the conservation law (2.7) of the approximated equation (2.1). In particular, it has been shown in f7] and [S] that the
44
conservation law for the total energy is violated in the case of the limiting solutions v* of non-conservative gas dynamic schemes in Lagrangian coordinates. It is interesting to note that the non-divergent operator A, can have as high an order of smallness with respect to /I as may be desired on smooth functions and, in spite of this, as follows from Theorem 2, it will lead, in the general case, to the violation of the conservation laws for the limiting solutions v* as h - 0. Remark 2. The result obtained may be considered as a direct extension of the well-known results [S, 61 on the convergence of difference schemes in a class of discontinuous coefftcients to the case of the convergence of continuous calculation schemes, which are invariant with respect to similarity transformations, on discontinuous solutions of conservative hyperbolic systems. 3. Let us now consider in greater detail the question of the convergence, as h - 0, of the solutions (2.5) and (2.6) to the discontinuous solution v*. In order to do this, we introduce the following notation: @I) is a monotonically decreasing positive function such that Iims(h)==O,
lii[k/e(k)
b-.*
h-O
R--R\[-~81,
3-0, (2.15)
llfll~-=M~) L.1
1,
where R is the set of real numbers. The functions (2.5) and (2.6) converge, as h - 0, to the discontinuous function v* on the whole of the R axis in a certain integral norm. However, by far the greatest interest lies in the convergence of these functions in the norm C. It can be readily shown that, for all t > 0, these functions converge in the space C[RJ. However, they do not converge in the whole of the space C[R\O] in which the function v* is defined. The following lemmas show in what sense it is nevertheless feasible to speak of the convergence of these functions in the norm C on the set R\O. Lemmu 4. Let there exist he > 0 and k 2 0 such that Then:
Proof. Assertion (1) is obvious. Let us now prove assertion (2). From the condition of the lemma, it follows that there exists h, such that, for every ei > 0, there exists N(h,,e,) such that, for every e>N(h,, E,) , we have lu,,,(E)-u’l<.s,/]~l’. Using (2.5) and replacing the variable f = h f/h,, we obtain Iv~(%‘)-v’I<(k/k,)‘el/l%‘l*.
(2.16)
This inequality is satisfied subject to the condition that 1r / >hN( h.,, eJ/h,. It fohows from (2.15) that a S > 0 Whence, aRowing for (2.16), we obtain that, for each exists such that, for every h 4 6, we have e(h) >hN(h,, et)&. e2 = q’ht, a Sexists such that, for every R < S and every ffr/>s(h) , }~~(%‘)-v*l<[kle(k)]*e~. The lemma is proved. It follows from this iemma that the functions (2.5) and (2.6), which satisfy these conditions, converge as h - 0 to the function v* in the following cases: (1) for all e > 0, t = const, in the space C(RJ and the order of this convergence is greater than !I’. (2) in the single parameter family of spaces C[RJ which pass into the space C[R\O] as h - 0. In this sense one may speak of the convergence of these functions in a norm C on the set R\O. The order of this convergence is greater than [h/e(h) I”. The following lemma is proved in a similar manner. Lemmo 5. Let there exist ho > 0, i L I and k 2 0 such that the following are satisfied: 1) for all
PER there exists v~(E)=d’v&Q’;
2) lim %%l’ (%)-0. I-r Then: (1) for a11k > 0 and
2) 2
&ER , there exists vhio (g) such that Iii ply;‘) (g) -0; r-r-
~8k~b)iV-‘JUYA(f~c%,ila,,,,~o.
It follows from this lemma that, when k > i and h - 0, the functions vb”‘(n converge to a function which is identically equal to zero in the following cases: (1) for all t > 0 in the space C[R,] and the order of this convergence is greater than hX-‘. (2) in the single-parameter family of spaces C[Rdr,] subject to the condition that the quantity /P/e’(h) is bounded as h - 0. The order of this convergence is greater than k”-‘le*(k). As an example of the use of Lemmas 4 and 5 we shall cite the solutions obtained in [4, p.631 of the differential equations of gas dynamics in Lagrangian coordinates with a linear artificial viscosity which, after the formal replacement of the coefftcient of viscosity Y by h, assume the form of (2.5) and (2.6). It can be shown that these solutions satisfy the conditions of Lemmas 4 and 5 for all k L 0 and i 1: 1. These solutions as we11as their derivatives
45
of all orders therefore converge, as h - 0, in the family of spaces C[R,,,] and the order of this convergence is equal to h’ (greater than hk for all k). In concluding this section we note the following important fact. If a certain solution v,(C) from the family of solutions (2.5) and (2.6) has oscillations on the shock wave front (in the neighbourhood of the point 5 = 0) then all the solutions from this family have oscillations of the same amplitude on the wave front. This is associated with the fact that the solutions ~~(6)when h = h, are obtained from the solution u&) as a result of compression of the latter by a factor of h,,/h along the f axis towards the origin of coordinates 6 = 0. As h - 0, all oscillation of the solutions (2.5) and (2.6) are thereby contracted towards the origin of coordinates while their amplitude remains constant. This means that, in the calculation of a shock wave using the corresponding difference scheme, it is impossible to attain any substantial reduction in the amplitude of the oscillations of the numerical solution solely by means of a reduction in the step size of the scheme h.
3. A method for the theoretical estimation of the unbalances of a non-conservative scheme on a shock wave 1. The violation of the conservation law (2.7) of the approximated equation (2.1) on the shock wave in the case of the non-conservative difference scheme (2.4) leads to the development of unbalance 6v,-vl-; in the determination of the variables behind the wave front and 6LB--6-a in the determination of the velocity of propagation of the wave front (; and @ are quantities associated with the exact relationship (2.7)). Let us assume that Z, Se!, and 6v, are quantities of the same order of magnitude which are small compared with v, and @. By varying Eq. (2.8), allowing for this, we obtain
This relationship, which is written in the form 6v &iv,(M, I), is conveniently considered as an unusual Hugoniot condition for the unbalances on a shock wave. For the approximate determination of the integral Z, let us suppose that Eq. (2.1) occurs in the system of equations
of the form of (2.1) which is closed, possibly by taking account of certain algebraic relationships G(u)--0
(3.3)
of the type of equations of state in gas dynamics, while the difference equation (2.4) occurs in the difference scheme which approximates this system and allows of the similarity transformation (2.3). With the help of the first differential approximation of this scheme [lo], we find the viscosity w which is characteristic for it which, in the general case, is written in the self-modelling variable in the form
where 6 is a certain sufficiently smooth function of two arguments such that e(u, 0) -0. By introducing this viscosity as an artificial viscosity into system (3.2) written in terms of the self-modelling variable f and integrating it after this with respect to 6 from +w to the current value off with the boundary conditions
lh
E-+0
u 6) -vo,
lhl
lim - = 0, re+r dt
(3.4)
we obtain
f&h;) Then, since a k L 1 exists such that the operator Z will be approximately defined by the formula
-Idu)-mmll.:. An/Z? approximates
(3.5) a certain differential operatorf,
the integral
where ;I(0 is the solution of (3.5),(3.3) with the boundary conditions (3.4) which correctly transmits the profile of a smeared out running shock wave [4,p.62]. Further simplifications are possible if the system (3.5),(3.3) taking (3.4) into account, like the analogous system of equations in gas dynamics [4], is solvable in the form
u,-u,(un,93),
t=4,2,.
(3.7)
. . ,N-i,
where, for all @ > 0, the function ~(uN, S?) is continuous and does not vanish in the interval (USA,%;I).The integral (3.6) can then be rewritten approximately in the form (3.8)
where the functionf, is obtained fromfby replacing all the fun~ions u,(z) and their derivatives in the latter by functions of uHand @ with the help of (3.7). Formula (3.8) will subsequently be the basis for determining the unbaiances. Remark 2. The method described above is based on a purely heuristic assumption concerning the “closeness” of the solution $0 of system (3.Q (3.3) to the corresponding solution of the difference scheme. Nevertheless, the unbalances in actual differences schemes of gas dynamics are determined using this method with a sufficiently high accuracy [7,8,12]. 2. As an example let us determine the unbalances which appear during the calculation of a shock wave using a difference scheme which is non-conservative when g # 0.5 Vt+V,,,t%=vVS, where Y-v(r,
t), YC==I-*(T~-E)*V, u,-~-~_‘(T,--E)~u; q=~-,*u., r--Ah,
which approximates
v-&h,
(3.9) u~~,-@E+(~-~)Z’_,~V,
A, c-const,
(3.10)
the non-linear transport equation
~+u$-0. Allowing for (3.10), this scheme admits of the similarity transformation Rewriting scheme (3.9) in the form -vvs+v,+(u’/2)r=-h(p-0.5)
and applying Theorem 2 to it, we obtain that its limiting ~~continuous relationship
(3.11)
(2.3).
(3.12)
(I&)*
solution U* (when u, = 0) satisfies the
At the same time, since here there is a fixed sign function, which is identically different from zero, under the integral sign, 1st 0 when j3 # 0.5. This means that the non-conservative scheme (3.9) on the shock wave does not satisfy the conservation law of equation (3.11). With the help of the first differential approximation we find that the characteristic viscosity in the case of scheme (3.12), which is the sum of an artificial and an approximate viscosity, has the form
By writing the viscosity in the self-modelling variable and transforming it after this taking account of (3.10) and also of the fact that, by virtue of the differential approximation of scheme (3.9), the equality
holds, we obtain
By introducing this viscosity as an artificial viscosity into the right hand side of Eq. (3.11) written in terms of the variable t and integrating it after this with the boundary conditions (3.4), WCobtain
au
-P
4
1
--u(rn-u). 2&h
Now, by applying formula (3.8), we find I=_
[email protected], 2Cih
5
&se
U (2s -z&u)&=
2 fP- 0.5)P 3Cl
Whence, taking account of (3. l), we have 2(~-0.5)SY ?r
8v,--2&B-
(3.13)
*
The quantity @ on the right-hand side of this formula is replaced by the quantity 3 which is close to it. 3. As an example of the use of this formula let us consider the numerical solution _ofthe problem of the propagation of shock waves which are specified by two different boundary conditions: ; =I 2 * @ = 0.5 (Fig. I, lines l-3) and G = 1 * 0 = 0.5 (lines 4-6) with respect to a uniform background v,, = 0. In this case, since u, = i and, consequently, au, = 0, we obtain from (3.13) the formula for determining 6 @ : (3.14)
32.6
3625
w
w.t5
,5p
3x1
.z
FIG. 1 b-0.1,S-4& linesZ-3- for A- b.2,C-0.5, T-40. 4-S- for A= 0.4,T-80; lines I-when@dt.~, (conscrvativc scheme), j 2 - when @i, 3 - when B-0, 4 -when b-s/,, C-V,, 5 -when B-S/s, C-.x0, 6 -when fi--I/,, C==‘/,. The unbalance (3.14) in the velocity of propagation
of the shock wave leads to the emergence of an unbalance 85~S-S--693T
(3.15)
in the distance traversed by the shock wave after a time ?‘(where S-62’ is the distance traversed by the wave front during the numerical calculation and S-&’ is the distance traversed by the shock wave when there is no unbalance). TABLE 1
The results of a comparison of the theoretical unbalances Ss,, mkuiated using formulae (3.14) and (3.15). and the unbalances &Sc which arise during the numerical caiculation are shown in Fig. 1 and Table I. Retnurk 4. In Fig. 1 the configuration of the shock wave front is shown by means of a vertical line. The broken lines are the positions of the fronts computed taking account of the theoretically calculated unbalance SS, (formulae (3.14) and (3.15)). The dotted and dashed lines are the shock waves which are obtained in the numerical calculation using scheme (3.9) with the corresponding parameters. The coordinate of the front of the wave in the numerical calculation is determined by the point of intersection of the corresponding graph with the line p = 0. The quantity A’S is defined bv the formula
4. The unbalances of non-conservative gas-dynamic schemes in Eulerian coordinates 1. Let us consider the system of differential equations of gas dynamics in Eulerian coordinates: (4.la)
9=-P&
e-Pe,
T-_(r-lke
(4.lb)
48
with a linear artificial viscosity w--vpi3u/as.
g=pfw,
(4.2)
By integrating this system in the self-modelling variabie t with the boundary conditions (3.2), we obtain, after transformations and allowing for the fact that r+,= 0 [4],
where
We approximate system (4.1),(4.2) by the differential difference scheme (4.4a)
(4.4b) Taking account of example (1.6), we obtain from this scheme the difference equation for the total energy:
where the unbalance of the Cst order with respect to T, which vanishes with the weights has the form 6E,=-r
[ (o,-a,)gl
G
5,=5,,
ape;,
5,-5,~0.5,
+ (0&.5)g@l’ $
1, while the unbalance of the order of rr has the form &E, - _$F
(&)‘.
By applying Lemma 1, it can be shown that the sum 6E, + 6E, is non-divergent for any choice of the time weights u,. The following theorem follows from this when account is taken of the fact that the divergence of its differentialdifference analogues is a necessary condition for the divergence of a difference scheme. Rkeorem 3. Among the two-layer difference schemes for which the schemes (4.4) are the differential~iffe~n~ analogues, there is not one scheme which is conservative. It follows from Theorem 3 that, among the above-mentioned schemes, there is also not one which is completely conservative. This result can be considered as a generalization of the analogous result in [9] to the case of an arbitrary spatial approximation. 2. Let us approximate system (4.1) by a difference scheme for which scheme (4.4) is the differential-difference analogue:
et--
(4Sb)
(~(51310,1)t_g(.,)~~~,
where f-Z’_& &=Tsa~a, g,--(2h)-‘( T,-T-J og. In this difference scheme not only the time approximation is non-conservative but also the spatial approximation (that is, the differential with respect to time and the spatial difference analogue of thii scheme is also nonconservative). A scheme with such a spatial approximation arises, for example, in the Euler direction of the combined Euler-Lagrange method proposed in [ 131. Under the conditions (3.10), the difference scheme (4.5) is invariant with respect to the similarity transformation (2.3) and, by virtue of this, Theorem 1 is true for this scheme. By transforming formulas (3.1), which correspond to scheme (4.5), we obtain
47P* 6% + (r+,)a’m
- -
By applying formula (3.6), we find, after some transformations,
2(7-f) D’(7+-i) (V-9) that
(4.6b) z
49
(4.7)
I=I,-kls,
where
(4.8a)
h--o,-o~-o,+0.5, +a
xP=u‘--o*,
h,-u‘-o,-u,+0.5,
(4.8b)
We shall refer to the quantity L, as the unbalance in the time approximation of the difference scheme (4.5) since this quantity is completely defined by the differential-difference scheme (4.4) which corresponds to this scheme. Similarly, we shall refer to the quantity Zs as the unbalance in the spatial approximation since it is completely determined by the differential with respect to time and spatial difference analogue of the difference scheme (4.5). It thereby follows from (4.7) that, in the approximation of the method being considered, the total unbalance of the difference scheme is equal to the sum of the unbalances of its time and spatial approximation. By supposing now that its artificial viscosity (4.2) is the characteristic viscosity in the case of scheme (4.5), we obtain from (4.8), making use of (4.3) and (3.8), that n zT_
(r+i)Tn% j 2v
- - (Yy;vrTnO
a(n) (Q-‘l)
(?-Pt)drl
;h,l,+LEi.+h,l,),
where
(4.10) 3. As an example of the use of the formulae which have been obtained, we shall give a calculation, using an explicit scheme (4.5) with the weights a ,-c~~=a~==cr~=u~=O,~,==-a,-I, u, -0.5*~,--lr-0, h-,-O.5 , of the problem of the decay of a discontinuity which arises when a diaphragm which separates cold gas (X) from a warmer gas (H) is instantaneously removed, where nOx-nOn-l, pox-O, p,,%- 2, r-2. When this is done, a shock wave will propagate through the cold gas with a velocity d = 1.13755 while a rarefaction wave will propagate through the warmer gas. The unbalances Sf,‘, 6D on the shock wave are related by Eqs. (4.6) where, when account is taken of (4.7), (4.9), and (4. lO), z=
-
(r+OD 24v
(h-r@)
[ (2no’+Snof)-4’) (no-q) -6~9
inn+] .
The unbalances on the rarefaction wave are related by the equations [ 121
where
c,=(yp,qo)”
is the initial velocity of sound. They are obtained as a result of a variation of the equations
which the gas-dynamic parameters on the rarefaction wave satisfy [3, p.1881. From this, by equating the unbalances of the pressures and velocities of the gases on the boundary of contact, we obtain a closed system for determining the unbalances 6~~-6u,~=6u~~, 6ps=6p+8p,‘, 6D, 6qix, 6qt’. After these have been found, the unbalances in the internal energy are obtained from the equation aa,-( ~-1)-1(@6~,+~6p,),which follows from the equation of state (4. lb). USER26:6-D
TABLE 2
1.1
0.0147 -0.0158
:i
0:s 9.4 0.8
-Fzz -0:0380 -0.Oli8
The results of a comparison of the theoretically computed unbalances with the unbalances which arise in the numerical calculation (h-0.1, z-0.02, v-0.05) are shown in Table 2, where
The author thanks S. M. Shugrin for useful remarks. Tronsiotedby E. L. S.
REFERENCES 1. SAMARSKII A. A., Theory ofDifference Schemes (Teoriya raznostnykh skhcm), Nauka. Moscow, 1983. 2. GODUNOV S. K. and RYABEN’KII V. S., Dijference Schemes (Raznosmye Skhemy), Nauka, Moscow, 1973. 3. ROWDESTVENSKII B. L. and YANENKO N. N., Systems of Quasi-Lineor Equations (Sistcmy kvasilineinykh uravnenii). Nauka, Moscow, 1978. 4. SAMARSKII A. A. and POPOV YU. P., Difference Methodsfor Solving Problems in Gas Dynamics (Raznwtnyc metody resheniya aadach gazovoi dinamiki), Nauka, Moscow, 1980. 5. TIKHONOV A. N. and SAMARSKII A. A., On the convergence of difference schemes in a class of discontinuous coefficients, Dokl. Akod. Nauk SSSR, 124.3, X9-532.1959. 6. TIKHONOV A. N. and SAMARSKII A. A., On uniform difference schemes, Zh. vychisl. Mot. mut. Fiz., 1, I, 5-63, 1961. 7. OSTAPENKO V. V. and SAPOZHNIKOV G. A., On non-conservative difference schemes of gas dynamics, in: Numerical Mefhods of the Mechunics of u Conrinvous Medium (Chisl. metody mekhan. sploshnoi sredy), Vol. 14, No. 2, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1983, pp. 103-l 13. 8. OSTAPENKO V. V., Estimates of the unbalances in the gas-dynamic parameters which arise in the calculation of a running shock wave using non-conservative difference schemes, in: Numerical Methods of the Mechanics of o Conrinuws Me&m (Chisl. metody mckhan. sploshnoi sredy), 14,6, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1983, 130-W. 9. KUZ’MIN A. V. andMAKAROV V. L., On an algorithm for constructing completely conservative dilfcrcncc schemes. Zh. vychisl. Mar. mot. Fir, 22, 1. 123-132, 1982. 10. SHOKIN YU. I., The Method of Diff erential Approximation (Meted differcntsial’nogo priblizheniya), Nauka, Novosibirsk, 1979. 11. YANENKO N. N. ET AL., Classification of the difference schemes of one-dimensional gas dynamics by the method of differential approximations, in Numerical Methods of the Mechanics of u Continuous Medium (Chisl. mctody mekhaniki sploshnoi sredy). Vol. 12, No. 2, Inst. Twr. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1980. pp. 123-159. 12. OSTAPENKO V. V., Estimates of the unbalances in the gas-dynamic parameters which arise in the calculation of the simplest discontinuous solutions using non-conservative difference schemes, in: Numerical Me:hodr of the Mechanb of o Continuotu Medium (Chisl. mctcdy mckhaniki sploshnoi srcdy), Vol. 15, No. 6, Inst. Tear. i Prikl. Mekhaniki Sib. Otd. Akad. Nauk SSSR, Novosibirsk, 1984, 110-124. 13. FRANK R. M. and LAZARUS R. B., A mixed method using the variable of Euler and Lagrange, in: Computationof Merhods in Hydrodynamics (Vychisl. melody v gidrodinamikc), Mir, Moscow, 1967.55-75.
U.S.S.R.
Cornput. Maths
Mafh.
Phys.. Vol. 26, No. 6,
Printedin Great Britain
pp. 50-59, 1986
0041-5553/86 E10.00+0.00 0 1988 Pergamon Journals Ltd.
A METHOD OF COMPUTING IDEAL GAS FLOWS IN PLANE AND AXISYMMETRIC NOZZLES WITH CONTOUR BREAKS* A. N. KRAIKO, N. I. TILLYAEVA and S. A. SHCHERBAKOV A method for computing inviscid and thermally non-conducting gas flow in a smooth or unsmooth nozzle, including one with a sudden contraction, is described. It combines a difference scheme that preserves the approximation on any mesh (in the case of an unsmooth generator its lines condense into a “sonic” or “supersonic” break), the method of characteristics, and analytic solutions, valid in the neighbourhood of the break and on the axis of symmetry. When using the method of characteristics to compute the supersonic flow
lZh.
vychisl. Mat.
mat Fiz.. 26,
11, 1679-1694, 1986.