THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Department of Aeronautical and Astronautical Engineering Aircraft Propulsion Laboratory THE EFFECT OF LONGITUDINAL VISCOSITY ON THE FLOW AT A NOZZLE THROAT Martin Sichel ORA Projects 05143 and 05800 under contract with: THE ARMY RESEARCH ORGANIZATION DURHAM, NORTH CAROLINA GRANT NO. DA-ARO(D)-31-124-G385 administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR June 1965

TABLE OF CONTENTS Page LIST OF FIGURES iii SUMMARY v ACKNOWLEDGMENT vi LIST OF SYMBOLS vii 1. INTRODUCTION 1 2. THE SOLUTION OF TOMOTIKA AND TAMADA 5 3. THE VISCOUS TRANSONIC EQUATION 10 4. SIMILARITY SOLUTION OF THE VISCOUS-TRANSONIC EQUATION 16 5. CONSTRUCTION OF NOZZLE FLOW FIELDS 24 6. DISCUSSION 32 REFERENCES 38 FIGURES 40 APPENDIX: ANALYSIS OF SINGULARITIES IN THE (q - p) PLANE A-1 ii

LIST OF FIGURES Figure l(a) Taylor type nozzle flow. Figure l(b) Meyer type subsonic-supersonic nozzle flow. Figure 2 Transition from Taylor to Meyer type of nozzle flow as postulated by Emmons (1946). Figure 3(a) Z vs. S from solution of Tomotika and Tamada (1950). Figure 3(b) Nozzle flows corresponding to the solution of Tomotika and Tamada (1950). Figure 4 Phase plane behavior of Tomotika and Tamada solution. Figure 5(a) Phase space showing the surface __ and the solutions q = 0, p = - a and q = 09 p = 2a for a = 1. 0. Figure 5(b) Contour map of the surface> _, for a 1. 0. Figure 6(a) Crossing trajectory singularities for Z = 3. 0, a = 1. 0. Figure 6(b) Crossing trajectory singularity for Z = 1. 0, a = 1. 0. Figure 6(c) Crossing trajectory singularity for Z = - 1. 0, a = 1. 0. Figure 6(d) Crossing trajectory singularity for Z = - 3. 0, a = 1. 0. Figure 7 Crossing trajectories for Z = - 1. 0. 2 Figure 8 Numerical solutions of the equation Z"' - 2(ZZ')' + 2aZ' + 4cr = 0 for a = 1. 0. Figure 9 Projection of solution trajectories in the Z'- Z" plane. Figure 10(a) Isotachs and streamlines corresponding to Z = 2S (MeyerVs solution). iii

Figure 10(b) Isotachs and streamlines corresponding to curve A in Figo 8. Figure 10(c) Isotachs and streamlines corresponding to curve B in Fig, 8. Figure 10(d) Isotachs and streamlines corresponding to curve C in Fig. 8. Figure 11 Detailed plot of nozzle contours. f = 0, 25, e = 0. 10, c = 1. 0O iv

SUMMARY An inviscid transonic theory appears to be inadequate to describe the flow near the throat of a converging-diverging nozzle during the transition from the symmetrical Taylor (1930) type of flow to the subsonic-supersonic Meyer (1908) flow. A viscous-transonic equation taking account of heat conduction and longitudinal viscosity has been developed previously (Cole 1949, Sichel 1963, Szaniawski 1963). An exact nozzle type of similarity solution of the viscoustransonic equation, similar to the inviscid solution of Tomotika and Tamada (1950), has been found. This solution does provide a description of the gradual transition from the Taylor to the Meyer flow and shows the initial stages in the development of a shock wave downstream of the nozzle throat. The solution provides a viscous shock like transition from an inviscid supersonic accelerating flow to an inviscid subsonic decelerating flow. v

ACKNOWLEDGMENT The author would like to express his appreciation to the Army Research Organization in Durham, N. C., who supported this work under Grant DA-ARO(D)31-124-G385. The author also would like to thank Mr. Thomas David who carried out many of the computations. The contents of this report have been submitted to the Journal of Fluid Mechanics for publication. vi

LIST OF SYMBOLS A r[ 1 + (y - 1)/Pr]1 a speed of sound a coefficient of thermal expansion, integration constant / h/Rt C specific heat at constant pressure y ratio of specific heats 0 lV*' i/ Eca* E expansion parameter, proportional to max deviation from sonic velocity r Iv*?/Ea*, of the order of weak shock thickness h nozzle half height X characteristic dimension p. shear viscosity p2 bulk viscosity z" u longitudinal viscosity vu ~ kinematic viscosity p pressure; Z' - meaning clear from the text Pr" longitudinal Prandtl number q ZZ' or p' Rt nozzle wall radius of curvature at the throat vii

p density S X + aY, or entropy. Meaning clear from text. u x velocity component u(1) coefficient in expansion for u U alternate notation for u(1) v y velocity component v(1) coefficient in expansion for v V -l/2 V(1) x, y space coordinates X,Y A(x/7r), Ari1/2 E1/2 Z centerline velocity distribution (-) dimensional quantities are indicated by a bar ( )* asterisk refers to value at the critical point where the velocity equals the speed of sound ( ) value on a streamline viii

1. INTRODUCTION As the back pressure decreases the flow through a converging-diverging nozzle changes from one which is symmetrical with respect to the throat to an asymmetrical flow with subsonic flow upstream and supersonic flow downstream of the throat. The two types of flow are illustrated in Fig. l(a) and Fig. l(b). The transition between these two classes of nozzle flow has formed the subject of many investigations, Many important features of such transitional flows are adequately explained by a simple one dimensional or hydraulic theory with normal shocks in the supersonic portion of the nozzle located to satisfy the downstream pressure boundary condition. However, to resolve the details of the flow near the nozzle throat, which is intimately related to the nozzle wall curvature, solutions of the two dimensional or axisymmetric gas dynamic equations must be investigated. The initial phases of the transition during which shock wave formation first starts in the neighborhood of the throat are of particular interest; however, the results of a number of investigators indicate that the inviscid equations alone are unable to provide an adequate explanation of such a transitional flow. It appears that to properly explain the initial stages in the shock formation, or development of shock wave structure near the throat, at least the effect of longitudinal viscosity must be included in the conservation equations. Development of a viscous theory for such transitional flows forms the subject of the present paper, 1

The two dimensional asymmetrical flow at a nozzle throat was first calculated by Meyer (1908) using a truncated series solution of the exact potential equation of an inviscid perfect gas. The calculations are straightforward and the solution appears to give a quite reasonable description of subsonic-supersonic nozzle flow. A similar series approach was applied by Taylor (1930) to the two dimensional symmetrical nozzle throat flow. As the maximum velocity on the nozzle axis increases pockets of supersonic flow begin to develop near the nozzle surfaces as shown in Fig. l(a). Taylor (1930) found that carrying terms up to the fourth degree in the double series in x and y symmetrical solutions no longer exist when the peak velocity on the nozzle axis exceeds some maximum value. For a ratio of nozzle half height h to wall radius of curvature R, (h/R = 1/4) there are no solutions for maximum velocities exceeding 0. 855 a, where a is the speed of sound. Gortler (1939) showed that the series employed by Taylor tends to diverge as the velocity near the throat approaches sonic velocity, and suggested that the difficulty in Taylor's solution may be due to the neglect of higher order terms which are cut off by the truncation process. Gortler (1939) attempted to extend Taylor's solution to the case of transitional flow by relaxing the requirement of symmetry with respect to the nozzle throat; however, a number of artificial assumptions regarding the series coefficients were required making the convergence of his solution suspect 2

Emmons (1946) used the method of relaxation to obtain numerical nozzle flow solutions of the inviscid gas dynamic equations. Emmons (1946) postulated that the transition from the symmetrical Taylor to the asymmetrical Meyer type of nozzle flow starts with the formation of shock waves within the pockets of supersonic flow near the wall as shown in Fig. 2. This postulate was borne out by the calculations. Below a peak centerline Mach number M of 0. 812 the compressible solutions were very much like the flow through a venturi, However, as the maximum centerline Mach number increased beyond this value shock waves had to be placed in the pockets of supersonic flow in order to eliminate residuals in the relaxation calculations, and for sufficiently large M the shock waves within the two supersonic pockets joined at the nozzle center, On the other hand the numerical results contained several inconsistencies The appearance of shock waves is sudden; that is, rather than gradually growing outward from some point in the flow the shock wave when it first appears is of finite length, A second difficulty is that there is a discontinuous rise in velocity immediately behind the shock waves, Emmons (1946) points out that this effect is caused by a discontinuity in the streamline curvature which occurs when a weak normal shock wave is adjacent to a curved wall, A similar effect was observed experimentally by Ackeret, Feldman and Rott (1946) and has also been discussed by Oswatitsch and Zierep (1960) and by Pearcy (1962), Since gradients behind weak shock waves adjacent to curved surfaces must be of the same order as gradients 3

within the shock structure the assumptions which permit the use of Hugoniot jump conditions across the shock are clearly violated, As Emmons (1946) has observed, a perfect fluid theory including shock discontinuities across which the Hugoniot conditions hold appears to be inadequate to describe the nature the transitional flow when shock waves first appear within the nozzle, rather a theory which includes the effects of fluid friction is required. All the solutions described above are in some sense approximate solutions of the inviscid gas dynamic equations. Tomotika and Tamada (1950), on the other hand, found a mathematically exact nozzle type of solution of the transonic equation, which is an approximate equation valid only for regions of inviscid flow with Mach numbers near one. This approach avoids questions of the convergence of either the series or numerical methods. Tomotika and Tamada (1950) obtained an exact similarity solution of the transonic equation describing both the Taylor and Meyer type of flow; however, they concluded that the flow of Meyer's type cannot be approached in a continuous manner from the group of solutions for the flow of Taylor's type, at least on the basis of the inviscid equations. It appears that an adequate explanation of the transition from the Taylor to the Meyer type flow requires consideration of an equation which includes viscous terms such that the formation of shock waves is inherent in the equations themselvesc It has been shown (Cole 194, 1949 Sihel 1963 and Szaniawski 1963) that in regions of transonic flow in which the longitudinal or compressive 4

viscosity is dominant such as in the interior of a weak shock, the flow can be described an equation which is identical to the transonic equation except for an additional viscous term, and which has sometimes been called the viscoustransonic or V-T equation, A nozzle type similarity solution, similar to that of Tomotika and Tamada (1950), has been found for the V-T equation and appears to provide a reasonable picture of the gradual transition from the Taylor to the Meyer type of flow. This viscous-transonic solution forms the subject of this paper; however5 because of the close relation to the work of Tomotika and Tamada (1950) their solution will first be discussed in detail, 2o THE SOLUTION OF TOMOTIKA AND TAMADA Approximate equations for inviscid two dimensional transonic flow have been derived by Guderley (1962) by introducing the series expansions -U -uJ+ ( ) (1) a* + u + V 1/2 ((1) 2 (2) (2) -^=v=e (ev + v +.o ~) (2) and the coordinate stretching transformation X= AX y E1/2 (3 xit ' y=e (3) into the equations for adiabatic inviscid flow of a perfect gas. In the above barred quantities are dimensional, A is a characteristic dimension of the 5

flow, a* the critical speed of sound while E is a small parameter proportional to the deviation of u from the sonic value of unity. To first order, for a perfect gas, the equation v 1) - (+ 1) u(1) u (1) = 0 (4) y x describing flow with only small deviations from the sonic velocity is obtained. If the undisturbed or upstream flow is irrotational, the transonic flow will also be irrotational so that u (1 v (1) (5) y x the nonlinear convective term in Eq. (4) is what distinguishes transonic small perturbation flow from linearized subsonic and super sonic flows, and of course makes the transonic equation very difficult to solve, Eliminating v the transonic equation can be written in the form U - (U2)XX (6) YY xx0 with u(1) replaced by U to simplify notation, and where 1 3/2 X= 2 (y +l )x, Y=[(y+ 1)/2] y (7) Upon introducing the similarity transformation U = Z(S) + 2a 2y2 (8) S = X + (Y2 6

Tomotika and Tamada (1950) found that Eq. (6), collapses to the nonlinear ordinary differential equation (ZZ)'? - z- 22 = 0 (9) Since the flow described by Eq, (8) is symmetrical with respect to the X axis it can be considered to represent a nozzle flow. Tomotika and Tamada (1950) obtained the implicit analytical solution 2 2a3 (Z - 2aS) (Z + aS) 3= (10) for Eq. (9) where the constant of integration a determines the nature of the solution, The arbitrary constant a determines the slope of the two special solutions Z - 2aS, and Z = - aS corresponding to a = 0. In their original paper Tomotika and Tamada (1950) used a = 1, The more general transformation above was introduced in a later paper by Tomotika and Hasimoto (1950)o From Eq, (5), (6) and (9) it follows that the y component of velocity is given by 2 2 2 V =2aY(Z+ 2YaX + ) (11) where i = / 1 )1/2 v(1) 7

To first order in e the streamline slope will be (d7y/ds) = 3/2 v(1) (12) s or (dY/dX) = e V (12a) Integration of Eq. (12) or (12a) then yields the streamline shape, and any particular streamline may be chosen to represent the nozzle wall. Tomotika and Tamada (1950) in particular chose that streamline for which the ratio of the nozzle half height h to the streamline radius of curvature at the throat was 1/4 to coincide with the calculations of Taylor (1930) and Gortler (1939)o The behavior of the function Z(S), which is equal to the velocity U on the nozzle axis, depends upon the constant of integration a. Figure 3(a) reproduces the curves of Z vs. S obtained by Tomotika and Tamada (1950) for various values of a, and for a = 1 0. The solution curves have four branches separated by the special solution curves Z = 2aS and Z = - cS corresponding to a = 0. Branches A and A9 correspond to a < 0 while branches B and B7 correspond to a > 0c From Fig, 3(b) it can be seen that nozzle flows constructed from branch A type solutions correspond to the Taylor type, asymmetrical, nozzle flow. As a- 0 curves of branch A asymptotically approach (1)-P-(4), which has a discontinuous slope at point P the sonic point, and represents the limiting Taylor flow such that the maximum velocity on the 8

nozzle axis is just sonic. The special solution Z = 2aS yields the Meyer type asymmetrical nozzle flow as shown in Fig. 3(b)o It is shown later that the solution Z = 2aS actually is identical with the first few terms of the Meyer solution, Solution curves from branches A' and B have infinite slope at the sonic point, and so are not physically meaningful. Branch B' represents entirely supersonic solutions which are not of interest here. It should be remarked that Tomotika and Tamada's solution has been presented here somewhat differently than in the original paper. Also the direction of the flow as indicated by arrows in Fig. 3(a) and 3(b) seems to have been reversed in the original paper, The behavior of the solution curves follows from the fact that the sonic point Z = 0 is a singularity of Eq, (9) for Z. If Eq. (9) is written in the form ZZ" + (Z' - 2a)(Z' + a) = 0 (9a) it becomes clear that only the singular solutions with Z' = 2a or Z' = - a can pass through the sonic point with Z"p or curvature finite as borne out by Fig. 3(a). Tomotika and Tamada (1950) suggest that the limiting Taylor flow (1)-P-(4) will change discontinuously to the Meyer type (1)-P-(2) provided the nozzle exit conditions change sufficiently. However, their solution does not permit a continuous transition from the limiting Taylor solution to the Meyer solution as becomes particularly apparent when the solution is plotted in the phase plane as in Figo 4. In this plane the sonic line Z = 0 acts as a barrier such 9

that subsonic solutions can never become supersonic and vise versa except for the two singular solutions Z = 2aS and Z - - aS, The question now to be examined is whether taking account of the viscosity in the formulation of the flow equations can resolve this difficulty. 3, THE VISCOUS TRANSONIC EQUATION Within the structure of shock waves the terms of the Navier Stokes equation due to compressive or longitudinal viscosity, and due to heat conduction are of the same order of magnitude as the nonlinear convective terms, for it is the balance between the steepening convective terms and the smoothing dissipative terms which leads to the existence of steady state shock wave structures. One dimensional Navier Stokes solutions of shock wave structure are well known (Hayes, 1958); however, there are regions of flow, which might aptly be called non Hugoniot Shock waves, where the main effect is still a balance between convection and dissipation but where the flow is not necessarily one dimensional In the transonic case, approximate equations describing the flow within such a non-Hugoniot shock layer may be derived from the full Navier Stokes equations (Sichel 1963, Szaniawski 1963) by using Guderley's approach of expansion in a small parameter coupled with coordinate stretching as described above for the inviscid transonic equation. 10

All flow parameters are made dimensionless using sonic point conditions as reference quantities except for the dimensionless pressure p which is defined as P- P(13) p*a* All parameters may then be expanded in the same manner as u in Eq. (1) except for p, and v. Because of the reference quantity used in (13) p must be expanded as p (1) + (2) 14) p-+ p +.14) 3/2 From the theory of characteristics it can be shown that v ~ O(e3/ ) in inviscid transonic flow (Guderley 1962) and this result is assumed to hold in the viscous case as well so that v is expanded according to Eq. (2). Substituting the expansion and coordinate stretching above in the full Navier Stokes equations and equating coefficients of like power of E then yields first, second, and higher order equations for the expansion coefficients. Details of this procedure are described by Sichel (1963) and Szaniawski (1963)o In the present case it is useful to introduce the parameter 0 *- Ep*a*X c a*X where ~p, the longitudinal viscosity, is related to the bulk and shear viscosities a? and gjby ~?' - 4/3 p4+ /, and the asterisk refers to conditions 11

at the sonic point, Lighthill (1956) has shown that the thickness of weak shock waves will be O(v*I7/ea*) so that&- 0(1) implies that the characteristic dimension X and the weak shock thickness are of the same order. It appears logical that the above might be the case within a non-Hugoniot shock layer, and in the derivation below it has been assumed, at least to begin with, that C- 0(1). Equating the coefficients of the lowest power of e in the continuity, momentum, energy, and state equations respectively yeilds the following set of first order equations: (1)+ u (1)= 0 (15a) x x u (l) + P() - 0 (15b) X PX v(1)+ P (1)_ 0 (15c) Y T -1) (a*a*2/C *) p )= 0 (15d) p ( p (1)+ a*T*T (15e) X X X where the subscript denotes partial differentiation and a* is the coefficient of thermal expansion at the sonic point, (a* = 1/T* for a perfect gas), With the help of (15a) and (15b), Eq. (15e) reduces to Eq. (15d) and so is redundant; hence, the Eqs. (15) cannot be solved for the first order coefficientso Use of 12

the expansion and the thermodynamic relation - aT TdS =CdT --- dp (16) P P shows that to first order the entropy remains constant along streamlines, Thus, if the flow originates from a region of constant entropy, as might be (Q ~~ () ()(1) the case in a nozzle, then p p) and T) will simultaneously be zero when u) - 0. The Eq. (15), then, do yield the relations (1) + u(1) 0 (17a) u'' + pM = 0 (17b) T(1) a*a (1) (17c) p which are identical to the relations between the pressure, density, temperature, and velocity perturbations within an acoustic wave. From Eq. (15c) and (17b) it now follows that v 1) = u1) (18) x y so that the first order flow is irrotationalo Equating the coefficients of the next higher power of e yields a set of (1) (1) (1) second order equations containing the first order terms u, p 9) v ) p () and T ) The second order equations are again redundant so that the second order quantities u(2), p (2) T 2, and p 2) can be eliminated. After 13

algebraic reduction, and use of relations (17) it is found from the second order equations that if0c- 0(1) the first order velocities u(1) and v1) must satisfy (1+ Y- r)ul) u - 2Fu1) u v = 0 19) Pr"/ xx x y where Pr" is the Prandtl number based on the longitudinal viscosity while r = (l/a) [ a(pa)/ap ] For a perfect gas r = (1/2)(y + 1). Without the viscous u ( term Eq. (19) would be the same as Eq, (4) xx (1) for inviscid transonic flow. On the other hand if the term v is left out y Eq. (19) becomes the steady Burgers equation (Hayes 1958) which has Taylor's (1910) weak shock structure as a solution. These properties suggested the name viscous transonic or V-T equation. It has been assumed that2 - O(1). On the other hand - O(E) would imply that the characteristic dimension of the region of interest is much larger than the shock thickness. Consideration of Eq. (19) asfy - 0 leads to a singular perturbation problem such as has been considered by Ludford (1952) and Szaniawski (1964a) for example, As jf- 0 the inviscid transonic equation provides a valid description of the flow everywhere except within certain thin viscous or shock layers. With - O0(1/) viscous terms already appear in the first order equations but are no longer balanced by the nonlinear 14

convective terms, The resultant equations, which can also be derived by a straight-forward linearization of the one dimensional Navier Stokes equation (Reissner and Meyeroff 1948) do not have steady one dimensional shock structure solutions which remain finite as x -+ co. The dissipative effects can be much larger than the convective effects only during the unsteady portion of shock structure development; therefore, cQ V 0(1/E) is a physically unrealistic assumption here. Basically the choice ofag!depends upon the region of interest. If the effect of two dimensionality upon the jump condition across nonHugoniot shocks, or the effect of shock structure development on nozzle flow transition is of interest, as in the present case, then-< 0(1) is the appropriate choice. Introduction of the transformation Y = (A/c)Pr/2y Ar12 E1/2 (7/r) X = (A/IJ)x- A(i/7/) (20) U - a) v r -1/2 v) where A= F(1 +2); v*/Ea* reduces the viscous transonic equation and the condition of irrotationality to the universal dimensionless form 15

UX - 2UU +Vy =0 (19a) Uy= VX (18a) or if V is eliminated between (19a) and (18a) UX - (U2)xx + yy 0 (21) Thus using r7, which is of the order of the thickness of a weak shock, as the characteristic length of the flow yields a viscous transonic similitude, for as discussed below, each solution of the dimensionless Eqs. (19a) and (18a) or (21) corresponds to a family of similar solutions in the physical plane, Except for the viscous UXXX term Eq. (21) is identical to Eq. (6), the inviscid transonic equation expressed in terms of U. A nozzle type similarity solution of (21) will now be considered. 4. SIMILARITY SOLUTION OF THE VISCOUS-TRANSONIC EQUATION With the transformation U = Z(S) + 2a2y2 (8) S=X+ aY2 Eq. (21) for U reduces to the ordinary differential equation ZTM - 2(ZZT)' + 2aZv + 42 - 0 (22) 16

Thus Tomotika and Tamada's similarity transformation also works for the viscous transonic equation! The resultant solution again represents a nozzle type flow symmetrical with respect to the X axis. As before the function Z(S) represents the centerline velocity distribution. Equation (22) is identical to Eq. (9) except for the viscous Z'" term. If Eq. (22) is written in the form Zv' - 2ZZ" - 2(Z' - 2c)(Z' + c) = 0 (22a) it is evident that the special inviscid solutions Z = 2aS (23) Z = - aS are also solutions in the viscous case, With the presence of the viscous Z'" term in (22a) the sonic point Z = 0 is no longer a singularity so that solution curves passing through the sonic point are not restricted to the two special solutions (23). So far it has not been possible to obtain any other analytical solutions of Eq. (22); however, it should be possible to obtain solutions Z(S) numerically. Equation (22) is such that for finite Z, Z' and Z" choice of initial conditions Z(S0), Z (S0), and Z"(S0) at some point SO will determine a unique solution (Coddington 1955); however, the question of what initial conditions to choose is certainly not a trivial one. In what might be termed the direct nozzle problem specification of the nozzle contour and conditions upstream and 17

downstream of the throat lead to a boundary value problem for the viscous transonic equation (21), Sichel (1963) has discussed the question of properly set boundary conditions and given a uniqueness proof for the viscous transonic equation valid for subsonic flows while Szaniawski (1964a, 1964b) has investigated the direct viscous-transonic nozzle problem using series expansion methods The present problem, on the other hand, is indirect in that the question asked is whether any of the flow fields corresponding to solutions Z(S) of Eq. (22) satisfy boundary conditions representative of flow through a nozzle throat, while also representing the transition from the Taylor to the Meyer type of flow. In some similarity analyses, such as in the case of the Blasius flat plate boundary layer solution, the boundary conditions which the ordinary differential equation obtained from the partial differential equation must satisfy are precisely specified; however, this is not the case here. All that is known is that the transitional solutions being sought should start where Z(S) < 0 (subsonic flow) and Z'(S) > 0 (velocity increasing), must pass through a maximum which may occur at either a subsonic or supersonic velocity, and then must decrease; however, at this point it is not even known whether such solutions of Eq, (22) exist. Consequently the general properties of Eq, (22) must first be studied to provide a guide for the evaluation of numerical solutions. For this purpose the qualitative behavior of solution trajectories in the phase space will be investigated. 18

Since Eq. (22) is of third order it becomes necessary to consider the behavior of solution trajectories in the three dimensional Z, ZV, Z" phase space, a more difficult problem than the more usual phase plane analysis of second order systems. Equation (22) can be integrated once to yield the second order equation Z" - 2ZZ' + 2aZ + 4a2S = C (24) where C is a constant of integration: however, since (24) contains the independent variable S it is no longer autonomous so that a separate phase plane is needed for each S, Hence, it is more straightforward to deal with the original third order equation, and with the three dimensional phase space, Letting p = Z7, q = p' = Z, the solution trajectories in the (Z, p, q) space satisfy the equations dZ _ dp_ dq pq 2 pp - 2[ )(p + or)+ Zq] ) Thus dZ and dp are zero on the planes p = 0; and q = 0 while dq = 0 on the surfaceX defined by (p - 2a)(p + a) + Zq = 0 (26) If r is the magnitude of the radius vector in any direction in the (Z, p, q) space and rZ is the magnitude of the radius vector in the same direction but to some 19

point on, then it is readily shown that dq > 0 when r > r (27) dq < O when r < ra The surface, also bears an interesting relation to the inviscid solution of Tomotika and Tamada for if their solutions are plotted in the (Z, p, q) space it follows from Eq. (9a) that the inviscid trajectories are constrained to lie on the surface_. There are no singular points where dp = dq = dZ = 0 simultaneously; however, the trajectories p = 2a, q = 0, and p = - o, q = 0 corresponding to the two inviscid solutions Z = 2aS and Z = - aS are singular lines in the sense that dp = dq = 0 on each of them. The phase space with the two singular trajectories and a portion of the surface j is shown in Fig. 5a for the special case a = 1, while contour lines of >j for constant values of Z are shown in Fig. 5b. The points P correspond to the sonic points of the two singular solutions Z = 2aS and Z = - aS. The plane q = 0 is quite similar to the inviscid phase plane of Fig. 4; however, the third order viscous term adds another degree of freedom to the problem and of course the p axis no longer acts as a barrier to solution trajectories. Comparison of the inviscid phase plane and viscous phase space shows how the viscous term drastically alters the nature of the transonic solutions. 20

Certain properties of the solution trajectories can be established immediately. In the region p > 0, Z will be increasing along all trajectories, and vise versa in the half space p < 0O The point at which trajectories cross the plane p = 0 will be a minimum, an inflection point, or a maximum of Z(S) accordingly as q > 0, q = 0, or q < 0. The transitional solutions sought here must thus cross the plane p = 0 where q < 0 and near the sonic plane Z = 0 Additional information can be obtained by studying the trajectories obtained when Z in Eq. (25) is held constant. These curves in the planes Z = const are not solution trajectories but are tangent to the projections of these trajectories at the point where they cross the Z = const plane. Consideration of such "crossing trajectories? for planes corresponding to various values of Z then yields a composite picture of the phase space behavior. In these (q, p) planes the points (0, - a) and (0, 2a) are now singularities where dp = dq = 0. Using well established methods for studying the singularities of second order systems (Minorsky 1962) it is shown in the Appendix that the point (0, 2a) behaves as a saddle point for all values of Z; however, the directions of the two separatrixes of the saddle point and of thek, contour line, which passes through both singularities, do vary with Z. The singularity (0, - a), on the other hand changes in character with Z from an unstable node to an unstable focus to a stable focus and finally to a stable node corresponding respectively to the ranges Z > fo 4 > Z > 0, > Z> - V6, and Z < - V The 21

qualitative behavior of the above singularities is shown in Fig. 6 for various values of Z, and for a = 1, 0, while Fig. 7 shows a more complete picture of crossing trajectories plotted by the method of isoclines for the particular value Z = - 1. The parabolas in Fig. 7 are lines of constant slope. In assessing the significance of the above results it is extremely important to recognize that the crossing trajectories are not solution trajectories and that there are no true signularities in the phase space as for example in the case of one dimensional shock wave (Ludford 1951) or detonation structure (Wood 1961). The phase space behavior is largely determined by the behavior of the solution trajectories near the singular solutions Z = 2uS, Z = - aS. Solutions starting near the Z - 2aS trajectory, no matter how close, will ultimately deviate from this trajectory as S increases. From the crossing trajectory diagrams and the behavior of the crossing trajectories near the singularity at q = 0, p = - 1, it appears that there may be solutions starting infinitesimally near Z = 2aS and with p < 2a and q < 0, which will pass through a maximum in Z on the plane p = 0 and will then, as Z decreases, asymptotically approach the solution Z = - aS. Numerical integration of Eq. (22) indicate that such solutions do indeed exist. Figure 8 shows a number of numerical solutions Z(S) obtained by starting the integration infinitesimally near the point (0, 2a) for different initial values of Z, Starting from an essentially subsonic Taylor type centerline velocity profile this sequence of solutions shows the gradual development of what 22

appears to be a shock wave. All these numerical solutions asymptotically approach the inviscid solution p = - a, q - 0 as S -oo. The initial conditions were adjusted so that the integration constant C1= 0 for then, as can be seen from Eq, (24), the phase of Z(S) will be such that the solutions will be asymptotic to Z - 2aS and Z = - arS As the maximum supersonic value of Z(S) increases the slope ZI(S) in the transition region becomes progressively steeper. If the dimensionless velocity upstream of a weak normal shock is 1 + E U then the downstream velocity will be 1 - E U provided the Hugoniot conditions hold. In Fig. 8 the downstream velocity at first overshoots the Hugoniot value; however, as Z increases the jump conditions more closely approach those of a max normal shock. As Z increases the large values of Z" and Z' in the max transition region make the terms Z" and ZZ' dominant in Eq. (24); however, the equation Z" - 2ZZ' = 0 formed keeping these terms alone is just the one describing the Taylor structure of a weak shock wave, Thus as Z increases, the supersonic-subsonic max transition on the axis of the nozzle seems to approach the structure of a weak normal shock. These results further suggest that with a < 1 solutions will be obtained such that there is essentially a weak normal shock near the nozzle axis which is modified by non-Hugoniot effects only for sufficiently large Y, 23

The expansion scheme used here is based on the assumption that U- 0(1); therefore, when Z is large the solutions are at the limit of what might max be considered consistent. However, in view of the close relation between the nozzle transition and the Taylor shock structure, which remains reasonably accurate for remarkably large upstream Mach numbers (Sichel 1960), these solutions may, nevertheless, be meaningful. Several (q, p) plane projections of trajectories corresponding to the solutions of Fig. 8 are shown in Fig. 9 and support the results of the "crossing trajectory" singular point analysis. Because of the unstable nature of the special solution p = 2a, q = 0, numerical solutions, though started very close to this solution do not, in general, correspond exactly to one of the solutions which asymptotically approaches Z = 2rS, as is evident from the plots of the detailed behavior near the two special inviscid solutions. The situation is similar to that encountered in plane shock structure problems where numerical integration must be started near the downstream saddle point. 5. CONSTRUCTION OF NOZZLE FLOW FIELDS A complete evaluation of the similarity solution described above requires the computation of the corresponding nozzle flow fields. For this purpose isotachs and streamlines must be determined and it is also necessary to relate the dimensionless solution in the X-Y plane to the physical x-y plane. 24

Since the dimensionless speed q = q/a* is given by 2 32 1/2 2 q = [(1 + U)2 + 3Vr]12 = 1+ eU 0(2) (28) it follows that isotachs correspond to contours of constant U to the present order of approximation. The streamline slope (dy/dx) is given by (dy/dx) = (i/u) = 3/2/2v/(1 + eU + ) = 3/2r1/2 + (E5/2) (29) or in terms of the stretched coordinates X and Y (dY/dX) = E rv (29a) The velocity V must now be determined from the similarity solution. From the condition of irrotationality (18a) and the similarity transformation (8) it follows that V = 2aYZ + 4a2XY + f(Y) (30) where f(Y) is a function of Y. From the V-T equation (19a), Eq. (30) above, and Eq. (24) for Z it follows that 4 3 3 (31) ^ ^laV -C^ C^ (31) f(Y) = aY - C (31) Since V(X 0) - 0 for nozzle flow, C2 = 0, The constant C1 depends upon the initial conditions used in evaluating Z and upon the origin SO chosen for S, and from Eq. (24) is readily seen to have the value 25

C1 = Z i(S) - 2Z(S0)Z'(So) + 2aZ(S0) + 4'S0 (32) C1 does not affect the functional form of Z(S) but merely determines the phase of the solution with respect to the S coordinate, It now follows that 2 2 32 C V = 2Y (Z + 2a2 X + 2 y - ) (33) 3 2 which is identical to the result of Tomotika and Tamada (Eq. 12) except for the constant C1 which they set equal to zero but which has been retained here. Streamlines can now be determined by integrating the differential equation (29a) for different initial conditions using V as given by Eqo (33). Any streamline can be considered as the wall of a nozzle; however, it is of particular interest to choose a streamline with a predetermined ratio of nozzle half height to radius of curvature at the throat in order to compare the viscous-transonic results to the inviscid calculations of Taylor (1930) and Tomotika and Tamada (1950). At the nozzle wall throat V = 0 so that aZ(Xt + cYt 2) + 2 X + Y - (C1/2) = 0 (34) where the subscript t refers to the throat coordinates. Letting h be the nozzle half height at the throat and Rt the throat radius of curvature of the nozzle wall, it follows from Eq. (29a) and (33) that (h/Rt) = h(d/d2)t = 2eYt2 [aZ' (Xt + aYt) + 22] (35) 26

Equations (34) and (35) are sufficient to determine the throat coordinates (Xt Yt) once (h/Rt) and e are specified, and with (Xt, Yt) known the wall streamline is obtained by integrating Eq. (29a). It now is necessary to establish the connection between the nozzle solution above, which is expressed in terms of the dimensionless coordinates X and Y, and the physical plane. The viscous transonic equation when expressed in the dimensionless form (21) provides the basis of a viscous-transonic similitude. For each solution of Eq. (21) there exists a family of physical flows corresponding to different values of the parameters e, r; v,*, a* and Pr'. The members of this family are similar, and each point (X, Y) of the dimensionless solution defines a set of corresponding points in the family of similar solutions. This similitude is closely related to the more conventional transonic similitude, discussed, for example, by Guderley (1962); the main difference between the two being the nature of the characteristic length X. In the conventional transonic similitude X is some characteristic dimension of the flow such as the chord length of an airfoil, for example, but in the viscous-transonic case X = (v*?/Ea*) as is evident from Eq. (20), i. e., X is of the order of the thickness of a weak shock wave. With decreasing E, X can remain fixed in the inviscid case but must increase in the viscous transonic case. 27

The definition of the parameter E, which characterizes the maximum deviation of the fluid velocity from the sonic value a*, is arbitrary, but usually related to the particular problem under investigation. In the study of nonHugoniot Shock Structure (Sichel 1963) it was convenient to let e = (ii1/a*) - 1 where u, is the velocity of the undisturbed flow upstream of the shock wave, while in flow about bodies the choice e = (M - 1), where M is the Mach number of the undisturbed flow, is frequently made. In the present case, since neither of the above definitions is suitable, E represents the value of (u/a*) - 1 corresponding to points where U = 1o As a consequence of this definition only that domain of the solutions of Eq. (22) for Z such that U - 0(1) is consistent with the expansion scheme used here. From Eq. (20) it follows that for given X and Y corresponding values of x and y are given by x= (r7/A)X; y = (7/Ari /2E12 (20a) -1 - 3/2 so that for fixed fluid properties and a*, x - e and y e/ at corresponding points, The streamlines in the x ~y plane corresponding to the streamlines passing through a particular point (X1, Y1) in the X, Y plane will be called corresponding streamlines. This reference point (1'yl) transforms accord3/2 ing to Eq. (20a) but since v O(3/2 ) the streamline slope, dy/dx, must also 3/2 be O(e )o On the other hand if all points on corresponding streamlines 28

tranformed according to (20a) the result would be (dy/dx) - ( 1/2 ) This strange behavior, which was also noted by Guderley (1962), is responsible 2 for the appearance of e in Eq, (30a) for streamlines in the X-Y planeo The relation between the dimensionless and physical nozzle solutions is now established, In the x direction the length of the region of interest will be O(r), Since (v*"/a*) is of the order of a mean free path, r = v*"/a*e will be very small unless the density is low or e is very small. The nozzle half height h is the other significant dimension of the flow, From Eq. (35) it follows that streamlines with (h/Rt) fixed will not at the same time be corresponding streamlines. Thus assuming that [ aZI(X + aYt ) + 2a ] 0(1) it follows from Eq. (35) that t -O[(h/Rt)l/2 (v*?/a*E2 (36) 1/2? and if (h/Rt) / 0(1) it follows that h/77 (1/e) (37) so that the nozzle half height will be much greater than the thickness of a weak shock provided that e << 1. For the particularly simple inviscid solution, Z = 2o(S - S0) explicit expressions for the isotachs, the vertical velocity, and the streamlines can be found. Equation (8) for U becomes U = 2c(X - S0) + 4c2Y2 (38) 29

so that isotachs will be parabolas in the X-Y plane. It can be seen that SO is equal to the value of X for which U(X, 0) = 0 i. e., it is the location of the sonic point on the axis of the nozzle. From Eq. (32) C1 = 4a SO and from Eq. (33) it then follows that V = 8Y[ a2 - - S) + (2/3) o3y2] (39) Letting f = h/Rt Eq. (35) yields the result 1/2 Yt 23//2 (40) from which it is again clear that as E varies the streamline with fixed 3 will not be a "corresponding streamline" in the sense considered above. From Eq. (21) relating X and Y to x and y it follows that h fi1/2 71 23/2r/2 (41) With Eq. (41) the dimensionless velocity perturbations v/a* and (u/a* - 1) can be expressed in terms of the physically more meaningful coordinates x/h and y/h so that - 1x y 2 1 - eU = - a*-e 2r (42) a* _V= h [h )h + 6 \h/ 30

where x and S are related by o o S -A0o r In the case of a perfect gas with F = (y + 1)/2 Eq. (42) is identical to the nozzle velocity distribution obtained from the first three terms of Meyer's double expansion for the velocity potential (Meyer 1908, Hall and Sutton 1964), except that in accordance with the discussion above the choice of the characteristic nozzle dimension is no longer arbitrary. From Eq. (34) it follows that (it/h) - (i/h) - 1/2 I /6 (43) Integration of Eq. (29a) for the streamlines is straightforward and will not be reproduced here. Isotachs and streamlines corresponding to the Meyer type flow with a = 1, are shown in Fig. 10(a) plotted in the (i/h, y/h) plane. The wall streamline has been chosen so that h/R = 0. 25 as in the calculations of Taylor and of Tomotika and Tamada. In the more general case streamlines can only be obtained by numerical integration of Eq. (29a), and it is no longer possible to obtain simple expressions for U and V in terms of x/h and y/h, Figures 10(b), (c), and (d) show the streamlines and isotachs corresponding to solutions Z(s) shown in Fig. 8. Figure 10(b) represents a Taylor type nozzle flow with subsonic velocities throughout but regions of high velocity flow near the nozzle wall. In Fig. 10(c) the maximum centerline velocity is just sonic, while there are pockets 31

of supersonic flow near the nozzle wall, Upstream of the throat the flow is similar to the Meyer type flow of Fig. 10(a), but downstream a certain crowding of the isotachs as compared to the Taylor flow of Fig. 10(b) is evident. In Fig. 10(d) the velocity along the axis becomes supersonic beyond the throat but this supersonic region is followed by a rapid deceleration to subsonic flow, Figure 10(d) appears to indicate an initial stage in the development of a shock wave downstream of the throat, There are undulations in the portion of the nozzle wall downstream of the throat which follow from the rapid changes in ZF(S). This boundary condition goes with the similarity solution for, as mentioned previously, there is no freedom to choose the streamline shape in the present case. Nevertheless, Figo 10(d) depicts flow through a nozzle with a throat or section of minimum area. The wall streamlines of Fig, 10 are shown in greater detail in Figc 1lo The nozzle contour of Fig. 10(d) coincides with the Meyer flow contour upstream of the throat while far downstream this contour (Curve C) approaches the shallower Taylor type contour. This result is not surprising since the flow in Figo 10(d) does change from a Meyer to a Taylor type of flow downstream of the throato 6. DISCUSSION By including the effect of longitudinal viscosity in the equation for plane transonic flow it has been possible to obtain solutions which provide a smooth 32

transition from the Taylor to the Meyer type of nozzle flow, and which show what happens in the initial stages of shock formation downstream of the nozzle throat. The difficulties near the sonic point, so characteristic of the inviscid analysis, disappear when the nozzle problem is formulated in terms of the viscous-transonic equation. The solution of the viscous-transonic equation found above is exact but is purchased at the penalty of not being able to specify an arbitrary nozzle wall shape. Hopefully this exact solution will provide a guide to future studies of viscous transonic flow. The results of this investigation indicate that the viscous-transonic equation may provide a key to the solution of problems where the use of shock waves with one dimensional Hugoniot jump conditions leads to contradictory results. The solution obtained here is in some sense related to Taylor's weak shock solution (Taylor 1910) describing the viscous transition between uniform, supersonic and subsonic flows. The viscous-transonic solution yields a viscous transition between a supersonic flow of increasing and a subsonic flow of decreasing velocity. In each case viscous effects vanish upstream and downstream of the transition. As the maximum value of the centerline velocity U(X, 0) increases the viscous transonic transition appears to approach the Taylor weak shock structure. The horizontal scale of the portion of the nozzle flow under consideration here is of the order of r7, the thickness of a weak shock wave, while the 33

vertical scale h is of the order of 1r/co Unless << 1, and the flow is one of low density9 the throat width of the nozzles under consideration here will be extremely small. For example for nitrogen at 27. 40C, with p = 1 atm, and = 0. 1, r O(, 001 mm) while h O(, 01 mm)o On the other hand with p = 001 atm9 and e = 0. 019 r - 0(1. 0 mm) while h 0(100 mm) which is certainly of a more reasonable magnitude. The applicability of the viscous-transonic nozzle solution to practical flows thus appears limited; however, it is possible to attach a broader interpretation to the results obtained here. The nozzle height h or characteristic y dimension is essentially the vertical distance over which non-Hugoniot effects will be important, and the nozzles considered above are such that non-Hugoniot effects are important over the entire nozzle widtho Under ambient conditions for nozzles with half heights several orders of magnitude greater than the shock thickness 77 the regions of non-Hugoniot flow are probably confined to relatively small regions near the nozzle walls as is to some extent borne out by the calculations of Emmons (1946)o It is possible to consider half of a viscous-transonic nozzle as representing such a region of non-Hugoniot flow near the wall, and perhaps not too unreasonable to suppose the flow beyond this wall region continued by means of more conventional Rankine-Hugoniot shock waves. It has already been noted that the viscous transition at the centerline, io eo Z(S) closely resembles the Taylor weak shock structure as the velocity U(X, 0) upstream of the transition increases. 34

As mentioned above it may also be possible to construct nozzle flows with normal shocks near the center by constructing solutions with r ~< 1. Since actually h - r7/oe such a procedure might also lead to more reasonable throat heights. Of course what is really needed is a solution approaching the one dimensional shock structure as the distance from the nozzle wall tends to infinity. It will be difficult to obtain a precise experimental verification of the results obtained here. The presence of the nozzle wall boundary layer will make it very hard to reproduce bounding streamlines or nozzle contours which agree exactly with the streamlines obtained from the similarity solution, and of course, slight shifts of the boundary can cause relatively large changes in transonic flow. The region of interest will be extremely small unless the density is low, and only slight deviations from the sonic velocity are considered. Under such conditions it is difficult to make accurate velocity and density measurements. Clearly the present investigation touches on the problem of whether it is possible to have regions of supersonic flow imbedded in a subsonic flow without the existence of shock waves. Extensive investigations of this problem based on the inviscid transonic equation have been made and are, for example, discussed by Manwell (1958, 1963) who, with others, concludes that it is not in general possible to obtain smooth inviscid solutions for the transonic flow in such regions. A detailed discussion of this problem in the light of the 35

viscous-transonic equation is beyond the scope of the present paper; however, the existence of supersonic regions within regions of subsonic flow does not appear to result in any difficulties when viscous effects are taken into accounts This is not a surprising result for the viscous transonic equation inherently contains the possibility for formation of steady shock structures where required by the conditions of the flow while the inviscid equations do not. Whether, in general, the proper inclusion of viscous effects can eliminate the difficulties encountered by Manwell and others in constructing transonic solutions is certainly a worthwhile subject for future investigation. There are two basic differences between Szaniawskivs (1964a, 1964b) studies of viscous-transonic nozzle flow and the present work. While Szaniawski permits an arbitrary nozzle contour his solutions are approximate rather than exact as in the present case. Also, the expansion scheme used by Szaniawski is different. The nozzle half height is used as the characteristic flow dimension with the result that while (y/L) - 0(1) the dimensionless coordinate corresponding to Y in the present paper has a maximum value of O( 1/2) and v O(e )E As a consequence the velocity U is a function of X only in the first order series solution, and to obtain details of the flow field, in particular a non-trivial result for the shape of the sonic lines, it is necessary to compute (2) the second order coefficient u In the expansion scheme used here Y 0(1) and all details of the flow are recovered from the first order solution. 36

Szaniawski finds, as in the present paper, that in the Taylor flow viscous effects become crucial near the nozzle throat as the maximum velocity approaches the sonic value. 37

REFERENCES Ackeret, J., Feldman, J., and Rott, N., 1946, NACA TM No. 1113. Coddington, E. A. and Levinson, N., 1952, Theory of Ordinary Differential Equations, McGraw Hill, New York. Cole, J., 1949, Problems in Transonic Flow, Ph. D. Thesis, Calif. Inst. of Tech. Emmons, H. W., 1946, NACA TN No. 1003. Gortler, H., 1939, Z. angew. Math. Mech., 19, 325. Guderley, K. G., 1962, The Theory of Transonic Flow, Addison Wesley, Reading, Mass. Hall, I. M. and Sutton, E. P., 1962, Symposium Transsonicum, (Editor, K. Oswatitsch), Springer, Berlin. Hayes, W. D., 1958, Fundamentals of Gas Dynamics, Sec. D, Princeton University Press, Princeton, N. J. Lighthill, M. J., 1956, Surveys in Mechanics, (Editors, G. K. Batchelor and R. M. Davies), Cambridge University Press, New York. Ludford, G. S. S., 1951, J. Aeronaut. Sci., 18, 830. Ludford, G. S. S., 1952, Quart. Appl. Math., 10. Manwell, A. R., 1958, Proc. Roy. Soc., A, 245, 481. Manwell, A.R., 1963, Arch. Ratl. Mech. Anal., 12, 249. Meyer, T., 1908, Uber Zweidimensionale Bewegungsvorgainge in Einen Gas, Das Mit Uberschallgeschwindigkeit Stromt, Ph. D. Dissertation, Gottingen (see also Carrier, G., 1951, Foundations of High Speed Aerodynamics, Dover, New York). Minorsky, N., 1962, Nonlinear Oscillations, Van Nostrand, Princeton, N. J. Oswatitsch, K. and Zierep, J., 1960, Z. angew. Math. Mech., 40T, 143. 38

Pearcy, H. H., 1962, Symposium Transonicum, (Editor, K. Oswatitsch), Springer, Berlin. Reissner, H. J. and Meyeroff, L., 1948, PIBAL Rept. No. 138, Polytechnic Inst. of Brooklyn. Sichel, M., 1960, JO Aerospace Sci., 27, 635. Sichel, Mo, 1963, Phys. Fluids, 6, 653. Szaniawski, A., 1963, Archiwum Mechaniki Stosowanej, 15, 904. Szaniawski, A., 1964a, Importance Des Effets De Dissipation En Ecoulement Transsonique, Paper No. 64-587, Internatl. Council of Aeronaut. Sci. Congress, Paris, Aug. 24-28. Szaniawski, A., 1964b, Archiwum Mechaniki Stosowanej, 16, 643. Taylor, G. I., 1910, Proc. Roy. Soc., A, 84, 371. Taylor, G.I., 1930, A. R. C. R. and M. No. 1381. Tomotika, S. and Tamada, K., 1950, Quart. Appl. Math., 7, 381. Tomotika, S. and Hasimoto, Z., 1950, J. Math. Phys., 29, 105. Wood, W. W., 1961, Phys. Fluids, 4, 46. 39

M>1 M<I\ M<1 - M<1 M>1 Figure l(a). Taylor type Figure l(b). Meyer type subnozzle flow. sonic-supersonic nozzle flow. >//// M>1 Figure 2. Transition from Taylor to Meyer type of nozzle flow as postulated by Emmons (1946).

z Ii? a>o 1, Z2S B' /// 0.2 /0 ua<O -B 0 -Y q B h A A' /1/ IFlow -o.F2 /(/ s 1so n of A aa< 0 -0.2- 0 0. 2 Figure 3(a). Z vs. S from solution Tomotika and Tamada (1950). q>a q>a A 9 <a,q >a Branch A (D- P-( (D-P-) Taylor Type Flow Limiting Taylor Meyer Type Flow Flow Figure 3(b). Nozzle flows corresponding to the solution of Tomotika and Tamada (1950).

dZ dS A' / B, A \ - I D2_/ \\ Al Figure 4. Phase plane behavior of Tomotika and Tamada solution.

z a- =1.0 I q-o P q Op, p - +2 Figure 5(a). Phase space showing the surface Z and the solutions q = O, p = -a and q = 0O p = 2a for a = 1,0,

,= -O. 5Z 0.0,I q -1Figure 5(b). Contour map of the surface, for a = 1.0 0.5 Figure 5(b), Contour map of the surface Z, for a = 1,0.

Z 3.0 ____\ _ I 1 2 -— 1- /Figure 6(a). Crossing trajectory singularities for Z = 3.0, a = 1.0. N ' ORZ — Lo q 1 2 3 / / 1 - Figure 6(b). Crossing trajectory singularity for Z = 1.0, a = 1.0.

or- 1.0 _ _ 1_ _ I I, q ^) a1C1 -2 -1 1 2 Figure 6(c). Crossing trajectory singularity for Z = -1.0, a = 1.0. fe>>^- 0 —1.0

-3. 0 -2:: - f — 2 -O F r. C g O fo2.0 3.0 Figure 7. Crossing trajectories for Z = -1.0.

8.0 Z(s) 6.0 -v 4.0 / 2,0 -Curve C~ x\ S Ah I t -2. 0 -LO 1 0\ 3.0O 4.0 5.0 6.0 Curve B -4,0 Curve A - o- =1.0 _Solutions Are Asymptotic To Z = 2S and Z = -S Figure 8. Numerical solutions of the equation Z'" - 2(ZZ')' + 2aZ' + 4a2 = 0 for a = 1.0.

Numbers near circled points Subsonic - Supersonic indicate corresponding values Subsonic of Z.18 -0.99 0.67 _ T^ " > -1.50 -4. 0 Z-1.04 -- -3.0 27Z 1. 28 -20 1 9 —8 -Z 2. 9 I I0 _! ~. I I I I -'- ~Z" -16.0 -12. 0 1. -1. 4.0 0 12.0 16.0 Detailed 1.2 2.2 -46 Q IBehavior Near 0.87 Detail ed Behavior \ -1.8 Z", Z'-2.0 3.0 -1.1 Near Z"0, Z' —1 0.16 _. -o: 8.8 4.0 — 1. 57 1 6 0.2 I Z" -5.0 -0. -0.2 -0.1 0.1 0.2 0 -4 1 -2.1 -3. 5 -0. 8 2.2 ZI ZI Zp Figure 9. Projection of solution trajectories in the Z'-Z" plane.

h ~ ~ ~ ~ ~~~~~~~~O.~~~~~25 E~O~j, )fml*4u/a l+E U Figure loa) Is~t~~ach.strai to Z 2S (Mleyer Is solution temie )re-odn *N I \^v^^y// Throat L.00 01 0 Li~ 0 () I Ot'c s an teajf0 to \e\0h A8r '' ~~^\ \j\ \\ ^W AC I *\' U\ \o! /\ /I L ~ I *^ / \I e-ai. y.1 4 *</^ 1.0 h E~ '0.1 u 1 'i/a*~i+EU Figure l o(bj 7,,,stcs adsra lie orsod to ~- = 2S (Mye' souin. ~1oo,^

- \ x \ \-.,Th o at:~ ~ \ o\ o\ \\ I I I ' / I: \ \ x \ \ o \ / [ I 1 I ijY~I I \ \ I | | I E=0.1, y=1.4, u/a*=1+cU C=h/Rt =0.25 Figure 10(c). Isotachs and streamlines corresponding to curve B in Figure 8. Th roat ~~1.0- ~ ~ ~ ~~ - 1. 0 =0.1, y=l.4, ula=1+EU, =.hlR, =0.25 Figure 10 Isotachs and streamlines corresponding to curve C in Figure 8. to curve C in Figure 8.

1.40 1.30 Meyer Y/ Curve A SoL tion hJi~-S o /x 1.10 Curve C Curve C and Meyer —. — Curve B le 0 Solution -1.0 0. 01.0 2o \h tJ Figure 11. Detailed plot of nozzle contours. P = 0.25, e = 0.10, a = 1.0.

APPENDIX ANALYSIS OF SINGULARITIES IN THE (q - p) PLANE The equations of the crossing trajectories in the q-p plane are (dp/dS) = q (A-l) (dq/dS) = 2[(p - 2a)(p + c) + Zq] (A-2) where Z in this case is treated as a parameter. The singularity at (03 2a) will be considered first. Choosing this point as origin and letting 4T = (p - 2a) Eq. (A-l) remains unchanged while (A-2) becomes (dq/dS) = 2[ (7 + 3a) + Zq] (A-3) To study the singularity, which is now at (0, 0) it is sufficient to consider the linearized form of Eq. (A-3) valid only in the immediate neighborhood of (0, 0). The problem then is reduced to an investigation of the linear system (da/dS) = q (A-4) (dq/dS) = 6 c + 2Zq The book by Minorski (1962) is one of many describing the procedure which is followed below. A transformation aq + + 7 (A-5) A r = yq + 65 A-1

reducing Eqo (A-4) to the canonical form (dd/dS) = D1 (A-6) (drl/dS) = D2r7 is sought. The nature of the constants D1i and D29 which are also called the characteristic values, then determines the properties of the singularity. In order that a non-trivial solution exist for a, i, y, and 6 in A-5 the constants D1 and D2 must for the singularity (0, 2r) be solutions of the characteristic equation D2- 2ZD - 6r= 0 (A-7) so that D, D= Z + (Z 2 + )1/2 (A-8) From Eq. (A-8) D1 > 0, D2 < 0 for all Z so that the singularity at (0, 2a) is a saddle point for all Z withdirectrices ~ = 0, along which solution trajectories approach the singularity and 77 = 0, along which trajectories leave the singularity, In the (q, ' ) plane the slopes of the directrices are given by (dTr/dq) 0 [(Z2 + 6r)1/2 + Z]/6c (A-9) (d/q -~(2 1/2 (d/dq)r = 0[(Z + 6a)1 - Z]/6a A=2

Near the singularity at (0, -a) with T = p + a the linearized equations are d7/dS-= q (A- 10) (dq/dS) = - 6o~ + 2Zq The characteristic equation is then D - 2ZD + 6a = 0 (A-1) with the solutions D1D2 Z + (Z2- 6(o)12/ (A-12) 1 22 For Z > V/, D1 > D2 > 0 so that the singularity is an unstable node with solution trajectories tangent to the line - = 0o For Z < - /6J, D2 < D1 < 0 and the singularity is a stable node with trajectories tangent to the line ri = 0, In the (q, "r) plane the slopes of the lines 5 = 0 and r7 = 0 are given by (diT/dq) =0 = [Z + (Z2 - 6)1/2]/6a (A-13) (d7?/dq) =[Z-(Z - 6)1/2]/6a In the range V- > Z > - V, D1 and D2 are complex conjugates. Since Re D1 > 0 and Re D2 > 0 for 4VJ > Z > 0 the singularity is an unstable focus 1 2 for this range of Z. In the range 0 > Z > - V, Re D1 < 0, and Re D2 < 0 so that the singularity is a stable focus, The results of the analysis above are borne out by the numerical calculations, A-3