UM-013610-2-T TECHNICAL REPORT VISCOUS TRANSONIC FLOW IN RELAXING GASES by M. Sichel and Y. K. Yin prepared for Army Research Organization Durham, North Carolina September 1971 Department of Aerospace Engineering The University of Michigan Ann Arbor, Michigan This document is subject to special export controls and each transmittal to foreign governments or foreign nationals may be made only with prior approval of the U.S. Army Research Office-Durham.

ABSTRACT A general differential equation governing the flow of a viscous, relaxing gas valid in the transonic range is derived by using a procedure similar to that used by Clarke and McChesney (1964) in obtaining the inviscid relaxing equation. Using this equation, the structure of a shock wave with relaxation is studied in detail in order to clarify the effect of viscosity in relaxing flows. Asymptotic solutions, valid in three different regimes of flow, depending on the order of magnitude of the free stream frozen Mach number Mf with respect to unity, have been obtained. The results show that the viscous solutions successfully eliminate some of the difficulties encountered in the inviscid theory such as the discontinuous velocity profile which occurs in the partially dispersed wave and the upstream corner which arises in the fully dispersed wave solution for Mf = 1. It is also shown that the inviscid relaxing flow equation is valid as long as Mf <K 1 and nowhere close to one. Furthermore, a bulk viscosity can be used to account for the relaxation when the equilibrium Mach number Me is close to one. The analytic results have been verified by the numerical solution of the general equation. Two dimensional versions of viscous transonic relaxing equations valid for frozen (Mf0 - 1) and equilibrium (Me0 - 1) regimes have ii

also been obtained. No attempt was made to solve these equations. However, the equilibrium transonic relaxing equation is found to be identical to the viscous transonic equation in an inert gas except for the fact that the dissipative term 1 + ((yf - 1)/Pr") v" due to compressive viscosity and heat conduction has been replaced by the term (af - 1) which is, essentially, a bulk viscosity. 00 iii

ACKNOWLEDGMENT The support from the Army Research Organization under contract DAHC 04 68 C 0008 is gratefully acknowledged. iv

TABLE OF CONTENTS Page LIST OF FIGURES vi LIST OF SYMBOLS viii I. INTRODUCTION 1 II. DERIVATION OF THE VISCOUS EQUATION FOR GASES WITH ONE INTERNAL DEGREE OF FREEDOM 10 III. ANALYTIC STUDY OF ONE DIMENSIONAL WAVE STRUCTURE 30 A. Inviscid Solutions 31 1. Fully Dispersed Relaxation Wave 31 2. Partially Dispersed Shock Wave 34 B. Viscous Solution 37 1. Equilibrium Transonic Approximation - the Bulk Visco 39 2. Fully Dispersed Wave 44 3. Partially Dispersed Wave 57 4. Partially Dispersed Wave Alternate Approach 71 IV. PHASE PLANE STUDY AND NUMERICAL INTEGRATION 83 V. TWO DIMENSIONAL VISCOUS TRANSONIC EQUATIONS 89 VI. CONCLUSIONS AND DISCUSSION 94 REFERENCES 103 v

LIST OF FIGURES Figure Page 1 Velocity Profiles through a Fully Dispersed Wave from the Inviscid Solution (3. 5) 33 2 Velocity Profiles through a Partially Dispersed Shock Wave from the Inviscid Solution (3. 9) 36 3 Phase Plane Plot of Inner Equation (3. 45) 50 4 Qualitative Sketch of the Matching of Inner and Outer Solutions for a Fully Dispersed Wave with af = uc 56 5 Phase Plane Plot of Inner Equation (3. 75) 64 6 Qualitative Sketch of the Matching of Inner and Outer Solutions for a Partially Dispersed Shock Wave 70 7 Comparsion of the Ideal and Actual Behavior of the Internal Energy Mode e2 through a Viscous Shock Wave 77 8 Qualitative Sketch of the Matching of Inner and Outer Solutions of a Partially Dispersed Shock Wave from an Alternate Method 81 9 Summary of Various Regimes and Their Solutions Depending on the Relative Order of Magnitude of uoo- a c, uca- a and af - a 82 IQO ecO oo eoo 10 Matching of Inner and Outer Solutions of a Fully Dispersed Shock Wave with af = uo 1. 0 96 11 Comparison between the Solution Obtained from the Method of Matched Asymptotic Expansions and the Result from Numerical Integration of Equation (4. 1) 97 12 Inner, Outer and Composite Solutions for a Partially Dispersed ShockWave with Mfa =1. 2 98 vi

Figure Page 13 Comparison between the Solution Obtained from the Idea that the Viscous Shock Drives the Relaxation and the Result from Numerical Integration of Equation (4. 1) 99 14 Velocity Profile of a Partially Dispersed Shock Wave with M = 1. 2 from Numerical Integration of Equation (4. 1) 100 15 Velocity Profile of a Partially Dispersed Shock Wave with a Very Weak Viscous Shock with Mf = 1. 05 from Numerical Integration of Equation (4. 1) 101 16 Velocity Profile of a Fully Dispersed Wave with Mf = 1. 0 from Numerical Integration of Equation (4. 1) 102 vii

LIST OF SYMBOLS A = 1 + (f- 1)/ Pr" =0(1) a sound speed B = 1 + (yf - 1)/Pr (Cvf/Cve) = 0(1) c specific heat at constant pressure c specific heat at constant volume v C2 specific heat associated with internal energy mode D coefficient of self-diffusion E = (U - us)/(uo- u ) e specific energy = el + e2 el energy content from translational mode e2 energy content from the internal energy mode h specific enthalpy K heat conductivity K modified heat conductivity defined as K(1 + Le) C Le Lewis number based on c2 M Mach number n number density Pr" Prandtl number based on p" p thermodynamic pressure q. heat flux vector qj viii * * Vlll

(int) q() the contribution to the heat flux due to the internal qj energy mode s specific entropy T1 translational temperature T2 internal energy mode temperature u. velocity vector u velocity at upstream infinity u velocity at downstream infinity s v. diffusion velocity vector r =( y +1)/2 e e rf =( y + 1)/2 y ratio of specific heats 6 small parameter associated with non-dimensional velocity differences e small parameter associated with the ratio of two characteristic lengths = AVX/c ~71 ~ small parameter for the intermediate region x a viscous length defined as v"/u or v"/a X thickness of a weak viscous shock defined as Xvl6 A. relaxation length defined as r- (Cvt/cve) a Tik viscous stress tensor r relaxation time defined as De/Dt = (e2 - e2)/ e ix

1. p f (1) shear viscosity longitudinal viscosity /L"/P gas density perturbed velocity potential viscous dissipation Subscripts 1 2 f e s b (-) translational energy mode internal energy mode frozen condition equilibrium conditions upstream infinity downstream infinity state immediately after the frozen viscous shock barred quantities are dimensional Superscripts (1) (2) first order perturbation quantity second order perturbation quantity x

I. INTRODUCTION Non-equilibrium gasdynamics has long been an area of extensive study because real gas effects such as, chemical reaction, vibrational relaxation, dissociation, etc., are often present in high speed aerodynamics. Much work has been published in the past. Most of it is based on the linearized equations, however, much important information is nevertheless obtained. A linearized inviscid equation governing the steady, non-equilibrium or relaxing flow about a slender object in a uniform parallel stream was obtained by Vincenti (1959 ), and Moore and Gibson (1960). The equation takes the following form -+- 9 /,. 2 a2 a2 a21 2 a 29 T- u (M l) -- -2+(Me 00 a 1)f 2.2 2 2 1xL a9x ay a9z oo ax 2 a2 0- — = 0 (1. 1) -2 -2 ay az where 4 is the perturbed velocity potential, Mf and Me are the free stream Mach numbers based on the frozen and equilibrium sound speeds, and f+ is the relaxation time evaluated at free stream 00 condition. The bars here denote dimensional quantities. 1

2 It is apparent that Eq. (1. 1) will no longer be valid, just as in the case of the equation for an inert gas,when Mfc - 1, as the first term may then be of the same order of magnitude as the terms already being neglected. Because of the dissipative character of the relaxation process, Vincenti's solution for the Ackeret wavy wall problem exhibits a drag force even in subsonic flow whereas the inviscid solution for the inert gas does not. Clarke (1961) and later Li and Wang (1962) worked out the steady three dimensional relaxing flow over a pointed body. Supersonic flow past a sharp corner was studied by Clarke (1960) later by Der (1961) and extended by Vincenti (1962). The analysis of all the above stated work made use of the integral transform method and Eq. (1. 1). Chu (1958) was the first to attack the one dimensional, unsteady, relaxation problem by solving the following linearized equation which can be found, in standard textbooks of gasdynamics such as Clarke and McChesney (1964) or Vincenti and Kruger (1965) T(r2tt - x + - -tt - 0 (1. 2) 00 xx 1_ f ae 00 00 In Chu's analysis, the Laplace transformation was used to solve the partial differential equation. Chu's solution sihows that the signal

3 or the discontinuous wave front, which is apparently the result of using an inviscid theory, runs ahead at the frozen sound speed independent of the value of relaxation time, and at the same time is attenuated by relaxation. The bulk of the wave, except perhaps for a very short time after the piston starts moving, essentially travels at the equilibrium sound speed. However, at large times, this linear theory predicts a wave which diffuses gradually with time and never reaches a steady state. This decay is caused by the neglect of the cumulative effects of the non-linear convective terms. This result again suggests that Eq. (1. 1) and Eq. (1. 2) are not valid in the transonic range. The one dimensional, non-linear, steady problem was discussed by Broer and Van Der Bergen (1954) and Lighthill (1956). Lighthill (1956) divided the problem into two parts depending on whether the Mach number based on the frozen sound speed Mf =o d/af is greater or less than unity. When the free stream velocity is in between the two sound speeds or afc > u >ae a fully dispersed wave is formed, according to Lighthill, and the convective steepening effect is balanced completely by relaxation. On the other hand, when Mf > 1 i. e. u > af a partially dispersed wave is formed. However, the relaxation alone is no longer sufficient to resist the convective steepening

4 effect and the viscosity must then play a role in determining the wave profile. Furthermore, afc and ae, were shown to be the upper and lower bounds for a fully dispersed wave to exist. Spence (1961), Ockenden and Spence (1969) and Blythe (1969) worked out the one-dimensional non-linear unsteady problem. Broer and Van Den Bergen (1954) were the first to study the relaxation wave structure by including also the viscous effects. Starting with the conservation equations and the linearized rate equation together with some very ingenious series expansions for velocity, translational temperature and rotational temperature, Broer and Van Den Bergen were able to reduce the conservation equations to the following relation dh I3h dh 1 Ph —- ^ --- ~~(1.3) dg 2 dg = (1 - g ) - ah where h and g represent the perturbed rotational temperature and velocity field respectively. The solution was then obtained through numerical integration. Furthermore two special analytic solutions for 3 >> 1 and 3 << 1 were also given where 3 is a parameter characterizing the ratio of a viscous length to a relaxation length while a1 apart from a factor, measures the viscous shock strength in terms of the critical strength. Broer and Van Den Bergen's approximate solution for the case j >>1 or when the viscous length is much longer than the relaxation length shows a pure viscous shock solution plus an extra term due

5 to relaxation which acts to broaden the velocity profile at the downstream side somewhat more than at the upstream side. For the case (3< 1, i. e. a very slow relaxation process, Broer and Van Den Bergen' s solution which is valid only for a~>>1 also consists of two terms, one due to relaxation, modified by a higher order term due to viscous effects. However, because of the various linearizations and approximations used in the analysis, both the analytic and numerical solutions give only qualitative results. In Broer's work, no attempt was made to compare the two approximate solutionsto the numerical solution. Furthermore the choice of variables and the series expansions used by Broer made it difficult to study the case when u -af (01 = 2, in Broer's notation), which will be an essential part of the present study. It is also necessary to keep in mind that for most gases it usually takes more molecular collisions to establish equilibrium for the rotational mode than the translational mode, although the characteristic time for translational and rotational modes of dense gases can hardly be distinguished and one has TT Tr (Losev and Osipov, 1962). For other internal degree of freedom such as those due to molecular vibrational and chemical reaction, many more molecular collisions are always needed to establish the equilibrium state. Therefore the case ( >> 1 is physically of questionable significance for rotational relaxation and definitely not applicable for other internal modes.

6 The most practical problem in relaxing gas dynamics is, perhaps, the nozzle problem when the flow is being accelerated from subsonic to supersonic velocity by passing through a converging-diverging passage. Most work on this problem is again based on the linearized equation (1. 1). Rhyming (1963) and Tirumalesa (1966) studied the non-linear transonic case. However Rhyming did not solve the nonlinear problem directly, instead the Oswatitsch (1956) transonic approximation was used to replace the non-linear convective term uu by cu where c is a constant corresponding to u evalux x ated at the sonic point. Tirumalesa (1966) derived, following the procedure used by Vincenti (1959), a non-linear differential equation for inviscid transonic flow. By assuming a series expansion for the 2 perturbed velocity potential in terms of y, where y is the vertical distance measured from the nozzle axis, Tirumalesa was able to obtain solutions up to second order. The frozen sonic line was found to be parabolic and downstream of the geometric throat. The difficulty arising from the use of the inviscid theory is readily seen in the one dimensional wave problem. When the free stream velocity is less than but very close to the frozen sound speed, Lighthill's solution for a fully dispersed wave is not accurate near the upstream portion of the wave (Clarke and McChesney, 1964), and

7 breaks down when u = af as a corner appears in the upstream portion of the wave profile. The upstream boundary condition then cannot be satisfied continuously. On the other hand, the inviscid solution for u > af is not unique and two velocities occur at each o00 00 location. A viscous discontinuity has to be introduced to solve the difficulty. It is to be expected that these difficulties will also exist in the case of two dimensional flow. From the success of viscous transonic theory (Sichel, 1963) in solving some of the difficulties encountered in inviscid transonic theory of inert gas flow, it is to be expected that viscosity and heat conduction will also play an important role in inviscid transonic relaxing gas flow. However, it is to be expected that the resultant governing differential equation will not only be very complex but will also be highly non-linear, which is essentially the feature of transonic flow. An approximate equation or equations will have to be found for analysis to be possible. Napolitano (1966) formulated the approximate equations of an inviscid relaxing gas valid in different transonic regimes. He showed that, depending on the parameter (X/Xc) IMe - 1 /IMf2 - 11, three different simplified differential equations can be obtained by stretching the variables as indicated below, where X is a characteristic macroscopic length and c is a characteristic relaxation length.

8 (a) IMe 2 _ ll=O(l);(Mf 2 _ l=-O(E) C0 0 in the frozen transonic regime,and proper independent variables are IM 2 e -ii 0o -- o (1.4) IM 2-f c IM -1 oo 0o (b) IMe 2 1IO(c);IMf 2 - 1=0(1) 00 00 in the equilibrium transonic regime, and proper independent variables are IM 2-11 e oo 3 c iMf 2-1 00 (1.5) IM 2-11 0e 2 IM -1ly 03 1. 2ilj e c M 2 1 oo 00

9 (c) IMe 2 1 - =0(); IMf 2 1- =O(e) 00 0 in the proper transonic regime and the relevant independent variables are a2 and r12. The first two regimes can be attained in any medium, while the last one can be attained only in those media for which the difference between the two speeds of sound is small. The present work will first attempt to derive, by following closely the procedure used by Clarke and McChesney (1964), the differential equation governing a viscous, relaxing gas with only one active internal energy mode. Then the one dimensional shock wave structure will be studied in detail with emphasis on the resolution of the discontinuous wave profile of a partially dispersed wave and the corner solution of the fully dispersed wave as u0 = af encountered in Lighthill s inviscid solution. It is hoped that this systematic study of the one dimensional case will provide a guide to the solution of two dimensional problems such as flow inside a nozzle or past a thin airfoil.

II. DERIVATION OF THE VISCOUS EQUATION FOR GASES WITH ONE INTERNAL DEGREE FREEDOM The procedure here used is essentially the one used by Clarke and McChesney (1964) in deriving the inviscid relaxing gas equation. The steady flow of a hypothetical pure gas is considered whose communicable energy states are those of molecular translation with energy content specified by a translational temperature T1 plus a contribution from the internal energy mode with an energy content specified by an internal energy mode temperature T2. It is also assumed that the energy states of the internal mode are volume-independent. This model is chosen for analytical simplicity. Then the following thermodynamic relations hold: T dS de +pd (1/p) (2.1) T2 dS2 =de2 (2.2) =e e1 + e2 (2.3) S S1 + S2 (2.4) where S1, 2 and el, e2 are the specific entropies and energies associated with the translational and internal mode. It is further assumed that the flow is only slightly out of equilibrium so that the linearized relaxation equation can be used 10

11 e2 T Uj aX- e2 ~ e2 (2.5) ] e where 1 is the relaxation time and subscript e means equilibrium value. The above equation can also be regarded as the definition of the relaxation time F. In general, F is a function of pressure and temperature. Furthermore, Eq. (2. 5) is exact for the harmonic oscillator (Vincenti and Kruger, 1965). If u. is the velocity, p the gas density, p the thermodynamic pressure and h the enthalpy, then the equations of continuity, motion and energy in Cartesian tensor notation are respectively au. UJ ap-+ P a = (2. 6) ] 3 au. aa F i ap ik pj - + - (2.7) ax. a x. a8. P a U = + 4i (2. 8) ] ] 3 where Tik' the viscous stress tensor, is given by au~ au. a-k a 1 ik b ) ik + ( (2.9) ik-' lik a a+ qj, the heat flux vector is 3

12 aTI (int) - -q (in (2. 10) and 4, the viscous dissipation, is given by au. aui -)ik( a +' -)/2 (2.11) k i In these equations, q (int)is the contribution to the heat flux vector due to the internal mode, K is the thermal conductivity, i is the viscosity and -i is the bulk viscosity which can, under special conditions, be used to account for the relaxation of an internal mode in addition to the single mode considered here such as rotational relaxation for example. However ib will be taken as zero in the present case as only one internal mode is assumed. The complex system of Eq. (2. 5)-(2. 11) will now be reduced to an equation with a single independent variable. Adding (2. 1) and (2. 2), and using Eq. (2. 3) yields de = T1 dS1 + T2 dS2 - pd (1/u) (2. 12) Using (2.4), Eq. (2. 12) becomes de -T1 dS + (T2 - T1) dS2 - d (1/n) (2.13)

13 or ~~~- ^^~~~as T U. - as a +p Ua( /p- ) t T (2.14) 1 j 8.= uj a. + p u a] x 2 1 a3'x J.J ] 3 Rewriting the energy Eq. (2. 8) as a 3 a(l/p) 1 a_ j uj -. + p u J = -( ) (2.15) and combining with Eq. (2.14) yields as 1 1 aq s U -=T [- p-,(% -aq) - T j ] (2.16) j ax. T= T px(. 2 axk Equation (2.16) is an entropy equation. Since the internal mode is assumed volume-independent, the density p is a function, only, of p and S1, consequently dp =(a P) dp+()3 dS1 (2.17) 1 p where ap aS ap =(aP)s T 2 -f (2.18) p Q p'2 2 S SIIT2 af (by definition). af is then the frozen sound speed in contrast to the equilibrium sound speed e which is defined as e

14 ap 1 (2. 19) F_ _: -2 S, T2- e a The two sound speed defined above are thermodynamic state variables. The calculation of af and ae for a gas with an arbitrary equation of state can be quite tedious. In general, Vincenti and Kruger (1965) show that af and ae can be written as 122- 1/ af-(@) =/ P (2. 20) af:a s, q = a~1/2 / h_+h- q 1/2 a P() ( P q Pe (2. 21) ae (S,a Fe- + _ q (1 where q is the nonequilibrium variable and is simply T2 in the present case. Equation (2. 17) can be written as I ajp -22 a(1/p) ap as - U a P 2 u a(l/ ( ) -U a-F (2. 22) -. 2 ] jx. x9. 2 = -p u. af and then using Eq. (2. 16) yields a- _2 2 a(/2 r +1 f2 (2.23) j a. =- p u. a- a'. aE - (2. 23) j 3f aR. ]a x i J J J

15 where l —- (AN= f (2. 24) 1 P Pf Cpf is the frozen translation specific heat at constant pressure and f the frozen thermal-expansion coefficient. By using the continuity equation which is in a slightly different form, Eq. (2. 23) becomes au. < a 3e ap 2 j 2 [ a2 1( j_ ( u. P+ p a + Ef — f = E0 (2.25) J x. f ax. f a. p aR. J 3 ] J Now, using the momentum equation (2. 7) to eliminate the pressure term in Eq. (2. 25) yields, after some arrangement, essentially the thermodynamic equation, a e2 1 a ( 2/2) aj.2 f aqj U* ae2 fafP- a "j^-'rryF 4%. f aR. c- aR. ip af P f 3 a ik - au.i (2. 26) The momentum equation (2. 7) can also be used to eliminate the pressure p in the energy equation (2. 8) which then becomes u =1+ k u+ 1 (2. 27), [j a3 J'a a 3 J J~~~~~~~~

16 but h = c T + 2 T2 (2. 28) pf 2 where c2 is the specific heat associated with the internal energy mode. All the specific heats will be assumed constant, a good approximation in transonic flow. The specific heat c2 is assumed to be volume independent so that c =c + c (2.29) Ve Vf 2 c = c + c (2.30) Pe Pf 2 e2 = 2T 2 (2. 31) and consequently 2 =C2 T1 (2.32) e Hence differentiating Eq. (2. 28) yields iaT1 aT2 ua.( -Ul( + C ) (2. 33) j x. ]I pf a(x. 2 ay. jgJ pfX.J ] Substituting Eq. (2. 33) into Eq. (2. 27), the energy equation can be expressed in terms of temperature.

17 U. C J Pf aT1 aT2 1 r a — =- U. C + 3x. 2 3a. p ] a (i2/2) "ax. + dX, - ik a j Ui ak-Sk - +] ( U iJ (2. 34) Using Eq. (2. 30) and (2. 31) and taking the substantial derivative of the relaxation equation (2. 5) yields a- i3U._ U ax I aT2 T U. -R ). aT1 ax. 3 aT2 - U 2 -j a (2. 35) Combining Eq. (2. 26), (2. 34) and (2. 35) yields __ a T C U. _pf i aix a (u2/2) p u. -_a 2 J -Paf ax. 3 - 2 p i3f - f Pf a3. ax. a* k j a c Pe 1 c2 pa 2E -- a(U/2) PU aj [ 3 i - a2 3 -paf axR. J -2Pjf aPq af c ax Pf 3 _ a -k I1 -_a (u_2/2) p i I3x a- k aq ja3rk, a3j j a x. a3 (2. 36) -2 f It can be shown that af e' = yf - 1 (Clarke and McChesney, 1964, p. 216). Hence Eq. (2. 36) can be rearranged to yield

18 V K Cvf 3(U /2) 2 J Yf k] a _ -- U. - u. a. c lax. ij ax. f ax. p ax. P p ax, v 1 ] ] ] e 2 a (U /2) yf- 1 vf =- a _ u + __(_q (2.37) ie - -U au.'c p P U x 3 J k Equation (2. 37) is the governing differential equation of the flow of a viscous, relaxing gas in its most general form. By discarding the heat conduction and viscous terms, the inviscid relaxing equation is recovered in the form -vf au2. ~ ja(2( /2) 22 a' Vf u a [a 2 a (uC / 2) a (2. 38) T - U a. U. a - (2. 38) c i ~x. f ax. a J. ax. e ax. v 1 ] ] ] 3 e In the equilibrium case corresponding to T =0, Eq. (2. 38) reduces to the classic gasdynamics equation u a(2/2)-a =0 (2.39) 3 ax. e ax. 3 3 e For one dimensional flow i. = U;. =x, and Eq. (2. 38) becomes j 3 - - -- 2 2- du Vf - d r;, - 2):U a( (2.40) c dR e d-x e which is the one dimensional inviscid relaxing equation considered by Clarke and McChesney (1964).

19 The three dimensional linear equation (1. 1) can be recovered from Eq. (2. 38) by considering small perturbations of a uniform velocity field u = i such that ooo u. = u + u 1 00 u2 v (2.41) u3 =w where u, V, w, are perturbation velocities which are small compared with i such that oO u -, W ~1. (2.42) u u iu 00 00 00 Then, substituting Eq. (2. 41) into (2. 38) and letting af = afd ae = eO and neglecting higher order terms yields - a (f2 au av 2 1au + U (Mf - 1) i - + (Me2 -1) co ccax c ax ay dZ ( co OX x a- o=0 (2.43) ay az with' =O (cvf /Cv )(af /ae ). The flow can also be shown to cc co"co f0 eco ccc be irrotational to first order; thus, Eq. (1. 1), or the Vincenti equation is recovered if a velocity potential is introduced. Equation (2. 37), although quite general, is very complicated; however, the viscous dissipation t can be shown to be negligible in transonic flow, and the heat conduction term q can then also be expressed in terms of the velocity as is shown below.

20 In a gas having only one internal degree of freedom, the heat flux vector is T1 (int) q. -K - +q. (2.44) where q., the contribution to the heat flux due to the internal energy mode is (int) e (int) 2 (2.45) q ii V ie n D (2. 45) q a a. a ax. a 3 i (int) (int) and =e2 = (l/n) En e and e = a energy state with V. being the diffusion velocity, i the number density and D, the coefficient of self-diffusion (Clarke and McChesney, 1964). If the flow is assumed only slightly out of equilibrium and T2 is nowhere different appreciably from T1 or T2 - T1, then Eq. (2.45) becomes q =nD dT2 dxo aT D - c2 ax. (2. 46) Equation (2. 44) can thus be written as Equation (2. 44) can thus be written as aT q - K(1 + Le) a (2.47) ]

21 or aT q cK 1 (2.48) j c ax. where i D c2 K =K(1 +Le) and Le = — (2.49) If Le < 1, the effect of diffusion upon heat flux can be completely neglected. Otherwise Kc, the modified conductivity has to be used. In the transonic flow of an inert gas, the entropy variation is of higher order; therefore, the relation dT d du can be shown to be valid to the first order so that dT can be replaced by du in the analysis (Sichel, 1962). With this result in mind, the entropy variation of a viscous, relaxing gas in transonic flow will be considered, to see whether the temperature gradient can also be replaced by a velocity gradient in this case. The entropy equation (2. 16) can be rewritten in terms of the total entropy S, and is j a~- T [ a+ -= D ]+ (T - T)] (2 50) U+ p (1 1 a but

22 ae2 i ax = 2 (Ti T- ) jaR. 2- iT2 dS2 = de2 = - 2dT2 Hence, Eq. (2.50), upon rearranging the heat flux term, becomes q. _ >2 AT) aT (T - T2)2 s j] 1 1 1 2 12 The first term of the right hand side is the entropy transport term, or a reversible entropy source. The rest of the terms are the irreversible entropy production terms due to viscous dissipation, heat conduction and relaxation respectively. The magnitude of the various terms in Eq. (2. 52) will first be considered in the interior of a weak viscous shock. The characteristic length is then the thickness of a weak viscous shock Xs= ioo /af s 00 f 00 with e = Mf - 1 (Lighthill, 1956), while the characteristic velocity is af. 00 For simplicity two dimensional flow is considered. All quantities are non-dimensionalized using the undisturbed free stream values except that u and v, the velocity components in the x and y direction are non-dimensionalized by af, while x and y are non-dimension00 alized by XA. However, using an argument due to Cole and Messiter

23 1/2 (1956) y must be further stretched by a factor E, if it is assumed that the flow is irrotational in the transonic regime. Consequently x,-; y 1/2 y (2. 53) As As S S The assumption of irrotationality cannot be rigorously justified a priori. In general, the assumption implies that only the compressive viscous stress due to the gradient in volume dilation is important in the first order flow, the viscous stress due to the shearing of the fluid being negligible (e.g. Rae, 1960). Obviously, the approximations made here are not valid in a boundary layer. The non-dimensionalized quantities u, T1, A", g, C and K are now expanded in p c series of the form LL + c L( + (2. 54) while v is expressed as v =a (c3/2 v +..) (2.55) 00 to be consistent with the irrotationality assumption. Furthermore the entropy s is non-dimensionalized by C so that o0 S = C (2. 56) Poo

24 The slightly out of equilibrium assumption which is in general true in transonic flow implies T2 nowhere differs appreciably from T1 and consequently T2 can be written as T2= T1 + (2. 57) so that the entropy equation (2. 52) becomes aS 2 1 a 2T 1 3 1 32T(2) a2T() () (1) axe Prf a2 + Prf. 2 2 +T ax ax ay ox Yf- 1 u(l) s 2 2 + (f - 1)(1 + )(u ) + s 2 0 (2. 58) J 2 Pr" is the Prandtl number based on the longitudinal viscosity C f P 0 Pr- — = k while X is a relaxation length defined as X =' oaf (2. 59) 00 Throughout the present analysis, the ratio between the viscous shock thickness Xs and the relaxation length X will be assumed to be of the order of s which can be considered either as an assumption or the necessary order such that flow is irrotational; hence Xs/X = 0(E) (2. 60) 5 C

25 Although this assumption is arbitrary it is certainly true for the vibrational energy mode and for dissociation that Xs/X ~< 1 as discussed in the introduction. For the purpose of analysis it is necessary to specify the degree of smallness of X /Xc. The actual value of Xs/X depends on the detailed physics of the relaxation process. The net change of entropy across the shock contributed by the entropy transport in Eq. (2. 58) is zero as uniform conditions are assumed at infinity. Therefore the entropy change across the viscous shock is third order in E while entropy variation within the viscous shock will be of second order in c. And the contribution of relaxation is of third order or higher depending on the order of magnitude of c2. Next, the entropy change in a relaxation wave with characteristic length Xc = F af will be considered. Then the entropy equation (2. 52) becomes ~S 1 /ai f 2 (1) F /a f 2 2T(2) (2T(1) () 1 as 1 0o/ to a T + 00 foo 2 [a( + a T + TMa ax Pr?? X 2 Xc 2[y2 c ax c = rx ay. 1 Yf -f1 au(l) 2 2 + (Yf - 1 )I I +E c2 02 + (Pr? )('xf-) )+e c2/ (2.61) In general, such as in a fully dispersed wave, it is irrelevant to define a viscous shock length, and a viscous length X which has V

26 the meaning of the order of magnitude of molecular mean free path is 2 used instead, XJ/X will be assumed to be the order of e. It can be considered, from Eq. (2. 61) as the necessary order of magnitude such that the entropy production is of higher order; hence XA =O (2 ) (2. 62) c This length was used by Moran and Shen (1966) in solving the viscous shock propagation problem. The entropy change due to relaxation is now of second order. Entropy variation due to transport is of third order while entropy changes across the wave due to viscosity and heat conduction are of fourth order. From the above discussion, it follows that if a solution of Eq. (2. 37), accurate to first order is desired, the flow can be considered as irrotational. The approximation dT-diu is obtained by considering Eq. (2. 34) 2 aT1 3 u T1 - 1 2 u ik lj I I u.C - u —-c --— +c — (2. 34) U CPf ax - ax. 2 p axk pax. p Inside a weak viscous shock, the characteristic length is Xs= X / where e is defined as IMf - 1 I. Nondimensionalizing Eq. (2. 34) according to the upstream values yields

27 2 aT 9 U- X C u 1 -(y-1)U. - c.. (T -T) pf x. x Pf^ a Xf 1/a, ax 2 X 1 2 auk aT1 a X a u. Dx. x. + + + higher order terms (2.63).... Xk... + ~k Pr"p 2 As X A is assumed to be of order E, then to the first order, only c the left hand side and the first term on the right hand side of the equation remain. Rewriting these lowest order terms in dimensional form yields: -2 aT a1 ~1 _ _2 u. C -u. — (2.64) j pf ax. ] ajx On the other hand, it can also be shown that in a relaxing zone where the characteristic length is X, the same approximation can also be reached if c2 is assumed small compared with Cpf or Cpe or c2/C = e. The nondimensional form of Eq. (2. 34) in a relaxing region is then

28 2 aT 1- c2 Pf a = -(if-1)uj (T1 - T2) 1f JI J P auk aT1 A u. x. X ax. v I 1 + -- 1 3 + higher order terms (2.65) X p axk X Pr" p ax. c k c ] To the first order, Eq. (2.64) is again obtained. Since the upstream flow is uniform, Eq. (2. 64) can be written as -2 dT 1 d( -) (2.66) Pf Substituting Eq. (2.66) into (2.48) yields 2 aT. a,,aT a ax. ax. c ax. Pr" aR. ax. 3 J J 1 3 However by using the vector identity u-Vu = V U- + (Vxu) xu (2. 68) 2 ~ ~ and recalling that in transonic theory, Vxu, (Vu-) are of higher order, Eq. (2.67) can also be written as 2_ Qqg u - a u. a Pr U a. ax. + higher order terms (2.69) J J. J

29 Then substituting Eq. (2. 69) into (2. 37) yields -2 CVf a a- ) 2 a TC Ui ax Uj ax. af ax. v i I 3 e Yf - 1 P Pr' j 2_ a U. a - -- U axikca p j axk ak I _ 2 e au. _1 ]1 _2 a(u ) - u. - 3 dX. i (Yf- ) Vf,,t p C5 l PrU, j e 2_ a a 21o p j ax k (2.70) The dissipation 4) above can be neglected since it is of the order 3 of E.

III. ANALYTIC STUDY OF ONE DIMENSIONAL WAVE STRUCTURE In order to understand the nature of flow and to establish where viscous terms are important, the one-dimensional shock wave structure will be studied first. There are several small parameters aside from X /xC or Xs/X in the problem and they are related to various velocity differences such as af - a, u - af and u - ae. In the wave structure probf e ao fno co so lem, only the relative order of magnitude of any one of the velocity differences to the wave amplitude u - US is important. Although u - US is, in general an order higher than u or af as required by transonic theory. Therefore, u - u is always one order higher than u, and af - a, u - af or u - ae is either at the same order or one order higher than u - u. The one dimensional version of Eq. (2. 70) is C f U d- 2 U2 d f - 1+ u 2 eU _ 2 ( +(1 +Pr-V -di (.1 c x dx Prcc oo -2 v dxv e 2f Vf 2_ v dxf 30

31 Without viscous and heat conduction terms, Eq. (3. 1) reduces to the familiar inviscid equation for a relaxing gas. c c - u a- -(a a) (3. 2) cZ d- d e dx v e Lighthill's inviscid solutions for fully and partially dispersed waves were reproduced by Clarke and McChesney (1964) by using the above equation (3. 2) and will be given here. The information obtained from the inviscid solution will then serve as a guide to the study of the viscous equation (3. 1). Some paradoxes in the inviscid solution will also be examined carefully. A. INVISCID SOLUTIONS 1. Fully Dispersed Relaxation Wave A fully dispersed relaxation wave is a steady compression wave in which the dissipative effect of relaxation completely balances the convective steepening effect,and occurs when the upstream velocity lies in between the frozen and equilibrium sound speeds. The problem is then to solve Eq. (3. 2) with the boundary conditions dU X - - = 0 x- - oo u=u; du 0 oo dx (3. 3) dd X - + 0 - -= d s dx

32 u is the velocity upstream of the wave while us, the velocity downstream of the wave, is found from the conservation relation (Clarke and McChesney, 1964) across the wave: oo e e s U s (u = -= (T +l){ (3. 4) uO s Y (y u e +1) s The solution of the above system is readily found and is (Clarke and McChesney, 1964). 2 2 2 22 af - u a - u f oo f s T (u — ) - I in (u u- ) u - u o u - u s 00 S oo S c ve (ye+ 1) C x + T (yf + 1) u + constart (3. 5) vf The constant is arbitrary and determines the location of the origin. Figure 1 shows three solution curves with different upstream velocities with respect to the frozen sound speed af. The upstream portion of the wave becomes narrower and steeper as u approaches af and in the limit, u -af a corner appears, and 00 00 the upstream boundary condition can no longer be satisfied continuously On the other hand, as u approaches aeod the wave profile becomes more symmetrical while the wave amplitude decreases. and in this

33 00 1.0 u o00 u S U 00I X Fig. 1 c case an approximate but useful form of Eq. (3. 5) can be obtained by using the fact that 2 _2. 2 2 af - u = af f x0 f e so co co (3.6) 2 2 2 2 af - u - af - a af f e s 00 00 If uo - Us is indeed very small compared with af2 a2, then Eq. (3. 5) becomes

34 u - u (u - U )(Ye + 1) In l -s, (3. 7) U -U 2 - x0 Vf 2 2 - (af - e V o0 00 e Equation (3. 7) has the same form as the solution for the structure of a weak viscous shock with the quantity IB given by c Vf (af 2a 2) (3.8) V 0O 0O e playing the same role as the compressive viscosity. Consequently B is sometimes called the bulk viscosity. Thus the relaxation effect of the internal mode is equivalent but different in order of magnitude to the action of compressive viscosity. Fully dispersed waves are at least an order of magnitude thicker than the equivalent viscous shock wave. 2. Partially Dispersed Shock Wave When the upstream velocity is greater than the frozen sound speed, the solution (3. 5) is no longer valid as it then possesses an infinite slope and is double valued at each location (Becker, 1968; Polyakova, Soluyan and Khokhlov, 1962). In the framework of the inviscid theory, this difficulty is resolved by inserting a viscous shock discontinuity upstream of the relaxation

35 zone within which the internal mode is assumed frozen. The detailed structure of this viscous region is ignored. The solution downstream of the shock discontinuity is (Lighthill, 1956) u =u + (Ub - se) e (3.9) where ub is the velocity downstream of the viscous shock as determined from the Rankine-Hugoniot relation: 2 2 2 (" -a ) 00 u -Ub (y+ l) U (3.10) Furthermore, it should be pointed out that the viscous shock should be sufficiently strong so that the majority of the velocity drop occurs across the viscous shock if Eq. (3. 9) is to be valid. This type of wave has been called a partially dispersed wave by Lighthill. Therefore, the solution for a partially dispersed wave does not hold for the case when u is greater than but extremely close to af. The complete partially dispersed shock wave solution is shown in Fig. 2 below. The discontinuity in the solution for a partially dispersed shock wave is clearly the result of neglecting viscosity in the theory. Mathematically speaking, the proper characteristic length in the upstream portion of the wave should be the viscous length X = v /U or the viscous shock thickness s rather than the relaxation length S

36 u 00 Viscous Shock Discontinuity Relaxation Zone Downstream of the Viscous Shock u1.. x Fig. 2 A = T u. The existence of two characteristic lengths as above C 00 0o suggests a singular perturbation problem. In view of the difficulties arising in the inviscid solution, that is the corner solution for a fully dispersed wave when u = -f and the discontinuities in the velocity and velocity gradient in a partially dispersed wave, the viscous equation (3. 1) will now be used to study the wave structure. It is also hoped that the one dimensional analysis will shed some light on more practical problems such as flow inside a nozzle or over an airfoil.

37 B. VISCOUS SOLUTION The problem is to solve Eq. (3. 1) c f - d 2 2 du r c u d- (a= - u ) d-7 c I/ f 2 v L e d2u _2 _ 2 du + A v" u 2 = (u - ae ) d-= =2 e d dx - + B i u dx" (3.1) with the boundary conditions X — - 0 U =U oc - + oo00 U =U 5 diu; d=0 )? ~dx (3.11) du 9 CIR where f - 1 A= 1 + =0(1) (3.12) yf - 1 Vf B = + Pr c (1) V e As mentioned before, there are two characteristic lengths in the problem. One is the relaxation length defined as Xc =F u, 0-af, or T ae depending on the nature of the problem. The other length is associated with the viscosity and is defined as

38 X = v t/u" or V "/af. The characteristic velocity will be u V 00 o0 00 0 00 for the partially dispersed wave, and for the fully dispersed wave either af or ae Odepending on whether the flow is frozen or equilibrium transonic (Napolitano, 1966). The distinction between frozen and equilibrium transonic flow can best be illustrated by rewriting Eq. (3. 1) in the form Kf u =K Ku (3.13) f e where Kf and K are the frozen and equilibrium transonic operators respectively defined as K -vf - d 2 2 d1 dxcf v e e (3.14) -2 2 di! K = (U -a ) d_ e e dx in the inviscid theory, and as c vaf vf - d2 d F d Kf = u- 2 _- u2) x -+ A'" ud2' f c dxK[ dx _2 v e (3.15) -2 2 d d2 K (U _ a 2) + B e e de dx ~+ B P" u~2

39 in the viscous theory. In the frozen transonic regime where 2 tMf - 1 = O(E), an approximate equation with the non-linear terms from Kf u plus a linear term from Ke u can be obtained by e stretching the x-coordinate with a factor iMf2 - 11. Similarly, an approximate equation valid in the equilibrium regime where iMe 2 - 11 =O() with a linear term from Kf u and non-linear terms from K u can be obtained by stretching x-coordinate with a factor e 1 / IM 2 - 1. However, since the sonic point is always a singular Oo point in the inviscid equation, Napolitano's stretching is not valid for 2 the case Mf = 1. Therefore, the stretching of the x-coordinate 00 will be different in the viscous theory where X /X = O(e2) is used V c as a stretching factor when Mf2 = 1. It is also expected that the viscous equation thus obtained should also be true for cases where IMf2 -1 = 0(). 1. Equilibrium Transonic Approximation-The Bulk Viscosity Case The equilibrium transonic case will be studied first. Nondimensionalizing Eq. (3. 1) by the characteristic speed 5e and length X ='ae where - =T (Cvf/Cv ) yields ~ c cO cO oo o v e 2 A 2 u [(af 2 u2)u 2 du V d u (2 a v d -U ) +A (_ + B U d 2 (af Vu 2... - u ax.. (3. 16) dx c dx (3.16) where X = v /ae <<. v oo oC c

40 The precise order of magnitude of X v/ is unimportant in the equilibrium transonic regime so long as X /X < 1. Apparently the waves are thick and gradients are sufficiently small that viscous 2 effects are negligible. However X /Xc =O(2 ) will be assumed here in order to be consistent with other cases. u is now expanded in the parameter 6 which is of the order of magnitude (ui - aeo/ae i. e. represents the wave amplitude, and 6 << 1. Then u= 1 +6 +62 u2 +... (3.17) _) 62 (2) af af +6 af (1) a +. 00 00 (3.18) a =1 + 6 a (1) 62a (2) + e e e u - a oo e 6= IM -11- 00 (3.19) e a so e oo The coordinate x is then stretched, following Napolitano, by =x6= 6 = (3.20) The fact that Xv/6 is a thickness of a weak viscous shock suggests that X /6 might be considered as a thickness of a weak fully dispersed wave. Keeping only the largest terms, Eq. (3.16) becomes

41 d 2 6 - (af 1) 2 du() - 1)6 dc x +A C 3 d2u(1) dx = 2 (u(1) - a(1)) e 3 du(l) 6 adx d1x xu 3 + B C C d2 (1) du2 A2 dx (3.21) 2 Since X /X~ O(e ) here Eq. (3.21) further simplifies to Sn c 2 (af 2 oo d2u(1) (1) - 1) 2 (u dx a(1) du(1) e dx (3.22) In transonic flow, it is readily shown that, to first order, u(1) - (1) e +1 _e (1) _2 u (3.23) Thus following Clarke and McChesney (1964, p. 218) 2 a e = IQ - U) U (3.24) where 2 F =P + po 00 00 Q = p u *o 00 At x -- oo a =a e e oo u= u oo and then

42 F Q a e c 2 -2 + u o o Y / e 1 +Y =I (3.25) y u e oo Introducing the expansions (3.17) and (3.18) yields (1) 6 u 1 +28.a = 1 -' 6 u1) e e + 6 u() + 0(62) and finally (1) _ e (1) e 2 to first order (3.23) Similarly, it can be shown (1) (1) f (1) u - af = -2 u to first order (3.26) Equations (3. 23) and (3. 26) will frequently be used below. Using Eq. (3. 23), Eq. (3. 22) can be written as 2 (af o d2 (1) -I d ( e dx +1) u(1) du() + 1) u dfic (3.27) Integrating once yields the equation (a 2 o0 du(1) +1 2 _Ye (1)2 - 2 =c (3.28) From boundary conditions

43 X - - 0 U=U 00 u(1) 1 U -a oo e oo as = 6 a e 00 so that c = - (e + 1)/2. The solution for u(i) is then readily found and is (y + 1) e A () (af - 1) (1) = I+u oo 1 - (3.29) Equation (3. 29) can be re-expressed in terms of the original physical variables A X x= - 6 T a oo e oo u - a u-a (1) eo 1 u = a o e 00 (3. 30) u - a oo e 00 a e 00 so that to first order in 6. (7e +1)(u0 - u) x U-U S u - u s e 0a - oo 2 - 2 2 (c /c ) f 2 - a 2) f e oo oo (3. 31)

44 This solution is identical to that obtained by Clarke and McChesney (1964). Thus Clarke's inviscid solution for the bulk viscosity case is completely recovered in Napolitano's equilibrium transonic approximation. Furthermore, the physical nature of the solution (3. 29) can best be shown by rewriting it in the form (ye + 1)(G -U ) u() - tanh e (3.32) 4MB This solution can now be seen to be identical with Taylor's weak viscous shock solution with the bulk viscosity replacing the longitudinal viscosity. Since this classic solution for the bulk viscosity approximation has been derived from the viscous equation it appears that the effects of viscosity can be neglected in the equilibrium transonic regime. 2. Fully Dispersed Wave Lighthill's solution, Eq. (3. 5), for a fully dispersed wave is good for almost the entire velocity range ae < u < af - However, a corner appears in the solution as the free stream velocity approaches the frozen sound speed from below and there is then a discontinuity in the velocity gradient at the upstream boundary. The viscous equation (3. 1) will now be used to study this problem. The governing equation is again Eq. (3. 1)

45 Vf _ d 2 _2 du d T - (af - u +A v't U 2 v dx e 2, 2 2 dii d21 u-2 a 2) du + B d u d (3. 1) dx with boundary conditions d5 X~ —oo U=; = 0 o o dx (3. 33) di! X~- _- - dux —+ou; u~=0 + =Us dx Lighthill's solution for a fully dispersed wave, obtained by simply solving Eq. (3.1), (3. 2) but without the viscous terms, has already been shown to be -22 _2 -2 af -u a - f 0o f s T 00- _ In (u - u-) - r s- In ( - u ) u u oo u u s o0 S 00 S 1 Ve _ _ (e +1) c x + (Yf + 1) u + constant (3. 5) vf For the case when af 2 2- a 2 becomes small as indicated by oo ^ fO0 2 2 Clarke and McChesney (1964), the scale parameter i (af - u )/(f - - s) for departure of u from ui in the early parts of the wave will be very small and the variation of flow properties correspondingly rapid.

46 Therefore, it is doubtful that the inviscid solution will remain valid. Furthermore when af = u-I a corner appears upstream. When afo = u Eq. (3. 5) becomes 2 - 2 af -u c s v e - T. l-n In (u- )= (y +1) x u -u s 2 e c oO s Vf { (yf + 1) u + constant (3. 34) and boundary condition (3. 33a) cannot be satisfied, as u = u occurs at a finite value of x and a corner appears. The slope of the velocity profile at the corner can be determined by rewriting Eq. (3. 2) in the following form. c (af d- 2 (e + 1)(u - - us) v e and is di- (1 e +1 )(U -Us) (3. 35) dx 4 cv rv e The failure of Eq. (3. 34) to satisfy the upstream boundary condition is apparently due to the neglect of the compressive viscosity terms near the upstream region. The compressive viscosity terms

47 in Eq. ( 3. 1) should be of the same importance as the relaxation terms in a small region near the upstream portion of the wave. Equation (3. 1) is first non-dimensionalized by the relaxation length X =T I af, where T T =TOo(vf/Ve), and by the characteristic speed af so that oo d, 2 2 du 2 d 2 u f 2 dxu = - u2 + B (3.36) L c d2 J c dx Again it is assumed that X V/ O(e2) so that, keeping only the largest v C terms, Eq. (3. 36) becomes u d [(af2 2) du] (u2 a 2) du (3. u- L a -u -a (3.37) dx- f dx e dx This equation will be called the outer equation and is just the inviscid relaxation equation. The solution in terms of the outer variables is then 2 2 f s 00 s This solution will not be valid near the upstream part of the wave where another equation is needed. In order to keep both the viscous and the relaxation terms in the equation, the first three terms in Eq. (3. 36) will have to be of the same order. If 6 is defined as 0(1 - Us), it follows that the x-coordinate is to be streched by 5,

48 /Vlk/X (1 - u ) = /61/2 and the perturbation in u is (c61/2). V C S However, for definiteness, the particular case e/6/20(6) will be considered, but this choice is not essential. Thus x and u in the inner region are = x/6 (3.38). = 1 +6 (1 -u ) u.( +.. (3. 39) i s i Substituting Eq. (3. 38) and (3. 39) into Eq. (3. 36) and collecting the lowest order terms, then yields the inner equation d[d2 2U.()d u. (1) d u. (1) u (I) 1 A,A 2w,3du i du.2.... A - 2r! u.....= r (34O ) dx d2 f I d j e dr where Yf +1 y +1 r; J and r e(3.41) f 2 e 2(341) Equation (3.41) can be readily integrated once to yield d (1) du (1) d u. A - 2f u( = r u.U. (3.42) fi d e i

49 The constant of integration is zero by applying the upstream boundary condition. In the absence of viscosity, Eq. (3.42) has the solutions u.(1) = (3.43) or re 1 r x + 1 (3.44) Clearly, Eq. (3. 43) is a trivial solution but solution (3. 44) describes the upstream behavior of the fully dispersed wave for the case af = oo 00 Unfortunately, Eq. (3.42) cannot be further integrated analytically; however by letting f = - u., Eq. (3.42) becomes Ad2+f df A —+ 2r f- r f=o (3. 45) dx which, except for the numerical coefficients, is exactly the same as the equation studied in detail by Cole (1968) in connection with what he refers to as a corner solution. The phase plane of Eq. (3.45) has been plotted by Cole and various possibilities for matching have been outlined. The following phase plane plot is reproduced from Cole (1968)o The asymptotic behavior of the inner solution can be readily

50 re/2rf Integral curve for inner solution f Fig. 3 determined. Letting p = df/dx, Eq. (3.45) becomes dp Ap + 2 f p - r f=0 (3.46) or Ap dp f r -2p r df Integrating once then yields A [-p - ln (e - 2rf )]= 2f + c (3.47) f~~~~ f and since

51 A r x —oo f=0, p=, c, = n r 2 e 4rf e Then Eq. (3.47) becomes r 2 r p- ln2 l - ( Ip)] 2 (3.48) ~f f e Upstream u —u oand both f and p are expected to be small. Than the logarithmic term in Eq. (3.48) can be expanded in terms of a power series of p to give r 2F 4r2 A pe( f 1 f 2 + O f2 2r1 -P 2r rP2 2(P fL f e r e F e so that A2 2 p = f + higher order terms (3.49) e Taking the square root of Eq. (3. 49) yields df p=~/r7A f (3. 50) The solution of Eq. (3. 50) is readily seen to be f=c/A e 2

52 From the boundary condition x- - cc f = 0 and since f = - u. the final result is u. =- c2 e (3. 51) The constant must be determined from matching with the outer solution. Thus the asymptotic solution of the inner equation approaches the upstream velocity exponentially and also satisfies the boundary conditions upstream. Equation (3. 51) can also be expressed in terms of the physical variables to give /(f e )/af 00 00 00 XK u exp 1+,, ~-Ta 1 Pr"' a ~oo f /Praf co fvce o00 e and after some rearrangement yields (1) X u =exp Yf iI (1)=e L (......1) J (3.52) Equation (3. 52) clearly shows the combined role of compressive viscosity and relaxation in the upstream part of the wave. It is interesting to note that the compressive kinematic viscosity is unrelated to the viscous shock amplitude; whereas the bulk viscosity

53 depends on the maximum strength of a fully dispersed wave. This difference explains why an increase in IB steepens the wave and an increase in v" smooths the wave. At the downstream end of the inner region as f -o, a different asymptotic representation is also available. Cole's phase plane plot suggests that p is a monotonic increasing function as long as p < r /2rf and as f - o p — e/2rf. Near the downstream end of the inner region, the solution is then found by solving the equation df r p =- (3.53) dx 2 f f r e f-2 r +c3 or.(1) e 1. u 2r( x 3 (3. 54) Equation (3. 54) is the asymptotic form of the inner solution near the downstream end and is to be matched with the outer solution. It is interesting to note that Eq. (3. 54) contains no parameter related to the compressive viscosity.

54 The idea of matching, according to Van Dyke (1964) or Cole (1968) is that the behavior of the outer expansion as the outer variable approaches the inner limit is in agreement with the inner expansion as the inner variable approaches the outer limit, i. e. there is an overlap domain where both expansions are valid. Therefore it is desirable to express both the inner and outer representations in an intermediate region where the intermediate variable is defined as u - af 65(1) 00 D u -u. Da-u) as D —O; -- 1 oD a l-i) a 00 (3.55) 6 ~ x =DX It is more convenient here to consider u rather than x to be the independent variable. Then the inner limit become u — O while outer limit becomes u -- o. The inner asymptotic representation in the intermediate region is readily found from Eq. (3. 54) and written in terms of the intermediate variables becomes re Ux =- x + c 3 (3. 56) 7 3 where c3 must be obtained from matching.

55 The easiest way to obtain the asymptotic representation of the outer solution in the intermediate region is from the outer equation rather than from the outer solution. Equation (3. 37) can be integrated once to give 2 2 du 1 (af - u) - = (ye + 1) (u - us) (u- 1) (3. 57) In the intermediate region u- af oo u- 1 7r D a D (1 - ) co af af= (1 + D(1 -u) af ) f ra then u7 - af = r u. Equation (3. 57) becomes du f dx - r u (3.58) This equation is readily integrated yielding r e u Inj 2r X7 + c (3. 59) The constant c in Eq. (3. 59) is arbitrary and determines only the location of the origin. If the origin is located so that u = u at x = 0 ~then c =-=~ O~ and ~0o then c = 0 and

56 F e u = -- 2r xv (3.60) Comparing Eq. (3.60) and Eq. (3. 56), the result c3= 0 follows. Thus the inner and outer solutions are matched to the first order; however due to complexities of the equations higher order terms in the asymptotic expansions were not evaluated. In general, since both the inner and outer equations are autonomous, each solution will have a constant which establishes the origin of the coordinates. However only one of these constants is arbitrary. In the present case, the constant in the outer solution was fixed by requiring u = u at x = 0. Then the constant in the inner solution was found from matching. A qualitative sketch of the solutions is given below. Mathematically speaking, the theory is now complete. U f -0(6)0 - 0(1) -I Uco afoo Inner Solution Outer Solution Common Part s______ |" - )/(u s) = (re/2rf)(x/xc)' —-----------' x —---— f c Fig. 4 A v

57 However the complete inner solution still can be found only by numerical integration. Equation (3. 52) provides the initial values needed to start the integration. An uniformly valid composite solution can now be constructed by adding the inner and outer solutions and subtracting the common part which from Eq. (3. 54) is in the present case 2r (-f ) x - r.f(1-u) (3.61) e af(us) 00 The final result is - 2 _2 a -u_ _ af S U-U 5 S x- = r In _ _ + 1) (ue u - u -u o00 S 00 S 2r (U - a ) f 00 r( -- i (3.62) e ( oo-Us where u. is the inner solution. 1 3. Partially Dispersed Wave The partially dispersed wave corresponds to the case when the upstream velocity is greater than the frozen sound speed afoo By considering that the relaxation due to the internal mode is a much slower process than the translational diffusion which accounts for

58 viscous effects, it can be assumed that the internal mode is frozen inside the viscous shock and then relaxes slowly to the final equilibrium state required by the conservation laws. The viscous shock is usually represented by a discontinuity in the inviscid theory (Lighthill, 1956; Clarke and McChesney, 1964). For the viscous problem, the governing equation is again (3. 1) _ f - d 2 _2 du- d u C-u = (af - u ) id + A v" u2- c dx - 2 v dx e 22 -2 dii d (3.1) = (U2 _ a + B v" U 2u (3.1) dx with boundary conditions d2oo dx' 2 dx (3.11) du du x- = d+ =0 d; u dx s d3R 2 where ui, the velocity at downstream infinity, is related to u by S o0 the equilibrium Hugoniot condition - 2 2(u2 e 2) u - u (- +1).... (3.4) 0o s (Ye + 1)u e c~~

59 Utilizing Eq. (3. 11) after integrating Eq. (3. 1) once yields T C-A V" U 2 +(af 2- u) dX (ye + l)(u - u)(u - ) (3. 63) v dx e Equation (3.63) is nondimensionalized with characteristic velocity uii and length X = T (C /f ) u oand takes the following form c Vf v e X 2 d u(a2 2 du I A u + (a - + 1)(1 -u)(u - us) (3.64) c dx where X iy "/j A v fu _ Q(E2) c Vf e by assumption. Then the outer equation is obtained by letting e - 0 in Eq. (3.64) and is again an inviscid one 2 2 du 1 (af -u2 - ( + 1)(1 -u)(u - ) (3.65) This equation is expected to satisfy only the downstream boundary condition and an approximate solution valid downstream of the wave is easily found, after integrating the approximate equation, 2 2 du (af -u ) r (1 )(u -u) f s dx e s) s

60 to be (1 - U) S -re 2 2 af - u u = s+ (ub - us)e S s (3. 66) where ub, the velocity immediately downstream of the frozen viscous shock, is obtained from the frozen Hugoniot condition. Apparently, solution (3. 66) is not valid near upstream infinity where a different equation is sought. The inner equation for a partially dispersed wave with a very weak viscous shock or the case u - u.... b = 0(6), U -u 00 S will be derived first. This can be done by considering a small region around ub of order 6 compared with the total wave amplitude (1 - us). Velocity u and sound speeds will then be expanded in the series u = b + 6(1 - s)Ui (1+.. (3.67) a a + 6(1 - u)a( +.. Substituting Eq. (3. 67) into Eq. (3. 64) yields

61 2 2 X du.s 2- 2 2 du1) A h(1 U2 + I (1 - 2 fus i c s dx s u- afbb afb ub l-U~b -0~ D~ =~ 2 b (3.69) S f thermore, using the frozen Hugoniot relassumption - 2 2f( 2 2 a a( 1 - u - lU0: ub ub f (3.70) 1r (1- u.(1)) D -US =b u (-2 r S e 1 t-u The x-coordinate will be stretched as in the case of a fully dispersed Eq. (3.68) can be rewritten as = by (3. 70):r ( i ( 1 - u s wave when U = -a by af oo

62 6rd = EV 6 c(1 - u ) 2/3 2 or 5 = 2/3 if (X/A) = e and (1 - u) =6. Then x = x/6. V C S Equation (3. 70) becomes d2 (1) ( du. 1) u dx (1) d'(1) Ub s A -+ (1 - 2() = - u>) (3. 71) d2 2i 1 di e lu where (ub - us)/(l - u ) is, in this case, of order of one. In order to give more physical insight, A / c(l - u ) will be written in a slightly different form, \ v \? V o0 0 0 0 (3. 72) A(1-u) c 2 c S V u - u V U - f_ 0f 5 -i o T _U 0- TV - - - Cv o u c co v u r e e oo e If indeed u is close to a as assumed then 00 I0 A IJ "/L V 0_ f010 er (3.73) c S) af 00 e Therefore, XA / (1 - us) can now be interpreted as the ratio of the two characteristic lengths associated with the longitudinal and bulk viscosity aside from a constant r. e

63 The procedure used to solve for the fully dispersed wave with U = =a can also be applied here by letting co =foo 1 - 2u() = f (3. 74) Equation (3. 71) becomes df df A + rf-= rE(f + 1) (3.75) dx2 fd e where E has been used to represent (ub-us)/(1 - u ). This equation is again very similar to Cole's corner equation except for a constant term r E. Its phase plane equation can be e easily written, by letting df =d p (3. 76) to yield dp f-ffp+ rE(1 + f) Idp f -- (3. 77) df p A qualitative sketch of the phase plane plot with three different E is shown below

E = 3 E = 1.5 _E = 1/4 Fig. 5 Upstream infinity, u = u, (du/dx) = 0, or f = -1, p = 0 is again a saddle point as in the fully dispersed wave. From these plots, the role of Eq. (3. 75) is now clear. The terms on the left hand alone is just the equation describing a weak viscous shock while the terms on the right hand side is due to relaxation of the internal mode. Near upstream, the terms on the left hand side dominate the equation unless E is large,while near the end of the inner region, the viscous term is no longer important and the inviscid terms dominate the equation-an indication of the matching solutions of Eq. (3. 71) and the inviscid outer equation.

65 Although, Eq. (3. 71) is arrived at by assuming u -U 00 S and U- Us u -b = E = 0(1) 00 ub it can be extended to solve for general cases of partially dispersed waves by simply considering that E is a parameter for the equation and allowing its order of magnitude to vary. The phase plane plots indeed indicates that is true. The extension of Eq. (3. 71) is apparently due to the fact that varying the order of magnitude of E does not affect other terms. For the case of a reasonably strong viscous shock or Ub - u _- -s = 0(6) U -U 00 S and - u- =0(1) U - u 00 S the term on the right hand side of Eq. (3. 71) can, in general, be neglected because it is of higher order. However, dropping of this

66 term makes it impossible to match the inner and outer solutions since the relaxation effect has then been completely left out in the inner equation. Therefore E should be retained as a parameter as discussed above even if it is of higher order if the present method is to be used. However Eq. (3. 71) will have to be derived in a slightly different way if it is to be used for the case E =0(6). The expansions of variables will be different u= b + 6u(l) + (3.78) a = ab+ 6a.(1 +. It is clear that the velocity in the inner region for the present case is of order 6 (the same as the total wave amplitude) contrary to the previous case when the velocity variation in the inner region is of order 6(1 - u ) or 0(6 ). s The stretching factor 6 is defined as A A -u) = O(E) (3. 79) Because of the reasonable strength of the viscous shock, the thickness of the viscous shock s = (X /(1 - ub)) is used as the characteristic length in the inner region. The the inner independent variable is

67 x = x with 6 = E (3. 80) E Substituting (3. 78) and (3. 80) into Eq. (3. 64) and neglecting all the higher order terms except the one expressing the relaxation effect, i.e., Ub Us U -ub oo D which will then be treated as a parameter as discussed before, yields d2 u(1) du (1) (1) (1) A 1( + r(1 - 2u.) - = 1(1 - ) E (3. 81) dx2 f 1 dx e 1 dx Equation (3. 81) is identical to Eq. (3. 70). It is also interesting to note that if X /A is fixed to be of order e, then the length of the inner region for the case when (tulo/ub)/(u - us) = 0(1) will be of 2/3 order E / but of order E for the case when (i - ub)/(u - u ) = 0(6). This is indeed true as the stronger the viscous shock strength, the thinner the shock. Solution of Eq. (3. 71) or (3. 75) by following the same procedure used in the fully dispersed wave with uoo = afoo will now be considered. Near upstream infinity where f = -1 an approximate form of Eq. (3. 75) can be obtained by letting f + 1 = g clearly g << 1. Equation (3. 75) becomes

68 A - 2 rf d = r gE 3(3. 82) dx2 fdR- e dx The solution of Eq. (3. 82) is rf v r. +4rFE 2A — x g = c1e (3, 83) c1 is the constant of integration and has to be determined by matching of the full inner solution and the outer solution. It can be easily seen from Eq. (3. 83) that a combined effect of viscosity and relaxation determines the structure of weak viscous shocks in a relaxing gas contrary to the fact that viscous effects alone determine the structure in an inert gas. However, as the viscous shock strength increases and E decreases, the relaxation effect in the upstream portion of the shock diminishes and vanishes in the limit of a strong viscous shock. At the downstream end of the inner region where the viscous effect is negligible as found in the phase plane study, as f becomes large compared with unity, or f - o, Eq. (3. 75) can be approximated to yield r df e d eE (3.84) f

69 r f =-EEx~+c asf-+oo (3. 85) r U(1) =-2- E +E +2+ c 2]? 2e E2 x+ +c(3. 86) The asymptotic form of the outer solution needed to match Eq. (3. 85) or (3. 86) is best obtained by solving the outer equation in an intermediate region where x 5 ~ D X= D-;- -00 c x D D c (3. 87) - b 1 6 (1) 77 ub D D The outer equation is then, in terms of the intermediate variable, du 2r u dx = - r3 E ur (3. 88) f 7 dx 3 ru The solution is r e u = - Ex (3. 89) 2rf Matching Eq. (3. 89) to Eq. (3. 86) which is written in terms of the intermediate variable F e 1 u = 2-Ex + +c) (3.90) 1 C —' 9 2

70 A qualitative sketch of the solutions is shown below u 00 0(6) U S Fig. 6 The exact matching of inner and outer solutions of a fully dispersed shock wave with Mf = 1 is shown on Fig. 10. The inner solution is obtained from the numerical integration of inner equation (3. 42) while the outer solution is from equation (3. 34). All calculations are based on yf= 7/5 and ye= 9/7. This particular solution has been compared with the result from the numerical integration of the full equation (4. 1) and is in very good agreement which is shown on Fig. 11.

71 A different method which gives more physical insight is now employed to solve for the partially dispersed wave with reasonable viscous shock strength. 4. Partially Dispersed Wave-Alternate Approach A method which shows more physical meaning will be given next for the solution of the structure of a partially dispersed wave with reasonable viscous shock strength. The disadvantage in the previous method in solving this problem lies in that there is no analytical composite solution available. However, the method which will now be described- a physical approach, although lacking mathematical justification, does give an analytical composite solution not only with a very simple form but also shows clearly the relation between the viscosity and the relaxation effects. Since ideas from the method of matched asymptotic expansions will again be used to solve the viscous equation, it is useful to rewrite the discontinuous inviscid solution in a slightly different form. Thus the inner solution, or shock discontinuity can be written u. =u + (ub - u)H(x) (3. 91) 1 00 00 where H(x) is the Heaviside function defined as H(x) =0 xO < 0 (3.92) H(i) =1 x > 0

72 The outer solution is then in the present form x 7T (1-e) u 0 -I (3.93) The zero order composite solution is 1 = Db (3.94) or x u - u( 1 - e (-ub - u)H - H) - Us)(1 - e) s U=U + 00 ] (3.95) The viscous problem is to solve Eq. (3. 1) c V c v e - d u - dx [(if 2 - 2 du d2u - u )a+ A'u -2 (3.1) ( 2 - 2) e _- d u2 du BPu,, - dK+B ft'2 d~i with the boundary conditions X -0 - 00 U = U 00 X - + 00 U = U S d= * = 0 dx Ca and u > af 0oo f9

73 Non-dimensionalizing Eq. (3. 1) by the characteristic velocity afoo and length X = T'o with' = (Cvf/e) yields d [(af2 2) du A d u d dx c dx - (3. 96) X 2 = 2 du v d u =(u- ae) + B -u 2 c dx 2 where X /A is of order of E as assumed before. By letting E -0; (X/ /X) -0, the inviscid equation is recovered and will be called the outer equation u (af- u ) dx =(u2 ae) du (3. 97) ud (aj fe2u2) ] du However, the outer solution cannot satisfy the upstream boundary condition as the outer equation will not be valid in a region where the viscous length X, rather than the relaxation length X is the proper characteristic length. In this inner region where the characteristic length is AX, the appropriate variables are X/f c X X" - 2 (3.98) = E U(X) = u(x)

74 Substituting (3. 98) into Eq. (3. 1) and collecting the lowest order terms yields the equation 2 2u d d) A d2 = 0 (3.99) dK2 with the upstream boundary condition d2 — ~o = u -0 0 (3. 100) dX - Integrating Eq. (3. 99) once yields 2 2,2di d d (af - u + A u 2= e and with boundary condition Eq. (3. 100), c = 0, to yield 2,d 2 2d ii d (a - u2) d+ A U 2 0 (3. 101) In the transonic regime, this equation is often called "Taylor's weak shock equation" and its well known solution is u -(u (f+l)(U0-ub)x (3. 102) =e u-ub or in a different form u = OO - 1 (u - ub) [1 + tan h4(yf + 1)(u - ub)xJ (3. 103)

75 The outer equation d (a2 du ](u2 a2 du (3. 104) will now be examined in detail. This equation has exactly the same form as that for a fully dispersed wave, and expresses the balance between the convective steepening and the diffusive effect of the relaxation. However, the nature of the solutions for the fully and partially dispersed wave are quite different. For a fully dispersed wave, the solution of Eq. (3. 104) is either a shock like profile, i. e., Eq. (3. 5) or the trivial solution, u = u. For a partially dispersed wave, the translational temperature rises from T1 = T to T = Tb through the viscous shock, while the temperature T2 associated with the internal energy mode still lags behind and is almost equal to the upstream value due to the assumption that o fc _ 2 7 I= - O(E. co foo In the limit, the slow mode is frozen inside the inner region, so that the internal mode temperature T2 remains constant at T. T2 still equals T immediately downstream of the inner region where the equilibrium value should be Tb with Tb > T. It is this difference between the equilibrium and frozen value of the internal mode

76 temperature, which depends on the viscous shock strength, that drives the relaxation wave downstream. In other words, T1 - T2 is zero ahead of the inner region and suddenly jumps to a finite value equal to (Tb - T ) behind the inner region. The inner equation as formulated here is independent of the outer equation. This contradicts the familiar feature of the method of matched asymptotic expansions in which the determination of an unknown constant in the inner solution depends on matching with the outer solution. In the present case the inner solution serves as a forcing function which acts on the outer equation over a very small but finite region. Equation (3. 104) can also be written in the following form since both the upstream and downstream flows are assumed to be in equilibrium, 2 2)du 1 (af - u dx= - 2 (e + l)(u - us)(u - u) (3. 105) This equation is not expected to be valid anywhere near the singular point u = af~ without modification. Furthermore, it does not have the information of how it deviates away from u which has to be supoo plied by the inner solution. Physically the relaxing flow considered here is similar to one dimensional flow with heat removal as can easily be recognized by examining Eq. (2. 26) in one dimensional form

77 but without the viscous term: de de2 1 [ 2 2) du1 oru dx=-1 f -u or (3. 106) de2 2 du (f - 1) dx dx 2 2 (af- u ) (3. 107) de2/dx can be considered as a heat sink and is zero upstream of the inner region but suddently jumps to a finite value which then drives the wave downstream of the viscous shock discontinuity. A qualitative sketch is shown below x of Shock Lighthill Solution Viscous Shock Actual Behavior Fig. 7

78 It is the deviation of u from ui within the viscous shock that 00 initiates the relaxation process and will determine the initial behavior of (de2/dx), as shown in the sketch of Fig. 7. Thus the inner solution initially drives the outer equation suggesting that an improved outer equation could be obtained by replacing (u - u) in Eq. (3. 105) by u - u = (u - Ub) [ + tanh b + 1)(u - ub) (3. 108) corresponding to the variation within the viscous shock. Then the modified outer equation becomes 2 2 (af -u) du (u% -ub) dx (e +1) 1 =- e (u - u)[1 + tanh((7f+ l)(u - u) (3. 109) If indeed, the viscous shock is not too weak, it is possible to use the approximation a ae; ub -u (3. 110) s s so that Eq. (3. 109) can be further simplified to yield du 1 1 d x2u (u -u - 1>[ + tan+ 1 - 111' (3. 1) s ~L 0

79 with boundary conditions x —0 u = ub (a) (3.112) x - +oo u = u (b) boundary condition (3. 112a) can be considered the result from matching with the inner solution. Then the solution is readily found to be u-u 2 2 S E = exp x +, ub s s (f +l)(u - us) 00 S (3. 113) (Yf+ 1)(uu) \In cosh - x 4E2 Letting 2 -1 - E E (3. 114) 4(7f+ l)(u -us) Equation (3. 113) can be written in the simplified form E 2x 2us u-u E s (l+eE b T - u- 2 ) +c (3. 115)

80 from boundary condition (3. 112a), it follows that E 2u S c C = (3. 116) The final version of the modified outer solution is then E 2x s u -s / E \ = 1+ e X 2x =u = (ub-u)(l + e (3. 117) or E 2u s (3. 118) The composite solution, obtained by adding Eq. (3. 118)- and (3. 103) and subtracting the common part ub, is F 11l U = Ui + u ub = (uoo ) [ - { 1 + tan h (f + 1)(u -lab) (3. 119) E 2u 2x s +u + (%u s) +eE ) $

81 The qualitative sketch of the inner, outer and the composite solutions is shown below S u 0(6)5 I ~~iL__Inner Solution 0(6) Composite Solution 0(E) Outer Solution 5 l I i = _IX us - X c Fig. 8 The composite solution and the exact matching of inner and outer solutions of a partiallydispersed wave with Mf = 1. 2 is shown on oo Fig. 12. The comparison between the analytic solution (3. 119) and the result from the numerical integration of the Eq. (4. 1) is shown on Fig. 13. The two solutions agree well except in the middle part of the wave where a small discrepancy exists. This is apparently due to the fact that the data used in getting the analytic solution (3. 119) do not quite satisy the approximation Ub-u for in the present case (u- Us)/( u - us) 1/4.

82 A diagram which shows various regimes depending upon the relative order of magnitude of ui - and u - and sumaco foo oo eoo marizes the various solutions is shown below. U0O Composite Solut idea that the vis relaxation proc< ion is (3. 119) with the cous shock drives the lss downstream Frozen Transonic Regime Partially Dispersed Shock Wave Inner Solution is (3. 76) or (3. 81) af, foo -. Corner type equation in the Inner Region (3. 42) _ - --------------------— I Inviscid Flow Regime Eq. (3. 2) Inviscid Solution is Eq. (3. 5) Equilibrium Transonic Regime, Bulk Viscosity Case, Solution is (3. 52) r a ecc I. All Outer Equations are Inviscid Fig. 9

IV. PHASE PLANE STUDY AND NUMERICAL INTEGRATION Since it is hoped that the method of matched asymptotic expansions used for the one-dimensional case could eventually be applied to more practical problems such as flow over an airfoil or inside a passage, the solutions for 1 - D wave structure obtained analytically will be compared with solutions obtained by numerical integration of the full Eq. (3. 1). In general, the full equation for some two or three dimensional problems may be too complicated to be solved even numerically. The phase plane of Eq. (3. 1) T- -U -= a - u V (dx f f dx u v e (3.1) (- ae2) + B vdu d =2 2du du d2 will be studied in detail in order to provide the information needed to start the numerical integration. However, Eq. (3. 1) can be integrated once to yield f d u 2 2 d2u T p, A u + ( - u )___= c [2 vf dx duB T (ye + l)( - u )(-us) (4. 1) 83

84 where boundary conditions dx,- * ~- 2 U~ dx have been used. Next, non-dimensionalizing Eq. (4. 1) by the characteristic velocity u and length A and also considering as before, (Xs/Xc) = 0o0 c s c O(e2), yields 2 d2u 2 2 du Au - +(af - u)dx dx2 (4. 2) 2du 1 =BE dx 2(Ye + l)(u - u )(1 - u) Furthermore, by letting (du/dx)-= p, Eq. (4. 2) becomes dp - (a - u2) + B E2p - (Ye+ 1)(1- u)(u - u) du=......... s —------ (4. 3) du 2 AE pu It is clear from this equation, that p = 0, u = 1 and p = 0, u = u are the two singular points of Eq. (3. 1) in the(p, u) plane where p = 0, u = 1 corresponds to upstream infinity and p = 0, u = uS corresponds to downstream infinity.

85 The singularity at upstream infinity will be examined first. It is, then, convenient to choose this point in the phase plane as the origin by simply letting dr 1u - u=, q dx = (4.4) In order to study the singular behavior, a very small neighborhood around - = 0, q = 0 is considered, and it is reasonable to linearize Eq. (4. 3) to yield d7T dxhi AEth 2 o t 2 2 2 while the higher order terms such as, irq, q, Tr etc. are neglected. The book by Minorsky (1960) is one of many describing the procedure which can be used to investigate the system (4. 5) and is followed below. A transformation = aq + PT (4. 6) r = yq + 6fT reducing Eq. (4. 5) to the canonical form

86 -=D dx 1 (4. 7) dX - dx' 2 is sought, where D1 and D2 are the two characteristic values which determine the nature of the singularity. In order that a non-trivial solution exist for a, j, y, 6 in Eq. (4. 6), the constants D1 and D2 must for the singularity (0, 1) be solutions of the characteristic equation 2 2 (af - 1) e + 1.D +. D + A2 D - 2 (l-u) =0 (4.8) As whose solutions are 1) (af - + 2(e + 1)(1 - us)AE -^ ^afoo - s D1, D2 2A (4 9) 2A2 Since the two roots are always real and opposite in sign, no matter whether u > af or ifco > u the singularity at (0, 1), or upstream infinity is always a saddle point. This is in agreement with Broer and Van den Bergen's results (1954); however, it is just the reverse of the viscous shock plane for which the saddle point occurs at downstream infinity. Therefore the numerical integration has to start from upstream infinity.

87 On the other hand, near the singularity (0, u ) or downstream infinity, the linearized equations are d7 dx q (4. 10) dq (ye + 1) 1 2 2 -- (1 - u7T - (a U )q dx 2 (-u) (afs s 2AE AE where do 7= U -U d q (4.11) The characteristic equation is then 2 1 (2 2'e +1 D + (af -u )D+ 2 (1 -u) =0 (4. 12) AE s s 2A2 s AE 5 2AE Since E is, by assumption small, the two roots are real and have the same sign; the singularity at downstream infinity or (0, u ) is a nodal point. The numerical computations are performed on the IBM System 360 computers. The Runge-Kutta method is used. The integration step size in the x-coordinate is chosen small enough to insure the stability of the integration. For example, the step size used is 0. 02 compared with a value for the thickness of the wave )c which C

88 is roughly 7 for a fully dispersed wave with Mf = 1. The integration stops when the velocity gradient du/dx reaches a value of the order -3 (10 ). The integration converges well, however, the solution could diverge if the step size chosen is not small enough. Because no difficulty was experienced in the present case, no attempt was made to establish the limit of the step size. The integration starts infinitesimally near the saddle point (ir = 0, q = 0) on the directrice or the axis which passing through the saddle. In the (vT, q) plane, the slope of the axis on which the integral curve leaves the singularity is given by (a 2 _2 + r (1-) A2 (a 2 _ 1) dq 0o o (4.13) (4.13) d7 r rF (1 - us) A e s The results are shown on Figs. 14, 15 and 16 for cases Mf = 1. 2, 1.05, and 1. 0 respectively. In the case for a partially dispersed wave, the velocity profile shows a sharp but continuous viscous shock-like front in the upstream portion of the wave. For smaller frozen Mach numbers the curves become less steep. The com-,parisons between the results obtained from the analytical study and the numerical integration of Eq. (4.1) for Mf = 1. 2 and 1. 0 are shown on Fig. 13 and Fig. 11 respectively and are in good agreement.

V. TWO DIMENSIONAL VISCOUS TRANSONIC EQUATIONS Because of the complexity of the general equation for the two dimensional flow of a V-T relaxing gas, approximate equations valid in equilibrium and frozen transonic regimes will be derived. The general equation is, 2 cv U Yi'L; - - 2 a-u ax p j k u + (2.70 ) ax. i ax. a j a~ yVe e 2,-_ 2 a flv au a U ci 1 ax2 i k a a xk P " UaR e... aR. Pr kk v e where i, j, k = 1, 2. 1. Equilibrium Transonic Approximation In this regime, the characteristic speed is a and the characeoo teristic length is a length related to the relaxation length or more precisely, the length of a fully dispersed wave, X /6, where 6 = IMe - 1 c e00 and c = T (v/C v e) a.If perturbations from the free stream are small, substitution of the expansions small, substitution of the expansions 89

90: u(1) u — =1+bu +.. aeoo v = 63/2 (1) V = — =86 / V +... eoo p = = 1 + Fp^+(.) (5. 1) POO a =-=a +a(1)+ af - +'f f a foo f eoo and coordinate stretching x X c (5.2) Y= y 1/2 c into Eq. (2. 66) yields the following equation 32 ( 1) u(1) v(1) 2 a u.(au av- (1 (af - 1) 2 = (y + 1)u ax (5.3) ax Equation (5. 3) is identical to the viscous transonic equation for an inert gas derived and studied extensively by Sichel (1962, 1963 ) except for the fact that the dissipative term (1 + (y - 1)/P " ) V" due to compressive viscosity and heat conduction has been replaced by the term which is, essentially, abulk viscosity

91 due to the internal energy mode. Thus a conclusion is that for velocities very close to the equilibrium sound speed the flow will behave like a viscous transonic inert flow except that the compressive viscosity is replaced by a bulk viscosity. And in this region of flow, the various existing viscous transonic solutions will, thus, be applicable. As it is well known, measurement and observation of the structure of weak curved shock waves is extremely difficult because such curves are very thin at ordinary temperatures and pressures. On the other hand at densities low enough for observation of the structure, it is difficult to apply optical measuring techniques such as interferometry or Schlieren photography. However, as reported by Strehlow and Maxwell (1969), very thick dispersed relaxation waves have been observed. It is, therefore, hoped that Eq. (5. 3) together with experiments in relaxing gases may provide a tool for the study of two dimensional shock structure and possibly of the interaction of shocks, with the boundary layer. (2) Frozen Transonic Approximation As the free stream velocity approaches the frozen sound speed af or (ii - afo/(afo - aeo = 6 < 1, the situation becomes more difficult since now the viscous terms of Eq. (2.66) can no longer be neglected as suggested by the study of one dimensional wave structure. The proper stretched coordinates are then

92 6 X X C (5.4) -y. 61/2 y 61/2 C where 6 is defined, as in the one dimensional case, as 6 = (5. 5) and 6 = /3 if X/X =2 and (af - a)/f = 6 00 00 00 af is now being used as the characteristic speed. Substitution of Eq. (5. 9) and the expansions u =- =1 +5 (1 - a ) u(I) +.. af oo 00 u = v =63/2 (1-a ) (1) +... af e 0000 (5.6) (5.7) a a =- a e a e 00 C + 6 (1-a ) a(1) +... e e 00 (5.8) c Vf X = TF - a c 0oo f V 00 e (5.9)

93 into Eq. (2. 66) and keeping only the lowest order terms yields ^a(l+p - )2 fU(1)a +_ au (5. 10) l+ rr, 2- 2Iu(1).u(1) (1) au(1) a_(I + f a a-u- a- + (X. -0) No known solution of Eq. (5.6) for two-dimensional frozen transonic flow are available. However, it is expected that this equation will play an important role in transonic flow problem when viscous and relaxing effects are both important. Clearly, more work is needed to complete the study of Eq. (5. 10).

VI. CONCLUSIONS AND DISCUSSION 1. The viscous transonic differential equation for a relaxing gas has been derived and used to study the one dimensional shock wave structure problem. Results show that viscous effects cannot be neglected whenever the velocity is close to the frozen sound speed. And the viscous solutions successfully eliminate some of the difficulties encountered in the inviscid theory. 2. Approximate one dimensional equations have been obtained in equilibrium and frozen transonic regimes. The Method of matched asymptotic expansions has been successfully applied to solve these equations. It is expected that the same technique can also be used to solve more complex problems. 3. Approximate two dimensional equations valid in the equilibrium and frozen transonic regimes have also been obtained. No attempt has yet been made to solve these equations and it is clear that more work should be done. However the identical form of the equilibrium transonic relaxing equation to the viscous transonic equation for an inert gas except for a coefficient suggests that this equation might be used to interpret the experimental results obtained by Strehlow and Maxwell (1968). 94

95 4. In a relaxing flow problem with reasonably strong viscous Shock, the idea that the viscous shock drives the much slower relaxation process can be used.

Matching of Inner and Outer Solutions of a Fully Dispersed Wave with af = u0 co af = UO -oo Numerical Solution of Inner Eq. (3. 40 ) u s Outer Solution (3. 34) Common Part x Fig. 10

Comparison between the Solution Obtained from the Method of Matched Asymptotic Expansions and the Result from Numerical Integration of Eq. (4. 1) u = a u o o Numerical Solution Analytical Solution Fig. 11 Velocity Profile of a Fully Dispersed Shock when u o = af co

Uc0 Inner Solution Ub 00 Composite Solution Modified Outer Solution Fig. 12. Inner, Outer and Composite Solutions for a Partially Dispersed Shock Wave with Mf = 1. 2. u S

Comparison between the Solution Obtained from the Idea that the Viscous Shock Drives the Relaxation and the Result from Numerical Integration of Eq. (4. 1) Numerical Solution Analytical Solution — +-+-+-+ — u s Fig. 13 Velocity Profile of a Partially Dispersed Shock Wave with Mf = 1.2 0o

Numerical Solution of One Dimensional Viscous Transonic Relaxing Equation (4. 1) u IP 0 0 TO TO u S 0 x Fig. 14. Velocity Profile of a Partially Dispersed Shock Wave With M = 1. 2 uoo u = 1.2 o00 us = 0.819 af =1.0

Numerical Solution of One Dimensional Viscous Transonic Relaxing Equation (4. 1) U u CO u s 0 x Fig. 15. Velocity Profile of a Fully Dispersed Wave with Very Weak Viscous Shock u = 1.05 00 u = 0.8965 S af =1.0 soO

Numerical Solution of One Dimensional Viscous Transonic Relaxing Equation (4.1) af =u0 so 00 u Co 02 u s x Fig. 16. Velocity Profile of a Fully Dispersed Shock Wave when af = u u = 1.0 00 u = 0.92857 S af. =1.0 0o

REFERENCES Becker, E., Gasdynamics, Academic Press, New York, London, 1968. Blythe, P. A., t Non-linear Wave Propagation in a Relaxing Gas, " J. F.M., V. 30, pt. 1, 1969. Broer, L. J. F., " On the Influence of Acoustic Relaxation on Compressible Flow," Appl. Sci. Res., V. A2, p. 447, 1950. Broer, L. J. F. and Van den Bergen, A. C., " On the Theory of Shock Structure II;'Appl. Sci. Res., Sec. A, V. 4, p. 157, 1954. Chu, B. T., "Wave Propagation in a Reacting Misture, " Heat Transfer and Fluid Mechanics Institute, Stanford, California, 1958. Clarke, J. F., " The Linearised Flow of a Dissociating Gas, " J. F.M.,V. 7, p. 577, 1960. Clarke, J. F., " Relaxation Effects on the Flow over Slender Bodies, " J. F. M.,V. 11, p. 577, 1961. Clarke, J. F. and McChesney, M., The Dynamics of Real Gases, Butterworths, London, 1964. Cole, J. D., Perturbation Methods in Applied Mathematics, Blaisddl, 1968. Cole, J. D. and Messiter, A. F., Guggenhiem Aeronautical Laboratory, California Institute of Technology, Rept. TN 65-1, 1956. Der, K. K., " Linearised Supersonic Non-equilibrium Flow Past an Arbitrary Boundary," Tech. Rept. NASA R-119, 1961. Li, T. Y. and Wang, K. C., " Linearised Dissociating Gas Flow Past Slender Bodies, " Rensselaer Polytechnic Institute, TR AE 6201, 1962. Lighthill, M. J., "Viscosity Effects in Sound Waves of Finite Amplitude, " in Surveys in Mechanics, p. 250, Cambridge, 1956. 103

104 Losev, S. A. and Osipov, A. I., " The study of Non-equilibrium phenomena in Shock Waves, " Soviet Physics, Uspekhi, V. 74, No. 3-4, p. 525, 1962. Minorsky, N., Nonlinear Oscillations, Van Nostrand, Princeton, New Jersey, 1962. Moore, F. K. and Gibson, W. E., "Propagation of Weak Disturbances in a Gas Subject to Relaxation Effects," J. Aero. Sci., V. 27, p. 117, 1960. Moran, J. P. and Shan, S. F., " On the Formation of Weak Plane Shock Waves by Impulsive motion of a Piston," J. F. M. V. 25, pt. 4, 1966. Napolitano, L. G., " Transonic Approximation for Reacting Mixtures, " Israel J. of Tech., V. 4, No. 1, p. 159-171, 1966. Ockendon, H. and Spence, D. A., " Non-linear Wave Propagation in a Relaxing Gas," J. F. M., V. 39, pt. 2, 1969. Polyakova, A. L., Soluyan, S. I. and Khokhlov, R. V., "Propagation of Finite Disturbances in Relaxing Medium, "Soviet Physics - Acoustics, V. 8, No. 1, 1962. Rae, W. J., Cornell University Rept. No. TN60-409, 1960. Rhyming, I. L., " Non-equilibrium Flow inside a Wavy Cylinder," J.F.M., V. 17, pt. 4, p. 551, 1963. Sichel, M., " leading Edge of a Shock-Induced Boundary Layer," Physics of Fluids, V. 5, No. 10, 1962. Sichel, M., " Structure of Weak Non-Hugoniot Shocks, " Physics of Fluids, V. 6, No. 5, 1963. Spence, D. A., " Unsteady Shock propagation in a Relaxing Gas," Proc. of Roy. Soc., Ser. A, V. 264, p. 221, 1961. Strehlow, R. A. and Maxwell, K. R., " Two-Dimensional Diffuse Shock Waves in Carbon Dioxide," AIAA J., V. 6, No. 12, 1968.

105 Tirumalesa, D., " Sonic Line in Non-equilibrium Flows, " UTIAS Tech. Note, No. 94, 1966. Van Dyke, M. D., Perturbation Methods in Fluid Mechanics, Academic Press, 1964. Vincenti, W. G., " Non-equilibrium Flow over a Wavy Wall, " J. F. M., V., p. 481-496, 1959. Vincenti, W. G., " Linearised Flow over a Wedge in a Non-equilibr - ium on Coming Stream," J. de Mecan 1, 193, 19o2. Vincenti, W. G. and Krugger, C. H., Introduction to Physical Gas Dynamics, Wiley, New York, 1965.

UNIVERSITY OF MICHIGAN 3 9015 03524 5128