Propagation of weak discontinuities along the bicharacteristic curves in a vibrationally relaxing gas

hf. 1. Engng Sci Vol. 20, No. 7, pp. t?Jl-%l, Printed in Great Britain.

1982

PROPAGATION OF WEAK DISCONTINUITIESALONG THE BICHARACTERISTIC CURVES IN A VIBRATIONALLY RELAXING GAS V. D. SHARMA and RADHE SHYAM Applied Mathematics Section, Institute of Technology, Banaras Hindu University, Varanasi-221005,India

Abstract-The transport equations representing the rate of change of discontinuities in the normal derivatives of flow parameters are obtained along bicharacteristic curves in the characteristic manifold of the differential equations governing the flow of a vibrationally relaxing gas, an explicit criterion for the growth and decay of weak discontinuities along bicharacteristics is given and the specialization is made to the waves of plane, cylindrical and spherical geometry. It is investigated as to how the effects of vibrational relaxation, the magnitude of the initial discontinuity, the initial curvature of the wave front, and the upstream flow Mach number will influence the growth and decay properties of the wave front. I. INTRODUCTION

IN THE theory of hyperbolic systems of partial differential equations in four independent variables (x, y, z, t), the characteristic surfaces play an important role. One of the important properties of a characteristic surface is that discontinuities in the derivatives of the dependent variables can exist only along the characteristic surface. At any point of a characteristic surface there is a privileged direction given by the tangent of the bicharacteristic curve passing through that point. Making use of the results of differential geometry, Coburn[l] derived an equation containing derivatives only along the bicharacteristic curve from the system of equations representing steady supersonic flow of an inviscid and non-conducting gas. Whitman[2], while studying the propagation of weak waves, derived an equation giving the rate of change of the wave amplitude along the rays (bicharacteristics). His theory has been used very successfully in computing the intensity of the sonic bang produced by a supersonic plane. Gubkin[3] making assumption of short waves, derived an approximate quasilinear equation giving the rate of change of the amplitude of a weak disturbance in the bicharacteristic direction. This work was further generalized by Prasad[4] for a general system of partial differential equations. Thomas[5], using the singular surface theory, derived the transport equations for the discontinuities in the normal derivatives of the dependent variables along orthogonal trajectories for an isentropic flow with a state of rest on the rear side of the singularity surface. When the singularity surface is adjacent to a region of rest, it is natural to consider the transport of discontinuities along orthogonal trajectories, as has been done by Thomas[5]. But, if the discontinuity surface is adjacent to a region of uniform flow in that case it has been shown by Elcrat [6] that it is natural to consider the transport of discontinuities along bicharacteristic curves in the characteristic manifold of the differential equations. In this paper, we have derived transport equations for the discontinuities in the normal derivatives of flow parameters along the bicharacteristic curves in the characteristic manifold of the differential equations governing an unsteady flow of a vibrationally relaxing gas. It is observed that the transport equations along bicharacteristics reduce to the corresponding equations along orthogonal trajectories when the singularity surface is adjacent to a region of rest. An explicit criterion for the growth and decay of weak discontinuities along bicharacteristics is given and specialization is made to the waves of plane, cylindrical and spherical geometry. In these cases the bicharacteristic curves are, of course, straight lines emanating from the centre of the wave propagation. 2. BASIC EQUATIONS

The differential equations, governing an unsteady flow of a vibrationally relaxing gas, are

$+

Uip,i + PUii

851

=

0,

(1)

852

V.D.SHARMAandRADHESHYAM ++

pUjl4i.j +p,i= 0,

(2)

(3) where the summation convention on repeated indices is employed and a comma followed by an index (say i) denotes a partial differentiation with respect to a space variable xfi The range of Latin indices is taken to be 1,2,3. The symbols appearing in (l)-(3) are as follows: p is the gas density, Uiare the gas velocity components, p is the gas pressure, h is the enthalpy and t is the time. The equations of state, when the gas is capable of being excited in its vibrational mode of internal energy, are p=p%T h=

LP.,” V-lp

(4)

(5)

where 9, T and v, respectively, denote the gas constant, the temperature and the frozen specific heat ratio, and e, is the vibrational energy per unit mass given by

2 +ui e,i

= o(ef( T) - e,) = x

where e*,(T) is the equilibrium value of the vibrational energy at the local temperature given by eS( T) =

SMJ{exp(&I T) - 1)

where 0, is the characteristic temperature of the molecular vibration. The quantity o in (6) is the relaxation frequency given by [7] (p. 67) w = A p exp(-BT-1’3) where A and B are positive constants depending on the physical properties of the gas only. Equation (3) with the help of (l), (5) and (6) is conveniently transformed into

(7) where uf = (VP/~)“’ is the frozen speed of sound.

3. VELOCITY OF PROPAGATION

Equations (l), (2), (6) and (7), using matrix notation, can be written as U,, +A’ U,i + B=O

(8)

where

u,t = au/at. The set U is a column matrix with six components p, Ui,e, and p. A' are the square matrices of order six and B is a column matrix which can be read off by inspection of (l), (2), (6) and (7). Let us now introduce a characteristic wave surface Z, +(x1, x2, x3, t) = 0, across which the set of dependent variables U and the interior derivatives of U along 8 are continuous but the exterior derivatives of U may have a jump across it.

Propagation of weak discontinuities along the bicharacteristic curves

853

The jump in a quantity f is denoted by [f]; the quantity f may represent any of the parameters p, Ui,e, and p. The square bracket stands for the value of the quantity enclosed immediately behind the wave surface minus its value just ahead of the wave surface. Since the parameters p, Ui,e, and p are essentially continuous across 8, we infer that the quantities ab o and x are also continuous across Z and that they will have their subscript-O values at the wave front, where the subscript-O indicates a value in the region just ahead of the wave surface 2. The medium ahead of the wave is assumed to be uniform and in equilibrium, we thus have x0

=

0.

(9)

Now introduce new coordinates p, t’, s’, 6’ to be chosen such that ,$‘,[‘, [’ determine a point on X which is itself defined by 4 = 0, and set p = 4. In terms of these new coordinates, the system (8) can be written as (Z4,, + A’4.i) U,, + A’($

u,zi=0

(10)

where Z is a unit matrix of order six. Forming the jumps in the usual way, we obtain from (10) (Z4.t+ A’4,iNUs+1= 0

(11)

This follows because the interior derivatives U,J are continuous across 2; the discontinuities existing only in the exterior derivative U,. Normalizing eq (11) by dividing by (4.i 4,i)“’ and setting 4,J(4,i 4,i)1’2= ni (the unit normal to the wave surface 2) and 4,J(4.i 4,i)1’2= -G (the speed of propagation of 8 along ni), the system (11) will have a non-trivial solution if det (A’ni - GZ) = 0, which yields V4(p - aL) = 0, where V = G - uflois the speed of wave front X relative to the fluid and 40 =

llio Ui*

Since, the case V = 0 in which the surface moves with the fluid is of no physical interest, we shall consider the case where V = + ufo.Here we shall consider only the advancing disturbance and, therefore, we shall take V= +a,,, or

G = u,,, + ufo.

(12)

The order of the characteristic matrix of the system (11) is equal to six, its rank is equal to five when (12) holds. Thus, for V= q,, the characteristic matrix has an eigen vector, the components of which are (1, uf, ndpo, 0, u;~). This means that the jumps in the normal derivatives of p, Ui,e, and p denoted, respectively, by 5, Ai, E and 5 are given by f =.poh/afo= du;,, = u

and

l

=0

(13)

where A = Aini and u is a non-zero scalar.

4. WAVE GEOMETRY

AND BICHARACTERISTICS

Let us recall that a bicharacteristic curve Xi= Xi(t) for the system (l), (2), (6) and (7) satisfies the equation [8] (14) where d/dt is the time derivative moving with the wave front along the bicaracteristic curve

854

V.D.SHARMAandRADHESHYAM

defined as

aa a ;i?=pia, (vi are the bicharacteristic velocity components) and F = (14,1+ A’+,i( = 0 is the characteristic equation for the system, which for the present case (i.e. V = a,J takes the form F

3

4.t + Ui4.i +

%(d,A,i)“‘.

(16)

Equation (14) together with (16) yields Vi = dt

=

Ui +

flf&

(17)

We define a/& to be the time rate of change moving with the wave front along the normal ni as follows

S=i+&.A St

at

’ aXi

(18)

which is the 6- derivative of Thomas [lo]. Equations (15), (17) and (18) can be combined to yield (19) where q, = %x~ay” are the components of the projective tensor of the discontinuity surface Z, Xi= xi(y’, y’, t); ~“(a = 1,2) being the Gaussian coordinates on the surface Z and ua = UiXka represent the velocity components along the discontinuity surface S. Thus for any tIow variable f defined on X, we have the relation

(20) since Xi,, f,i = f,=. It is evident from (20) that the two derivatives are identical if the discontinuity surface is adjacent to a region of rest. 5.DERIVATION OFTHE GROWTH EQUATION In this section we shall first derive transport equations for the quantities A and 6 along orthogonal trajectories and then we shall pass on to the bicharacteristics. If we differentiate the eqs (2) and (7) with respect to x~, take jumps across S, and multiply the resulting equations by nk, we find on using the relations (4), (9), (12) and (13) and the second order compatibility conditions due to Thomas [9] that +5-poarni=0 Lb($)+ +u”h,)

(21) (22)

where I= [ui,jklf4 nj nk,I= b,ijl 4 5 & = V(v - 1)2wet* exp

Propagationof weak discontinuitiesalong the bicharacteristiccurves

855

and

n =f g”@b,, is the mean curvature of Z with guB and b,, being the first and second fundamental forms of 2, respectively. Eliminating the term (I- poa,,i) between (21) and (22), we get

PO%(~+U.A.“)+(~+u”L.)+2(lb-o,n)z +u+4 pears

2

(23)

=a

II

In view of the relations (13), eqn (23) yields a differential equation for A and, therefore, one for 5 and one for J along the orthogonal trajectories of X. In eqn (23), the terms involving surface derivatives cause some difhculty in its interpretation, but if we transform (23) to a differential equation along bicharacteristics, this difficulty disappears. This can be done easily if we make use of the relation (20) in (23). Thus, eliminating from (23) the term ,$ in favour of A by using (13), the transport equation along bicharacteristics for the quantity A can be written as

g

+

(A,, - q&-W + t!+d A2 =

0.

The transport equations for the quantities ,$ and [ can be obtained from (24) by making use of the relations (13). Thus, the eqn (24) is sufficient to predict the growth or decay of weak discontinuities associated with Z along the bicharacteristic curves. For the situation envisaged here, i.e. for a wave surface travelling into a region which before the arrival of the wave is in a uniform state, the eqn (24) can be integrated to yield

(25) where A0 is the value of A at the wave front at t = 0. In order to predict the physical aspect more clearly, we discuss particular cases of plane, cylindrical and spherical waves. In these cases the bicharacteristic curves are, of course, straight lines emanating from the centre of the wave propagation.

6. DISCUSSION

Case (i). Plane waues For such a wave front n = 0 and the eqn (25) reduces to A=

AOexp(-A,

1 +?{l c

t)

-exp(-&

(26) t)}

where AC= 2A,,/(l+ v). Equation (26) shows that if A0 > 0 (i.e. an expansive wave front) then A + 0 as t + 01,and the wave decays. Also if A0 < 0 (i.e. a compressive wave front) and if it has the magnitude less than An then A + 0 as t + 00,i.e. a compressive wave decays and ultimately disappears. Further, if A0 = - A, then A = A0 and the wave attains a stable wave form. But if A0 is negative and has a magnitude greater UES Vol.

20, No. 7-E

856

V.D.SHARMA and RADHESHYAM

than A,, then A increases beyond all bounds for a finite time t approaching the value -I

t,=fQlog

( > 1 +k

.

(27)

Thus, at the instant t,, the velocity gradient at the wave front becomes infinite and this signifies the appearance of a shock wave. The initial value of the discontinuity A,,= -A, plays an important role in determining whether the wave shall grow or decay. Accordingly, we designate the value A,,= -A, as the critical A, for a plane wave. From (27) it follows that (&jaA,) > 0, i.e. an increase in A,, causes an increase in t,. Thus, the vibrational relation has a delaying effect on the formation of a shock wave. The behaviour of A/A, vs 7(=IAo/t) for different values of pz(=Ao/lAo)is exhibited in Fig. 1 for the adiabatic exponent v = 5/3. It is found that the effect of vibrational relaxation causes an expansive wave of given initial discontinuity to decay at a faster rate relative to what it would be in the absence of vibrational relaxation. Case (ii) Cylindrical waves

Consider an outward travelling wave front which at t = o is a cylinder of radius RO.Then, at any later time the wave front is a cylinder of radius R = & + G t. For such a wave R = -1/2R, uno= u. the radial velocity ahead of the wave and thus G = a&U0 + l), where MO= (uo/af,J is the upstream flow Mach number. Then (25) becomes

Fig. 1. Effect of vibrational relaxation (p2) on the growth and decay of plane waves for v = 5/3.

Propagation of weak discontinuities along the bicharacteristic curves

857

where = 2n&&Jaf0, P = 2n&oRlaf,, CLO and n = l/2@&+ 1) are positive quantities. I(n, x) is an improper integral, given by I(n, x) = J,”t-” exp(-t) dt, which is convergent for positive values of n and x considered here and ,iC= (((1 + v)/2ho)& exp(& Z(n, h)}-’ is a positive critical value of the initial discontinuity in the sense discussed below. The term inside the bracket in the denominator of (28) increases monotonically from 0 to 1 as R increases from R. to 03.Hence if A0> 0, or if A0< 0 and IA01 < iC then it follows from (28) that A+0 as R-03, the wave decays. However, if A,<0 and IA,,/= & then at any time t, the discontinuity (A)is given by

which approaches the value 2&&l + V) (critical Alofor a plane wave) as R-m, i.e. the wave takes a stable wave form. But if A,<0 and lAoI> A, then there exists a finite time iCgiven by

Zh

AL +1.4)= (l

-6,Zh PO)

(29)

such that IAI+ m as t + &.This signifies the appearance of a shock wave at a finite time &.Equation (29) shows tha,f the shock formation time iCdepends on Ao,Ao, MOand Ro. The dependence on A0 is such that (&/alA,-,()< 0 which means that an increase in the initial value of the discontinuity leads to a rapid shock formation. Also it follows from (29) that (&/JR,) < 0 which means that an increase in the initial curvature will cause the shock formation time fCto increase. The differential eqn (24) can be written in the following non-dimensional form

(1+ g+(P2-PI)s+-~1

2

ho

s2 = o

(30)

boi

where 6 =

pI

=

-

AlAo,

(a/2&’

T =

/Aok

+ (1 + Mo)$’

~2

=

Ao/lAol,

and

p3 = (~Ao~Ro/af,)-‘.

In writing (30), we have used the relation 2Cl= - a/R, where LY= 0, 1 and 2 for plane, cylindrical and spherical waves, respectively. Some of the solution curves of (30) for a = 1, and v = 5/3 obtained by using Runge-Kutta method of order four, are exhibited in Figs. 2-4. Figures 2 and 3 indicate that the effects of vibrational relaxation and that of the wave front curvature are to reduce the steepening rate of a compressive wave and to increase the flattening rate of an expansive wave; Fig. 4 shows that the upstream flow Mach number has an exactly a reverse effect on the growth and decay properties.

Case (iii) Spherical ivan3 In this case also the growth and decay phenomenon is very much similar to those of plane and cylindrical waves. Now, consider an expanding wave front which at t = 0 is a sphere of radius Ro. Then at any time t >O, the wave front is a sphere of radius R = R. + Gt. For such a wave CI= -l/R and (25) reduces to A

=

AO(~RO)-mexp(~0

-P)

(31)

V. D. SHARMA and RADHE SHYAM

858

0

“I

3NM

3hlSS3HdVf03

MOJ

OY

-

3AVM

7

3AVM

3hlSS3Yd)(103

HOd

f

-

3 L

b

-

3AVM

3hlSNVdX3

3hlSNVdX3

9 MO4 F

MO4

OY F

6

2

Propagationof weak discontinuitiesalong the bicharacteristiccurves

0

Y)

ul

b 3AVM

0

3hlSS3MdW03

804

OY yj-

n

OY 3AVM 3hlSS3MdW03 YOj y -

-

3hVM

0

3AiSNVdX3

MO4

ox F

m*

b

-

MVM

3hlSNVdX3

dO4 ‘+

V. D. SHARMA and RADHE SHYAM

860

T

n

0

3AQM

MlSS3YdNO3

MO4

$-

-

3AVM

3hlSNVdX3

HO4

Ti-

Propagation of weak discontinuities along the bicharacteristic curves

where

861

m = 2n, and ic=

K

!$

>

j.~$exp(~~)Z(m, P&‘.

(32)

It follows from (31) that for A,,< 0, there is a positive critical value h’, given by (32) such that waves whose initial discontinuity is less than ii, damp to zero and waves whose initial discontinuity is greater than h’, grow without bound in a finite time &, given by the solution of Z(mAoi+ti)=(l-&Z(m,&. For A0< 0 and [A,,[= n’,,it follows from (31) that the wave can neither terminate into a shock nor it can ever completely damp out. Using L’Hospital and Leibnitz rules, it follows from (31) that JAI+2ho/(l + u) (critical ho for a plane wave) as R+m. Thus, a compressive wave with [A,,[= h’, takes ultimately a stable wave form. A few solution curves of (30) for (Y= 2 and v = 5/3 are exhibited in Figs. 5-7. It is ivident from Figs. 5 to 7 that the behaviour of solution curves of (30) for a spherical motion is very much similar to that as for cylindrical motion. It is interesting to note that the time taken for the shock formation by a spherical wave is greater than that for a cylindrical wave (Figs. 2,3,5,6). Acknowledgement-One of us (RX) is thankful to C.S.I.R. (India) for financial support. REFERENCES [l] N. COBURN, Q. A&. Math. 15, 237 (1957). [2] G. B. WHITHAM, .I Fluid Me&. 1, 2% (1956). [3] K. E. GUBKIN, J. Appi. Mar/r. Mech. 22,787 (1958). [4] P. PRASAD, J. Mat/t. Anal. Appl. 50,470 (1975). [5] T. Y. THOMAS, J. Math Mach. 6,455 (1957). [6] A. R. ELCRAT, Int. .I Engng Sci 15,29(1977). [7] K. N. C. BRAY, In Non-Equilibrium Flows (Edited by P. P. Wegner), Part II, Marcel Dekker, New York (1970). [8] R. COURANT and D. HILBERT, Mefhods of Mathemafical Physics, Vol, II. Wiley, New York (1%2). [9] T. Y. THOMAS, Plastic How and Fracture in Solids. Academic Press, New York (l%l). (Received 15 August 1981)