THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Department of Aerospace Engineering Gas Dynamics Laboratories AN AXISYMMETRIC SIMILARITY SOLUTION FOR VISCOUS TRANSONIC NOZZLE FLOW M. Sichel and Y. K. Yin ORA Project 07146 under contract with: THE ARMY RESEARCH ORGANIZATION DURHAM, NORTH CAROLINA CONTRACT NO. DA-31-124-ARO-D-276 administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July 1966

TABLE OF CONTENTS Page SUMMARY iii ACKNOWLEDGMENT iv LIST OF FIGURES v LIST OF SYMBOLS vi 1. INTRODUCTION 1 2. THE AXISYMMETRIC NOZZLE SOLUTION 1 3, ASYMPTOTIC BEHAVIOR 7 4. NATURE OF THE FLOW 12 5. DISCUSSION 17 REFERENCES 18 ii

SUMMARY An axisymmetric viscous-transonic equation is presented. A nozzle type similarity solution of this equation has been found, which describes the initial stages in the development of shock-waves downstream of a converging-diverging nozzle throat. This solution is an extension of a two dimensional solution found previously (Sichel 1966). By an appropriate choice of an arbitrary scaling constant solutions were found such that there is essentially a weak normal shock near the axis with effects of wall and shock wave curvature occurring only at a sufficiently large radius. The upstream and downstream asymptotic behavior of these viscous-transonic nozzle solutions has been investigated. iii

ACKNOWLEDGMENT The authors would like to express their appreciation to the Army Research Organization in Durham, N. C. for their support of this work under Contract DA-31-124-ARO-D-276. The contents of this report have been submitted to the Journal of Fluid Mechanics for publication. iv

LIST OF FIGURES Figure l(a), (b), and (c). Numerical Solutions of the Equation Z'" - 2ZZ" - 2(Z' - w1c) (Z' + wo2) = 0 Figure 2(a). Isotachs and Streamlines Corresponding to Curve A in Fig. l(a) Figure 2(b). Isotachs and Streamlines Corresponding to Curve B in Fig. l(c) v

LIST OF SYMBOLS A (1/2)(y+ 1) [ 1+(y6 - 1)/Pr"]1 a speed of sound p rt/Lt y ratio of specific heats E expansion parameter, proportional to maximum deviation from sonic velocity 71 r/rt 71 Tv*"/ Ea*, of the order of weak shock thickness;f' longitudinal viscosity v kinematic viscosity Pr" longitudinal Prandtl Number r radial coordinate (dimensional) R /(1 y/2) ( + 1 ) E/ A/ a arbitrary scaling factor S X + aR u x velocity component U E [ (u/a*) - 1], also a standard solution of Weber's equation V E 3/2 [(1/2)(y + 1)] 1/2 v/a*); also a standard solution of Weber s equation v r velocity component x axial coordinate X A(xA/) Z centerline velocity vi

LIST OF SYMBOLS (cont.) C1'%2 Z - 1aS, Z + caS 0 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 vii

1. INTRODUCTION A similarity solution of the inviscid transonic equation describing flow near the throat of a converging-diverging nozzle was found by Tomotika and Tamada (1950) in the two dimensional case and by Tomotika and Hasimoto (1950) in the axisymmetric case. These solutions describe both the symmetrical Taylor flow with subsonic velocities upstream and downstream of the throat and the asymmetrical subsonic-supersonic Meyer flow, but do not permit a smooth transition between the two types of flow. Since this transition is accompanied by the formation of shock waves downstream of the nozzle throat this difficulty appears due to the neglect of viscous effects. If the longitudinal or compressive viscosity and the thermal conductivity are taken into account the inviscid transonic equation should be replaced by a "viscous-transonic" equation (Cole 1949, Sichel 1963, Szaniawski 1963). Sichel (1966) found nozzle type similarity solutions of the two dimensional viscous-transonic equation, that do permit the smooth transition from the Taylor to the Meyer type of flow and display the initial stages in shock wave formation downstream of the nozzle throat. An axisymmetric viscous transonic nozzle solution has also been found and is the main subject of the present paper. 2. THE AXISYMMETRIC NOZZLE SOLUTION The viscous-transonic equation for the axisymmetric flow of a perfect 1

gas can be shown to be UXXX XX U2RR riving (1) UR= - (2) In these equation X and R are dimensionless axial and radial coordinates and U and V are the corresponding dimensionless velocities. Equations (1) and (2) are derived from the full Navier-Stokes equations by a simultaneous coordinate stretching and series expansion, and since the derivation is almost identical to that of the two dimensional viscous-transonic equation (Sichel 1963) the details will not be reproduced here. Without the third order viscous term equation (1) is the same as the inviscid axisymmetric transonic equation (Guderley 1962); without the term UR/R and with R replaced by Y equation (1) becomes the two dimensional viscous-transonic equation. The stretched dimensionless coordinates X and R, and the corresponding velocities U and V are related to dimensional coordinates x, r, and velocities u, v by X A (x/i); R= - 72)(Y + 1) 2 A (r/) (3) a*U; /a*) 3/2

where A= (1/2)(y + 1) [1 + (y - 1)/Pr"] In equation (3) a* is the critical speed of sound while E is a small parameter proportional to the maximum deviation of (u/ a*) from the sonic value. The characteristic length r7 is equal to (,* "/ E p* a*), which is of the order of the thickness of a weak shock. ii" is the compressive or longitudinal voscosity (Hayes 1958) while p is the density, and the asterisk refers to conditions at the sonic point. The Prandtl number Pr" is based on the viscosity p1" and is assumed constant. The relations between the deviations of the pressure, density, and temperature from their critical values and the velocity perturbation U are identical to those within an acoustic wave (Sichel 1963). The transformation U= Z(S) + 2a2R2 (4) S = X + cR2 which was also used by Tomotika and Hasimoto (1950) reduces the axisymmetric viscous transonic equation to the ordinary differential equation Z"' - 2ZZ'" - 2(Z' - wla)(Z' + w2r) = 0 (5) where 3

w1 =V5 + 1 2, 2 = /5 - 1 The flow described by equation (4) can be considered to be a nozzle flow by choosing one of the streamtubes as the nozzle wall. Z(S) will be the value of U on the nozzle axis R = 0. The arbitrary constant a is related to the streamline curvature and will be discussed further below. Except for the value of the constants wl and w2 equation (5) is identical to the ordinary differential equation considered in the two dimensional case; therefore, the properties of equation (5) are similar to those of the two dimensional equation, which has been discussed in detail by Sichel (1966). As before the inviscid solutions Z = wl (S - b) (6a) Z = - w2 (S - b) (6b) also satisfy the viscous equation, and (6a) represents the inviscid Meyer type subsonic-supersonic accelerating flow. The arbitrary constant b locates the sonic point Z = 0. Again the behavior of solutions of (5) in the Z", Z', Z phase space can be established by studying the two dimensional trajectories obtained when Z is held constant, and there will be singularities where the inviscid solutions pierce the Z = constant planes. The point Z' = wla, Z" = 0 will be a saddle-point for all Z while the point Z' = - w2a, Z" = 0 will be an unstable node, an unstable spiral point, a stable spiral point, and 4

a stable node respectively for Z corresponding to the ranges > 2a( + 2); 2a(w + 2) > Z > 0; 0> Z >- V2a(w + w); - V2c(1 + w2) > Z Thus any solution starting near the inviscid accelerating solution will diverge from it for all Z; however, some of these solutions pass through a maximum and then asymptotically approach the inviscid decelerating solution. Numerical solutions of equation (5) representing stages in the transition from the Taylor to the Meyer type of flow are shown in figures 1. a, b, and c for a = 1. 0, 0. 5, and 0. 1, and were obtained by choosing initial values very close to the accelerating inviscid solution and lying on the directrix of the saddle-point in the corresponding Z = constant plane. Integrating equation (5) once yields z - 2(zz')+ 2o(w1 - w2) Z + 2ww2 crS = C, (7) and initial conditions Z(So), Z'(S ) and Z"(So) were chosen so that the constant C1 = 0 for then it follows from equation (7) that the transitional solutions will be asymptotic to Z = - w2oS as S - + oo and to Z = wlCS as S - - ao. In figure 1, Z represents the nozzle centerline velocity distribution so that Z = 0 corresponds to the sonic point. As in the two dimensional case, figure 1 shows the gradual development of what appears to be a shock wave 5

as the maximum of Z increases beyond the sonic value. With increasing Z max the velocity gradient steepens in the region of transition from supersonic to subsonic flow. The expansion scheme which provides the basis for the derivation of the viscous-transonic equation will be valid only if U, and hence Z are 0(1); therefore, the solutions with Z greater than 2. 0 to max 3. 0, while of interest with regard to the overall behavior of equation (5), cannot accurately represent the transition from the Taylor to the Meyer flow. As the parameter a decreases the supersonic-subsonic transition shifts to larger values of S for a given value of Zax With a = 0.1 these transitions max appear to closely approximate a normal shock wave with almost uniform upstream and downstream flow. Weak normal shock waves are to order E symmetrical with respect to the sonic point for if the stream velocity u1/a* = 1 + EU1 the velocity u2/a* downstream of the shock will be 1 - EU1. Supposing Zmax to be U1 it can be seen that, as in two dimensions, nozzle flow transitions overshoot the corresponding downstream Hugoniot value U = - Z. For subsonic Taylor z I~nmax type flows with S < 0, Z diverges from the inviscid solution Z = claS very slowly, but for large positive values of S the solution Z(S) very rapidly deviates from the inviscid solution. On the other hand, even for S >> 1 the solutions approach the decelerating inviscid solution Z = - w2oS very gradually. This behavior can be verified analytically by studying the asymptotic behavior of Z near the two inviscid solutions as discussed below. 6

3. ASYMPTOTIC BEHAVIOR Although equation (5) could only be solved numerically it is possible to analytically determine the asymptotic behavior of Z(S) where it lies near the inviscid solutions. Thus from Z = owlas + 1 (8) Z = - 2aS + 2 it follows that the perturbations,' 2 << 1 respectively where Z asymptotically approaches the inviscid solutions Z = wlaS and Z = - w2aS. Substituting equations (8) into equation (7) for Z and dropping terms of 0(1 ) and O(C2 ) then yields the following linear differential equations for 1 and C2 1 - 2co-S S1 - 2w2ua1 = 0 (9) 52 + 2w2aS C2' + 2cla(2 = 0 Equations (9) are readily reduced to Weber's equation and have the solutions 1 S2 2=e 1 2 1 2 1 (10) 1 2 W2 CraS _ 11 2-= e 5-U +2 ) +-C6V - ' S) 7

where U(a, S), V(a, S) are parabolic cylinder functions of S with parameter a as defined and tabulated by Miller (1964); and C3, C4, C5, and C6 are arbitrary constants. The numerical results for Z(S) (Figure 1) indicate a difference in the behavior of C1 and C2 for large S and this can be verified from equations (10). Thus, using the asymptotic expansion for U and V as S - + oo (Miller, 1964) and keeping only the largest term in each expansion it follows that w2 21 (S) CS 1 + 0(1/S) + C4S e' 1 + 0(1/Sl] -waS (11) 1 21 2 - 2 21 21 2 C (S) C 5S 2e [1+ O(1/S) + C6S 1 + 0(1/S) For a given SO the choice of (So), 2 (S), (S),and (S) determines the constants in equations (10). Asymptotic expressions for the derivatives 1' (S) and I2' (S) are thus required and can be determined by differentiating equations (9) and determining the asymptotic behavior of the solution of the resultant Weber's equation for 1' and 2' with the result that as S - oo 8

12 W2 (1' (S) - 1 + O(1/S2) + 2w1aC 4S e 1 + 0(1/S2) (12) w 1 ( ~1 w l 3 2 s2 ~o C ' (S) - 2wC S 2 C 2 [ + 0(1/S2) - CS 2 1 + 0(1/S2) The need to differentiate the asymptotic expansions (11) directly, a procedure generally not valid (De Bruijn, 1958), is thereby avoided. The form of (12) is such that C3, C4, C5, and C6 are the same as the constants in equation (11), and equation (11) can be recovered by integrating the asymptotic expressions (12). The constants in equation (10) can now be evaluated with the result that "2 -2 S 0 W, co aO4 (SS2 - S 2) (S) = D E S -e (13a) w1 / w2 2 - a(s2-So l2 w2 - (S)= F Se + G (S (13b) (s)= s2 - _2_ (s2W-So \D 1o 2 OaSOj+2 2 2 2 9

_ I_ ^10 '2,10 2 f 2 (a S o2 S 0 2W, 2aS 2 2 o2+ a1S... 2| 2 C I O 0 G 22So 2a 22S02 where the subscript zero indicates values at S = S0. Equation (13a) shows that, except for the special case E = 0, C1(S) - 21/ 1 increases very rapidly for S > S as do the numerical solutions in Figure 1. The term (S/S) dies out with increasing S and so will have little effect upon 1. Equation (13a) is consistent with the saddle-point behavior near Z = couS discussed previously, for w1(S) will tend to zero with increasing S only in the special case i10 = - (w1/w2) So 10' for which E = 0, otherwise 1 increases with increasing So As S decreases, i. e., S < S, the exponential part of (13a) decays rapidly but the - c/o 1 (S/Se) term increases unless D = 0. The special solutions D = 0 and E = 0 thus correspond to solutions lying on the directrices of the saddlepoint, and the solution which approaches Z = cOlCS as S decreases is clearly 10

the one with D = 0. The numerical integrations were started near Z = w1aS and at points lying on the directrix of the solution diverging from the point Z' = w1a, Z" = 0, in the appropriate Z = constant plane, a procedure equivalent to choosing D = 0. Nevertheless, because of round-off errors the solutions always diverge from the inviscid solution when integrating backward from the starting point. Both terms of the expansion for (2(S), equation (13b), tend;to zero with increasing S in accord with the fact that the point Z' = w2a, Z" = 0 1/2 ___ is a stable node in the Z = constant plane for Z < - 2a12 2w + 1. The exponential term decays rapidly when S > SO so that %2(S) will be dominated - W1/ 02 by (S/S0) and so approaches the inviscid solution Z = - w2oS relatively slowly as do the numerical solutions in figure 1. The asymptotic expansions for U(a, S) and V(a, S) presented by Miller (1964) are not valid for S < 0; however, by using the expansion for U(a, -S) (Whittaker and Watson, 1952) together with appropriate recursion formulas for U and V and redefining the standard solution for Weber's equation in the cases S < 0 equations (13a) and (13b) for 1 and C2 can be shown to be 2 2 valid for S - - oo. Now, however, for S > SO, S < SO so that the exponential term in equation (13a) decays rapidly while the power term (S/S7) increases. Again in accord with the numerical curves, Z deviates from the inviscid solution very slowly with increasing S as SO - -- oo, particularly since c2/o1 = 0. 382 in the axisymmetric case. 11

From the inviscid solution of Tomotika and Hasimoto (1950) for Z it 2/- //- / 2 can be shown that C1 ~ il and C2 Isl as S - oo The exponential terms in equations (13a) and (13b), which account for the difference in the behavior of the subsonic and supersonic solutions, thus reflect the effects of viscosity. The discussion of asymptotic behavior given above is fully applicable to the two dimensional case (Sichel, 1966) provided the two dimensional values c, a2 = 2. 0, 1. 0 are used in place of the axisymmetric values w' w2 = (/-5+ 1), (5 — 1). 4. NATURE OF THE FLOW To compute the streamlines corresponding to the solutions in figure 1 the vertical velocity, V, must be known, and can be determined from the irrotationality condition and equation (7) with the result 2 + 42'12 3R3 31 c2 V = 2cRZ + 4aRX+ 2 R -2-2 R + - (14) The arbitrary constant C2 is taken to be zero since only solutions with finite V on the axis R = 0 are of interest. The constant C1 merely translates the origin of S and so, for convenience, will be taken as zero. Unlike extended to the two dimensional case. To the present order of approximation 12

dr- / 2 2 dx s/ V (15) dx 2 where rS(R) is the variation of the radius r along any given streamline. In terms of Rs(x) equation (15) becomes dR o 2 i dRS 2 )V (16) The significance of E in (16) and the relation between solutions in the X-R and physical x-r plane are the same as in the two dimensional case (Sichel, 1966). Generally it is desirable to specify the ratio of throat radius rt to the radius of curvature of the nozzle wall streamline. To the present order of — 1 approximation the streamline curvature, L, at any point in the flow is given by 1 -S 3/2 (Y + 1/ dV 23/2 y+ 11/2 A L dx2 2 dx 2 where L is the streamline radius of curvature. For most cases of interest the transition from supersonic to subsonic flow occurs well beyond the throat; hence, the flow in the immediate vicinity of the throat follows the inviscid solution, Z = o1oS, so that (L/I ) is given by 13

(L/7) = 2E32 (1/2)(y + 1) Ra (w1 + 2)1 (18) With the ratio rt/ Lt of throat radius to wall radius of curvature fixed at t t3, the throat radius rt will then be A 1/2 -tr // (19) t Aea + 1)(w1 +2) 1/2 The axial throat coordinate xt is derived using the condition V = 0 at the throat as in the two dimensional case (Sichel 1966). Isotachs, or lines of constant speed in the r-x plane correspond to curves along which U is constant to the present order of approximation, and so can be determined from equation (4) and the numerical solutions for Z. In the region of inviscid flow the velocity perturbations EU and /3/22"/2 EU - 1U - +- 7 (20) [(y + 1) (w + 2)] 2 3/_2 -2(w1 + 1) 1p3/2 [(1/2)(y + 1)]1/2 3/2(1i/2)(, + 1) V = v T1 + [ 2)3/2 (21) + 2)_3/2 when rt is chosen as a reference length so that (/rt); 71 - (/r=( t) 14

Integration of equation (15) using (21) for v/a* then yields the following expression for the wall streamline in the portion of the flow where Z = wC uS - 1/2 r _ q/~\np2K 4(wb 1) 13/2r1/2 r (1/2)2 et (Deo 22 rt e 3- +2)] dX (22) = _to 2 [2( +2)]3/2 J With t, the throat coordinate, given by 21/2 (w) + 1) r1/2 = (xt/r) =- 3/- (23) (t ( t/r t [2(w + 2)]3/ Equation (22) shows that r] - co for sufficiently large, but for those solutions with Z, and U O(1) transition to subsonic flow occurs long before 71 diverges to infinity. Once Z deviates from Z = oloS equation (15) can only be integrated numerically. The relation between the arbitrary constant a and the nozzle flow field can be seen from equations (18), and (19). With E, 7r, and y fixed the streamline radius of curvature L varies inversely with a for a given fixed radius R, and the velocity gradient of the inviscid solutions ahead of and behind the viscous transition decreases as is evident from figure 1. As 15

discussed previously with this decrease in velocity gradient the viscous transition, at least on the nozzle centerline, approaches the Taylor (1910) shock transition. For fixed 3 and E the nozzle throat radius rt varies inversely with a, (equation (19)). Typical isotach contours and streamlines for a = 1. 0 and 0. 1, corresponding to curves A and B in figures l(a) and l(c), are shown in figures 2(a) and 2(b) in the 5-7j plane. Figures 2(a) and 2(b) have been drawn using numerical solutions with the same peak value of Z, and with wall streamlines corresponding to i3 = 0. 21 as in the paper by Tomotika and Hasimoto (1950). Since the flow is axisymmetric the isotachs are really the intersections of constant speed surfaces with planes through the nozzle axis. The strange shape of the isotachs downstream of the region of rapid deceleration in figure l(b) results from the slight increase in Z(S) immediately behind the shock like transition. The a = 0. 1 wall streamline has a second minimum some distance downstream of the supersonicsubsonic transition; however, for streamlines with I8 sufficiently small this second throat disappears. An inherent property of similarity solutions, such as presented here, is, of course, the inability to specify streamline shapes a-prioriL The shock like nature of the supersonic-subsonic transition is clearly indicated particularly in the case of a = 0. 1. 16

5. DISCUSSION The axisymmetric nozzle similarity solutions are quite similar to the two dimensional solutions found previously (Sichel 1966). In the present case the asymptotic behavior of the solutions for the centerline velocity Z has been investigated, and the difference in the behavior of the numerical solutions in regions of subsonic and supersonic flow has been verified analytically. The effect of the parameter a upon the nozzle solutions has been examined. Figure 2(b) shows that when a << 1 solutions are obtained such that there is essentially a weak normal shock near the axis with effects of the wall and shock curvature occurring only for sufficiently large r as was anticipated previously. 17

REFERENCES Cole, J. 1949 Problems in Transonic Flow. Ph.D. Thesis, California Institute of Technology. DeBruijn, N.G. 1961 Asymptotic Methods in Analysis. Amsterdam: North-Holland Publishing Co. Guderley, K.G. 1962 The Theory of Transonic Flow. Reading: Addison Wesley. Hayes, W.D. 1958 Fundamentals of Gas Dynamics, Sec.D. Princeton: Princeton University Press. Miller, J. C.P. 1964 Parabolic Cylinder Functions. Chapt. 19. Handbook of Mathematical Functions (ed. M. Abramowitz and I.A. Stegun). U.S. National Bureau of Standards, Washington, D. C. Sichel, M. 1963 Phys. Fluids, 6, 653. Sichel, M. 1966 Journ. of Fluid Mech. (in press). Szaniawski, A. 1963 Archiwum Mechaniki Stosowanej, 15, 904. Taylor, G.I. 1910 Proc. Roy. Soc. A, 84, 371. Tomotika, S. and Hasimoto, Z. 1950 J. Math. Phys. 29, 105. Tomotika, S. and Tamada, K. 1950 Quart. Appl. Math. 7, 381. Whittaker, E.T. and Watson, G.N. 1952 A course of modern analysis. 4th ed. Cambridge, England: Cambridge Univ. Press. 18

Z(s), Z (s) 6.0- i 6.0- - /2 - ' - 0 4.0 o-= 1.0 4.0 -1 o=0.5 2.0 o- / \ Pi -2.0 2.0 4.0 6.0 -2.0 /2.0.0 6.0 8.0 Curve A /..' (a) (b) Figure l(a), (b), and (c). Numerical Solutions of the Equation IZ' - 2ZZ" - 2(Z' - w o) (Z' + o 2o) = 0

Z(s) 6.0 — 4.0- Curve B 2.0 _ \\ S.-4o — - 0- - -4.0 - ao. 0 8.0 00. 12.0 14.0 16. 18.0............... (c) Figure 1 (Concluded)

0.5 0.0 -0.5 1.0 i A — 1.0 0.0 1.0 2.0 -0.14 y|.l4 U/'t ItU I-t/Lt 0.21 (-= I Figure 2(a). Isotachs and Streamlines Corresponding to Curve A in Fig. 1 (a) 00 Throat 0 'T 1 '?S0.5 0.0 1.0 2.0 2.0 - 1.5 ~\ \ \\1.0 777,.oq~s~p^^\(..?/// / /-'-' Fig. l(c) OO I.0. 2.0 Figz.,(c)