THE STRUCTURE OF WEAK NON-HUGONIOT SHOCK WAVES by M. Sichel, Y.K. Yin, T.S. David Final Report on Grant DA-ARO(D)-31-124-G385 and Contract DA-31-124-ARO-D-276 with The Army Research Organization Durham, North Carolina Report 07146-3-F The University of Michigan Department of Aerospace Engineering Gas Dynamics Laboratories OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR January 1968 Distribution of this document is unlimited.

FOREWORD This report is a combined Final Report for work sponsored by the Army Research Organization in Durham, N. C., both under Grant DA-ARO(D)-31124-G385 and the succeeding Contract DA-31-124-ARO-D-276. This support is gratefully acknowledged. The object of the research program summarized here was to examine the effect of viscosity on transonic flow, and was directed by Martin Sichel, the chief investigator. Following the Introduction, Sections II andIIIof the report describe the personnel who have been active in this program, and summarize the research which has been completed. Unpublished studies of the mathematical properties of linearized forms of the viscous transonic potential equation, xx -x + yy = 0, are described in detail in Section IV. Solutions of the equation 6 + 0y = 0, xxx yy on a finite domain by separation of variables are described. The asymptotic behavior of solutions of the equation ~ + P2x + 0 =0 for large values xxx xx yy of y has been developed. Finally difference and exact solutions of 0 + 0 = 0 xxx yy on a semi infinite domain have been compared, as preliminary step in the finite difference solution of the full viscous transonic equation. ii

CONTENTS Page FOREWORD ii I. INTRODUCTION 1 II; PERSONNEL 4 III. PUBLISHED WORK 4 A. Papers 4 B. Reports 5 C. Papers presented at annual meeting of the Fluid Dynamics Division of the American Physical Society 5 D. Seminars 5 IV. ANALYTICAL AND NUMERICAL SOLUTIONS OF THIRD ORDER PARTIAL DIFFERENTIAL EQUATIONS 6 A. Solution of xxx + yy = 0 on a Finite Domain by Separation of Variables yy 6 1. Formulation of an Eigenvalue Problem 6 2. Properties of the Eigen Expansion 14 3. Numerical Example 18 B. A Study of the Linearized Viscous Transonic Equation + 2 0 +- ( = 0 22 xxx x yy C. Numerical Investigation of xxx+ = 0 on the SemiInfinite Domain y > 0 x YY 28 BIBLIOGRAPHY 38 iii

I. INTRODUCTION In the study of transonic flows the use of the inviscid transonic equation coupled with the treatment of shock waves as Rankine-Hugoniot discontinuities sometimes leads to paradoxial results or solutions in serious disagreement with experiment. These difficulties usually arise when the application of the above theory leads to large velocity gradients or discontinuities in streamline curvature, and it appears that a viscous transonic theory must be applied in such cases (Cole 1949, Sichel 1962 and1963 Szaniawski 1963, Ryzhov and Shefter 1964). It has been shown (Sichel 1963) that regions of transonic flow where the effect of compressive or longitudinal viscosity is of the same magnitude as convection are governed by a viscous-transonic equation, which, in normalized form (Sichel 1966), is U - 2UU +V =0 sxx x y U =V y x or written in terms of a normalized potential, q q -f p +A =0 xxx- xx xx yy In this viscous transonic or V-T equation U and V are normalized x and y components of the perturbation to the sonic velocity. The undisturbed velocity is presumed to be near sonic and in the x direction and x and y are normalized coordinates The ivestigation of the above viscous-transonic equation was the subject of the investigation carried out under ARO Grant DA-ARO(D)-31-124-G385, and Contract DA-31-124-ARO= 1

D-276. The question to be answered was whether physically meaningful solutions of the V-T equation resolving some of the difficulties of the inviscid theory could be found. The investigations, summarized in this final report, have shown that the V-T equation does indeed provide a reasonable description of those transonic flows in which compressive viscosity is important. One problem in which the inviscid transonic theory fails is in explaining the nature of the transition from purely subsonic to subsonic-supersonic flow in a converging diverging nozzle. Exact similarity solutions of the V-T equation, which do show the nature of this transition have been obtained (Sichel 1966, Sichel and Yin 1967). In the case of radial and spiral flows the velocity gradient becomes infinite near the sonic radius; however, exact spiral-radial flow solutions of the V-T equation have been found which show that this difficulty can in part be eliminated if viscous effects are taken into account (Sichel and Yin 1967). External transonic flows about two-dimensional airfoils are, perhaps of greater practical interest because of their relations to the performance of transonic airfoils. The question of how the viscosity affects the far field has been investigated by Ryzhov and Shefter (1964) and Szaniawski (1966). Near the airfoil surface there is often a region of supersonic flow imbedded in an external subsonic flow. While potential flow solutions for such supersonic regions, in which the flow decelerates smoothly to subsonic velocities, have been found, experiments almost always indicate that the supersonic region is terminated by a shock wave. This conflict between theory and experiment led to the well known "transonic 2

controversy", extensively discussed by Nieuwland (1967), and Manwell (1963), as to whether shock free pockets of supersonic flow embedded in a subsonic flow can exist. Since it can be shown that the inviscid solutions for such flows are not continuously dependent on the boundary data, it was concluded that such shock free pockets of supersonic flow could not exist. As indicated by Nieuwland (1967) this controversy has been reopened by the experiments of Holder (1964), and Pearcy (1962) which showed that shock free supersonic pockets could be attained under special conditions. Viscous effects may play a role in this controversy, particularly since experimentally it is impossible to say whether the flow under observation is a "pure potential" flow or a viscous flow. The discontinuous dependence on boundary data, on which nonexistence of supersonic pockets is based, may disappear if viscosity is properly taken into account. When the supersonic pocket is terminated by a shock it is found that the pressure rise across the shock may be appreciably less than the Rankine-Hugoniot value (Holder 1964). This is not surprising for application of the Rankine Hugoniot conditions to a weak shock normal to a curved surface leads to infinite streamline curvature where the shock touches the surface (Emmons 1946)o Emmons suggests that the shock structure is no longer one dimensional near the wall and that a viscous theory is required to properly treat this portion of the shock. It would appear that solution of the V-T equation for two dimensional external flow past bodies of arbitrary shape might shed some light on the difficulties discussed above. As a preliminary step in this direction approximate solutions for 3

V-T flow past a wavy wall have been obtained (Sichel and Yin 1967). In the viscous case it is found that supersonic pockets are terminated by smooth, shock like, transitions from supersonic to subsonic flow. Several mathematical aspects of the V-T equation have also been explored. This work is described in Section IV below. II. PERSONNEL M. Sichel, Principal Investigator, Associate Professor of Aerospace Engineering T. S. David, Graduate student, active 1962-1964, after which he left the University B. J. Florsheim, Graduate student, active summer, 1965 Y. K. Yin, Graduate student, active 1965-presento (Since this work is being continued under a new contract I anticipate that work completed so far may contribute toward Mr. Yin's Ph.D. degree.) R. S. B. Ong, Professor of Aerospace Engineering. Professor Ong has acted as a faculty consultant on several occasions. III. PUBLISHED WORK A. Papers "The effect of longitudinal viscosity on the flow at a nozzle throat, by M. Sichel, Jour. of Fluid Mech., Vol. 25, Pto 4, p. 769-786 (1966). "An axisymmetric similarity solution for viscous transonic nozzle flow, " by M. Sichel and Y. K. Yin, Jouro of Fluid Mech., Vol. 28, Pt. 3, p 513-522 (1967). 4

"Viscous transonic flow past a wavy wall, " by M. Sichel and Y. K. Yin, to appear in Fluid Dynamics Transactions Vol. 4, (Proceedings of biennial symposium on Advanced Problems and Methods in Fluid Dynamics; held in Tarda, Poland, September 18-23, 1967). B. Reports "Viscous transonic spiral and radial flow, " by M. Sichel and Y. K. Yin, University of Michigan ORA Report, August 1967. "The effect of longitudinal viscosity on the flow at a nozzle throat, " by M, Sichel, University of Michigan ORA Report 05800-1-F (1965). "An axisymmetric similarity solution for viscous transonic nozzle flow," by Mo Sichel and Y. K. Yin, University of Michigan ORA Report (1966). C. Papers presented at annual meeting of the Fluid Dynamics Division of the American Physical Society "The effect of longitudinal viscosity on the flow at a nozzle throat, " presented November 1964 at Cal. Inst. of Tech. by M. Sichel. "Viscous transonic spiral and radial flow," presented November 1966 at Stanford Univ. by M. Sichel. Do Seminars (presented by M. Sichel) "Viscous transonic nozzle flow, " Cornell University, April 1965. "Viscous transonic flow, " University of Michigan, December 1966; University of Pennsylvania, April 1967; Aberdeen Proving Ground, May 1967. 5

IV. ANALYTICAL AND NUMERICAL SOLUTIONS OF THIRD ORDER PARTIAL DIFFERENTIAL EQUATIONS The viscous transonic potential equation Lx -^ ^ +cf =0 (1) oxxx - x xx yy 0 is non-linear and of third order. Various approximate and iterative methods of solving Eq. (1) require the solution either of the equation fxx + y = 0 (2) xxx yy or of the linearized V-T equations 6 i2 0 +0 = 0 (3)'xxx xx yy o +c + = 0 (4)' xxx -' x yy While the literature on second and fourth order partial differential equations is extensive references on third order partial differential equations are extremely limited. Hence various aspects of the behavior of equations (2), (3) and (4) have been investigated as described below. A. Solution of 0xx + kyy = 0 on a finite domain by separation of variables 1. Formulation of an eigenvalue problem The solution of Eq. (2) on the semi infinite domain, y > 0, has been discussed extensively by Sichel (1961). Here we are concerned with solution of (2) on a finite domain by separation of variables. A key objective of this work is to establish the behavior of the eigenfunctions which occur in this problem. 6

The problem considered is the following: ^k +pf =0 xxx yy (O, y) = <(7r, y) = Ox(0, y) + b 6x(Y, y) = 0 (5) 4(x, O) = h(x); f(x, ) = 0 i. e. homogeneous boundary conditions are chosen. Assuming a solution of the form k =X(x) Y(y) it follows that the functions X and Y must satisfy the ordinary differential equations X"' + p X = 0 (6) Y - p3Y = 0 (7) where p is a separation constant. The third order ordinary differential equation (6) has been discussed by Kamke (1944) who cites several other papers dealing with eigenvalueproblems related to this type of equation. The solution of Eq. (6) can be written in the form wlPX w2PX w3px X = Cle + C2e +C3e (8) where w1 = - 1, w2 = (1/2) + (V372)i, and w3 = (1/2) - (i/32)i are the three cube roots of - 1. C1, C2, and C3 are arbitrary constants. The boundary conditions (5) require that X(O) = X(7) = X'(O) + bX'(7r) = 0 (9) 7

and from (9) it follows that unless 1 1 1 wl p w2pff w3 p7 e e e =0 (10) P'l~libbYPjj p4(1iw2P~) w3 (1+ beW3Pi pwl ( +be 1 ) pw2(1 +be ) pw +be the solution X will be the trivial one with C1 = C2 = C3 = 0. The characteristic equation (10) will determine the eigen values p. In general p will be complex so it is convenient to let p = 5 + ir? where 5 and 71 are of course real constants. Now the characteristic equation (10) can be written in the form F(5,Ti) + iG(,77) = 0 (11) which implies F(, qT) = G(, ri) = 0. A further-reduction is possible in that F(, T7) and G(,, 71) can be written in the form F(, 71) = 7f(5, T) - 7g(~, 7) = 0 (12a) G(5, 7) = 71f(, 7) + 5g(, 71) = 0 (12b) and from (12) it follows that f(, 1 0) = 0 (13) g(0,)7) =0 (14) provided (E + 11 ) t 0. 5 = TI = 0 is also a solution of (12) but this special case will be discussed separately below. The functions f(, Tr) and g(, ri) are given 8

as follows -7Th3 - (5 + V77) v3 1 1 f(, 7) = e sin - sin - e sin + - be sin - 2 -.7/ - ->tt-47), 3 1 b+ e sin -2 + 27 - 3 7 -+ bej( g(,7) = e- cos nrT - e cos | ^4 - -7 + - JU-ff ^ 1 1 - e cos 2 +2?7+ 3 -be 5 cos T/71 (16) -2( Q 3Z) 3 1 I + be- cos ( —^ - (6)7 + be cos 2 + 2 3 Equations (12a) and (12b) must be solved simultaneously to determine values of 5, and 77. This was accomplished graphically using the method of "grid plotting" described by Hartree (1958). 9

From the boundary condition (9) it follows that the constants C1, C2, and C3 must be in the ratio _ w~p7 2ff w3p l C1:C2 3 P (w -w3) pb +b w 2 e w3 e: W3pU w3pPf P(w3 - w1) +pb W3 e3 - wl e )]: (17) ~~r Wwlp7T w2P7 p(w -w2) + pb wl e W2 e To correspond with (17) we now define the eigenfunction X, corresponding to the particular eigenvalue pn as follows /=/ \I, ow2Pnn 3 n W l nPnX Xn IPn(W2 - w3) + Pnb w2 ew eW P e + Pn(W3 w1) + pnb w3 e w e nj e2 (18) F w 1p W2Pn7T W 3Pn + Pn(Wl - w2) + pbwl e -w2 e e Now Eq. (7) has the solution p3/2 -3/2 Y=K1 e +K2 e (19) It will be found later that one of the constants K1 and K2 must be chosen zero to obtain solutions that remain bounded as y -co. For reasons discussed later 10

the choice K = 0 was made. Then it follows that p3/2 = K e n X (20) 2 n In order to satisfy the condition f(x, 0) = h(x), K2 must represent the constants An in the expansion of h(x) in terms of the eigenfunctions X, such that 00 E W An X n (21) h(x) = A Xn (21) n=1 In view of the strange nature of the functions Xn, it is initially not evident that arbitrary functions, h(x), can be expanded according to (21). To determine the constants A it is necessary to use the property n X (x) V(x) dx =; m i n (22) n m 0 where V (x) is the adjoint eigenfunction corresponding to X (x) V (x) is a m m m solution of the adjoint equation V - p V= 0 (23) satisfying the boundary conditions V(0) = V(7T) = bV'(0) + V(7r) = 0 (24) The characteristic equation for determining the eigenvalues of the adjoint problem is 11

1 1 1 -Wl1PT -w2Pf -WT3PT e e e 0 -W1Peb -W2P7 -w3pl -Pw1(b + eWP - Pw2(b + e ) - pw3 b + e (25) wlplr Upon multiplying the first row of the above determinant by e, the second by w2P'7 e and so on it can be shown that Eq. (25) is identical to Eq. (10), so that the eigenvalues of the adjoint problem are identical to those obtained from Eq. (10). As in the case of X (x), it is readily shown that the eigenfunction V (x) can be defined as (x) r Pn W3pn7 W2P n7r - Wlpn X) = Pn(w3 - w2) + w3 -w2 e - wP e e Pn -W Pn 7T -W Pn7T -W2Pnx + Pn (w -w2) + bT1 w e - w3 e e (26) p 1P-W2+ Pn -w -we P +nW2 - W1) + w2 e 1w3PX It should be mentioned that the definitions of Xn (x) and Vn(x) are by no means unique, and that no attempt has been made to express these functions in normalized form. From Eq. (22) it now follows that 12

7f h(x) V (x) dx A = (27) n If X () V (x) dx 0 As mentioned above T = 71 = 0 is a solution of the characteristic equation (10). We shall let p1 = 0 be the corresponding eigenvalue. In this special case X1(x) cannot be obtained by simply substituting p1 in Eq. (18) for Xn(x); rather, it is necessary to return to the original differential equation (6). For p = p1 = 0 it then follows that 1 = C1 + C2x + C3x (28) Application of the boundary conditions (9) now yields C1 = 0 and requires that 1 T = 0 (29) b+ 1 27r if C2 and C3 are to be different from zero. Equation (29) can only be satisfied for b = 1 and then since C2: C3 = -7: + (1/7r) it is possible to define X1 as 2 X = —; b= 1 (30) X1=0; b 1 13

2. Properties of the Eigen Expansion Simultaneous solution of Eqs. (12a) and (121J for n and 77n, the real and imaginary n n parts of pn is a major source of difficulty in the boundary value problem of Eq. (5). However certain properties of the expansion (21) and of the eigenvalues pn can be deduced without actually solving Eqs. (12a) and (12b), and these are considered below. a. Eigenvalues pn occur as complex conjugates. Proof: Inspection of Eqs. (15) and (16) shows that f(, r-) = - f(, -r) (31) g(t,77) = g(5, -7r) From Eqs. (13) and (14) it then follows that if (n, 77 ) is a solution of the characteristic equation, ( n - n) must also be a solution so that the eigenvalues occur as pairs of complex conjugates. bo The expansion (21) is real provided g(x), the constant b, and the variable x in (9) are real Proof: Using the property w1 =w1; w3 = w2; w3 =w2 (32) which follows from the definition of wl, w2, and W3, with the bar denoting the complex conjugate, it is readily shown from Eq. (18) that Xn(x, n) = -Xn(x, (33) 14

Similarly from Eq. (26) it follows that Vn(x, ) - Vn(x, Pn) (34) Then if h(x) is real equation (27) for An (p) yields the result A (Pn) - A (Pn) (35) so that An(Pn) Xn(X pn) = An(Pn) Xn(X, p) (36) Since except for n = 1 the eigenvalues pn occur as pairs of complex conjugates, and the eigenfunction and the adjoint corresponding to p1 = 0 are real it follows from (36) that the expansion (21) for h(x) will be real. c. The characteristic equation (10) in symmetrical about rays from the origin of the (5,7r) plane making angles of 27r/3 and 4,7/3 with the 5 axis. Proof: From the following properties of the cube roots of - 1:. 7T 2 3. 27r w2 = - w3 = e (37) w3w2 = - 1 = 1 and the characteristic equation (10) it is readily shown upon writing the characteristic equation as 15

() = 0 (38) that (-w2P) = Up) (39) Then from (39) 4/(+w2 p) = (-W2P) = (p) (40) 4 2 8 2 Since - w2 = exp (i ), + w2 = exp (i- ) = exp (i - r) the above symmetry follows. It should be observed that not only 4(p), but its real and imaginary parts 4 a 2 F(5,7 ) and G(5,7 ) will be symmetrical with respect to exp (i 3 7r) and exp (i- 7r). d. F(S, T7) = Re[l(p)] = 0 on the rays: -1 2w 7 2S 27w 71 0 = tan -=0, 3, 7T, 3, - L-0 3 3 3 3 for all values of b. In the special case b = 1,F(Q, 71) is also zero on the rays =tan1 _ 7 5 X 7 5t7 0 - tan d =2 "6 2'6 6 2 6 Proof: From Eqs. (12a) and (15) it follows that F(~, 0) = 0 for < K 0. From the symmetry properties (c) it then follows that F(5, T7) is zero on rays at angles 7T 27T 2w7 7T 3 n 3 7 " 3 3' In the special case b = 1, F(O, T7) = 0 for all 77 < 0. Then using symmetry arguments similar to (c) it follows that F(5, T) will be zero on rays at angles 7T 7w 57 7T T 5w7 6' 2' 6' 6' 2' 6 16

e. Special properties when b = 1: The case b = 1 which implies the symmetry X'(O) = - X'(7r) in the boundary conditions results in considerable simplification and so merits special discussion. First, if the function h(x) is symmetrical about x = (1r/2) so that h(x) = h [(7r/2) - (x - 1r/2)] = h(7r - x), then b must have the value unity to be consistent with the boundary conditions. From property (d) and the grid plot it follows that with b = 1, the eigenvalues will lie on the rays 0 = tan (1/A) = + 7T/6, + r/2, + 57/6. If 0 is restricted to -7 < 0 < ir, then if the solution (20) is to remain bounded the eigenvalues lying on 0 = + 1r/2 and + 571/6 must be rejected, and only the values lying on 0 = + 7r/6 need to be used. If K2 were chosen as zero in (19) then the eigenvalues on 0 = + 7T/2, + 5?7/6 would have to be used; thus the choice K1 = 0 is the simplest to make. The fact that certain eigenvalues have been rejected raises the question whether the expansion (21) for h(x) is complete; this question remains to be investigated. Since the eigenvalues, pn, lie on the rays + = + 3r77 it is readily shown that the real parts of Xn(x) and Vn(x) are symmetrical about rT/2 while the imaginary parts of X (x), V (x) are asymmetrical about f/2. Thus (41) Re +[V(x)] + i Im V(x)] = Re [V ( - x)] - i Im [xn (T - x) Then if h(x) is symmetrical about 7T/2 it follows from 17

7T h(x) Vn(x) dx n i7f Xn() Vn(x) dx 0 that the expansion coefficients A will all be real. n 3. Numerical Example The expansion 00 h(x)= A X (21) n=1 was carried out in detail for the function h(x) = sin x; < x < 4 4 (42) h(x) =0; 0< < 4; 4 <x<7 Since this function is symmetrical about x = 7T/2 the choice b = 1 appears appropriate; however the function (42) is in some sense degenerate since b 1 will in no way conflict with the boundary conditions. However, as shown previously, when b I 1 the eigenfunction Xl(x) = 0, and then the series (21) does not seem to converge to h(x). In the present case the eigenvalues Pm and their complex conjugate Pm were found to be 18

PP =(m-j2)(+~i); 2m=2,3,4... i. e. two eigenvalues correspond to each value of the index m. Since A (pm) = - A(pm) and the A's are all real in this case Am(pn = - Am(pm) so that we have only one [Am I corresponding to each m. It is found that Am takes the form -3 (m - )7r A =a e; m = 2,3,... n m with the coefficients a as tabulated below: m a m a m ________m 2 0.17240 7 -0.01547 3 -0.02640 8 -0.02761 4 -0.06400 9 0.01217 5 0.02070 10 0.02151 6 0.03884 11 -0.01005 12 -0.01775 From Eq. (27) it also follows that A = 0. 9943. Values of the function X exp [- 3 (m - ) 7r] are given in Table I below. In terms of the index m m the series expansion (21) becomes 00 h(x) =AX1 + 2 Re [Am (x) (43) m=2 19

the factor two being necessary since two eigenvalues correspond to each index m. Counting the eigenfunction corresponding to p = 0 only once has not been justified rigorously, although this procedure appears to yield reasonable results. The sub-sums Si(x)=A1X1 +2 Re[A Xm(x)] m=2 are shown in Fig. 1 for various values of i and are tabulated in Table II below. The series appears to be converging toward h(x) except near x = 7r/4, 37T/4 where h(x) is discontinuous. Here the series appears to be approaching the mean value of the function. A number of questions remain to be answered, among them the following a) Are there restrictions on the type of functions h(x) which may be expanded according to (21). b) Can a more rigorous justification for neglecting certain eigenvalues be found? c) Is the expansion (21) complete? 3/2 d) How does b behave for y > 0? Since pn will be complex the complex antisymmetric part of X will affect the behavior of f for y >0. n It appears that f will cease to be symmetric with respect to 7r/2 for y >0. 20

Table It. Re{Xn exp [- (n - 3) x - 0. 0 7i/8 37r/16 7f/4 37i/8 7r/2 n 2 10 -1.000 -0.934 -0.583 0.430 0.950 3 -0.785 -0.062 -0.828 -0.249 1.000 4 -0.074 0. 934 0. 485 -0. 867 5 -0. 658 -0. 599 0. 870 -0. 967 6 0. 973 -0. 501 -0. 493 -0. 499 7 -0. 694 0. 991 -0. 870 0. 255 8 -0.005 -0.263 0.495 0.866 9 |0.708 -0.791 0.869 0.967 10 -1. 000 0.866 -0.493 0.500 11 0.707 0.134 -0.870 -0..253 12 0.000 -0.965 0.495 -0.865 tThe tabulated function is symmetrical about x = 71/2. i Table II. S(x) = AX1 +2 Re [A X (x)] m=2 x - O. 0.0 /8 37i/16 7r/4 37r/8 7r/2 S 0. 0 0.342 0.472 0.586 0.732 0.781 S2 -0.003 0.150 0.385 0.880 1.126 S -0. 035 0. 034 0. 367 1. 004 0. 945 S6. 0.014 0.010 0.365 0.926 1.064 S 0. 035 -0. 044 0. 364 0. 870 0. 978 S 0. 0.10 -0.026 0. 364 0. 915 1. 045 S12 -0. 005 0 006 0. 364 0. 951 0. 990 h(x) 0.0 0.0 0.0 0. 354 0. 925 1. 000 Si(x) is symmetrical with respect to n/2. 21 21

B. A Study of the Linearized Viscous Transonic Equation xx + 2 + = 0 xxx x yy If the viscous transonic equation (1) is linearized by replacing the velocity 2 ox in the convective non-linear term by a constant 3, then the VT equation becomes qx - 2+ + 2 + =0 (44) xxx xx yy with plus sign corresponding to subsonic and the negative sign to supersonic values of x. The Dirichlet boundary value problem O(x, 0) = f(x) (45) I O(x,l | I<M will be considered. If X(A, y) is the Fourier Transform of +(x, y) such that -00 -_00 then it is readily shown from the properties of Fourier Transforms (Sneddon, 1951) that b must satisfy the ordinary differential equation 2 (-i)3 4(X, y) ~ 2 (-iX)2 4~(X, y) + = 0 (47) dy so that /i3 2 3/ 22,. ) C e-iA + A y+C -iAk + A y (48) ~(X,y) =C1e - 2+Ce y (48) 22

Since Re [-iA3 A2 <0 it follows that C2 = 0 if, according to (45), 4 is to remain bounded as y -ca Arg ( -iX i 23 A ) is restricted to angles between 0 and 21r. Applying the boundary condition at y = 0 yields 00 C; =F(X) = V= jf() eid dF (49) -00 where F(h) is, thus, the Fourier transform of f(x). The solution q(x, y) is obtained by taking the inverse, Fourier transform of 4(X, y) such that (x, y) = (A, y) e dX -00 and in the present case this becomes co 00 (x, y) J f(2) e7-i i( -X dX dS (50) -00 -00 where it has been assumed that the order of integration between - and X can be reversed. In the present case it is more convenient to write Eq. (50) in the form 23

00 00 ~(x, y) = Re ( y) ei x) d d 0 -oo (51) i 3 22 =T Rell f( ) e -iy + Ye - i(-x)d d -o0 0 So far no way has been found to carry out the integration with respect to X in Eq. (51) in closed form; however, it is possible to determine the 3 22 asymptotic behavior of O(x, y) as y - o Since Re (-iA ~ 3 A ) < 0 it follows from the exponential term in (51) that the integrand will contribute to ~(x, y) only for X << 1 if y >> 1. In the subsonic case, with X << 1, it is 3 2 2 1/2 then possible to expand (-iX + 3 A ) y according to o 9 1/2 i 2 3 (-iX X2) +y-3 —+X 2 1. (52)_ /y -y 84 - + 16 6 It will be assumed throughout that P ~ 0(1). The real part of the integrand of (51) then can be written Re e /-iX 3 2 i( 2 3 - e cos A(t - x) - 2 sin (t - x) - cos XA - x) 2P 8 3 X4 + XY sin X( - x) + O(4y2 ) (53) 165J 24

In carrying out the expansion (53) it has been assumed that even with y >> 1, X2y << 1; thus the results will be valid only for extremely large y. Keeping only the two leading terms in (53) and carrying out the integration with respect to A in Eq. (51) now yields the result 00 f2y + (x - d) 00 2 2 2 1 f f( (x -)y [3y - (- -) ]d + (54) 2T 2,2, 3 -00 ~ [a3y + (x - )] Equation (54) represents the first terms in an asymptotic expansion with the first term of O(1/y) and the second term of O(l/y ). The first term of (54) is a solution of the inviscid equation 2 + 0 3' k +4 =0 x sxx yy with 4(x, 0) = f(x) The situation is different in the supersonic case when the expansion of 3 22 21/2 (-iA - 32X21/) y yields 3 2 2 1/2 i 2 (-iA. 22) /y i (55) Then the leading term of the integrand of Eq. (50) for the inverse transform becomes 25

21 2 -ix3 _- 2 y i[ y + (-x) 2 [ ] (56) e e h iX P ih4x=e e It is now possible to write 1 2 (A, y) F(X)e e 2 (57) Since 00 -1 F(A) eiY= F(X) e-(-Y dX -oo00 - i(x - 3y) (58) where - denotes the inverse Fourier transform, and since - 1 A x22 -1 2 2y (59) it follows from the convolution theorem (Sneddon 1951) that 2 co (x-1) 23 = / f~ f(i - y) e 2y (60) -00 With a simple change of variables Eq. (60) can also be written in the form I ( - y - e d (61) -00'26 26

Again the expression (61) for f is the first term of an asymptotic series, The nature of this asymptotic solution can best be indicated by considering several specific and simple forms of the function f(x) = O(x, 0). Considering first a pulse type disturbance f(x) = a; xl < x f(x) =0; 1xl > x (62) f(x) = (a/2); xl = xl it follows from (61) that -=L erfc (x-x - x y) / - erfc (x+x1- y) / (63) The solution (63) represents a decaying pulse propagating along the line x = 3y, which is one of the characteristics of the inviscid equation 2 - 2 xx + + = 0 As a second example the initial function'xx yy f(x) = a; x < f(x) = a/2; x = (64) f(x); x > 0 yields the solution = erfc (x - y) Y (65) 2 2y which represents a transition from k = O to q = a which propagates along x = 3y and broadens with increasing y. Of coursejf(x) does not satisfy the 27

o00 condition J If(x) Idx required for convergence of Fourier transforms -00~~~ ~~~~~~~~1 AX2 (Sneddon, 1951) however because of the expoential factor exp (- 2 - y) in (57), the inverse transform exists nevertheless. This latter point has been discussed by Sommerfeld (1949), for example. The asymptotic results obtained above in both the subsonic and supersonic cases appear reasonable from a physical point of view. However, the behavior of the transform (AX, y) in the complex A plane should be investigated in greater detail to lend more rigor to the results derived above, and higher order terms of the asymptotic series for y > 1 should be investigated. C. Numerical Investigation of + 4 = 0 on the Semi-Infinite Domain y -— >O —-----.xxx yy The complexity of the V-T equation (1) suggests that numerical methods may have to be used in the solution of certain problems. As a preliminary step in this direction the finite difference solution of the equation (2), xx + yy = 0 on the semi-infinite domain, was investigated. Since analytical *xxx yy solutions are available for this linear problem (Sichel, 1961), it can be used as a means of testing various finite difference methods. As the highest order terms of equation (2) are the same as those of the V-T equation, there is some hope that finite difference schemes which are applicable to equation (2) may also apply to the V-T equation. The work described below is still in a preliminary stage. 28

If central differences are used, then the partial differential equation + =0 (2) xxx yy can be replaced by the finite difference equation 1i+1,j 3i,j + 3i-1, j i-2,j + i,+1 ij + i, -1 0 (66) h k The grid points and grid spacing associated with this equation are shown in the sketch below. (i, j+1) 0 k 0 _ 0 Y Y (i-2, j) (i-1,j) (i, i) (i+l j) i y 0 (i,j-l) x This grid is consistent with the results of the uniqueness proof for Eq. (2) (Sichel, 1961) which shows that on a rectangular domain > must be specified on the closed rectangular contour surrounding this domain, except for the left hand boundary on which it is also necessary to specify x, Equation (66) can be rewritten as w 3w w 1 i,j 2 + 3w i+l,j 2 + 3w ~i-l,j 2 + 3w i-2, j 2 + 3w i, j+1 2+3w -= 0 (67) 2 + 3w 9i, j29

where w = (k /h ). It is also possible to consider (66) as a relation for i j in terms of the values of 0 at neighboring grid points so that 1 1.(w + 3w. + 0. + - w. (68) i, j 2 + 3w (w i+1, j 3 i-,, j+1, -- wi2, j The finite difference equation (67) is actually a system of simultaneous algebraic equations for the unknown nodal values 0i, and the finite difference solutions of (2) requires a solution of this simultaneous system of equations. This system of equations can also be expressed in matrix form. Suppose the finite difference grid contains.m rows and n columns of unknown nodal points. Let u be the n component vector 0s's2'. sn' i.e. the unknown s s1's2 sn nodal values of 0 in row s, Then the system of equations (67) can be written Az = k (69) where z is the vector Ul, u2,.. o u and k is a constant vector with mn components which reflects the boundary conditions. Since each u has n components, z is also an mn component vector. The mn x mn matrix A can be expressed in terms of submatrices A1 and B of rank n x n such that A1 B 0 0 A = | B A B 0 0 B AI B 0 0 B A B Al B 0 O 0 B A1 B 0 B A.1 30

In the case of Eq. (67) 1 -(2 3w) A = _3w 1._w 2+ 3w 2 + 3w |w_ 3w w O2+3w "2 + 3w 2 + 3w * w 3w w 2 + 3w 2 + 3w 2 + 3w(1 o (23) ~~(2 + 3w while -B -2 +3w I (72) where I is the n x n unit matrix. In principle Eq. (69) could be solved using Cramer's rule or, what is the same thing, by inverting A so that - -1 -- z =A k (73) In practice this procedure is impractical if the number of grid points, mn, is large, for then the inversion of A becomes very difficult and time consuming, even on a high speed computer. Thus for large mn Eq. (69) must be solved iterativelyo The Gauss-Seidel method (Varga 1962, Todd 1962) is an iterative method which is particularly well suited to the use of a computer. Therefore this 31

method was applied to a solution of the boundary value problem (x, o) = 1 I; xl.< 1'. O(x, o0) =; Ix > 1 (74) +(x, 0) = 0.5; Ixl= -0 as x-+oo, y-+oo on the semi-infinite domain y > 0 since an analytical solution is also available in this case (Sichel 1961). A key question is whether the Gauss Seidel iteration procedure, or for that matter any other, when applied to Eq. (69) will converge at all. The convergence properties of any iterative scheme depend on the properties of the matrix A. In particular Varga(1962). Theorem 3. 4, proved that if A is a strictly or irreducibly diagonally dominant n x n complex matrix then the Gauss Seidel method is convergent. From Varga, Definition 1. 7: "An n x n complex matrix A = (a.i) is diagonally dominant if n ai, i> la.1 (75) j=1 j/i for all 1 < i < n. An n x n matrix A is strictly diagonally dominant if strict inequality in (75) is valid for all 1 < i < n. Similarly, A is irreducibly diagonally dominant if A is irreducible and diagonally dominant with strict inequality in (75) for at least one i". The matrix A defined by Eqo (70) is neither 32

strictly or irreducably diagonally dominant and so does not satisfy the conditions of Varga's theorem. Since Varga's theorem only provides a sufficient condition, the failure of A to satisfy the conditions of Varga's theorem does not necessarily mean that the Gauss-Seidel method will diverge. Varga (1962) also presents a number of theorems applicable when A is Hermitian; however, A as defined in (70) is non-Hermitian. Clearly more work is required to prove convergence of iterative schemes in the present case. In spite of the above difficulties it was considered worth while to study the convergence of the Gauss-Seidel method by applying it to the boundary value problem (74). Although a semi-infinite domain is under consideration only a finite number of grid points covering a finite region can be used and this may be a source of considerable error. To eliminate this error the exact solution, which is available in the present case, was used to determine the boundary values at the first two columns, at the last column, and at the last row of the rectangular grid of size (n+2)h in the x direction and (m+l)k in the y direction. In the present case there are n columns and m rows of unknown grid points Two iterative finite difference calculations and the corresponding exact solutions are prescribed in Tables III, IVa, and IVb. These results have been used to construct the graphs of q(0, y) vs y and q(x, 1) vs x in Figso 2 and 3. Preliminary indications are that the Gauss-Seidel iteration converges, though apparently very slowly, but this result is not a proof of convergence~ 33

Additional work is clearly necessary in establishing finite difference methods for the V-T equation. A proof of the convergence of the GaussSeidel or other iterative methods should be found or developed in the present case. Results of the direct method of solving (69) should be compared to the exact solution for a small enough grid to make the direct method tractable. The possibility of stretching the y coordinate to take care of the boundary condition at y -ooshould be examined. Thus the transformation rt = e - changes Eq. (2) into 6o + T7 1 + 71 A = 0 (76) with boundary conditions applied at r] = 0 and?7 = 1. 0. Now of course 7U = 0 is a singularity of Eq. (76). 34

Table III. Comparison of Exact and Iterative Solutions; m = 7, n = 8, h = 1. 0, w = 0. 5 Exact Solution y 1 -5.0 -4.0 -3.0 -2.0 -1.0 0 1.0 2.0 3.0 4.0 5.0 x0.000.000.000.000.000 500 1.000 0. 500.000.000.000.000 0.707.011.019.036.096.272.605.651.262.016.000.000 1.414.020.036.061.114.226.393.504.391.161.032.001 2.120.031.044.067.118.196.301.389.381.246.105.031 2.830.034.052.075.115.174.249.316.336.278.170.070 3.540.039.053.077.111.158.214.267.294.271.200.115 4.240.039.055.077.106.145.190.233.259.254.212.142 4.950.042.056.076.102.134.171.207.232.236.211.162 5.660.043.057.075.098.126.157.187.210.218.204.170 c,1 Gauss-Seidel; 75 Iterations* y i -5.0 -4.0 -3.0 -2.0 -1.0 0 1.0 2.0 3.0 4.0 5.0 x0. 00.000.000.000.000.500 1.000.500.000.000.000.000 1. 71.011.019.038.099.344.616.518.232.072.014.000 1.41.020.036.062.124.268.431.444.307.155.057.001 2.12.031.044.073.128.223.331.371.313.206.106.031 2.83.034 o052.078 124.193.269; 313.295.228.145.070 3 54.039.053. 07 8 o 18 71.229.269.271.234.174.115 4.24.039.055.078.111.158.199.235.248.230.190.142 4.95.042 o056..077.104 138.176.209.227.224.199.162 5.66.043.057 075.098.126.157.187.210.218.204.170 *Values inside rectangle determined by iteration. Values outside rectangle are boundary conditions taken from the exact solution.

Table IVbo Comparison of Exact and Iterative Solutions; m = 13, n = 8, h = 0. 5, w = 0. 5 Gauss-Seidel; 75 Iterations* y -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 x0.00.000.000.000.500 1.000 1.000 1.000.500.000.000.000 0.25.026.047.115.381.713.861.847.588.245.068.000 0.50.049.081.155.326.553.710.735.598.364.165.032 0.75.062.098.168.293.455.591.640.573.419.247.107 1.00.073.109.171.268.391.503.560.536.438.306.187 1.25.080.116.170.248.344.437.495.496.437.343.249 1.50.086.116.165.231.310.387.442.457.426.361.291 1.75.085.118.161.217.283.348.399.422.409.366.308 2.00.090.118.157.205.261.318.364.390.389.364.325 2.25.092.118.152.195.243.293.335.363.370.356.329 2.50.092.117.148.186.229.272.311.339.351.345.326 2.75.093.116.144.178.216.255.291.318.333.333.319 3.00.093.114.140.171.205.240.273.300.316.320.310 3.25.093.113.137.165.195.227.258.284.301.307.301 3.50.092.111.133.159.187.216.244.269.287.296.291 *Values inside rectangle determined by iteration. Values outside rectangle are boundary conditions taken from the exact solution.

Table IVa. Comparison of Exact and Iterative Solutions; m = 13, n = 8, h = 0. 5, w = 0. 5 Exact Solution y -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 x0.00.000.000.000 0.500 1.000 1.000 1.000 0.500.000.000.000 0.25.026.047.115.309.701.923.867.667.262.016.000 0.50.049.081.150.287.506.730.783.665.423.163.032 0.75.062.098.162.263.419.585.682.635.485.278.107 1.00.073.109.159.249.364.490.584.594.490.348.187 1.25.080.116.164.235.325.425.508.539.494.387.249 1.50.086.116.162.222.296.378.449.487.471.396.291 1.75.085.118.158.210.273.341.403.443.443.397.308 2.00.090.118.155.200.255.313.367.405.414.387.325 2.25.092.118.151.192.239.289.337.373.387.373.329 2. 50.092.117.147.184.225.270.312.346.363.357.326 2.75.093.116.144.177.214.253.291.322.341.341.319 3.00.093.114.140.170.204.239.273.302.321.325.310 3.25.093.3.137.164.195.227.258.285.303.310.301 3.50.092.111.133.159.187.216.244.269.287.296.291

Bibliography Cole, J. 1949, Problems in Transonic Flow, Ph.D. Thesis, California Institute of Technology. Emmons, H.W. 1946, NACA TN 1003. Hartree, D. Ro 1958, Numerical Analysis, Oxford, at the Clarendon Press. Holder, D.Wo 1964, Jour. of Royal Aero. Soc. 68, 501. Kamke, E. 1944, Differentialgleichungen Losungsmethoden und Losungen, Leipzig, Akademische Verlagsgesellschaft. Manwell, A.Ro 1963, Arch. Ratl. Mech. Anal. 12, 249. Nieuwland, Go Yo 1967, Transonic Potential Flow Around a Family of Quasi Elliptical Aerofoil Sections, Ph. D. Thesis, University of Groningen, Netherlands. Also NLR Technical Report T. 172. Pearcy, H. Ho 1962, The Aerodynamic Design of Section Shapes for Swept Wings, Advances in Aeronautical Sciences, Vol. 3, London. Ryzhov, O.S. and Shefter, G.M. 1964, Prikl. Mat. i Mekh. 28, 996. Sichel, M. 1961, Properties of the Equation - + vyy = 0- Princeton University, Dept. of Aero Engrg., Report 1166. 1962, Phys. Fluids, 5, 1168, 1963, Physo Fluids, 6, 653. 1966, J. Fluid Mech., 25, 769. Sichel, M. and Yin, Y K. 1967a, J. Fluid Mech., 28, 513. 1967b, Viscous Transonic Spiral and Radial Flow, University of Michigan ORA Report. 1967c, Viscous Transonic Flow Past a Wavy Wall, to appear in Fluid Dynamics Transactions Vol. 4. (Proceedings of biennial symposium on Advanced Problems and Methods in Fluid Dynamics. Held in Tarda Poland, Sept. 18-23, 1967. 38

Sneddon, I. N. 1951, Fourier Transforms, New York, McGraw Hill. Sommerfeld, A. 1949, Partial Differential Equations, Academic Press. Szaniawski, Ao 1963, Archiwum Mechaniki stosowanej, 15, 904. Szaniawski, A. 1966, "The Asymptotic Structure of Weak Shock Waves in Flows over Symmetrical Bodies at Mach Number Unity, " Royal Institute of Technology, Stockholm, Sweden, Report KTH No. 70. Todd, J. 1962, Survey of Numerical Analysis, New York, McGraw Hill. Varga, R. S. 1962, Matrix Iterative Analysis, Englewood Cliffs, N. J., Prentice Hall Inc. 39

Y-N |/ \ Plot of the Subsums I.o i — ~. y-'".' Si (x) 2A A X(x) /... nn ~/ //^^"^^^\\ ^v ---— ~~S,(x) old i) _.6-.........S x) It1'.4 " 27' \I T7 7r 3nr 7r 57r 37r 77r 7r 8 4 8 2 8 4 8 Figure 1. Plot of Eigen-Expansion Sub-Sums.

COMPARISON OF EXACT AND NUMERICAL SOLUTIONS GAUSS-SEIDEL 75 ITERATIONS BEHAVIOR AT X=0 -.0 — Exact Solution +++++ m=7, n =8,h =.O,w=0.5 00000 m=13,n = 8,h =0.5,w=0.5 0.8 BOUNDARY CONDITION 0.6- (xo) =;Ix I > I 0.4 +. \^ q(x,o) =0.5;1x= 0.4 - 0.2 O0 1.0 2.0 3.0 4.0 5.0 6.0 y Figure 2. Comparisons of Exact and Numerical Solutions.

COMPARISON OF EXACT AND NUMERICAL SOLUTIONS GAUSS-SEIDEL, 75 ITERATIONS BEHAVIOR AT Y=1.0 --- Exact Solution 00000 m=13,n=8,h=0.5,w=0.5 0.6 0 0.5 0.4 0.3 / BOUNDARY CONDITION 0.2 /(xvo); Ixl<l f(xo) =O;Ixl>l (C). I }</ ^ (x,o) =0.5; x > =l 0II - (xo)I I.I I 5 1 1 I -3.0 -2.0 -1.0 0 1.0 2.0 3.0 x Figure 3. Comparison of Exact and Numerical Solutions.

Unclassified Security Classification DOCUMENT CONTROL, DATA - R & D (Security classification of title, body of abstract and indexing annotation must be entered when the-overall report is classified) 1. ORIGINATING ACTIVITY (Corporate author) 2a., REPORT SECURITY CLASSIFICATION The University of Michigan Unclassified Ann Arbor, Michigan 2b. GROUP 3. REPORT TITLE The Structure of Weak Non-Hugoniot Shock Waves 4. DESCRIPTIVE NOTES (Type of report and inclusive datee) Final Report 5. AU THOR(S) (First name, middle initial, last name) Martin Sichel Thomas S. David Y.K. Yin 6. REPORT DATE 7 TOTAL NO. OF PAGES OF REF January 1968 42 22;. CONTRACT ORF GRANT NO. go. ORI'GNATOROS REPORT NUBI ER(S) DA-ARO(D) —31-124-G385/ b. PROJECT NO. DA 31-124-ARO-D-276 07146-3-F c. 9b. OTHER REPORT NO(tS (Any other numbers that may be aaaigned this report) cL 10. DISTRIBUTION STATEMENT Distribution of this document is unlimited. 1. SUPPL.EMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY The Army Research Organization Durham, North Carolina 13. ABSTRACT The object of the research program summarized here was to examine the effect of viscosity on transonic flow, and was directed by Martin Sichel, the chief investigator. Following the Introduction, Sections II and II of the report describe the personnel who have been active in this program, and summarize the research which has been completed. Unpublished studies of the mathematical properties of linearized forms of the viscous transonic potential equation, kxxx - ox oxx + kyy = 0, are described in detail in Section IV. Solutions of the equation 4xxx + kyy = 0, on a finite domain by separation of variables are described. The asymptotic behavior of solutions of the equation Oxxx ~ 32 4xx + 0yy = 0 for large values of y has been developed. Finally difference and exact solutions of oxxx + Pyy = 0 on a semi infinite domain have been compared, as preliminary step in the finite difference solution of the full viscous transonic equation. DD lPOr 473 REPLACLS DO FORM 147. 1 JAN 6,. WHICH IS H U MOV.1 41 / OBSOLETE FOR ARMY USE. Security Classification