LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 I£6r no. 188-199 cop.SL Digitized by the Internet Archive in 2013 http://archive.org/details/numericalsolutio190shao Report No. 190 ID. «PV J. lor o * 9° NUMERICAL SOLUTION OF PLANE VISCOUS SHOCK REFLECTIONS by Tzu-Sien Shao September 30, 1965 li M " JAN 07 t%6 Report No. 190 NUMERICAL SOLUTION OF PLANE VISCOUS SHOCK REFLECTIONS* by Tzu-Sien Shao September 30, 1965 Department of Computer Science University of Illinois Urbana. Illinois This work was submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Mathematics, August, 19&5, and was supported in part by the AEC under contract US AEC AT(ll-l)lU69 and by the NSF under contract NSF-GP hG^G . ACKNOWLEDGMENTS The author wishes to thank Professor C. W. Gear of the Department of Computer Science for his guidance in the research and help in the prepara- tion of this work and his wife, Jin-Ing, for her patience and understanding. Thanks are also extended to the members of the drafting department of the Computer Science Department for their assistance in providing the final drawings and illustrations for this work, to the operational personnel of the ILLIAC II computer for their help in obtaining computer output and to Mrs. Anna Rita Ferris who typed the manuscript in its final form. -in- TABLE OF CONTENTS Chapter p age I INTRODUCTION 1 II DEVELOPMENT OF SHOCK REFLECTION THEORIES 6 2.1. Failure of Local Theories 6 2.2. A Semi -Global Shock Reflection Theory 8 2.3. Global Theory with True Viscosity 17 III DIFFERENTIAL EQUATIONS SUITABLE FOR NUMERICAL SOLUTION. . . 19 3.1. Time -Dependent Equations Versus Stationary Equations. 19 3.2. The Cell Method and the Point Method. 20 3.3. Differential Equations in Lagrangian Form 23 3.^-. Equations in Streamline Coordinate System 25 IV NUMERICAL SOLUTION OF VISCOUS SHOCK REFLECTION PROBLEMS . . 35 k.l. Stability of One -Dimensional Difference Equations . . 35 h.2. Numerical Solution of One-Dimensional Viscous Shock Problem 39 h.J). Stability of Two-Dimensional Difference Equations . . hk h.k. Numerical Solution of Regular and Mach Reflections . 50 V CONCLUSION 65 BIBLIOGRAPHY. 69 VITA. . 71 -IV- NUMERICAL SOLUTION OF PLANE VISCOUS SHOCK REFLECTIONS Tzu-Sien Shao, Ph.D. Department of Mathematics University of Illinois, 19^5 ABSTRACT A numerical method is presented for the solution of two-dimensional, viscous fluid flow problems . Emphasis is placed on the formulation and solution of weak shock reflections. The basic equations used are the time -dependent Navier-Stokes equations without heat conduction. A version of the shock re- flection theory dividing flow into various regions is discussed. A "behavior similar to experimental results is displayed "by the numerical results of this theory, although they do not quite agree with the experimental data. To show that the Navier-Stokes description of the fluid flow is an adequate model for the weak shock reflection problems, equations of motion are numerically inte- grated by a difference scheme to obtain asymptotically stationary solution. A transformation of independent variables is first made to reduce computation time. Stability analysis is performed to give an upper bound for the time step. Numerical results are given for asymptotically stationary weak shock reflections. Reflected shock angles for both regular and Mach reflections obtained agree well with those given by shock tube experiments. CHAPTER I INTRODUCTION It is known that almost all hydrodynamical phenomena cause a development of "discontinuities/ 1 so-called shocks and boundary layers, sooner or later. The essential difference between a shock and a boundary layer may be picturesquely described as follows. A streamline that enters a shock comes out again downstream, whereas a streamline that enters a boundary layer remains within a region of shear. Throughout this work, we shall consider shocks only. Boundary layers will be ignored completely. Physically speaking, the shocks are regions in the flow where higher order derivatives of certain flow varia- bles containing small coefficients, e.g., coefficients of viscosity and heat conduction, play an important role. Due to mathematical difficulties, these derivatives are usually neglected to simplify the partial differential equations describing the flow. This causes the flow variables to become discontinuous across the shocks which now are curves of zero thickness. Using the integral form of the equations of motion, the Rankine-Hugoniot equations can be derived to determine the jump conditions. This is the classical approach. It works quite well for simple shock problems. When a plane shock strikes a smooth rigid wall, it gives rise to regular reflection for small angles of incidence and to an irregular type of reflection, known as the Mach reflection, for large angles of incidence. In regular reflection or the two-shock case, the reflected and the incident shock meet at a point on the wall (Fig. l.l). The flow direction changes through the incident and reflected shocks must be equal and opposite to satisfy -2- the conditions of flow parallel to the wall. The Mach reflection is charac- terized by the intersection of three shocks at a point detached from the wall (Fig. l.l). This is known as the three-shock case. The shocks are called incident, reflected and Mach shocks respectively. REGULAR REFLECTION MACH REFLECTION INCIDENT SHOCK ^REFLECTED SHOCK ////////J////////// WALL INCIDENT SHOCK REFLECTED SHOCK MACH SHOCK WALL FIGURE 1.1 REGULAR AMD MACH REFLECTIONS For this class of shock reflection problems, theories based on the classical approach do not work very well. In particular, they fail badly for Mach reflection with weak incident shocks. There is a wide discrepancy between the reflected shock angles predicted by the three-shock theory and Smith's [ 1^+] experimental results obtained in shock tubes. Another anomaly 'irst observed by Smith is the persistence of regular reflection somewhat nd the extreme angle where the two-shock equations cease to apply. Since three-shock theories deal with shocks in the neighborhood of the inter- Lnt, we shall refer to them as "local" theories. Presumably local -3- theories do not work because they do not consider the influence of the wall and the flow downstream of the intersection region. One then hopes to include these effects by enlarging the intersection region. Theories based on this idea will be referred to as "semi-global" theories. The flow field can be divided into various domains to be treated separately with matching conditions along the boundaries . A version of the semi-global theory based on the Navier- Stokes description of fluid flow will be presented in Section 2.2. It does not yet give satisfactory results, although it seems to work better than local theories. Therefore one must solve the whole problem including boundaries and downstream flow flied. This gives rise to a "global" approach. The partial differential equations of compressible, viscous and heat conductive fluid flow, namely the Navier-Stokes equations, with boundary conditions are too complex to be solved analytically. Numerical solution seems to be in order. Many authors [3>8] have numerically solved stationary fluid flow problems with artificial viscosity and no heat conduction from a global point of view. Due to numerical difficulties, they have integrated time-dependent equations of motion to obtain asymptotically stationary solution. For simple shock problems, their results agree quite well with classical theory and experimental data. However, no method seems to work for three-shock problems with weak incident shocks « Since the true viscosity effect in the flow is more important in this case, we shall include it in the differential equations describing the flow. A numerical method for solving asymptotically stationary, two-dimensional, viscous fluid flow problems from global point of view will be presented in this work. Emphasis will be placed on the formulation and solution of weak shock reflections. The purpose is to determine whether or not the Navier-Stokes equations are an adequate model in this case. Since the shocks are weak, the corresponding changes ■4- in flow variables are small throughout the flow region. Therefore we will neglect heat conduction terms in the equations of motion without effecting the physics of weak shock reflection problems too much. The coefficient of viscosity may be assumed constant for the same reasons. The viscous time- dependent equations are nonlinear and parabolic, and therefore very time- consuming to solve numerically even for simple shock configurations. Some efforts have been made to make our method more efficient and easier to apply than the existing ones. The main features of the present methods are: [l) it gives a faster technique for numerical solution of shock reflection problems and it is directly applicable to stationary viscous fluid flow problems in general; (2) the results for weak shock reflections agree very well with Smith's shock tube experimental data. The latter is a strong evidence for the adequacy of Navier-Stokes description of shock reflections., In Chapter II, a survey of previous works on shock reflection problems will be presented together with a discussion of their failure for weak three-shock cases. A new semi-global theory will also be given there. In Chapter III, the time-dependent, non-heat-conductive Navier-Stokes equations will be transformed into a streamline coordinate system. A more detailed explanation of using time-dependent equations to obtain asymptotically stationary solutions will be given. The equations of motion will be written a form suitable for numerical solution. Chapter IV forms the main text this work. Complete numerical procedures for solving one-dimensional and two-dimensional viscous fluid flow problems will be derived and stability for both cases will be performed. Numerical results of the one- nal viscous shock problem and plane shock reflections made on the iversity of Illinois ILLIAC II computer will be reported. Plots of flow -5- variables drawn "by a CalComp 67O digital plotter will also be shown. Graphs consisting of data obtained from local theory, semi-global theory and our numerical solutions as compared to experiments will be given. Good agreement for weak shock reflections between our results and Smith's experimental data will be demonstrated. CHAPTER II DEVELOPMENT OF SHOCK REFLECTION THEORIES 2.1. Failure of Local Theories In the absence of a complete shock reflection solution due to mathematical difficulties, von Neumann [18] suggested a local three-shock theory to apply to the neighborhood of the intersection of the incident , reflected and Mach shocks. The idea is that the streamline passing through the triple point divides the flow into two contiguous regions . (This dividing streamline is usually referred to as the slip stream.) Combinations of Mach and reflected shocks can be found so as to produce the same static pressure and flow direction in the two regions downstream of the intersection point. Essentially the same ideas give a solution for regular reflection (two-shock theory) which corresponds with the experimental observations. Unfortunately for weak incident shocks, there is an apparent wide discrepancy between the reflected shock angles predicted by the three-shock theory and Smith's [Ik] experimental results obtained in shock tubes. Stimulated by the discrepancy of the local three-shock theory, a number of asymptotic solutions under valid limiting conditions have been found [1,9]- However, none of these solutions is suitable for explaining the failure of the theory for reflected shocks with finite strength. As a result of the limiting conditions, the reflected shock is forced to have zero strength at the triple point. Later, it was felt [2] that perhaps the theory fails for curved shocks rather than weak shocks. Taub [16] worked out the effects of curvature and its -6- -7- derivatives in the neighborhood of the triple point. Clutterham [k] calculated reflected shock angles as a function of the angle of incidence and the incident shock strength based on the curved shock theory. Nevertheless this local theory of curved shocks still suffers from any weakness that the three-shock theory suffers. This is due to the fact that it contains the three-shock theory in the limit when the neighborhood under consideration shrinks to the triple point. In 1959* Sternberg [15] examined the effect of shock thickness on the boundary conditions at the triple point. He showed that at the intersection there must be a non-Rankine-Hugoniot shock zone separating the three R-H shocks. According to him, the three-shock theory of von Neumann and Guderley [7] fails because they used nonviscous models. With viscosity included, there is a smooth transition of flow variables inside the non-R-H zone. This non-R-H zone acts as a buffer zone separating R-H shocks in the upper and lower domains. It then completely eliminates the singularity at the triple point appearing in nonviscous models. However, the two-dimensional character of the non-R-H zone introduces a theoretical complication. The unknown downstream elliptic field has a significant influence on the structure of the shock zone. This makes the mathematical problem indeterminate from a local point of view. Later Gear [6] tried to introduce viscosity into the curved shock theory. For very weak shocks, the Mach number M ahead of the incident shock 1 2 is close to unity. Writing M = 1 + — e , he expanded flow variables in power series of the small parameter e. For the structure of shock regions, he used a linear approximation. Retaining only low order terms in e under consistent assumptions, a system of nonlinear algebraic equations was derived relating the reflected shock angles and shock curvatures. The only numerical -8- solution of the algebraic system found prescribed a negative shock in violation of the second law of thermodynamics. Sukurai [13] also tried to use power series expansion in terms of radial component of the independent variable to integrate the viscous differential equations. His choice of boundary conditions, however, remains to be justified. 2.2, A Semi-Global Shock Reflection Theory Since all local theories fail for weak Mach reflection problems, we propose in this section a semi-global theory to improve the results. This theory includes the effects of the wall and the slip stream behind the inter- section region. The shocks are treated as domains with finite thickness to avoid singularities. As in Gear's work [6], the flow field is divided into several regions. The intersection region is thought of as a mixture of incident, reflected and Mach shocks, outside of which viscosity can be ignored. The boundary condition at the wall to which the Mach shock is attached will be assumed to correspond to the condition that the wall exerts no viscous drag on the fluid. This condition can be realized when two identical shocks inter- sect each other symmetrically. The central dividing streamline is physically equivalent to the wall exerting no viscous drag. The idea is to integrate equations of motion along the boundaries of the intersection region. A system of nonlinear algebraic equations relating the flow variables and geometric variables will then be obtained. Appropriate solution of these equations will yield reflected shock angle and Mach shock height for given incident shock. For simplicity, we use the following two-dimensional, stationary and non-viscous equations of motion for a perfect fluid in rectangular coordinates [17]. (l) Continuity equation c* ( pu ) + "oy ^ = ° ; (2.2.1) (2) Momentum equations (pu ,+ p> + >x - (puv) = , bx % (2.2.2) 3x" (pUv) + ¥ (pV + P) = j (2.2.3) (3) Energy equation pu( 2 2 u +v + -2-£) 7-1 p-J fr 2 2^' /U +v _7 V\ pv( — - — + — — — ) w K 2 7-1 p\ = (2.2. k) where 7 is the ratio of specific heats . We shall discuss the semi-global theory in two cases using different geometric models for the intersection region. They are as follows. (A) 1-Region Model ■8 X — H«s— g INTERSECTION REGION /7777777777777777T7777 D rrrrrm WALL FIGURE 2.1 SEMI-GLOBAL 1-REGION MODEL -10- Here we consider a model (Fig. 2.l) where the intersection region is connected to the incident and the reflected shocks and extended to the wall. The velocity component normal to the wall at A' and D is taken to "be zero. The shocks are assumed to be stationary and oblique to the flow direction. The flow conditions in zone 1 are assumed uniform and given. By specifying the incident shock strength £ T and shock angle CC we can calcu- late the flow conditions in zone 2 by using the conditions for a stationary- oblique shock [10]. Similarly, flow conditions at B can be determined if the reflected shock angle a is assumed in computation. Thus the unknowns in K this model are the flow variables at D, the reflected shock angle a and the K height of the intersection region h. There are actually 5 unknowns, namely p , u n , p_, & and h since v is taken to be zero. If we integrate D D D t\ D Eqs. (2.2.1) - (2.2. k) over the region ABDA.' with respect to x and y, we can get k equations relating the above unknowns. We shall illustrate the inte- gration procedure by using Eq. (2.2.1) as follows. Integrating Eq. (2.2.1) over ABDA.' with respect to x and y, we have B (pu)dy B f B B + / [pv] dx - (2.2.5) A J A U In order to numerically handle the above integrals, we first observe that the intersection region can be thought of as having a mirror image with respect to the wall. (This is equivalent to assuming that the wall exerts no viscous drag.) The following modified Simpson's rule can be derived to take advantage of the symmetry x f f(x)dx = I [f(x Q ) + 2f(x 1 )] - |- f (lv) U) , (2.2.6) o -11- where h = x and x < i < x o o o 1 This formula can be used for integrations from D to B. Integrals from A to B are handled by the standard trapezoidal rule. Using the conditions v. = v . , = v ~ 0, Eq. (2.2.5) now becomes § [(pu) B + 2(pu) D ]^ + -± (pv) I + -f [(pv) 1 + (pv) B ] = , (2.2.7) where 5 and 6 are the oblique thicknesses of the incident and reflected I K shocks respectively. They can be calculated by using the solution of one- dimensional viscous shock problem [17]- Further simplification of Eq. (2.2.7) using the uniformity condition of the incoming flow in zone 1 leads to — (pu) D = h(pu) A - ^(pu) B - — (pv) B 6 +5 I R / x -g- (pvJj. • (2.2.8) Similarly, Eqs . (2.2.2) - (2.2. k) can be integrated and written in the form 2h / 2 s _ , 2 v h t 2 s & R t v I +5 R / % — (pu +p) D = h(pu +p) A - - (pu +p) B - — (puv) B - — - — (puv). I (2.2.9) 5 I +5 R h ( , 5 I +5 R , 2 s 5 R , 2 , & R -2~ P D - 3 (puv) B + -5— (pv + p) I + - (pv +P ) B - - p A (2.2.10) and 2h 3 /u 7 V\ = h D /U 7 P\ R 2 2 /U +v 7 P\ pv( ~F- + — p } j 7-1 p y J A pv( 2 2^ / U,+v 7 Pn J B 5 +S I R B 2 2^ u +v _7_ P\ 2 7-1 p J J I respectively (2.2.11) -12- We now have k equations in 5 unknowns. The extra equation can be supplied by assimung that the flow behind the intersection region is shear- less. That is to say, the horizontal velocity components at B and D have the same value. We thus have u D = u B . (2.2.12) Eqs. (2.2.7) _ (2.2.12) form a system of 5 nonlinear algebraic equations in 5 unknowns p_, u^ p , GC and h where Ot appears implicitly through the flow D D D K K conditions at B. Further manipulations of the equations can reduce them to one single equation involving a . A numerical solution was found using the n binary chopping technique, giving a for given % and a . The Mach shock height h can also be calculated., The results for £ = 0.9; 0.8, 0.7 and a = k6 (l )70 are given in Table 2.1. The shock angles are also plotted in Figures k.^ - ^■•9- Since the equations are very nonlinear, several extraneous solutions were discarded. Fortunately the extraneous solutions can easily be distinguished. They all fall into the following categories: (l) a corresponds to no reflected shock, i.e., u T = u_; (2) a corresponds to 1 B ' K negative value of h; and (3) 0! is less than a . n I (B) 2-Region Model. Encouraged by the 1-Region model, we consider a more complicated and presumably more realistic model. There is a slip stream extending from the intersection region through BC . (See Fig. 2.2). All assumptions and conditions are taken to be identical to case (A). The unknowns here are low variables at C and D, the reflected shock angle QL,, the height of slip tream s and the distance from slip stream to the wall h. With v taken to be zero, we have 10 unknowns, namely, p Q , u , v , p , p D , u , p D , a R , s and h. ■13- SLIP STREAM S A' JNJER SECTION REGION /T7777777T77777777777T777777777 WALL FIGURE 2.2 SEMI-GLOBAL 2-REGION MODEL As in case (A), we integrate Eqs . (2»2.l) - (2.2.^) over regions ABDA" and A'CDA" separately in x and y directions. Using the conditions v. = v. , = V = V D = °' Eq ° C 2 - 2 ' 1 ) S ives h 2h . (pu) c + — (pu) D + - (pu) c = h(pu) A + s(pu) A - - (pu) B - — (pv) I 2" [(P V ). B + (p*)^ (2.2.13) and — (pu) D = h(pu) A - - (pu) c 5 I +5 R r , -g— (pv) c (2.2.1^) respectively for integrations over regions ABDA" and A'CDA". Substituting Eq. (2.2.11+) into Eq. (2.2.13) and simplifying, we can eliminate p , u , and h and set -Ik- 5 I +S R 5 I +6 R | (pu) c - -—^ (pv) c = s( pU ) A - - (pu) B — (pv) I - — (pv) B . (2.2.15) Similarly, Eqs. (2.2.2) - (2.2.4) can "be integrated and simplified to yield s / 2 x o (P u +P), S I +6 R ~2— (puv) c s(pu +p) A - - (pu +pJ B 5 +S I R / v ~2— (puv) — (puv) B , (2.2.16) and I (puv) c -g- (pv +P) C & R 2 (pUv) B & I +5 R ( 2 , 5 R , 2 , — — (pv +p) I - — (pv +p) B (2.2.17) 2 2 /U +v 7 p-v pu(-g- + ^ -) 6 I +8 R J c 2 2 /U +v 7 P\ pv(-^— + — -) 7- = s J C /U 7 P\ pu( — + — — -) N 2 7-1 p ; pus 2 2 rU. +V 7 P" 7-1 p' V & R J B 2 2 / U +v 7 P\ pv(-g- + ^ ") J I pv( 2 2 u +v ,_?._ pn 7-1 p ; _ B (2.2.18) respectively. In Eqs . (2.2.15) - (2.2.18), unknown flow variables p , u , p and the distance h have "been eliminated. We now have k equations in 6 unknowns p , u , v , p , a, and s where a again appears implicitly through L U U K R flow conditions at B. To supply two more equations for the problem, we assume that flow variables across the slip stream satisfying the classical conditions. These conditions can be derived by integrating Eqs. (2.2.1) - (2.2.4) across BC and letting s -* 0. They are simply P C = P B 111 Ib u c (2.2.19) (2.2.20) -15- Eqs . (2.2.15) - (2.2.20) now form a system of 6 equations in 6 unknowns. Further manipulations and simplifications can again reduce them to a single equation in a . Numerical solutions are given in Table 2.1 for | = 0.9, 0.8, 0.7 and a = h6°(l O )70°. In Figures 4. 7 - k>9, we shall plot the results of 2-region model as compared to 1-region model, Smith's experiments, global theory and local theory. The extraneous solutions discarded in this case fall into the categories: (l) CC corresponds to no reflected shock; (2) QL, R R corresponds to negative: values of s or h; (3) oc is less than a, and (h) a K In corresponds to unreasonably large values of s and h. The calculated shock reflection angles in both cases do not quite agree with Smith's experimental data as we can see in Figures ^.7 - ^■•9* They seem to give better results than local theories for weak three -shock intersections, however. Strangely, examination of Table 2.1 and Figures ^.7 - K.y, shows that the 1-region model consistently gives better results than the 2-region model, although the heights of intersection region in both cases agree fairly well as shown in Table 2.1. We do not yet have a good explana- tion for this phenomenon. The assumption of shearless flow behind the inter- section region used in 1-region model does not seem to have much physical justification. This more or less arbitrary assumption might accidentally give better answers. In the mean time, conditions (2.2.19) - (2.2.20) used in 2-region model may not be proper when viscosity is present. A better understanding of the viscous effect along the boundaries of the intersection region is needed to improve the results. It is conceivable that some approach similar to the semi-global theory discussed in this section could be useful for shock intersection problems. The flow field could be divided into various regions where certain flow parameters may be more important than others in a given region. • 16- o II H 4JU1 H 0) O d o •H bD - O CM L/N CO HLT\ _;)- _4 -4 Lr\LrALf\Lr\Lr\Lr\VOVQVOVO t — t- — C— L^CO CO H QJ Ti O o •H bD (D K i H ^ VO h-00 ON H -4 VO O -4 OnVO KN OJ -4 CO CO N~\ na, ka r- ON H KN UA t> ON OJ LTNCO OJ -4 Lr\UAir\l/NirNLr\VOVOVOVOVO t^t^-L^-co CO o II H rH - OA-4 rH rH (ACAONOl ChH l"-C0 VQ ON H H OJ ua OJ ON OJ O O NA t-~ H ua ON KNCO -4 OAVO OJ ON [>- ua ir\ ua VO CO OJ t^VO CO t-00 ON H OJ KA UAVO CO ONH K>4 VOCO O OJ J" VO OA H -4- D— _4 .4 _4 Lr\Lr\uAUALr\uAUAVOvavovovo t>-r--l>-[^L"-cococo rH 0) O O •H bO « i H rC OJ KA r-ONO OJ NA UAVO CO O H KA ua l>- ON H KA ua CO H ua _4_4 ualTni/AUAUAUAVOVOVOVOVOVO [T- [— l"— O-CO CO ; 6 ii H •xn rH CD o O •H bD VD O4OV0iALr\CM^aiO on on 0N0>a\0Na\0Na\oooHHaiai^4iAh-o4 KA rH 4COK>OON.ONH4rfMM40000 t^OO IA00 ON ON ON KA |>-0N0JLr\t>-O-4t-HUA0N-4-C0-4-0NLrN0JV0[>-V0 C*- O OJKNLrNNOt^ONOHKN4 LA f-OO O H KNLAVflCO O -4" UA LT\UAir\UALrNLTNVOVOVOVOVOVDVOI>-L^-L^-f~-[>-f— CO H QJ o o •H bD i _c« OOOHH[M4lAC0 04COaiOOLA44l s -CJWhOaJO\l^ OOOOOOOOOHHHC\lC\JKN4lAVOCOOOI\DOiAKN HHHHHHHHHHHHHHHHHHHWWOI^KN4 OJH OJ -4 CO KA O ONO-4 ONCO ON KN H KA ON O L^-ONCOVO KA l^ O ONH KA UA t— O KA UA OA OJ ua OA KA CO KA CO KA O VO KA H O O HVO VOCOONOH KA -4 UA VO CO OA O OJ KA UAVO CO OH KA UA 1^- OA rH KA _4J-J- Lr\Lr\Lf>LrNLrNLr\LrNLrNVOvovovovovo t"- t>-t>-L N —I>-L v -cocO M -CO ONO H CM fA4 LAVO f-CO ONO H Ol KN4 LAND t-CO ONO -4-_4-_4--4 U^UAl/NLr\UAUAlJ^l/NU^ir\VOVOVOVOVOVOVOVOVOVO o- -17- Integrating equations of motions in these regions separately with matching conditions along the boundaries, we could obtain flow conditions behind the intersection region and geometric variables associated with the regions. With the lack of success of local theories and the present incomplete semi- global theory , we proceed to solve the shock intersection problem numerically from a global point of view. This forms the main text of Chapters III and IV. 2.3' Global Theory with True Viscosity As we have seen, all local three-shock theories as well as the proposed semi-global theory seem to fail in one way or the other. The question naturally arises whether or not there is a global theory with true viscosity which would work for all incident shocks. Since the effects of wall and the flow behind the intersection are believed to be important, any global theory- must consider them in the boundary conditions of the whole problem. Analytical solution of the viscous equations of motion with boundary conditions appears out of the question. A complete numerical solution by means of high speed digital computers seems to be in order. Numerically, the appearance of discontinuities makes it difficult to integrate the basic equations throughout the entire flow region under consideration. Von Neumann and Richtmeyer [19] proposed a method for numerical calculation of hydrodynamic shocks. They introduced the so-called "artificial viscosity" to smear out the shock regions. Thus the shock is not considered as an interior boundary, but as a part of the solution,, The idea is to replace the true viscous terms in the differential equations by simpler dissipative expressions. Later Harlow [8] proposed the particle-in-cell (PIC) method for numerical solution of fluid dynamics problems. Rich [11] recently used explicit Landshoff type viscous pressure as well as implicit dissipation terms -18- arising from the "repartitioning" of the energy and momentum in cells for numerical calculations of detached shocks with success. All these techniques were mainly used for relatively simple problems nevertheless. For shock reflection problems, they are still too time-consuming with the existing computing equipment. In 196^-, Bur stein [3] numerically integrated the time -dependent inviscid equations of motion by the Lax-Wendroff difference method. He obtained results for asymptotically stationary regular and Mach reflections in air. The Mach reflection calculations for very strong shocks agree with experimental photographic data obtained from wind-tunnel tests. This is to be expected as the local three-shock theory has already proven satisfactory for strong shocks. However, the results for Mach reflection of weaker incident shocks were quite disappointing as compared to Smith's experiments. It seems, after reviewing previous attempts at the problem, that the true viscous terms play a vital role for the solution of weak shock reflections. Numerical solution of the viscous equations can give a global view of the whole problem. Due to the slowness of the asymptotic development of the shock intersection region as well as the complexity of the true viscous terms in the equations, a more efficient numerical method than the existing ones is needed. We will, in the following chapters, present such a global approach with its numerical solution to show that the Navier-Stokes equations are an adequate description of the flow under these conditions. CHAPTER III DIFFERENTIAL EQUATIONS SUITABLE FOR NUMERICAL SOLUTION 5»1« Time -Dependent Equations Versus Stationary Equations The class of plane shock intersection phenomena under consideration is stationary. That is to say, the value of a flow variable at a given point in a fixed coordinate system does not change in time. It is then natural to use stationary viscous equations of motion for numerical integration. However, in applying global theory to shock reflection problems, there are numerical difficulties. A different approach is needed. In this work, we shall integrate time -dependent viscous equations in Lagrangian form with initial and boundary conditions to obtain asymptotically stationary solutions. It is known that stationary fluid flow equations are nonlinear even without true viscosity terms. There are no convenient numerical methods to handle the nonlinearities that can be guaranteed to converge. The situation becomes worse when viscosity is introduced. If one considers the time- dependent equations, the time derivatives of the flow variables appear linearly. One can therefore difference the equations to iterate in time. The iterative process can be continued until values of all flow variables between consecutive time cycles differ by a negligible amount. Since the time derivatives are then small, the flow is the solution of a set of equa- tions which are disturbed from the original stationary equations by a small amount. In the limit, the solution thus obtained is stationary. On the other hand, if the problem is initially ill-posed, one may not get convergence This fact makes it plausible to use time -dependent equations of motion. There are in general two basic ways of writing time-dependent equations; they are called Lagrangian and Eulerian. In the Lagrangian form, -19- -20- the coordinate system is moving with the fluid. The main advantage of Lagrangian techniques has to do with the characteristics and the Courant stability condition. The equations of motion are also simpler in form. The accuracy of Lagrangian methods is, however, doubtful when applied to a system after distortions of the fluid have occurred. In the Eulerian view point, the coordinate system is fixed relative to the observer and the fluid moves through the mesh of coordinates. Eulerian methods have the advantage of applicability to problems with arbitrary distortions or slippage of the fluid. They also suffer from several disadvantages such as the false diffusion arising from forcing the meshs to be homogeneous. It is then desirable to combine the features of the Lagrangian and Eulerian methods and thereby eliminate some of the disadvantages. The technique we shall adopt is to first write the time -dependent equations of motion in Lagrangian form. An iterative numerical procedure is then derived from these equations . The corresponding difference equations are iterated in time to obtain asymptotically stationary solution. However, we are interested in the stationary solution from Eulerian point of view. Any solution of the Lagrangian equations of motion alone will not be sufficient. This problem can be solved by "repartitioning" fluid cells or "interpolating" fluid mesh points depending on the basic numerical subdivision used in the process. The basic idea is to re-evaluate the flow variables in the fixed Eulerian mesh after each Lagrangian time cycle. The transport effect of flow variables can thus be obtained. A more detailed discussion of the "repartitioning" and "interpolating" processes will be given in the next section. 5.2. The Cell Method and the Point Method In most numerical techniques for solving fluid dynamics problems, the fluid region is subdivided into small cells or mesh points . Finite ■21- difference approximations to partial differential equations discribing the flow are derived and applied to each cell or mesh point. They will be referred to as the cell method and the point method respectively. In 1957; Harlow and his colleagues [8] developed the particle-in-cell method, abbreviated PIC method, for the numerical solution of problems involving the dynamics of compressible fluids. It combines the features of the Lagrangian and the Eulerian methods and thereby eliminates some of the disadvantages. There are two computing meshes; one is Eulerian, the other Lagrangian. The domain through which the fluid is to move is divided into a finite number of cells o This is the Eulerian mesh. In addition, the fluid itself is represented by a finite number of particles which move through the Eulerian mesh, repre- senting the motion of fluid. This is the Lagrangian mesh. Associated with the mesh points of each system are certain variables whose history the calcu- lation develops. Thus for each Eulerian cell there is kept the velocity, the internal energy, and the total mass of each kind of material. For the Lagrangian mesh of particles, individual masses and positions are kept. The differential equations of motion, with transport terms neglected, are written in finite difference form relative to the system of cells. The transport effect is obtained later by allowing the particles to move through the Eulerian mesh according to the velocities in the neighboring cells. Later, Rich [11] modified the PIC method for the case of a continuous fluid. He applied the integral form of the conservation laws to the whole volume of each cell. As in the PIC method, the calculation of the configura- tion of the fluid at time t , from that at time t may be thought of as n+1 n proceeding in three phases. In the first, the fluid element in each cell is assumed fixed and no flow across the cell boundaries is allowed . By applying conservation equations, intermediate values of the velocity and specific -22- energy in each cell are obtained. The second phase consists of the calculation of the mass flows across the cell boundaries. Mass crossing a boundary is assumed to carry with it the intermediate velocity and specific energy of the cell from which it came. In the third phase , the new cell densities are computed; and by application of conservation of energy and momentum, the new velocities and specific energies are subsequently obtained. Rich's results on the supersonic flow of air past a variety of cylindrically symmetric objects agree quite well with experiments. In 1964, Burstein [3] differenced the inviscid conservation equations at each Eulerian mesl joint by the Lax-Wendroff method. It is based on the Taylor expansion of the vector function w(x,y,t -1- At) so as to include the 1 2 second order term — At w, . Since the difference scheme is made conditionally stable, with the addition of this term, a linear analysis allows the computa- tion of the linear stability limit. The time increment At and space increments 2 2 l/2 Ax and Ay are related to the flow Mach number M = (u + v ) ' /c (c being the local sound speed) by At = At < (M^l)- 1 ( j Ax Ay - gl/2 c ^• J - LJ This numerical scheme was tested for various values of incident shock strengths For regular reflection, the results agree quite well with exact solutions. The Mach reflection calculations for weak incident shocks did not compare well with Smith's experiments. Now, we have seen that the PIC method treats each cell as composed T many individual particles, each of them is free to move from one cell to the other. Rich considers each cell as a unit with continuously distributed it. These will be referred to as the "cell method." The equations of -23- motion are integrated with respect to the cell first, and then the integrals are approximated by suitable quadrature formulas . As in Burstein's work, we shall use the "point method" for numerical integration. The flow region will be divided by a mesh. The Lagragian equations of motion will be differenced at each mesh point so that one time step serves to displace the Lagrangian mesh points relative to the Eulerian mesh points. The physical positions of the displaced mesh points relative to the Eulerian mesh points can be calcu- lated from the velocity components « Then the new flow variables on the fixed Eulerian mesh can be evaluated from the intermediate values of the flow variables computed in Lagrangian time cycle by interpolation. 3«3° Differe ntial Equa tions in Lagran gian Form The theory of fluid flow that we will use is based on the Newtonian mechanics of a small body. For continuously distributed fluid particles in the flow, the essential part of Newton's principles can be described by the conservation laws of mass, momentum and energy of a fluid element [17] • For simplicity, we assume the fluid is not heat conductive. From the result of one-dimensional flow [16], principle features of the problem can be found under the assumption of no heat conduction. For the viscosity terms, we assume that the stress tensor is a linear function of the rate of strain tensor. The coefficient of viscosity is taken to be constant and the bulk viscosity coefficient is assumed to be zero= These assumptions are particu- larly justifiable when weak shock intersection problems are considered. Under these assumptions, we can write the time-dependent, two-dimensional, viscous equations in Lagrangian form [10] as follows. (l) Continuity Equation p + pu. . = , i = 1,2 (3-2.1) 1,1 ■2k- (2) Momentum Equations pu. = u(u . .+u . . ) - — Liu. . o . . -p6 . . , j = 1,2 (30-2) where u. is the coefficient of viscosity and 5. . is the Kronecker delta, (3) Energy Equation p(e + g u 1 u i ) LLU.(u. .+U. . ) - — liU, ,U.-pu. J i,J J,i 3 k,k 1 ^ ij ■p(e + - u.u. )u. . (3.3.3) where e is the specific internal energy In the above, we have used a dot over a function to denote, not partial dif- ferentiation with respect to time t at constant (x ,x ), but rather differentiation for a given particle whose positions changes according to its own velocity vector. This is known as the Euler's rule of differentiation; 5t + u 1 bx ± + U 2 bx 2 + u. x~ 1 ox. 1 at i = 1,2 . (3O.10 We have also used the summation convention in Eqs . (3.3.1) - (3-3.10 • This convention will be used throughout this chapter. Furthermore, the notation of a comma followed by a subscript i means partial derivatives of functions with respect to x.. The Eqs. (3.3.1) - (3.3.3), which describe Newton's principles of motion for a viscous fluid, will be referred to as time-dependent viscous equations in Lagrangian form. The motion of the fluid thus described has the property of following the position of each single fluid particle for all values of t. However, there are only four scalar equations in five unknowns, namely, p, U,, u , p and e. We need one more equation to completely specify the -25- problem. There exists no general physical principle which would supply a fifth equation to hold in all cases. For perfect gas, we can use p - (7 - l)pe (30.5) where 7 is the ratio of specific heats. 3° ^° Equations in Streamline Coordinate System The solution of Eqs„ (3-3«l) - (3-3-5) in general describes a family of world lines in (x , x , t)-space, each line corresponding to one particle. The projection onto the (x, ,x )-space of the respective world line is a particle path. If, in particular, the motion is stationary, i.e., at each point of the fluid the same thing happens at all times, there can be only one family of streamlines in the fluid. The family of streamlines thus carries most of the flow information along with it in this case. It is natural to use the streamlines as one of the fundamental coordinates in the computation. For stationary and psuedo-stationary flows, Taub [16] introduced a pair of independent variables, s and y. For fixed time t, the curves s = constant are the loci of particles which have crossed a given curve, say y = 0, at times earlier than t. If y represents the arc length along the curve s = constant, measured from y = 0, then one may transform the functions f = f(x ,x ) into functions of the form f = f(s,y) by defining a x . ~ = \. > i = 1*2 , (3.^.1) and s 1 dx. u. oy "' v i = 1,2 (3.^.2) where X. are the components of the tangent vector to y = constant; -26- u. are the components of the velocity vector; and o v = u.u. . (3-^-3) 11 This defines the streamline coordinates . In addition to the natural reason of using streamline coordinates, there are numerical advantages . In the point method that we shall use for numerical integration, a fluid particle can move into any one of four quadrants if a rectangular coordinate system is used. The "interpolation" process for recovering transport terms in the equations of motion has to handle all four possibilities depending on the various combinations of signs and magnitudes of the velocity components of the particle. However, in streamline coordinate systems, a particle can only move along a streamline. Since the streamlines are one of our basic coordinates, the "interpolation" process only need consider two possibilities. The corresponding "interpolation" formulas will be simpler in this case. A similar situation exists in the cell method- The amount of fluid moved out of a given cell may enter any one of eight neighbor- ing cells in rectangular coordinates. Using streamline coordinates, the outgoing fluid from a cell can only join one of two neighboring cells along the same stream tube. The saving in computation is even more apparent in this case. The second advantage is the simplification of boundary conditions in the point method. In streamline coordinates, the streamlines are straight lines. Flows past inclined boundaries are forced to follow them. The stream- lines are parallel to the boundaries. Values of flow variables at a boundary ^amline can be obtained easily. In rectangular coordinates, the appearance of irregular boundary curves may require complicated interpolation formulas ■27- to compute boundary conditions. Similarly, in the cell method, if a square mesh is set up in a two-dimensional flow region, non-square cells may appear due to the presence of inclined boundaries. The computation becomes more complicated there „ Again the use of streamline coordinates can save time quite appreciably. For time-dependent problems, the particle paths vary in time. They can not be used to define coordinates. The curves s = constant defined earlier are actually flow lines which now depend on time. They shift positions from time t (solid lines) to t + At (dotted lines) as shown in Figure 3.1. The correct interpolation process for flow variables at (£',k') y=y o at t+At, P s=s at t o FIGURE 3.1 MOVEMENTS OF FLOW LINES FROM t TO t + At must average the intermediate flow variables at (£-l,k), U,k), (i-l,k-l) and (i,k-l), referred to as the U-point interpolation . However, for asymptotically -28- stationary flow, the time dependence of flow lines approaches zero. In the limit, (i',k ! ) coincides with (J,k). The corresponding interpolation process at (-T,k' ) will only involve the points (i-l,k) and (£,k), referred to as the 2-point interpolation. Since we are dealing with asymptotically stationary problems, we propose to ignore the movements of flow lines in time, i.e., use the 2-point interpolation in our numerical procedure. This simplification can "be justified by the following arguments. First of all, the asymptotic solution of our time-dependent difference scheme, if one exists, is a solution of stationary difference scheme and conversely. This is because the flow lines are stationary (independent of time) in the limit and then the schemes are identical. The stationary flow lines of the 2-point interpolation are either fixed or moving constantly in one direction. The latter is not possible because of the finite flow region. The converse is true for the simple reason that stationary solutions of 4-point interpolation are solu- t ions of 2-point interpolation. Secondly, the solution of our difference scheme is stable for small perturbations which implies that it locally converges to the solution of stationary scheme with respect to time. We shall prove the stability of the difference scheme with 2-point interpolation process in Chapter IV. In view of the limiting process discussed above, it is possible for us to use the coordinates (s,y) defined by the transformations (3.^.1) - (3 A. 2) in the proposed numerical method. We shall refer to them as the streamline coordinates, although streamlines do not exist in time -dependent problems in the strictest sense. Now, for transforming equations of motion into streamline coordinates, we need the identity transformation t = t (3.^*0 -29- to take care of the time -dependent terms. Furthermore, it is convenient to introduce two new variables X and OL, associated with X. by the equations [l6]: X 2 = X ± X ± , (3.4.5) and tan a = r— . (3-^.6) The angle OL is of course the inclination of the curves y = constant in the (x ,x )-plane. From [16], the following differential identities are also reproduced: u i i = if" chT ' ^ n not summed ) » (3.^-7) > n y df v . df U k ,_ , ov f • = 3T t"7~ £ m \ - Y~ £ M 77" > (3-4.8) ,J oy U n jk k ds jk U n cp = e..u. • =;f (I 1 " ^~) > (3.^.9) u.f .=y| , (3.i+.10) and X.f . =^f , (3.4.11) where U. = u.\. , (3-4.12) til is the velocity tangential to the curves y = constant, scaled by \; U =u.€.A. , (3A.13) n 1 ij J is the velocity normal to the curves y = constant, scaled by X, Snd £ 11 = € 22 = ^ ' £ 12 = "*21 = 1 ' a ^ COS p =r— ds > da Sy = sin P d£ a. as J -50- If the angle (3 is defined as the angle between the normal to the curve y = constant and the tangent to the curve s = constant, then U = \v sin P, (3.4. l4) and U = A.v cos p. (3.4.15) The angle P is taken to he positive when the direction from the normal to y = constant to the tangent to s = constant is counter-clockwise. It can also be shown that (3.^.16) (3.4.17) = v cos (a + p - — ) = v sin £ (3.4.18) and u p = v sin(a + p - — ) = -v cos£ (3.4.19) where % = a + P . (3.4.20) Thus the angle (a + P - — ) gives the inclination of the curves s = constant in the (x ,x ) -plane. Due to the presence of viscosity terms, we shall need: (u i,Aj ■ ( \,k>,i + Ki\,fu\i - K*\i + ^u\t (3 - u - 21) Using the Eqs. (3.4.1) - (3.4.21), we can transform the Eqs. (3-3.1) ■3) into the new coordinate system as follows: (l) Continuity Equation dU 11+^^ = , (n not summed) (3-4.22) n u -51- or p = -pi , (5.^.25) where du v n U dy ' n " (5-^.2^) (2) Momentum Equations "h^KiK^hK^A'^A i = 1,2 . (5^.25) Multiplying Eqs. (3.^.25) "by u. and summing, we have pu . u . = u.u . (cpe . . ) . + — u.u . (u n . ) . - u . p . , K 1 1 ^ x xljl 3 x k,k ,1 1 ,1 where we have used Eq„ (3.4-21). Using Eqs. (3«^<-7) _ (3-^-H) an d simplifying, the above equation becomes 2n_ ,vU. v 1 dU Or, Pw = u.(— ^ - -77- k^-J + — [IV ^— (77- xT~) - V AT ^ V U ds U dy 3 dy U dy dy n n J ' J n ^ v ,dcp t dou ^ dt di U v ds v n 5 » *y (5.^26) Multiplying Eqs. (3 A. 25) by \. and summing, we have pX.u. = ^(cpc.p^ + - ^(u^ . - ^. . Or, 4 d v U t dcp 2 dcp N p\(|v cos p + v sin p) = n ~ (- g - \ ^) + - u ^ ( n dU . y n\ dp 3 ^ ds V U oy~ ; ' ds ' n Finally, using Eqs. (3.K°lk) - (3.^.20) and Eq. (3.J+.26), we can write the above equation in the form h . r dt jt d\J/\ Jl M .. ^£ _i ^P • P V = 5 ^ <£ - T d7 } " ^ v dy ds v dy (3.^.27) -32- (3) Energy Equation >(e + 5 Y l - u(u. .U. . + U. .U. . - — U . .LL ) u .(u. . ) . + u . (u . . ) J 1,1 , 3 J J,i ,ij - — uu . (u. , ) . 3 i k,k ,i - pu. 1,1 1,1 u.p . - p(e + - v )u. . ,i 2 j, j Or, (pe) + - (pv ) k .2 2 , y^6 dv d£ OVs + H v xr + u.(\|r . + (cpe.J J - - uv 3- p^ _ v ^ _ p( e + - v )\|r . Or, (pe)' = 4 .2 2 ■ ,d£ dv d£ dvx If ot v ,dcp u, .3 v oy ; _ - v J* - p\|r - (- pv ) -p(e-+ - v )t . Substituting Eqs . (3.^.23) and (3.^.26) into the above equation and simplifying, we obtain pe = u k- . 2 2 , ,d£ dv as dvx - pV (3.^.28) Let us define and U V (3.^.29) (3.^.30) We then have from Eq. (3.I+.2I+) •33- , 1 a /, Q \ v hr\ av •ty _ _ __ t\ v C os PJ = — -r- 1 + t— t\ oy t) ay ay where we have used Eqs . (3 A. 14) - (3-4.17) and Eq. (3.4.20) (3.4.31) Similarly , we can write cp in terms of r\ and £ as follows 1 /dv a£ «. cVv\ The differential equations for t] and £ are simply, ay " as " *• dy ' (3.4.32) (3.4.33) and oy " n oy (3.4.34) Substituting Eqs „ (3.4.33) - (3*4.34) into Eqs. (3.4-31) - (3.4-32) respectively, we get ♦ -?&-■« !>♦!' 1 /OV u, GV-, * = n ( ^ " s ¥ } as Finally, the equations of motion in terms of r\ and £ are: p = -pt , pv = ^ (- acp \ 4 ,, M ap oy 3 ay oy pvl - 3 ti ^as Q ay ; *- ay n ay 1 ap T) 3s~ ' and = n(7-l) 4 .2 2 . ,ai av as av^ 3 * + * + Uv ( ¥ as" " ai ¥ } 7pt , where we have used Eq. (3.3=5) to obtain Eq. (3»4.4o) (3.4.35) (3.4.36) (3.4-37) (3.4.38) (3.4.39) (3.4.40) Eqs. (3.4.33) - (3.4.40) form a set of eight first order partial differential equations in eight unknowns p, v, £, p, t\, t>, cp and ^o They describe the time -dependent, two-dimensional viscous fluid flow phenomena in streamline coordinates when appropriate boundary conditions are specified. CHAPTER IV NUMERICAL SOLUTION OF VISCOUS SHOCK REFLECTION PROBLEMS 4.1. Stabi lity of One -Dimensional Difference Equations In this section, we are only concerned with the stability of one- dimensional equations , We thus have U n = v, U t = 0. (4.1.1) This leads to the conditions that r\ = 1, £ = , (4.1-2) by definition. Substituting Eqs . (4.1.1) - (4.1.2) into Eqs. (3.4.31) - (3°4c38), we get the following one -dimensional viscous equations of motion P = -pf > (4.1.3) - - M - 1 > <^-« 4 2 p = - n(7-l)t - 7pt , (4.1.5) and *■¥• ( ^ 6) For the sake of stability analysis, it is convenient to eliminate \|/ from the above equations by using Eq. (4.1.6). We then have P=-p|, (4.1-7) pv^.g§-|, ft.1.8) ana p = ; . r-ixfr) - ^p^ • (k - 1 ' 9) For evenly spaced, mesh points, we employ central differences to approximate -35- -36- space derivatives. For the time derivatives, we use the one-sided difference approximation. Writing Eqs. (4.1-7) - (4.1-9) i n difference form, we have n+1 n n n ■ p » n Vl ~ Vl p. ^~ > (4.i.ioj At H H 2Ay n+1 n n n n n n n V l - V l 4 Vl ~ 2Y t ~ Vl P l+1 - P l-1 ,.-_.* Pi At" = 5 ^ ^2~ 2Ay" ' ( ^ 1 ' 11) n+1 n n n 2 n n P l ' P l 4 . lN /i+l " V I-1 N „ n v i+rVi ,, » and ______ ,_ M .( 7 _ 1 )( ^ y ) - 7V i ^ ( + .1.12) Assuming that a solution of Eqs. (4.1.10) - (4.1.12) can be expressed in terms of Fourier series, let us examine the conditions under which a typical term from the series can be a solution. Take as the representative terms respectively n _ ikJ n n . __x Pjl = p + Spe r , (4.1.13) n _ iki n ,, _ _ , v v = v + Sve r , (4.1.14) n _ ikl n , , n , _ N and p = p + ope r , (4.1.15) where p, v and p are the true solution of Eqs. (4.1-7) - (4.1-9) • Substituting them into Eqs. (4.1.10) - (4.1.12) and simplifying, we obtain (r-l)Sp = -ap( el - e X )6v - (^ At)6p , (4.1.16) p(r-l)8v + (vAt)Bp = I in(e X -2+e" 1 )&v-a(— |-S — -)op (4.1.17) and (r-l)Bp = | n( 7 -l)o * (£_^|__) &v - 7pa (£_^|-_) 5 v - 7 ^ At6p , (4.1.18) ire == p^ , (4.1.19) and T = At Ay 2 (4.1.20) Assuming \x, At and Ay are small and of the same order of magnitude, we can hold a and [IT fixed and let At and [ia -» in Eqs. (4.1.16) - (4.1.18). In the limit we get the following system of linear algebraic equations in 5p, &v and op (r-l)bp + (lap sin k ; bv = , 8 p(r-l) + ~ u: 1 ■ cos k) 5 v -»- i a s 1 n k I op - C , and L7po sin k)8v +- (r-l'Sp = (4.1.21) (^. 1.22) (4 1.23) A necessary and sufficient condition that Eqs (4.1 2L - (4.1.23) have a non-trivial solution for 5p., ov and op is that r-1 iap sin k C I r-1,) f — |J.t{1 - cos k' i.a sin k i?pa sin k (r-1 = (4:1. 2*0 Or 1 1 > 2 8 J r-1) + - q - 5 P 5 2 2 L - ccs k)(r-lj ■•- c a sin k (4.1.25! where P I U.1,26] In obtaining Eq< (4.1.25), we have discarded tne asymptotically stable root r = I. The solution of Eq. (4.1 25' can be written in the form r = 1 - nil I g 2 2 2 ' cos k) - A/a (l-cosk) -p sin k , wl e re a. and p. 4 1 3 " e ca (4.1.2?) (4.1.28) (4.1.29) -38- In order that solution (4.1.13) - (4.1.15) does not grow with time, it is necessary that |r| < 1. This condition imposes restrictions upon the parameters At and Ay. We will discuss the stability conditions in two cases: CASE I . a 2 (l-cosk) 2 < p 2 sin 2 k Then r is complex and |r| 2 = 1 - 2a(l-cosk) + p 2 sin 2 k = 1 - (l-cosk)[2a- P 2 (l + cos k)] . (4.1-30) Thus, in order that |r| < 1, it is required that 2a - p 2 (l + cos k) > . Using Eqs. (4.1.19), (4.1.28) and (4.1.29), we finally get 4 - 7P where we have replaced (l+ cos k) by its maximum value 2. Case II. a 2 (l - cos k) 2 > p 2 sin 2 k At < -* (4.1.31) In this case, r is real and -vj- is proportional to sin k so that dr dk |r| attains its maximum value at k = 0, it. For k = 0, r = 1. No restriction is imposed on At and Ay. For k = jt, we have r = 1 - 4a . (4.1.32) For stability, we require a<| . (^.1.33) This gives At < td^dl (.k.i.ik) 5» -39- For sufficiently small Ay, condition (4.1„3^) dominates. We will use it to determine the time step for numerical integration, k . 2 , Nume r i cal Solution of _One -Dimen sio na l Viscous Sh o ck Problem Lmagine a long pipe with one end closed, A uniform fluid flow is "being supplied continuously from the other end. When the fluid particles meet the closed end, disturbances are created. A shock front is generated which moves away from the closed end as shown in Figure k (4.2.3) where n <4 n n v - v 1+1 1 Ay (4.2.4) and n n v - v 1 1-1 Ay (4.2.5) n 1 / ,n , n 1 / ,n n N *! -2 ( * 1 + * 1 } (4.2.6) !+— !• 2 2 The reason for using half mesh values of \|f is to keep uniform error in difference equations. Since p , v , and p ' are the intermediate flow variables of the 1th mesh point at ( n+1 ) -time cycle, we need to combine them with the of appropriate neighboring mesh point for the "interpolation" process. -41- The weighted linear interpolation is performed according to the following formulas : CASE I . v^ > n+1 .— n+1 — n+1 ,, x Pi = Ap i-i + Bp i ' ^- 2 -T) n+1 A — n+1 -n— n+1 ,, „ ON Y * = A Vi + Bv i - ^- 2 - 8 > and p^ +1 . Ai£ + Bp^ +1 , (4. 2. 9) v At where A = — — — , ( + .2.10) n n Ay - v^_ x At + v^ At n Ay - v At and B = ■ =^±- . (4.2.11) Ay - v^_ x At + v'J At For allowing mesh points to move in opposite direction, we need CASE II . v^ < n+1 — n+1 — n+1 /, _ __\ p i =A ^i +Bp i+1 ' ( ^ 2 - 12) n+1 . —n+1 t,— n+1 / , _..,-, \ and p^ +1 = A^ +1 + B?£ , (4.2.14) Ay + v" At where A = — , (4.2.15) Ay + v^ +1 At - v^ At and B - - — . (4.2.16) Ay + v° w At - v« At The reasons for choosing above formulas are quite clear geometrically as shown in Figure 4.2. 42- ^V-VJL,At JC?M^ v£>0 -v^ n +1 — _ - ~ A1 2 2 Ah prj J B T i Vb 2 + C .3 5 (4.3-18) 2 2> where B = (l - cos k)(£ + T) ) + (l - cos i - £ sin k sin £) , (4.3*19) 2 2 2 2 and C = T] sin k sin I - 4tj (cos k - l)(cos I - l) (4.3.20) For stability, we require At < (Ah) 2 pi] 2 f b ; i ^B 2 + c (4.3.21) But ,^2„2k.2i/ 2k 2 1 ,^ , „ C = lor) sm — sin — (cos — cos — - 1) < (4o3-22) Therefore, we can simplify (4»3°2l) and get At < WW < ^ |i Max(B) f Bf \ Vb 2 + C (4.3.23) The extrema of B occur at k = 0>* and i = 0,it. We have the following cases : -50- (1) k = , I = , (2) k = , i = « , (3) k = Tt , I = , (k) k = jt , £ = a , Condition (4.3-23) now becomes B = , B = 2 , (minimum) B = 2(£ 2 + T! 2 ) , 2 2 B = 2(£ +q + l) (maximum) At < / (Ah) 2 pn 2 - 15 772 2 _v (h-3.2k) Since we neglected the term C in deriving (4.3*23) > the stability condition (4-3.24) may be a little too restrictive. However, it does give us a reasonable upper bound for At in terms of Ah. Comparing to the corresponding one- dimensional stability condition (4.1-34), we can get £ = 0, T] = 1 and £ s 0. 2 2 The maximum of B now takes on the value 2(£ + t] ) =2. (4.3-23) then becomes At < p(Ah)' - 16 (4-3.25) I" This gives a smaller upper bound than (4.1.34) because of the simplification made earlier. For the proof of stability of interpolation process in this case, it is identical to that given in Section 4.2 since the interpolation formulas are the same. We therefore have the stability of the two-dimensional difference scheme with interpolation. 4.4. Numerical Solution of Regular and Mach Reflection^ In this section, we shall describe a numerical procedure of solving shock reflection problems. To simplify the computation, we consider a two- dimensional rectangular box (Fig. 4.4) where uniform fluid flow comes in from -51- the left and exits at the right. The top and bottom boundaries are assumed rigid and smooth. The boundary layer effects will be ignored completely. The upper wall can be thought of as the dividing streamline of two identical shocks incident on each other as explained in Section 2o2. UPPER WALL ///////////////////////////////////////////// ////// LOWER WALL FIGURE k.k TWO-DIMENSIONAL SHOCK REFLECTION CONFIGURATION One way of generating a steady incident shock in the flow is to place a straight half wedge with angle go on the bottom boundary of the box. Flow will then be disturbed starting at the tip of the wedge. If the incoming flow is supersonic and tne wedge angle is small, a fairly straight attached shock will appear as shown in Figure h.K. Using this simple arrangement, desired incident shock angle a and its strength ^ can easily be achieved by fixing the half wedge angle u and the incoming flow Mach number M T . -52- In the present calculations , a square mesh is set up in the rectangular box in streamline coordinates. The finite difference equations applied at each mesh point, derivable from Eqs . (3 •^•37) ~ (3-^-^0), can be written in the following form n+1 n n n (k.h.l) 9 9 n+1 n At v„ = v. + l,m I ,m n l,nn ,m-g n, Ah - 5 9 , i _cp , i n * + 2' m £ -2- m !,m " Ah £,m £,m + M> ,n n . 1 . 1 n n i+— ,m ^-7T>m P. , - P„ -, 2' 2* l+l,m l-l,m Ah 2Ah (4. +.2) , n+1 , n I. = £. + At £,m i.m n n P - v I- 1 , Tl n ,n + 1" + 1 k " Ah £,m £ ,m l ,m 9 ! -cp ,n ,n * 1 " + 1 ^£,m Ah ' " ^ ~ Ah .n n n n n ^l,m P l+l,m " P l-l,m 1 P l,m+1 " p l,m-l n 2Ah " n 2Ah 'i,m 'i,m (^•3) and n+1 ~* l,m C + * fa-v [j <*? fB )2 + KJ »" * n " ii liv 11 / *+l,m l-l,m V l,m+l~ V l,m-l V i,m 2Ah 2Ah »n n n n /,m+l~ S l,m-l v f+l,m" V *-l,m 2Ah 2Ah n i n P £,m £,mj . (+.+.+) =53- In order to keep uniform error in the above difference equations, we use half mesh values of cp and \|/ in the corresponding difference quotients „ These half mesh values can be computed from Eqs . (3. 4. 35) - (3 <^«36) as follows t n i+-,m 2> m n i l+-,m £+^,m -y . i n n v -v £+l,m I ,m Ah (^•5) and f £ ,m+ 1 £ ,m+- £,nn- £,m+2 m 1 / n (v V 1 + 1, m £ ,m (+.^•7) V £,m+; 1 / n n v - v +■ v ) . 2 £,m+l |, m ' ' (h.k.8] «/ i :+2» m 5 1 5 1 £+-ym+l £+— ,m-l ~2Ah~ n n n n l+Vm+l^^ _ i .igLt- 1 £+l,m r l _" J^m^l , i4Ah (k.k.9) (5 s> 1 £,m+- £^m+l_ b £,m Ah (lj-.lMO) (U n y . i ,n ? n '£+l,m _ £,m "Ah (lf.i^.ll) •5k- and , vn ' 2 ' 2 b l+l,m+l l+l,m £-l,m+l "i-l,m (4.J+.12) are already available from previous calculati Since \k and ty are already available from previous calculations we have ,n 1 / ,n ,n ,n ,n 1 + ^ 1 + ^ 1 + ^ --,m i-g.m 2,n>f- £,_ ,n 1 / ,n ,n ,n ,n i+^ m *-3> m £ >™4 ^ m ~^ ) • (^-^.13) Similar expressions for half mesh values of cp can also be derived. As for the half mesh and mesh point values of "Q and t,, we numerically integrate Eqs. (3.^.33) - (3.^.3^) along s = constant and get two simultaneous equations in t] and c. ~ n ^h /» N n *- n ^, n ,, /, \n Ah ^n /„ N n i+p,m J ' ^+p^m ^-p^m £ > m ^-pjm i,m ' and Ah /. \n n «,n uX\ Ah n , „ >n ■T Vl,m\ 1 + K 1 = ^ 1 + "\ 1 (? y } Mn ' i+7j,m U-,m i--,m i--,m ' f4.lj-.15) With the values of T) , and £ given by previous calculations, the solution of Eqs. (4.^.1^) - (k.k.Vi) can be written , g — 2 y «,"> (U.U.:6) ■55- and .. /„ \n n L r Ah ,„ >,n -.21 «,n Ah ,, N n ,„ N n 5 1 i--,m ^ ' ^ i«2,m '4- i + [ f < ty) ; 3 2 (W.ifl i /* N n i+l,m i-l,m /. , _, n \ where (£ ). = > OAV *— , (4.4.18) w y i ,m 2Ah v ' and ( } n %m+l ^m-1 (4 4.1 9 ) We also have and n n Slightly more complicated expressions for rj and £ can "be derived I ,m+2 * ,nH"2 by using the same method. The set of equations (4.4.1) - (4.4.21) completely describes the numerical iteration procedure for two-dimensional viscous fluid flow problems. In order to recover the transport terms of each flow variable for asymptotical- ly steady flow, we need perform the "interpolation" process. In our coordinates, this process is very simple as expected. Since the mesh points are all moving unidirectionally from left to right along curves s = constant in the shock reflection problem, it is only necessary to use formulas like - 5 6- (J+.2.7) - (^.2.1l) to interpolate between them. The boundary conditions for the problem shown in Figure h.k are simply as follows, thanks to our coordinate system, (1) Left side (uniform incoming flow) n n n n n n n n P 2,m " P l,m ' V 2,m = V l,m ' ^2,m '~~~ *l,m ' P 2,m == P l,m ' VJ = , cp^ =0 for m = 1,2,..., M ; (2) Right side (continuous outgoing flow) n n n n „n „n Pt = Pt -, ) V t = V t -, > It = It n i L,m L-l,m L,m L-l,m L,m -L-l,m n n * n o p = p for m = 1,2,. . . ,M ; L,m L-l,m (3) Top side (rigid smooth wall) n n n n n n P i,M == P i,M-l ' Y £,M = Y i,M-l ' *l,M '~'~ * *1,M-1 ' and n n P* m = P. m_i for * = 1,2,..., L ; (*0 Bottom side (rigid smooth wall with straight wedge) n n n n n ,n P i,l = P H,2 > V i,l = V *,2 ' *f,l = * ~ *|,2 > p £,i = p ?,2 for i = 1 ^ 2 ^--- > L ; |--=Jt+2u-| for mesh points where the wedge lies. : appears in above conditions because we want to maintain zero Lty component normal to the walls. In actual computation, the starting -57- half wedge angle w should be taken small enough so that the inclined wedge face does not intersect any of the initially horizontal streamlines. It can he increased gradually after every few time cycles to achieve the desired final angle. The initial flow values at all mesh points are uniform as in Section h.2. For simplicity, we take them to he ^° in ° -, n ,° k I = 1,2,. ...L p *,m = 1 '°> V Mn = 1 -°> h,^ = 2> for m = l, 2 ;...;M . The value for initial pressure depends upon the incident shock strength | and angle a specified by the problem. Using oblique shock conditions [10], we have 1 7 - 1* z 1 , \ "i = 7^~ ' : ( u -^> and ? vl. 7 The choice of v. =1.0 implies that the initial speed of sound C T = — . X. A III _!_ X 1— But M > 1 if | < 1 from Eq. (4.4.22). Thus we are guaranteed to have supersonic incoming flow. The desired wedge angle w for generating an incident shock with angle o^. can be calculated using the formula [10], _L 2 2 M_ sin a - 1 oa = tan" 1 [ 2 cot a g ] . (4.4.21+) M (7 + cos 2a ) + 2 Other constants used for all calculations in this work are 7 - 1.4, (j, = 0.01, L = 50, M = 25, Ah = 0.05 and At = 0.02, where At has been chosen to satisfy (4.3.2 ). -58- The main computer program is written in FORTRAN II language with inner iteration loop coded in ILLIAC II machine language. We test for stationary solution of a given shock reflection problem by comparing all flow variables between every 40 time cycles. If the relative differences are found to be negligible, we stop the computation. The average number of time cycles required for stationary solution is about 600 which takes about 15 minutes of machine time. This time includes starting computation from initially uniform flow with an inclined wedge until flow becoming stationary and writing flow data on a magnetic tape as shown in Figure 4.5. Figure 4.6 gives the asymptotic approach of stationary solution. The changes in p and v along streamlines are plotted by a CalComp 67O plotter. Main results of this work are summarized in Figures 4.7 - 4.9. For given incident shock strengths £ = 0.9> 0.8 and 0.7, we plot reflected shock angle cc-. versus given incident shock angle a x obtained by our numerical calculations. We compare the results with those given by Von Neumann's local theory, the semi-global theory of Section 2.3 and Smith's experiments. For £ = 0.9, we also give numerical results for the case of no explicit viscosity effect, i.e., u = 0. Numerical results that we obtain by the global theory are in much better agreement with Smith's data than any other results. For f = 0.9, the results given by the 1-Region and 2-Region models of the semi-global theory agree quite well with experimental data. However, the agreement becomes worse as | decreases, i.e., as the incident shock becomes stronger. The method of obtaining reflected shock angles from our computer is demonstrated in Figure 4.5 for the case | = 0.9, a = 56 and >0.8 . There we draw constant density lines on a density variati on -59- chart of streamlines (upper diagram of Fig. ^,5). Transplanting the same constant density lines onto a streamline chart in physical space in rectangular coordinates (lower diagram of Fig. ^.5)> they become curved. The average of the angles of constant density curves made with the wall measured inside the reflected shock region is the desired reflected shock angle. The same technique can he applied to velocity and pressure variation charts to verify the results. •6o- TIME CYCLE > 600 , INCIDENT ANGLE = 86*0 , 6 8 9 A B C E G H 1.02 1.03 3 9 4 9 5 O 6 J 8 8 9 9 A A B B c c D c E E G G H H 1 1 n 2 « 3 " 4 ii 5 .i 6 H 9 A B C D E G 1 1 "^ 2 9- 3 <»-4 *■ 5 9- 6 ; 8 9 A B C D E G H 1 I U 2 'l 3 U 4 l-l 5 '') 6 1 8 9 A B C D E G H 3/4 / 5 / 6 7 8 9 INCIDEN1 - A B C D E G H ^1 / 5 / 6 7 8 9 SHOCK A B C D E G H J7~-®-J_6 7 9 *M*j! A B C D E G H 3rf — s A 8 C D E G H 1 ■3r~~~~* JJMO'O A 8 C D E G H 1 1 2 S^~®-i^3~~ — ' -3Jt}f&&~—- - A B C D E G H 1 1 2 3 4 st^~5~—~--: c D E G H 1 1 2 3 3 4 5 ^K-^. 4 5 6 7 — ^ *9T— - -fijF^ 1 A r A " c ■ Tfrr: nrp G H % ® 1 1 2 3 4 5 6 8 9 A -©— " A 1 1 2 3 3 3 3 4 5 6 4 5 6 4 5 6 4 5 6 7 7 8 8 8 8 9 9 9 A A A B B c c D E —f MICH —6 u — SHOCK g> ^2 1 1 2 B c D V) h F -G — H -#— £ ^ 1 1 2 3 4 5 6 4 5 b_^- 1 B "Sfe- aV .915*— E -<&- -c — H — n — SA 1 1 2 •*^-&— J 6 ~l ♦ - 8- - •""® o ~9j»4-®C ^$L s^w^ -.0— --E-; G__ — G- — H — 1 T 2 Jl_^—-®~~ 6 7 -"9 ""[> T:a4^j^— - — I G H 1 1 2_ "4" l 5 8 _S——-®' A 8 E G H l 1 4 i 5 6 -"H -~~9 9^^_@- r? C E E G G H H 1 1 1) 5~~ 6 7 — -9 — 6 A , b c D E G H r~J_3 4 a 9 A b c D E G H 4 5 6 7 8 9 A / b c D E G H 1 1 f ? 3 4 5 6 7 a 9 REFLECTED c E G H % % 1 1 l '1 2 U 3 U 4 U 5 1 6 7 a 9 SHOCK c D E G H ? / 3 / 4 / 5 /6 7 8 9 A B c E G H 1 1 f- 2 / 3 2 ID 3 i : s ; s 6 . 7 7 8 8 9 9 A A B B c c D E E G G H H 2-3 = 4 - 5 O 6 7 8 9 A B c E G H 7 8 9 A B c D E G H 1 CL 2 Q. 3 CL 4 ft. 5 <*. 6 7 8 9 A 8 c n E G H 2 3 4 5 6 8 9 A B c n E G H 2 3 4 5 6 ' 8 9 A B c n E G H 2 3 4 5 6 7 a 9 A B c n E G H 20 STREAMLINE NUMBER FIGURE MEASUREMENT OF REFLECTED SHOCK ANGLE FROM COMPUTER OUTPUT -6i- £ =0.9 , a T = 59°, 0^=65.8° O 8 u _l o >- o UJ 1.174 Ilk A 5 i DENSITY % 1 1.0 >- 0.95 H O 3 UJ > 091 O O If) o Ul 2 1st STREAMLINE FROM WEDGE I L- 1.167 /fiSHi ^STREAMLINE CLOSEST -* TO WALL >- fm 1- V) Itt STREAMLINE FROM WEDGE 1.088 g a in i|J|l Jp 094 >- -J UJ > 0.89 O 9 I UJ 2 U/'*'^^ 1 w 1 w 1 V- tUB 1 ■< 1 1 J 1.118 1.036 20 30 1.0 10 20 30 40 50 10 CELL NO. IN Y-DIRECTI0N FIGURE k.6 ASYMPTOTIC APPROACH OF SOLUTION OF PLANE SHOCK REFLECTION 40 50 -62- 90 fl Cl- 0.9 / x / X / XX X 85° / / i / • -/ / / X • 80° mm /./ R E F L • / E \ • / t75° \ // E D . ;,/ s -• SEMIt- GLOBAL THEORY 1 x H ( 1- REGION MODEL) c70° _ SEMI -GLOBAL THEORY ■ // K A "■ (2-REGION MODEL) N D D NUMERICAL RESULTS WITH NO // G EXPLICIT VISCOSITY • • E 65° • • • NUMERICAL RESULTS WITH / x« / - TRUE VISCOSITY * <*r° - - X X X SMITH'S EXPERIMENTS xAv 60° ■" THREE- SHOCK THEORY /// A' IWU-oMUUK 1 ntUnl y 55° A «so° J_l 1 1 1 /£\ 1 J 1 1 40< 45* 50° 55° 60° INCIDENT SHOCK ANGLE a} 65< 70° FIGURE U.7 c^ VERSUS a FOR | =0.9 -63- 90°i- Cx* 0.8 X X / // / // / / / • /-•-SEMI-GLOBAL THEORY ( l-REGION MODEL) — SEMI-GLOBAL THEORY (2-REGION MODEL) /"/ • • • NUMERICAL RESULTS // / /- XXX SMITH'S EXPERIMENTS THREE-SHOCK THEORY TWO-SHOCK THEORY 1 1 I I I I I I I I I 50' 55° 60 c INCIDENT SHOCK ANGLE Qj* J 65° FIGURE K.Q QU VERSUS a FOR | =0.8 -64- e x = 0.7 / / / / / / X X / X / / ' / / / / / / / / / / i / / ■ / / / / ,'— •— SEMI-GLOBAL THEORY ( l-REQION MODEL) / / / / / / / / • / SEMI -GLOBAL THEORY / / ( 2-REGION MODEL ) • • • NUMERICAL RESULTS XXX smith's experiments THREE -SHOCK THEORY TWO- SHOCK THEORY i / l/ 1 I I i I 1 I i i J— L I ' I ' ±J 45° 50° 55° 60° INCIDENT SHOCK ANGLE Qj 65* FIGURE 4.9 CL VERSUS a FOR £ =0.7 CHAPTER V CONCLUSION The purpose of this work was to show that the Navier-Stokes equations are an adequate model of fluid flow in the case of the weak shock reflection problem. It is also attempted to give simpler models that would agree with experimental observations and be computational simple so that more complex problems can be tackled. Previous attempts to approximate the Navier-Stokes model in various ways have failed to agree with experimental results for weak Mach reflections, Many of the efforts were discussed in Chapter II. One of the first, Von Neumann's, which uses Rankine-Hugoniot equations locally in the neighbor- hood of the shock intersection point worked for strong shock reflections and regular reflections but not for weak Mach reflections. This is presumably because, in the presence of viscosity, the differential equations are always elliptic and that this small effect becomes more noticeable as the flow Mach number decreases to 1. Therefore with small Mach numbers, the effect of the downstream boundaries must be taken into account. Many authors have made modifications of this local theory to take account of viscosity locally. They also failed to work for weak Mach reflections, presumably for the same reason. In Section 2.2„ , we proposed two new models which take into account the shock thickness (due to viscosity) and the effects of the wall to which the shock is attached. Numerical results were calculated from these models for a range of incident shock angles and three incident shock strengths. -6 5 - ■66- While these were pleasing in that they displayed a behaviour similar to Smith's experimental results, they were not as close to those results as would he desired. Furthermore, the more complex and presumably realistic model gave worse results in this respect than the simpler model. However, this so-called semi-global approach appears to be promising in that an arbitrary large number of regions can be used to represent the shock configurations, but at this time a better understanding of the flow relations that should be used across such regions is needed. Since these models were based on simplifying assumptions made about the Navier-Stokes model, the question of whether or not this is an adequate model naturally rises. Previously workers in this field have numerically integrated the equations of motion for boundary value problems. They either ignored viscosity or replaced it by an artificial viscosity term in order to try to avoid instability that often results near shock regions. Even with such simplifying assumptions, the computational process is quite time-consuming. Burstein has specifically studied the shock reflection problems. His approach neglected viscosity (the difference equations serve to smear out the shock region over several mesh points much as viscosity does to avoid singularities). His results agreed well with experiments and with the local theory for strong shock reflections, but not for weak shock reflections. We therefore proceeded on the assump- .on that the true viscosity terms should be retained. Because weak shocks were to be investigated, it seemed reasonable to assume that the viscosity ; was constant (the temperature does not change much) and also heat conduction terms can be neglected. [Terential equations were numerically integrated in a mixed i.shioned after the PIC method of Harlow. -67- A transformation of the independent variables was first performed in order to greatly reduce the amount of computation involved in this two-stage process. Numerical results were obtained for a range of incident shock angles over a set of three incident shock strengths whose pressure ratios were 0„9> 0-8 and 0„7 respectively. This computation was performed on the ILLIAC II computer in the Department of Computer Science, University of Illinois. Each shock configuration took about 15 minutes in which time about 600 time steps had been made on a 50 x 25 spatial mesh. The inter- mediate and final flow diagrams were plotted by the computer and hand measurements made on the reflected shock angles. They agree well with Smith's experimental results (see Figs. 4.7 - 4.9). Some computations were also made with the viscosity terms neglected. The results did not agree well with experimental data (see Fig. 4.7). This leads us to suspect that true viscous effect is indeed important for weak Mach reflection problems. It also shows that the Navier-Stokes description of fluid flow is valid in this case. There are several obvious ways to improve the efficiency of the present numerical scheme. For shock reflection problems, a stationary incident shock could be maintained in the two-dimensional computational region. This is because the generation of incident shocks by placing a wedge in a uniform flow turned out to be quite time-consuming. As for the iteration process itself, the rate of convergence may be increased by using appropriate relaxation schemes. For example, in treating asymptotically stationary problems, the time is no longer physically meaningful because it is used in iteration to avoid numerical difficulties. We can consider variable optimum time step at each mesh point instead of a uniform time step -68- for the whole mesh. This local time step may serve as a type of relaxation parameter to improve the convergence rate. Also higher order difference approximations may be used to reduce the number of mesh points required in the flow region and hence speed up the computation. BIBLIOGRAPHY 1. Bargmann, V. and Montgomery, Do , "On Nearly Glancing Reflection of Shocks/' O.S.R.D. Report No. 5171, (19^5). 2. Bleakney, W. , "Review of Significant Observations on the Mach Reflection of Shock Waves," Proc of Symposia in Appl. Math., Vol. 5 (195*0, PP- ^1-^7* 3. Burstein, S. Z. , "Numerical Methods in Multidimensional Shocked Flows," AIAA Journal, Vol. 2, No. 12, (19&0, PP° 2111-2117. k. Clutterham, D. R. and Taub, A. H. , "Numerical Results on the Shock Configuration in Mach Reflections," Proc. of Symposia in Appl. Math., Vol. 6, (1956), pp. i+5-58. 5. Crocco, L. , "A Suggestion for the Numerical Solution of the Steady Navier-Stokes Equations," AIAA Paper, No. 65-I, (1965). 6. Gear, C. Wo, "Singular Shock Intersections in Plane Flow," University of Illinois, Digital Computer Lab., Report No. 100, (i960). 7. Guderley, K. G. , "Consideration of the Structure of Mixed Subsonic - Supersonic Flow Patterns," Wright Field Report F-TR-2168-ND, (19^7). 8. Harlow, F. H. , "The Particle-In-Cell Method for Numerical Solution of Problem in Fluid Dynamics," Proc. of Symposia in Appl. Math., Vol. 15, (1965), pp. 269-288. 9. Keller, J. B. and Blank, A., "Diffraction and Reflection of Pulses by Wedges and Corners," Coram. Pure Appl. Math., Vol. k, No. 1, (l95l), pp. 75-9^" 10. Liepmann, H. W. and Roshko, A., "Elements of Gasdynamics," John Wiley and Sons, Inc., New York, (1957). 11. Rich, M. , "A Method for Eulerian Fluid Dynamics," Los Alamos Scientific Lab., Report No. LAMS-2826, (1963). 12. Richtmeyer, R. D. , "Difference Methods for Initial -Value Problems," Inter-Science Publishers Inc., New York, (1957)« 13. Sakurai, A. , "An Attempt to Remedy Discrepancies in the Theory of Three-Shock-Interactions," Unpublished report, Courant Institute of Math. Sciences, NYU, (1961). I^o Smith, L. G. , "Photographic Investigations of the Reflection of Plane Shocks in Air," O.S.RoDo, Report No. 6271, (19^5). -69- -70- 15. Sternberg, J., "Triple-Shock-Wave Intersections," Phys. of Fluids, Vol. 2, No. 2, (1959), pp. 179-206. 16. Taub, A. H. , "Determination of Flows Behind Stationary and Pseudo- Stationary Shocks, " Annals of Math. , Vol. 62, No. 2, (1955), pp. 300-325. 17. Von Mises, R. , "Mathematical Theory of Compressible Fluid Flow," Academic Press Inc., New York, (1958). 18. Von Neumann, J., "Oblique Reflection of Shocks," Explosives Research Report No. 12, U.S. Navy Dept. , (19^3). 19. Von Neumann, J. and Richtmeyer, R. D. , "A Method for the Numerical Calculation of Hydrodynamic Shocks," J. of Appl. Phys., Vol. 21, (1950), pp. 232-237. VITA The author was born on August 29, 193^ in Nanking, China. He attended Chien-Kuo Middle School at Taipei, Taiwan from 19 J +9-1952. He enrolled in the National Taiwan University, Taipei, Taiwan in September, 1952, majoring in Mechanical Engineering. He received his B.S. degree in June, 1956. From 1956-1958, he served in the Chinese Air Force ROTC as a second leiutenant in Taiwan. After honorably discharged from the military service, he entered the graduate college of the University of Illinois. He studied Aeronautical Engineering and worked half time as a research assistant in the same department from 1958 - i960. He received his M.S. degree in February, i960. Beginning in September, i960, he transferred to the Department of Mathematics and taught half time as a teaching assistant. In June, 1961, he went to work full time for the International Business Machines Corporation in New York as a numerical analyst and computer designer. He published a paper, co-authored with T. C. Chen and R. M. Frank, on "Tables of Zeros and Gaussian Weights of Certain Associated Laguerre Polynomials and the Related Generalized Hermite Polynomials" in Mathematics of Computation, October, 196k, for his work performed while at IBM. He returned to the University of Illinois in September, 19^3 with a university fellowship to study applied mathematics. He also taught one-fourth time in the Department of Mathematics in the second semester of the academic year of 1963-196^. He worked for IBM again in the summer of 196^. Since September, 196^ till now, he worked for the Department of Computer Science as a half-time research assistant doing the research reported herein. -71- -nm!