T H E U N I V E R S I T Y OF M I C H I G A N COLLEGE OF ENGINEERING Department of Mechanical Engineering Heat Transfer Laboratory Technical Report No. 3 A FINITE DIFFERENCE METHOD FOR COMPUTING UNSTEADY, INCOMPRESSIBLE, LAMINAR BOUNDARY-LAYER FLOWS Charles Luh-Sun Farn Vedat S. Arpaci John A. Clark o.;,ORA.:Project [05065" under contract with: AERONAUTICAL RESEARCH LABORATORY, OAR AERONAUTICAL SYSTEMS DIVISION AIR FORCE SYSTEMS COMMAND CONTRACT NO. AF 33(657)-8368 WRIGHT-PATTERSON AIR FORCE BASE, OHIO administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July 1965

This report was also a dissertation submitted by the first author in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan, 1965.

ACKNOWLEDGMENTS The authors are indebted to the many individuals who have contributed to this work. In particular, they would like to express their appreciation to Professors Arnold M. Kuethe and Wen Jei Yang, who served as members of the thesis committee. The research reported herein is part of the research program on "Oscillating Boundary Layer Flows" of the Aeronautical Research Laboratories, Office of Aerospace Research, of the United States Air Force, whose support is acknowledged. The assistance of Dr. Max G. Scherberg of this laboratory is appreciated. ii

TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vi NOMENCLATURE viii CHAPTER I INTRODUCTION 1 II. REVIEW OF THE PREVIOUS WORKS 5 A. Similar Solutions 5 B. Semi-Similar Solutions 8 C. Series Solutions 9 D. Successive Approximations 10 E. Integral Method 12 F. Other Methods 13 III. FORMULATION OF THE PROBLEM 16 A. Differential Form 16 B. Finite Difference Form 22 1. System of Space Grid Points 22 2. Approximation of Derivatives by Finite Differences 25 3. Finite Difference Form of Boundary-Layer Equations 27 4. Numerical Computation Procedure 29 5. General Error Analysis 33 6. Stability Analysis and Test 36 7. Convergence Analysis and Test 44 8. Singular Errors 56 9. Local Rounding Errors 61 IV. NUMERICAL EXAMPLES 66 A. Oscillations in Blasius Flow 66 1. Statement of the Problem 66 2. Numerical Solution for High-Frequency Case 68 3. Numerical Solution for IntermediateFrequency Case 72 iii

TABLE OF CONTENTS (Concluded) Page B. Impulsive Start of Wedge Flow 82 1. Numerical Solution of the Transient Velocity Profile 82 2. Numerical Solution of the Transient Temperature Profile 84 V. CONCLUSION 91 APPENDIX I. DERIVATION OF INEQUALITY (3. 66) 93 II. COMPUTER PROGRAMS 95 A. Calculation of Velocity of Oscillation in Blasius Flow 96 B. Calculation of Velocity of Impulsive Start of Wedge Flow 99 C. Calculation of Temperature of Wedge Flow after Sudden Drop of Wall Temperature 102 REFERENCES 104 iv

LIST OF TABLES Table Page 1 Systems of space domain 44 2. Series of space interval AXi 45 3. Distribution of disturbances at Xa 96 4. Damp of disturbances introduced at X = 1.95 x 105 97 5. Damp of disturbances introduced at X = 4.00 x 105 97

LIST OF FIGURES Figure Page 1. Curvilinear orthogonal coordinates. 16 2. System of space grid points. 23 3. Effects of the size of AT on the instability errors. 43 4. Effects of the size of AXi on the convergence of solutions of FDAEs (3.35) and (3.37) for Blasius flow. 57 5. Effects of singular errors on numerical solutions of the local velocity profiles in Blasius flow. 58 6. Effects of singular errors on numerical solution of the local skin friction in Blasius flow. 60 7. Phase of high-frequency oscillations in Blasius flow. 70 8. Amplitude of high-frequency oscillations in Blasius flow. 71 9. Phase of intermediate-frequency oscillations in Blasius flow (XW = 5.06). 74 10. Amplitude of intermediate-frequency oscillations in Blasius flow (XW = 5.06). 75 11. Phase of intermediate-frequency oscillations in Blasius flow (XW = 2.45). 76 12o Amplitude of intermediate-frequency oscillations in Blasius flow (XW = 2.45). 77 13. Phase of intermediate-frequency oscillations in Blasius flow (XW = 1.40). 78 14. Maximum lagging phase angle of intermediate-frequency oscillations in Blasius flow. 79 15. Time-averaged velocity profiles of high and intermediatefrequency oscillations in Blasius flow. 80 vi

LIST OF FIGURES (Concluded) Figure Page 16. Local velocity fluctuations of intermediate-frequency oscillations in Blasius flow (W = 5.818 x 10-5). 81 17. Wedge flow. 83 18. Development of the local velocity profile in wedge flow. 85 19. Development of the local skin friction in wedge flow. 86 20. Numerical solutions of the fully developed velocity profile in wedge flow. 87 21. Development of the local temperature profile in wedge flow. 89 22. Development of the local Nusselt number in wedge flow. 90 23. Rectangular space domain. 94 24. Function F(Y) and dF/dY. 96 vii

NOMENCLATURE Roman a arbitrary constant A arbitrary positive constant, Equations (3.14) Al point at coordinates (Xi,Tn,U ij Vijij,,j,j ~n ~o-n n n an A2 point at coordinates (Xi,Tn,Ui'jVijii, ij ji'j) b arbitrary constant c arbitrary positive constant Co maximum value of C. Ci positive constant, Ci = AY/AXi Cp specific heat under constant pressure process Er rounding error, Equation (3.74.a) Et truncation error, Equation (3.30) Et truncation error, Equation (3.32) Et truncation error, Equation (3.34) f() arbitrary function F() arbitrary function g() arbitrary function i any integer greater than or equal to 1 j any integer greater than or equal to 1 Ki finite positive constant K'i finite positive constant Li finite positive constant viii

NOMENCLATURE (Continued) Li finite positive constant finite positive constant, Equations (3.64) m positive integer, Figure 2 il positive integer, Figure 2 m2 positive integer, Figure 2 Mn maximum value of i Unji n any integer greater than or equal to 0 n Nn maximum value of leij | -(69/6Y)y=0 Nu Nusselt number, Nu = - ( e-eo ) p static pressure n Uo n P maximum value of I( ~-) + (Uo x)ni Pr Prandtl number, Pr = pvC /K q heat flux KU3 Q dimensionless heat flux, Q = q/( c vCp R local radius of curvature, Figure 1 Re Reynolds number, Re = ucxc/v S positive integer, Figure 2 t time variable T temperature u x-component velocity, Figure 1 Au amplitude of oscillation of u, Equation (4.1) U dimensionless form of u, U = u/uc ix

NOMENCLATURE (Continued) AU dimensionless form of Au, AU = Au/uc v y-component velocity V dimensionless form of v, V = v/uc W dimensionless frequency, W = /(uc/v) x space variable along the surface of the obstacle, Figure 1 Xi any independent variable X dimensionless form of x, X = x/(v/uc) Xi value of X at space grid point (i,j), Equations (3.16.a) and (3.16.b) AXi space interval between grid points (i,l) and (i-l,l), Figure 2 y space variable perpendicular to the surface of the obstacle, Figure 1 yj any independent variable Y dimensionless form of y, Y = y/(v/uc) Yj value of Y at space grid point (i,j), Equations (3.16.a) and (3.16.c) AY space interval between grid point (i,j) and (i,j + 1), Figure 2 and Equations (3.14) z any independent variable Z quantity defined by Equation (3.49.a) Greek za difference between U and U, Equation (3.57.a) maximum value of Icn 1, Equation (3.61.a) 1fi difference between V and V, Equation (3.57.b) 1n maximum value of I|ijl, Equation (3.61.b) ln maximum value of Bk (k=O,l,...,n), Equation (3.62)

NOMENCLATURE (Continued) Yn quantity defined by Equation (3.69) Yij rounding error, Equation (3.75.a) 6 some constant between 0 and 1 E instability error, Equation (4.49.b) quantity defined by Equation (3.56.c) quantity defined by Equation (3.56.a) qi rounding error, Equation (3.75.b) a dimensionless form of T, 8 = T/(u2/Cp) 8 phase angle of local oscillations ~K thermal conductivity v kinematic viscosity quantity defined by Equation (3.56.b) p density a some constant between 0 and 1/2 T dimensionless form of t, T = t/(v/u2) Tn value of T at the nth time level, Equation (3.17) 0p arbitrary function (D arbitrary function function defined by Equation (3.50) w frequency of oscillation in free stream, Equation (4.1) Subscripts b upper bound c characteristic quantity xi

NOMENCLATURE (Continued) ij condition at coordinates (Xi,Yj) o condition at free stream (Y + co) ob condition of the moving obstacle s initial condition (T = 0) w condition at the wall of the obstacle (Y = 0) X partial differentiation with respect to X Y partial differentiation with respect to Y YY second partial differentiation with respect to Y Xo condition at main stream Superscripts ()(i) ith approximation on condition at the instant Tn oscillating quantity, Equations (2.14) condition at some point between Al and A2 (7) time-averaged quantity, Equations (2o14) () exact solution of PDE (in distinction from the exact solution of the corresponding FDAE) computer solution of FDAE (in distinction from the exact solution of the same FDAE) Abbreviations cps cycles per second FDA(s) finite difference approximation(s) xii

NOMENCLATURE (Concluded) FDAE(s) finite difference approximate equation(s) fl[(B arithematic operation'B' according to the floating point number system fps feet per second PDE(s) partial differential equation(s) JBI absolute value of'B' xiii

ABSTRACT A finite difference method for computing velocity and temperature profiles of an unsteady, incompressible, laminar boundary layer around a two-dimensional cylinder of arbitrary cross section is develcped, Blowing or suction may be present on the wall of the cylinder. The thermal boundary condition at the wall may be either specified wall temperature or specified heat flux at the wall.l The governing finite difference equations are explicit, the velocity and the temperature at the next time step can be directly computed in terms of those at the current time step. Various errors associated with the finite difference method are studiede Great effort has been made to derive the stability and convergence conditions of the method. The upper bound of the local rounding errors is estimated. Two examples', one oscillation in Blasius flow and the other impulsive start of wedge flow, are given; the numerical results are compared with the existing analytical and experimental results. Tt is concluded that the present method can be used fcr the aforementioned computation with high accuracy except at and near a singular point where the singalar errors become significant. xiv

CHAPTER I INTRODUCTION It is well known that the motion of a viscous, incompressible fluid can be described by the Navier-Stokes equations coupled with the equation of continuity. The general solution of these equations has not yet been available within the present knowledge of mathematics. However, by assuming the viscosity effect important in a thin fluid layer over boundaries, Prandtl (1) introduced his well-known boundary-layer equations which have simpler forms than full Navier-Stokes equations. After Prandtl, most research has been concentrated on the steady flow. It is only during the last two decades, the unsteady laminar boundary-layer flows have been treated extensively. The unsteady laminar boundary layer theory presents many challenging problems which may be classified as the periodic and starting-ending problems. In the former case, the periodic behavior may be introduced by oscillating the obstacle or the fluid. As a result of this, the response velocities of the fluid particles in the boundary layer are also periodic with respect to the time. Many practical important problems fall into this category, such as the flow past a turbine blade or past an oscillating airfoil. The latter case includes all the flow problems in which the obstacle or the fluid starts to move from the rest with a time dependent velocity. The flow generated by an impulsive motion of an obstacle in a stationary fluid field is a special example of this category.

2 Generally speaking, the method of solving an unsteady, laminar boundary-layer problem will fall into one of the following six ways: 1. Similar solutions. 2. Semi-similar solutions. 3. Series solutions. 4. Successive approximations. 5. Integral method. 6. Other methods. A detailed review of the literature about the above methods will be given in Chapter II, Here we shall briefly discuss these methods. The method of similar or semi-similar solutions reduces the number of independent variables of a two-dimensional unsteady laminar boundarylayer problem from three into one or two, respectively. Such similar or semi-similar solutions exist only when the free stream velocity of a given boundary-layer problem behaves in a certain specific way. Therefore, these methods are restricted to problems in which the obstacle has a particular geometry and moves with a particular time-dependent velocity. The method of series solution is based on the expansion of the local velocity into a power series of some parameter or coordinate. Because of the mathematical complexity, only the first few terms of the series may be calculated. The series solution deviates appreciably from the exact solution when the value of the expansion parameter or coordinate becomes large. The method of successive approximations is actually an iteration procedure. In this case, a first approximation is obtained by physics, and the successive approximations are then obtained by iteration. Therefore,

3 the convergence of the final approximation to the exact solution of the problem largely depends on the choice of the first approximation. Also the higher approximations of this method are rather complicated algebraically. The integral method is an approximate procedure, because of the fact that the profiles obtained by this method only satisfy the boundary conditions and the integral form of the governing equations. None of the foregoing methods meet the requirements of generality and high accuracy. In recent years, advances in high-speed digital computer technology have made it possible to develop some numerical methods applicable to complex fluid flow problems which were beyond reach by the foregoing methods of solution. For example, De Saute and Keller (2) recently studied the transition form laminar to turbulent flow over a flat plate by finite-difference technique. Der and Raetz (3) developed a numerical method for general three-dimensional laminar boundarylayer problems. Fromm (4) computed the development of the vortices behind a two-dimensional obstacle. Also Yang (5) was able to study the hydrodynamic stability of two-dimensional unsteady incompressible laminar boundary layers by means of the numerical iteration scheme and Barakat (55) studied the transient natural convection in closed container by finite difference method. It is the purpose of this work to develop a finite difference method to determine both unsteady velocity and temperature profiles in a laminar

boundary layer over a cylinder of arbitrary cross section moving with an arbitrary time-dependent velocity in an incompressible fluid. Suction or blowing on the surface of the solid cylinder may also be included. The wall temperature or the heat flux at the wall of the obstacle may vary arbitrarily with the time and the distance from the front stagnation point. In order to illustrate the applicability of the method, two numerical examples, one periodic and the other starting-ending type, are considered. The results are compared with those obtained by the analytical methods and, whenever possible, by experimental measurements.

CHAPTER II REVIEW OF THE PREVIOUS WORKS In this chapter, a brief survey of the existing literature about unsteady incompressible laminar boundary layers will be given based on the methods classified in Chapter I. A. SIMILAR SOLUTIONS The solutions of the velocity component u in an unsteady two-dimensional laminar boundary layer are called "similar" if the two velocity profiles at the different coordinates x and at the different instants t differ only by a scale factor. More precisely, ~uU(('2~~r~t) = {(~r) ~(2.1) U0(Zt) where, F 2(r- t) (2.2) The transformation defined in the equation (2.2) is called the "similarity transformation" and the function'g' is called the "scale function." Through this transformation, the original governing PDE* for u can be reduced to an ordinary differential equation. Consequently, a considerable mathematical simplification of the problem is achieved. *Partial Differential Equation.

6 Schuh (6) showed similar solutions may be obtained if the velocity at the outer edge of the boundary layer is proportional to ta, ebt or x/t, where'a' and'b' are arbitrary constants. Yang (7) studied the similar solutions of a stagnation flow of which the free stream velocity can be expressed as ax/(l-bt). Perhaps, the most complete work about the similar solutions of an unsteady boundary-layer flow was the one given by Hayasi (8) who studied the possible similar solutions in an unsteady quasi-two-dimensional incompressible flow. The definition of a "quasi-two-dimensional" flow can be found in Reference (9). The similar solutions found by Schuh (6) and Yang (7) are only the special cases of Hayasi's results when the flow becomes exactly two-dimensional. A systematic study of the topics by means of the "free parameter method" is described in a recently published book by Hansen (10). The rest of this section is devoted to a brief review of this method. Let cp be some physical unknown to be solveds Assume it satisfies a certain governing PDE and depends on some independent variables xl, x2,..., xi and z. Let the boundary conditions in z be (I(t,, ~X, —-, X )-,,, f= (2.3.a) LoAM, (y1pjI, y — s(1G g - L-t, _ ) (2.3.b)

where, +}8-~>(^Xi,:,^X,; I~(2.5)'r' is the unknown'free parameter" to be specified. By means of the transformations given by the foregoing equations, one hopes to transform the governing PDE for p into an ordinary differential equation for'f' under certain imposed conditions. If the given physical problem meets these conditions, then a similar solution of Cp exists and vice versa. It is also very important to examine the boundary conditions for'f'. According to equations (2.3) through (2.5), the boundary conditions for'f' are ( 2.6.a) J, +(F-= I (2.6.b) where, Zo -- 1-'(~,,~., 1,;, &-o) ( (2.7.a) XR- T lfi2,.X_~~~-l Xi a) (2.7.b) If a similar solution of cp exists, the values of Fo and Fi must be independent of xk (k=l,...,i). These requirements become difficult to be achieved when both zo and zl take finite values, In general, the problems of finite geometry do not have similar solutions.

B. SEMI-SIMILAR SOLUTIONS The solution of the velocity component u in an unsteady two-dimensional boundary layer is called "semi-similar" if it can be expressed as a function of two independent variables. An example of such flow is the incompressible boundary layer whose free stream velocity takes one of the following forms: atb, a/(l+bt), f(\)ta, where x = xt+l The first two cases in the above example were treated by Cheng (11) and the last one was studied by Hassan (12). Recently, Hayasi (13) gave an extensive treatment to the semisimilar solutions of an unsteady quasi-two-dimensional incompressible laminar boundary-layer flow; two simplest cases, time-dependent and space-dependent semi-similar boundary layers, were studied. The semi-similar solutions are obtained by the "group-theoretic method" (10). Because of the length of the background required, this method will not be reviewed here. In short, one tries to find a oneparameter transformation group under which the governing PDEs can be transformed into the different variables, Therefore, by doing so, one reduces the number of independent variables by one. As pointed out in the Reference (10), there are two disadvantages to employing the group theoretic method. The first is the boundary conditions which are not taken into account until the entire analysis is completed. The second is the uncertainty in choosing a proper transformation group.

C. SERIES SOLUTIONS For boundary-layer problems of starting-ending type, the stream function is usually expanded in ascending powers of some parameter related to t or x. Blasius (14) employed this method to study the impulsive motion of a cylinder and the uniform acceleration of a cylinder in an incompressible fluid. In both cases, the stream function has been expanded in ascending powers of A which is defined as tuc/xc, where uc and xc are the characteristic velocity and length, respectively. In the first case, only the first two terms of the series have been calculated by Blasius, the third term has been given later by Goldstein and Rosenhead (15). Watson (16) generalized the Blasius' problems by studying the boundary-layer growth on a cylinder starting from rest with the velocity uob proportional to ta or e lbt. For the case when Uob = elbl t Watson expanded stream function in ascending powers of the parameter'k' which is defined as J uOb(t)dt. The boundary-layer problem of a semi-infinite flat plate starting initially from rest in an incompressible fluid with an arbitrary velocity uob was studied by Cheng and Elliott (17) by means of the series t solution method also. A parameter Q, which is defined as x/ uob(t)dt, has been used in power series expansion of the stream function. For periodic boundary-layer problems, the local velocity is expanded in ascending powers of Aus/um which is the relative amplitude of oscillation of the free stream. For instance, Kestin, Maeder and Wang (18)

10 utilized this method to study the effect of a longitudinal sine wave on the boundary-layer along a flat plate. Hori (19) used the same method to investigate the boundary-layer around a cylinder of arbitrary cross section in a fluctuating main stream. Clearly, the mathematical complexity does not allow the computation of many terms in a series solution. It is doubtful, unless checked by experiments, that the series solution converges to the exact solution when the parameter involved in expansion becomes large. As pointed out by Rosenhead (21), for example, the series solutions employed in Reference (14) are not valid for values of A greater than 0.3. The foregoing method of series solution is actually a small perturbation technique. In starting-ending problems,,the time or space coordinate is used as perturbation quantity, hence, the method may be classified as the coordinate perturbation. In periodic problems, the relative amplitude of oscillation is used as the perturbation parameter, therefore, the method becomes a parameter perturbation. It should be remarked that the choice of an appropriate perturbation parameter is by no means an easy task. The details of the small perturbation method may be found in the book by Van Dyke (21). D. SUCCESSIVE APPROXIMATIONS The process of the successive approximations has been described in Reference (22). An example on the application of this method can be found in Reference (25). The method is based on the physical argument

11 that if a solid obstacle is initially at rest relative to the fluid, the boundary layer formed on the surface of the solid obstacle is thin right after the start. Consequently, the viscous term in unsteady boundarylayer equations is large, while the convective terms retain their normal values. Therefore, at early stage of the transient phenomenon, it may be assumed that the viscous term is balanced by the acceleration plus the unsteady portion of the pressure gradient. According to the foregoing argument, the first approximations u(l) and v(l) satisfy the following linear PDEs: dA" b - (2.8) bt b t (r) + tO (2.9) with the boundary conditions y = 0: u(l) = v(l) = 0; y: u(l) = uo(x,t). The second approximations u(2) and v(2) may then be calculated from the following equations: zt S aLA +; - u,, U a - ( 2. rt At +b =o (2.11) with the boundary conditions y = 0: u(2) = v(2) = 0; and y + o: u(2) uO(xt). The third and higher approximations are obtained by repeating the above iteration procedure. However, the complexity of algebra grows at

12 a staggering rate. It should be noted that the successive approximation method fails to distinguish an unsteady boundary-layer along a semiinfinite flat plate from that along an infinite flat plate. This is due to an oversimplification in the first approximation of the boundary-layer equations. E. INTEGRAL METHOD The methods discussed so far give local profiles which staisfy the differential form of boundary-layer equations. A high accuracy may be obtained by carrying out enough terms. When it is desired to predict approximately the distribution of the skin friction, the position of the flow separation, etc., an'approximate' method such as the integral method becomes quite useful. Schuh (24) studied the unsteady boundary-layer problems by using this method. It was shown that in the general case of any type of velocity distribution over outer edge of the boundarylayer, regarding space and time, the integral of boundary-layer equations could be reduced to two simultaneous differential equations, one being an equation of characteristics. Two numerical examples were given: first, the growth of the boundary-layer around a circular cylinder for impulsive motion, second, the generation of the boundary-layer along a flat plate for impulsive motion and the steady acceleration. In the first example, Pohlhansen's (25) velocity profile was used, while in the second example, Hartree's (26) velocity profile was assumed. Schuh' s procedure was later modified and extended by Yang (27) to problems includ

13 ing heat transfer. Other works related to integral method may be found in References (28,29). Recently, Hayasi (30) developed an approximate method for calculating unsteady quasi-two-dimensional laminar boundary-layers. The velocity profiles were assumed to be given by one-parameter family of curves. The profile parameter which was selected to be the non-dimensional shearing stress at the wall was determined by the momentum and energy-integral equations. It was found, after comparison with the exact solutions, that Pohlhansen's velocity profile gave good results for decelerating main flows, while Hartree's velocity profile was good for accelerating main flows. F. OTHER METHODS For a special unsteady boundary-layer problem, one may make some assumptions which are particularly suitable to the conditions of the problem under consideration. Thus the method to be developed based on these assumptions may not be used for problems of the different nature. To illustrate this point, a brief review on C. C. Lin's (31) theory of the harmonic oscillations is given next. Lin considered an unsteady two-dimensional boundary-layer with a main stream oscillating periodically about a mean value. The mathematical difficulty of the problem was simplified by taking time-averages of the terms in boundary-layer equations in a manner analogous to that used in turbulent flows. Denoting the oscillating components by primes,

14 Lin assumes ULYXt)- uo(X) (x,t), IA= (2.14.a) U(Lt(Y,.4 L4t)\(x &(xtfLt),t L'=O (2.14.b) U'(),it)= — (a, ) t-I(X,,-,t), l[' -=0 (2.14.c) V(Xz~-t) 1 t k'(X,t),' - (2.14.d) By restricting the oscillation of the main stream to very high frequency range, the unsteady portion of the boundary-layer equation is linearlized in the form: -C = 1f it- o U a(2.15) After calculating the longitudinal velocity component u' from the equation above and the corresponding boundary conditions, the transverse velocity component v' is obtained from the unsteady portion of the continuity equation 0)G th a =0 )(2.16) The steady portion of the boundary-layer equations reduces to JiX to ~ s _ a, D ax,) at P }- (2.17) BX3L -- (2.18)

15 From the above two equations, the time-averaged velocity components u and v may be calculated with the corresponding boundary conditions. Clearly, the foregoing method of solution can not be used for fluid flow problems other than those related to high frequency oscillations. Finally, it is worth noting that the same problem was later considered by Gibson (32) by the so called "splitting solution' method.

CHAPTER III FORMULATION OF THE PROBLEM A. DIFFERENTIAL FORM A curvilinear orthogonal coordinate system attached to the surface of a two-dimensional cylinder is selected as shown in Figure 1. Free Stream Velocity, Uo(x,t) Local Velocity, U (x,y, t) Main Flow Velocity, U()(t) R(x) Figure 1. Curvilinear orthogonal coordinates. The laminar boundary-layer equations with respect to these coordinates are 16

17 h e - FR (5.1.b) +- R all + * =0O (532) bT- DX (feW j (5 5) where R(x) is the local radius of curvature of the wall, reckoned positive when the wall is convex outward. The derivation of the above set of equations can be found in texts on boundary-layer theory, for example, in References (20, 22). The equations are valid within the following assumptions: (1) The fluid has constant properties. (2) The solid obstacle may be stationary or may move with a translational motion parallel to the motion of the main flow. (3) The Reynolds number Re based on the characteristic velocity and length of the problem is large enough so that the viscous effects of the flow are limited to a thin layer on the wall. The thickness of this layer is one order smaller than that of the characteristic length. Outside the velocity boundary layer, the flow is assumed to be inviscid. Also the variation of the velocity in the direction of the y-axis is one order greater than that in the direction of the x-axis. (4) The Prandtl number Pr must be order of unity or higher. This implies the thickness of the thermal boundary-layer is the same order or less than that of the velocity boundary-layer. Consequently, the

18 variation of the temperature in the direction of the y-axis is one order greater than that in the direction of the x-axis. (5) The viscous dissipation is negligibly small. (6) The magnitude of the local radius of curvature of the wall is one order greater than that of the local velocity boundary-layer thickness. (7) The variation of the local radius of curvature along the direction of the x-axis is at most of the order of unity. In general, the Navier-Stokes equations are derived relative to a space-fixed coordinate system. Consequently, the resultant boundarylayer equations are restricted to a coordinate system fixed in space. However, when the coordinate system moves with a translational motion parallel to the main stream, the uniform inertial force per unit volume due to the translational motion of the coordinate system appears in equation (3.1.a) in the form of an induced pressure gradient. The transform from the inertial force to the pressure gradient is achieved instantaneously by the acoustic waves in an incompressible medium (33). Therefore, as indicated by the second assumption, the boundary-layer equations (3.1) and (3.2) are valid for both stationary and moving obstacles. Based on the third assumption, equation (3.1.b), it may be concluded that the magnitude of the pressure variation across the boundary-layer is second order. Therefore, hereafter the pressure is assumed to be constant across the boundary-layer, and the pressure gradient in equation (3.1.a) is same as that given by the momentum equation of the inviscid

19 flow at the outer edge of the boundary-layer. Thus _ +LI0 -. I _ _ P (3_4) -ava"l~ X Uo XI where uo and vo are the streamwise and the normal components of the velocity at the outer edge of the boundary-layer, respectively. Because of the boundary layer thickness, the first approximation uo and v0 of equation (3.4) may be assumed to be the corresponding velocity components evaluated at boundaries. Since the normal component of the fluid velocity vanishes at boundaries, equation (3.4) reduces to ov1a - t U >u~ -- I 38 -(3.5) bt ( b (3O5) For boundary-layer problems with blowing or suction, equation (3.5) is approximately true provided that the magnitude of blowing or suction velocity is smaller in order than that of the free stream velocity uo. Accordingly, equations (3.1) become bt a X + V t_ +_ _ a A (3.6) bt + X u 7-~-+ ~ — t X "~ )~ Equations (3.6), (3.2) and (3.3) are the governing differential equations for dependent variables u, v and T in a laminar boundarylayer flow. Their boundary conditions are for x = O, y > 0 and t > O; u = uo(O,t) (3.7.a) T = To(O,t) (3.7.b) for x > O y = O and t > O; u = O (37.c) 0or x >(7

20 v = vw(x,t) (3.7.d) T = Tw(x,t) (3.7.e) for x > 0, y -+ o and t > 0; u = uo(x,t) (3.7.f) T = To(x,t) (3-7.g) where uo and To are the velocity and the temperature of the free stream, respectively, while vw is the blowing or suction velocity normal to the wall and Tw is the wall temperature. When the blowing or suction is absent, vw is equal to zero. The thermal boundary condition given by equation (3.7.e) specifies the wall temperature. In case when the heat flux at the wall is specified, this boundary condition becomes then - h (39 )g' a = t(x0t) (3-7.h) where qw(x,t) is the heat flux at the wall, reckoned positive when heat flows from the wall to the fluid. The initial condition of the problem is at t = 0; u = us(xy) (5.8.a) v = vS(x,y) (3.8.b) T = Ts(x,y) (3.8.c) where us, vs and Ts are the initial distributions of velocity components u, v and temperature T, respectively. However, only us and Ts need to be specified according to the physical situation of the problem; once us is given, vs can be readily obtained by integrating equation (3.2). Next, the above governing differential equations will be nondimension

21 alized. Defining the following dimensionless quantities: U —7 1 -tt- 0- (tC) t where uc denotes some constant characteristic velocity, equations (3.6), (3.3) and (3.2) take the following forms, respectively: ~13 t~)X = ~U o -4(3.9) bit + - a t y -- V -+ U 3x ay tt + ~U +V y P (3.10) a: + =O (3-.11) The dimensionless boundary conditions are for X = 0, Y > 0 and T > 0; U = Uo(O,T) (3.12.a) e = eo0(oT,) (3. 12.b) for X > 0, Y = 0 and T > 0; U = 0 (3.12.c) V = V(x,T) ( 3.12.d) 8 = GW(x,T) (3-12.e) for X > 0, Y + oo and T > 0; U = Uo(X,T) (3.12.f) e = &o(X,T) (3.12.g) If the heat flux at the wall is specified, the thermal boundary condition given by equation (3.12.e) becomes (- aY)y= ( X Q,(XJt) (3.12. h)

22 The dimensionless initial condition is at T = O; U = Us(XY) (13 o a) V = Vs(X,Y) ( 3.13 b) e = -s(X,Y) ( 3-53 c) B. FINITE DIFFERENCE FORM 1. System of Space Grid Points To obtain the solutions to a set of partial differential equations by the finite difference approximations, it is necessary to establish a system of space grid points in the domain of interest. The space grid used in this work is shown in Figure 2. The interval AY between two neighboring grid points along a line where X = Xi is assumed constant. However, this constant varies with the value of'i' as follows: AY = A, for 1 < i < m2 (314. a) AY = 2A, for m2 < i < m (3.14.b) where i, m and m2 are positive integers, and'A', any positive constant. Likewise, the interval AXi between two points (i-l,l) and(i,l) on the Xaxis also depends on the value of i. Thus AXi = fn(i), for 2 < i < m (3515) where fn(i) is any real function of i, its value is positive and increases monotonously with i. All grids are smaller in the neighborhood of the stagnation point

Main Flow (m,2S +) H —--- -- _ H —- - H —-- - I —-- -- -- Am? I I~~c- L (IS(m 2Sl- ~ ~ ) FigureM 2. Sytem o sc rpn Figure 2. System of space grid points.

and become larger in the downstream region. The reason of choosing such nonuniform grid size will be explained in Section (B-8) of this chapter. Figure 2 also shows a "stair" shaped domain of interest. This is specially designed for boundary-layer problems. The space grid point (1,1) corresponds to the origin of the coordinate system described in Part (A) of this chapter. The relations between the coordinate point (Xi,Yj) and the space grid point (i,j) are given as follows: X1 - 0, Y1 - 0 (3.16.a) Xi-X AXk, for 2 < i < m (3.16.b) k=2 2 < j < S+1 (if 1 < i < ml) Yj - (j-l)AY, for < j <2S+1 (if ml < i < m) (3>16.c) where i, j, m, ml and S are all positive integers. Special attention is paid to the following two points: firstly, AY in equation (3.16oc) takes two different values according to equations (3.14); secondly, for a given set of AXi and AY, the integers m, mli, m2 and S are so chosen that the corresponding "stair" shaped space domain covers the entire boundary layer flowo Besides the space variables X and Y, time variable T is also involved in the study of an unsteady boundary-layer problem. Hence, in addition to the space grid, a uniform positive time step AT will be specifiedo The relation between the time variable Tn and the time step AT is

25 [],: 1=: }[. /5'~"~~ (3.17) where n is a positive integer or zero. Based on the foregoing space and time intervals, the correspondence between the coordinate point (Xi, Yj, Tn) and the grid point (i,j,n) is one to one. Thereafter, any physical property cp at the coordinate (Xi, n Yj, Tn) is uniquely described by CPi j. 2. Approximation of Derivatives by Finite Differences Let cp(X,Y,T) be any function with continuous and bounded partial derivatives of sufficient high orders. By Taylor's series with remainder, it yields ~qi, -o -l"~~ ~~~, i (A + E(3.18) - 9 /Y(ay)P (+AY)'( )?: (3.19) where 6 is some number between 0 and 1, and it generally takes different values in the different terms unless specified. Dividing the foregoing equations by AY, after some arrangements, they become ( 4 <3)4__ _1 6 (3.20) (-vyr),, --- 1 (aAY I _ __ If Taylor's expansions of cp(X,Y,T) in equations (3.18) and (3.19)

26 are carried out including the fourth term, they become I, -NY + _- + (3.22) Subtracting equation (5.23) from equation (5.22), then dividing the result with AY, gives A ('Y) (as),i- z (tlY)- - I Z [L~]it -- (3~.24)_~] tivesuoflt with AY.,rgespect t Y tw The first terms on the right-hand side of equations (3.20), (1.21) and its corresponding FDA is known as the truncation error of FDA. The truncation error of the central FDA involves the higher ordered derivatives of p with respect to Y than that of the forward or backward FDA. In general, the truncation error of any FDA approaches zero as the interval AY approaches zero. In this limiting case, FDA becomes the exact representation of the corresponding derivative. Finally, it is easy to show that the second partial derivative of cp with respect to Y at (Xi, Yj, Tn) can be represented by the following formula: _ C ) I (3425) *Finite Difference Approximations.

27 3. Finite Difference Form of Boundary-Layer Equations As is well-known, the boundary-layer equation given by (3.9) is valid only if there is no reverse flow within the boundary layer. Otherwise, the boundary layer becomes so thick that the third assumption in Part (A) of this chapter no longer holds. Therefore, the streamwise velocity U will be restricted to zero or positive values. Consequently, problems involving instantaneous local separation are excluded from the present investigation. Now equations (3.9) and (3.10) are considered at the coordinate (Xi, Yj, Tn) and equation (3.11) at (Xi, Yj-1/2,Tn+l) as follows: Ibecom ~n L7 nes 7- VI VI' t""'"" +U q) "a T ) + (3a) (5.26) ( fi UT +j17.i, P (I (3.27) +~): (4. u o (3.28) Note the use of forward, backward and central exact finite difference representations for the derivatives in the first, second and third terms on the left-hand side of equation (3.26), respectively. Furthermore, the last term on the left-hand side of the equation is represented by the formula given by equation (3.25), then PDE stated by equation (3.26) becomes (5.29) \a7'li' ~U~x~L E

28 where {"t.1" I \ ('rt i( 3f-i A i) T43U. )x2] +;,1 Y) ___Y~'r -Y4 V4 (5.530) Similarly, PDe given by equation (329) remain in27) becomes Next equation (3.28) is expressed in finite difference form. An LT IJ AXI I U Y where It should be remarked that the first two terms on the right-hand side of equation (ij-1/29) remain in differential form. These terms can easily be calculated once the functional form of the free stream velocity U0 is known. In boundary-layer problems, U0 is available from the corresponding inviscid flow problems. Next equation (5.28) is expressed in finite difference form. An exact central finite difference representation will be used for the second term of the equation. The first term of the equation, namely, the derivative of U with respect to X at the grid point (i,,j-l/2, n+l) is first (i,j-l,n+l) and (i,j,n~l), then each of these derivatives is replaced by

the backward finite difference form. The complete equation becomes then z(Ax) A t (3.33) where ( )( 2h ~aYV3t (3 34) and 6 is some number between 0 and 1/2, its value may vary from term to term. The above treatment of the continuity equation was first used by n+l Wu (34). Its main advantage is that the normal velocity V.. for entire space domain can be directly calculated according to the order j = 2, n+l 3,4,... provided that the streamwise velocity Ui,j is known for entire n+l space domain and the boundary condition Vi.1 is given for all values of'i'. In addition, aV/6Y term in the continuity equation is replaced by a central FDA which is a better approximation than a forward or backward FDA. 4. Numerical Computation Procedure Equations (3.29), (3.31) and (3.33) are the exact finite difference representations of PDEs given by equations (3.26), (3.27) and (3.28), respectively. Hence, no truncation error has been introduced in these equations. Suppose now the space intervals AY, AXi and the time step AT in equations (3.29), (3.31) and (3.33) are chosen to be sufficiently

50 small such that the terms Et, and Et and Et in these equations become negligible Hence, they can be approximated by the following FDAEs*, re spectively' ~ U rC ~ XIi A _ " i — It' (L ff'1,+ (EU tm] Lt L Voi aLi;,. Li' L-1- 3_ --- 1 (tY) ) (3'35) The truncation errors of FDAEs (3.35), (3536) and (35.537) are AT-Et, AT.Et and AY-Et respectively, where Et, Et and Et are given in equations (5.530), (3532) and (3.34), correspondingly. Directly interpreted from equation (53.12), the boundary conditions are n+l n+l for S+1 > j > 2 and n > 0; U1 = (UO)l (5.58.a) n+l n+l Qe1 = (eo)l (3.8b) n+l for m > i > 1 and n > 0; Ui 1 = (538.c) n+-1 n+l Vil = (VwDei (A 38ad) *Finite Difference approximate Equations.

31 n+! n+l G~i9.. (o)Q (3.38.e) n+I n+l Ui. k = (Uo)i 3 38rf) n+l n+l ei,k = ()i (338g): 2(3 &4 (~) ( (338oh) where k = S + 1, for ml > i > 1; k = 2S + 1, for m > i > ml. The initial condition, according to equations (3i13), is 0 U = (U) j (3.39.a) Vi j = (Vs)i j (3.39.b) 9iij = (eS)i~j (3.39.C) where (i,j) covers all the space grid points in the domain of interest. FDAE given by equation (3o36) is called "explicit" because of the fact that the temperature variable 9 at time step (n+l) can be directly computed from this equation once the values of O, U and V at time step n are known and the boundary conditions of 9 are specified. Suppose, instead of a forward FDA, a backward FDA were used for the first term on the left-hand side of equation (3.27), then the corresponding FDAE of equation (3.27) would be tL,~~=r N.f "" A +'~:'-. i t,\~ ck. ": i:>-~" lI,_.I-t.i..L I J _.41i..:i::.: i'~hY i'~~~~~~ (nsi"~'

32 which is equivalent to In this case, even the values of U and V at time step (n+l) are known, In this case, even the values of U and V at time step (n+l) are known, a large set of simultaneous algebraic equations must be solved in order to compute the temperature variable e at time step (n+l) in terms of the values of G at time step n. This kind of FDAE is called "implicit". In general, an iteration procedure is necessary to solve a large set of simultaneous algebraic equations on a digital computer. Therefore, it is clear that one advantage to choose an explicit system of FDAEs is to avoid time consuming iteration procedure. Equation (3.35) coupled with equation (3.37) is also an explicit FDAE. Equations (3.35) through (3.39) are actually used in the numerical. computation. The procedure of the numerical computation may be outlined as follows: (1) Determine the numerical values of m, mln m2 and S. (2) Choose a proper set of AT, AY and AXi. n n n (3) Set n 0, and compute Uij, Vij and nQij according to the initial conditions given by. equations (3.39) for all possible grid points (it ). (4) Set i -= 2. n+l n+l (5) Compute Uij and &iij according to equations (3.35) and (3.36) respectively for all possible values of'j' except the two values j=l and j = S+1 (for mi > i > 2) or j = 2S + 1 (for m > i > ml).

n+l (6) Compute Vi: according to equation (3.37) and the boundary condition given by equation (3.38.d) for all the values of'j' (7) Increase the value of'i' by one. (8) Repeat steps (5) through (7) unitl i = m is included. (9) Increase the value of'n' by one. (10) Repeat steps (4) through (9) until'n' takes any desirable value. In the fifth step, the stability conditions of FDAEs (3.35) and n+l (3.36) should be examined before proceeding to the computation of Ui, n+l and ij. If they are satisfied, the computation can be continued as outlined above. Otherwise, equations are unstable and the computation should be stopped. In this case, a new set of AT, AY and AXi should be chosen such that the stability conditions are satisfied. A general definition of the "stability" of a FDAE is given in the next section. The stability conditions of FDAEs (3-35), (3.36) and (3.37) are discussed in Section(B-6) of this chapter. 5. General Error Analysis To justify the use of a finite difference method, various errors associated with the method must be examined first. For a detailed discussion on the errors involved in the approximation of finite difference method to steady, three-dimensional, compressible laminar boundary-layers, Raetz (35) may be referred to. Here these errors are classified as the truncation, instability, singular, rounding and machine errors.

34 The truncation error of a FDAE is defined as the difference between the FDAE and its corresponding PDE. The existance of this error introduces the problems of "consistency' and "'convergence". A FDAE is said to be consistent with its corresponding PDE if its truncation error vanishes when all the space intervals and the time step approach zero. In our case, FDAEs given by equations (3.35) through (3.37) satisfy this requirement. A FDAE is said to be convergent if its solution approaches the solution of its corresponding PDE as all the space intervals and the time step approach to zero. Note that the consistency of a FDAE does not necessarily imply its convergence; the solutions of both FDAE and its corresponding PDE may not be unique and even may not exist. The convergence conditions of FDAEs (3.35) through (3537) are studied in Section (B-7) of this chapter. The instability errors are those due to the growth or amplification of any local errors (such as rounding errors) as the computation advances. They are directly related to the "stability" of the FDAE. A FDAE is called stable if its solution remains uniformly bounded as the computation advances indefinitely under a set of fixed space intervals. Thus, for a stable FDAE, if any local error is introduced in its solution at sometime level'n' of the computation procedure, it must either remain bounded or decrease as the value of'n' advances. Hence, if the decay of the local errors in a stable FDAE is faster than the generation of the local rounding errors, the computer solution of this FDAE will closely approxi

35 mate its exact solution. On the other hand, when the instability exists, a small error tends to grow as the computation progresses and finally leads to a meaningless computer solution. The stability conditions of FDAEs (3.35) through (3.37) are discussed in Section (B-6) of this chapter. The singular errors are due to the existance of mathematical singularities. At and near a singular point, the velocity or temperature changes so rapidly that the otherwise valid FDAs of their derivatives become inadequate. An example of this kind of error will be given in Section (B-8) of this chapter. The rounding errors are those due to rounding the numbers. IBM 7090 digital computer, by which all numerical computations are carried out, has thirty-six binary bits in every single word. This gives an accuracy of eight decimal digits if a floating-point number system is used in the computer. Any number has decimal digits more than eight will therefore be rounded in the computer, that introduces the rounding error. Likewise, any arithmetic operation on the computer which leads to a number of more than eight decimal digits also introduces the rounding error. An estimate of the upper bounds of the local rounding errors introduced in the numerical n+l n+l n+l computations for Ui,, Gi,j and Vij according to FDAEs (3.35), (3.36) and (3.37), respectively, are given in Section (B-9) of this chapter. The machine errors are those due to the incorrect functioning of the computer machine. The probability that such errors occur is very low.

6. Stability Analysis and Test Here the stability theory is briefly reviewed first. In literature (36,37), Von Neumann's theory has been often used to study the stability conditions for the governing FDAEs of the nonlinear fluid flow problems. Strictly speaking, Von Neumann's stability theory, which may be found in References (38,39,40,41,42) is only valid for initial-value problems with linear, homogeneous governing differential equations. The theory predicts the necessary and sufficient conditions for the stability of a multilevel FDAE which approximates a differential equation of the above category. John (43) studied the stability conditions of the following explicit FDAE: N t=2xatd)~+ At-t (az) (.40) which corresponds to the differential equation, ba =q:)b _ o(/,t) + a,(iX 2a x a tap)x ( are where ao > 0 and the functions ao, aao/ax, a ao//x2 al, aal/ax, a2 are bounded and uniformly continuous in all their arguments. Furthermore, a2 satisfies a uniform Lipshitz condition in cp, that is I an2(,tP,)-av(xtAs)I ~ C Va nd for any values of Pa and tL, Next the stability conditions for FDAEs (3-35), (3.36) and (3-37) are derived. First, these equations are rearranged in the following forms

57 3,,, -- (AT) (T)U (l - (A - T(T)-" -iz, (AT) + (A Tt!, ( t I (A,, (AT = ~u, ) (A()2) o,- ~ [ V Vi z~-)- +t/- LU!',+- tU+/'3 ) Then, it is assumed that (1),A n ) Y and Xi are all positive quantities. (2) A finite time domain Tb is assumed. Thus, Tb > nAT >_ 0. Hence n + 0 implies the limit AT + 0. However, the converse is not true. (3) A finite space domain Xb by Yb is considered. Furthermore, the space intervals AY and AXi are fixed and take finite but non-zero values during the numerical computation. n (4) Flow separation is excluded. That is Ui > 0. (5) The initial values (U) and ()i are bounded by ad o No respectively. Thus, according to equations (A 039) Ui,: <M0 and + V + ~Y 1)H- n+l (6) The bound ary value s (V)i given by equ onthation. (1) ATe AY andarXi are all positive quantities. 0~~~~~~~~~~n+

38 n+l uniformly bounded by (Vw)b. Thus IVill < (Vw)b for all possible values of'i' and at any time level'n' n n (7) The terms (6Uo/6T)i + (Uo1Uo/6X)i on the right-hand side of equation (3.35) represent dimensionless pressure gradient in the boundary layer. It is assumed that this pressure gradient remains finite in the interested space and time domains. Let its maximum value be Po, then ('C + (- _/ - for all possible values of'i' and at any time level'n' Based on the above assumptions, the stability conditions for FDAEs (3.42), (5343) and (3.44) are derived in the following theorem. n+l Theorem: The sufficient conditions for the uniform boundedness of IUi9 n+l given by FDAE (5.42) and li1j1 given by FDAE (3.43) are respectively, — 14 U-N,; 2 ~ 1( 45.a) A >X A (A)245 and {! + 2T (3.46.a) LAX, Pr (Ay)2 anl nl. A (-346.b) n+l n+ l Furthermore~ if IUi,j is uniformly bounded, then Vi j9 given by FDAE (3.44) is also uniformly bounded.

39 n+l Proof: First, the uniform boundedness of IUi, j is proved. According to FDAE (3.42), I Ui~, I I I~b,(AT) _ ( - IW. I under consideration. Hence,by the assumptions (1), (4), (7) and stability conditions (3.45), the above inequality can be written as ori~~~~ + (a-P (UZ (y) (A 2() L2 (AY () n+l Since Mn+l is the maximum value of 1Ui,j l for all the space grid points (i,j) under consideration, therefore, M~,, <_ Mn t z~'l-Po Successive application of the foregoing inequality until n = 0 gives Mlh +l Vh Mot (n ()) Pt or IU |I< ~Mo+ (t lT)RP

40 The second assumption implies nAT < Tb, therefore, |U1j|.-.4 Po~r (3.47) Since Mo, Tb and Po are all finite and are independent of'n', n+l hence, it is concluded that I Ui is uniformly bounded. The uniform boundedness of @i,n+ Under the stability conditions given by equation (3.46) can be proved in a similar way. jn+l Next, the uniform boundedness of IVi jl is proved from FDAE (3.44) as follows:,IL d_' t2(A[', I U"'+' 1-U,-+0,+J UlJ wV4,1 A ( IX l' l A~i. i, I b~Y fit'~l'. U. ji The fourth assumption and inequality (3.47) imply - U"tlt I OT.l i< P0f ~ 7,n+ lift (J-Q j I ovlkx. of UJ 1-,,a U*dleJ' Mot -P, It follows then tv..714 1V,7 ~"'Y (M0ot',- ) Repeatedly applying the above inequality until j = 2 and also making use of the assumption (6), the following expression is obtained: |VtJ |C (~ibb l)4) \(~

41 Furthermore noting from Figure 2 that j < 2S + 1, the above inequality may be rearranged in the form l" \ti I (V)t2;O PO (3.48) If AY/AXi is denoted by Ci, the assumption (3) implies that Ci is a finite number. Then in terms of the maximum value Co of Ci corresponding to all possible values of'i', equation (4.48) becomes I5,'I, -(Vv)b+2' Co (Clo( * P-o) Since (Vw)b, S, Co, Mo, Tb and Po are all finite and are independent of n+l'n', therefore, IVi,1j is uniformly bounded and the theorem is proved. A test on the computer has been carried out for Blasius flow. A standard space domain shown in Figure 2 is used. Its dimensions are given by system no.l1 in Table 1. The result shows how the size of time step AT affects the amplification of the local errors. At first, Blasius solution is read into the computer as the initial condition. Then the numerical computation is carried out according to the procedure described in Section (B-4). Time step AT is chosen to be 200. After 600 time steps, the numerical solutions of the local velocity become steady. Actually, the criterion employed to determine the steadiness is e..,o t a, boo S.ox I for all the grid points (ij) in the above space domain. At this stage,'n' is set to be 0, a new quantity Zi~j is defined as follows,

42 zi =tCt; (3.49.a) O Then a -5% error is artifically introduced in U, 12 which takes a numerical value of 0.8115774 before the introduction of the error. Now denoting n the instantaneous local errors by Eci j it yields E'~ _ LN1<-Zz (3.49.b) n where the initial value of ci j are E,12 -= 0.05 x 0.8115774 = -4.0578870 x 10 and 0 C..= 0 i]j for all the space grid points (i,j) except the point (3,12). During the test, the stability condition given by equation (3.45.b) is assured, however, that given by equation (3.45.a) is not imposed. n First Ui j is computed according to FDAEs (3.35) and (3.37) for AT = 200, 250, 300 and 350. Then local error E3,12 is computed from equation (3.49.b) correspondingly. Its numerical results are shown in Figure 3. For cases AT = 200 and 250, E 3 121 becomes smaller and smaller, and finally vanishes as'n' increases indefinitely. This implies, according to the definition n n of e given by equation (3.49.b), that U is uniformly bounded and ij,12 approaches its steady value Z3 12 as'n' increases. When AT = 300, 1c121 first decays, and finally remains constant with a numerical value

43 5.0 Condition of Stability A\r< 266. 4.0 3.0 o = 350. 2.0 o300 AT- 250. 1.0 Ant- 200. 0 O -2.0 -4.0.Error Artificially Introduced at n=O -5.0 Figure 3. Effects of the size of AT on the instability errors.

44 TABLE 1 SYSTEMS OF SPACE DOMAIN sem |no.1 no.2 no.3 no.4 dimens ion ml 4 7 7 3 m2 8 15 15 7 m 8 33838 30 S 20 20 20 22 A* 25.0 25.0 30.0 16.4 AXi** series series series series no.l1 no.2 no.3 no.2 *AY = A for m2 > i > 1, AY = 2A for m > i > m2. ** Series of AXi are listed in Table 2. n of 0.5 x 10-p as'n' increases. Hence, U 12 is still uniformly bounded but deviates from its steady value Z3,12 by an amount of ~ 0.5 x 10-3. In the last case that AT = 350,,I121 is amplified indefinitely as'n' increases. In fact, it becomes so large that the accumulator of the n computer overflows before'n' increases to 50. Consequently, U312 becomes unbounded and is meaningless in the limit as n oo. n+l It is important to note that, for U3'12 to be uniformly bounded, the stability condition given by equation (3.45.a) requires AT < 266. This estimate of the limiting value of AT agrees very well with the results obtained from the above test on the computer. 7. Convergence Analysis and Test If a FDAE is consistent with its corresponding PDE, the proof of its convergence is reduced to the proofs of the uniqueness and the ex

45 TABLE 2 SERIES OF SPACE INTERVAL AXi series no.1 no.2 no.3 A i X AX X AX. X AX. 1 1 1 1- Xi 1 - O. - O. - O. 2 1514.67 1514.67 700 35 700.35 1200. 1200. 3 1965.42 3480.09 814.32 1514.67 1400. 2600. 4 2404.02 5884.11 927.06 2441.73 1600. 4200. 5 2827.83 8711.94 1038.36 3480.09 1800. 6000. 6 3289.12 12001.06 1148.04 4628.13 2000. 8000. 7 3620.64 15621.70 1255.98 5884.11 2200. 10200. 8 3982.35 19604.05 1561.97 7246.08 2400. 12600. 9 1465.86 8711.94 2600. 15200. 10 1622.41 10334.35 2800. 18000. 11 1666.71 12001.06 3000. 21000. 12 1763.37 13764.43 3000. 24000. 13 1857.27 15621.70 3000. 27000. 14 1948.35 17570.05 3000. 30000. 15 2034.00 19604.05 3000. 33000. 16 2121.33 21725.38 3000. 36000. 17 2202.96 23928.54 3000. 39000. 18 2281.23 26209.57 3000. 42000. 19 2355.96 28565.53 3000. 45000. 20 2427.06 30992.59 3000. 48000. 21 2494.41 33487.00 3000. 51000. 22 2557.92 36044.92 3000. 54000. 23 2617.50 38662.42 3000. 57000. 24 2673.03 41335.45 3000. 60000. 25 2724.42 44059.67 3000. 63000. 26 2771.64 46831.51 3000. 66000. 27 2814.57 49646.08 3000. 69000. 28 2853.18 52499.26 3000. 72000. 29 2887.35 55386.61 3000. 75000. 30 2917.11 58303.72 3000. 78000. 31 2942.37 61246.09 3000. 81000. 32 2963.07 64209.16 3000. 84000. 33 2979.21 67188.37 3000. 87000. 34 2990.76 70179.13 3000. 90000. 35 2997. 69 73176.82 3000. 93000. 36 3000.00 76176.82 3000. 96000. 37 3000.00 79176.82 3000. 99000oo. 38 3000.00 82176.82 3000. 102000.

istance of their solutions. For a general initial-value problem, such a theory has not been available yet. However, for a properly posed initial-value problem, if the governing PDEs are linear and homogeneous, it is known (39,41) that the corresponding FDAEs are convergent if and only if they are consistent and stable; here a problem is characterized as "properly posed' if a solution of governing PDE exists and depends uniquely and continuously on the initial values. John (43) has treated the convergence of the solutions of FDAE (3.40). The result that consistency and stability imply convergence also holds in this case. Douglas (44) has studied the sufficient convergence conditions for the explicit FDAE which approximates the following nonlinear parabolic PDE, The convergence conditions for the implicit DAE which corresponds to The convergence conditions for the implicit FDAE which corresponds to the following nonlinear PDE, were investigated by Douglas (45) and Lees (46). In the rest of the section, the sufficient conditions are derived under which the solutions of FDAEs (3o35) (3~36) and (3537) respectively converge to the solutions of equations (3.29), (3.31) and (3.33) as AT, AY and AXi approach zero. It is important to note that equations (3.29), (3.31) and (3.33) are the exact finite difference representations of PDEs (3.26), (3o27) and (3528), respectively.

47 ~n+l n+l Now let Ui,j, V1 ij be the solutions corresponding to the exact n+l n+l finite difference equations (3.29) and (3.533) and Uij and Vi j be the solutions of FDAEs (3.35) and (3.37). Then FDAEs (3.35) and (3.37) are called convergent if A7,-~ U-, -, - o and a't,'y,,~X.~:L\-, TVi-"'L, ---— o AIIAYAXf-> Cj Lpj for all space grid points (i,j) and at any time level'n'. It is important to note that U and V are also the solutions of PDEs (3.6) and (3.2). For a problem of physical significance, there is no loss of generality by the assumption that the solution is a well-behaved function. Therefore, it is assumed that (1) The solutions U and V are uniformly bounded and their partial derivatives with respect to T) X and Y exist and are also bounded uniformly, Particularly, "y4 I 24,Y I L2' (2) The following two limits exist and are uniformly bounded: )y3~ ~ ~ ~ ~~~~~~~h

48 axi- o AX I - -T 1: ~' = I il', 4Y-~c Z(AY) for all possible grid points (n,ij). Based on these assumptions and the assumptions made in the previous section (except that AY and AXi are now allowed to approach zero), the convergence conditions for FDAEs (3.35) and (3.37) are stated in the following theorem. Theorem: The stability conditions given by equations (3.45) are the sufficient conditions for the convergence of FDAEs (3.35) and (3.37). Proof: Rewriting PDE (3.6) in the form, axn (dx, i, r and tV-3E x X UX + U, —U (3 5~) where x_.u _r _ sU UX> yd y, Uyy- yV ) the equations (3.29) and (3.33), which are the exact finite difference representations of PDEs (3.6) and (3.2), can respectively be transformed into the following inequality forms:

49 L A i(V it -F — vlUtK ST + * + ( 2 A- AY il -t,i j) -- [,,I-I 1 2(d]~.) L'l!.I iX - t-, t - I j ~ AYtL ~4(AXL) + K AY(A) + 6 (AY)Y) (3552) whe re'-Y U' -1t (3.53. a)..... - 53+,i AX1 (t;,A=ZliYe- (5.535b) (Ay)_ I> j1-2%-it l~tJ-t (3-53-c) Here it should be remarked that the first assumption was used during the transformation from equations (3.29) and (3.33) to inequalities (3.51) and (3.52) respectively. Similarly, FDAEs (3.35) and (3.37) can be written as follows: where yl ^w- V ('6') Vln_ V., U:....... (3 565a)

50 —, -1 2 j i;(4Y)z (3-56-c) bcd — yl (5.56.b) n+l n+l Def ining Cij and Bij as Jt'+] f-.h+1 14(3 357*a) St,~ t (3 57 b) and subtracting equation (3.54) from inequality (3.51) yield ( L j-IT( + AZJ~41 (ATAA) t.( X K(t())() (.58) where *2 and V1 are the function 4 evaluated at the points A2,(XiTn, ~n ~'n en ~n ~~n n n n n Ui,j, Vi,j, fli,j~ ~ ij ~i,9j) and A1, (Xi, Tn, Ui jVi, ij' n -ij)y respectively. Using the mean values, the following equation is obtained: l t-t =(X)(gj X )1 ( (Vj U; (. l IbI+ f',, —.) ('i) ( V)(j) t- { c'.

51 where all the quantities in the brackets ( ) are evaluated at some point between A1 and A2. It follows, by the definitions given in equations (3.53), (3.56) and (3.57), that VK —4 = xr u- (Uy) (j (Vt) - -K ( 24 (j tF 3 ( 59) (359) Inserting equation (3.59) into (3.58) yields iL, (1 *W) (aJ ( (A(ti)'z + T) ) t + (4-y + Aui I )<I+K (AT+)AX) ((KK) (AT3 ( 6o) n n If Xn and Pn are the lowest upper bounds of a i and Pi j, respectively, then, Max. |,y j —- (3.61.a) and Max. | - Ph (3*61.b) for all the grid points (i,j) under consideration. Furthermore, let Max fn (3.62) for k = 0, 1, 2,..., n.

52 Now if ____-2 __,_ (3063.a) (AM) Y) ( V)| (3y.6.b) then, by the definitions given in equations (3.61) and (3.62), inequality (3.60) can be expressed in the following form: |j +j ((t)a t + K-t gx 0 t r)i ffi,)+(g3 k )(^ whe re at~ -Max./~j ) (3.64.a) X- MIA6 O(iy)n (3-64ab) Since [(Ux)[ takes a value between [hi, jl and.i jl, from equation (3.53oa) and the first assumption made in this section, it may easily be seen that the maximum of |qi;jl is uniformly bounded whether AXi approaches to zero or not. Similarly, from equation (3.56.a) and the second assumption, the maximum of lI jl remains uniformly bounded in the limit as AXi and AY approach zero. When AXi and AY are finite, the conditions (3.63) imply that [Uinjl is uniformly bounded, and Iqi,jl still remains uniformly bounded. Consequently, t-l is a finite positive constants The same is true with a2. Since la i|j has a lowest upper bound anl,hence,

53 Applying the above inequality successively until n 0 O gives it 1 3. l h-Kl'I.',t(.. 6 /. The detailed derivation of this inequality is given in Appendix I. Since'-....,'~ l, iX |'\ |-IG I - L — I- l, - | —'! -,,-~ in the limit as A/, nAX and AY approach zero, inequality (3.66) becomes i X <"-L ) lti -'-: iZ e ) (3.6*7).'r it t; i t -.4 0 - -6 n+i u Next the inequality for Bi,j which was defined in equation (3.57.b) is established. Subtracting equation (3.55) from inequality (3.52) results in rH-l 7' (A'( ** -, <> 0 _ (g.) _ t (' "~ -. + -l _,- dk I.-I. + dY i 1k (4L) + K. (4\f) t ( ) The:n repeated use of this inequality until j - 2 gives' n~~~~t~' ]-4/LA I -- _ K -+ (r-I)A('~ i ~j..t 4AKT K,,_(AY)+K(Ay) where'NAI+TI - 1I ~. s

54 for all the space grid points (i,j) under consideration. Since =V', -V< = NW)- V = 0 and (]-i)Ay Yj Yb therefore, in the limit as AT, AXi and AY approach zero, inequality (3.68) becomes or,lQkaY" o Pntl-YbY;,t, (3.70) CT,AXZ,&Y->O where a~t^b~;,daYioa.oit AbW _N (I371) Applying mathematical induction to inequalities (3.67), (3.70) and equation (3.71), and the use of a' - P-= ~,,.I [iP I= Ma. X!, Vt, I-,l=.(vS- m,j1= ~ readily yields Un+tl= 0 AT, AX~,~A Pn+= Q= at any time level'n' According to equations (3.61) and (3.57),it may therefore be concluded that (r =- IJ * I,LAXl;,AYo{ Ij7-Vc, j ~0

55 for any point (n,i,j). The above proof of the convergence of FDAEs (3.35) and (3X37) is based on the assumptions given by inequalities (3.63). Hence, these assumptions are the sufficient conditions for the convergence of FDAEs (3.35) and (3.37). It should be noted that (U) in inequality (3.63.a) ~ —n nn takes some value between Ui, and Uij, while (V) in inequality (3.63.b) -~n n takes some value between Vi,j and Vij. Since the solution of any convergent FDAE can be made considerably close to the solution of corresponding PDE by choosing AT, AXi and AY sufficiently small, it may be assumed then and Accordingly, the convergence conditions given by inequalities (3.63) may be approximately written as ___ ___ _ (3.72.a) A, (AY1z L |IV, j ~AY\ (3.72.b) which are identical to those given by inequalities (3.45). Therefore, the theorem is proved. Similarly, it may be proved that the stability conditions given by inequalities (3.46) are the approximate sufficient conditions for the convergence of FDAE (3.36). This proof is omitted here.

56 Figure 4 shows how the solutions of FDAEs (3.35) and (3.37) approach the solutions of corresponding PDEs as the space grid size is refined. The solid curve in Figure 4 represents the Blasius solution U versus the similarity parameter Y/NTX. The broken curve no. 1 is the numerical solution Z8,j (for j = 1,2,...,41) which was computed in the previous section and corresponds to the local velocity profile at X = 19604.05. The relatively large deviation of Z8,j from the Blasius solution is due to the local truncation errors. The broken curve no. 2 is the numerical. solution of the local velocity profile at the same location after the size of AXi has been refined. The deviation of this numerical solution from the Blasius solution is much smaller than that of the previous case. It is reasonable to assume that as the grid size is further reduced, the numerical solution of the Blasius flow converges to the Blasius solution. 8. Singular Errors During the computation of the numerical solution for the broken curve no. 2 described in the previous section, a standard space domain shown in Figure 2 is adopted. Its dimensions are specified by system no 2 in Table I. The time step AT is taken to be 200, and the initial values of U and V are assumed to be Blasius solution. The steady state has been reached after 600 time steps. The numerical results of the steady local velocity profiles at five different locations are shown in Figure 5. At X = 55386.61, the numerical solution agrees very well

1.2 1.0 0.8 0 (p 0.6 Blasius Solution 0.4Q Numerical Solution (Ar 3200,AY= 25 and o Ax?, given as series no.1 in Table 2). ~~~0. ~2 Numerical Solution (Ar =200., AY=25. and A Xigiven as series no.2 in Table 2). 0 I 2 3 4 5 6 7 8 9 10 Figure 4. Effects of the size of AXi on the convergence of solutions of FDAEs (3-35) and (3.37) for Blasius flow.

1.2 1.0 0.8 p 0.6 - Blasius Solution o I = 29, X 55,386.61 0.4 i= 15, X 19,604.05 Numerical Solution o i 7, X 5,884.11 0. 2 P v i=3, X =1,514.67 0 2 3 4 5 6 7 8 9 10 Figure 5. Effects of singular errors on numerical solutions of the local velocity profiles in Blasius flow.

59 with Blasius solution. However, as the values of X decreases, the deviation of the numerical solution from the Blasius solution increases. Figure 6 shows that the above conclusion also holds for numerical solution of the local skin friction. The numerical computation of the local skin friction, (6U/6Y)y=o, is based on the following three-point Lagrangian interpolation formula (47): (zy, ( ~ (-3t1n' U 4laJ-U,3) 0 rI 2(AY) Hereafter this formula is used for all numerical computations of the local skin friction and the local temperature gradient. There are two reasons for the deviations mentioned above. Firstly, Blasius solution is not the exact solution of PDEs (3.26) and (3.28) at and near the origin X = O. Secondly, under the current system of the grid size, FDAEs (3.35) and (3.37) are not valid near X = 0 where a singularity exists. From the second reason, it is clear that the existance of singularity causes an unacceptable local truncation error. Complete elimination of these singular errors is impossible. However, by choosing a system of nonuniform grid size as described in Section (B-l) of this chapter, the singular errors due to the singularity at X = 0 can best be compensated. A close examination of Figures 5 and 6 shows that, under the present space grid size the numerical solution of Blasius flow agrees quite well

6.0, _- Blosius Solution o o o NumericolSolution 5.0 - 4.O0 to ot \ \ 0~~~ 3.0 It >~ 2.0,o.O2~~~~~~~.0 65.0 4..0 4.0 ID~o 20 30 4~o4 60 70 8.0 9 0 0 I.O Xxts4 f inl e re 6. Effects of sinB~lar errors on nerical solution of the local skin figureicts Of s ic~nin Blasills flov. fr a.io

with Blasius solution for X > 1.96 x 104.* However, for X 1.9 x 104 the deviation from Blasius solution becomes considerable because of the presence of the singular errors. Therefore, in the examples to be given in Part (A) of Chapter IV, the numerical results in the domain that X < 1.96 x 104 are not used. 9. Local Rounding Errors The procedure of the numerical computation described in Section (B-4) of this chapter has been carried out on IBM 7090. Since IBM 7090 has only thirty-six binary bits in a single word, the aforementioned rounding errors are inevitable during the numerical computation. The only concern is that whether the rate of generation of the local rounding errors is higher than the rate of decay of the local errors or not. If it is, then the numerical solution to be obtained from the computer is meaningless even though FDAEs are stable and convergent. Thus becomes necessary to estimate the upper bounds of the local rounding n+l n+l n+l errors during the numerical computations of Ui j, i,j and Vij according to FDAEs (3-35), (3.36) and (3-3.7), respectively. However, before this estimate can be made, the maximum possible rounding error introduced in a single floating-point arithematic operation must first be studied. In the floating-point computation, each number z is represented *If the main stream velocity is assumed to be 20 fps and air to be the medium of flow, the point X = 1.96 x 104 corresponds to a location about 1.9 inches from the leading edge of the plate.

62 by an ordered pair'a' and'b' such that z = a.2, where'b' is called the exponent and takes an integer value, positive or negative;'a' is called the fractional part and satisfies 0.5 < al < 1. In IBM 7090, the word bits are so allocated that'b' has eight binary bits and'a' has twenty-seven binary bits. This gives an accuracy of eight decimal digits for'a' or, accordingly, for z. Meanwhile'b' is restricted by Ibl < 127. The above restriction on'b' implies that IBM 7090 is only capable to take a number whose absolute value is smaller than 2127 Let fl [B) denote the computer result of the arithematic operation'B' according to the floating-point number system described in the previous paragraph. It is not difficult then to derive the following identities (48): ~ X,t 4- (, 2j) O(l t Er) (373. a) +k t,) t-~ + Y L -X + - - ) t4v( tJl4 ) t XL( litti) (5 3 y7b)

63'/l.rl },-%2 X, (itEr) (,.73 c) (3 73.d) whe re -IErI 2-7 (5.74.a) tt-Z.2B L z7!+)1 -r, t? I (l i ) 9 H f, ~ i 1 ) (574 ) (I_ z"l t~i-< I ~' ( I tZ2-, -.to r - Since the number 2-27 is much smaller than unity, all the second and higher order power terms of 2 may be neglected:in all binomial expansions in inequalities (3.74) Consequently, Yij and Can be expressed as follows: txj = ( C +aj -l) i(3 *75.a) "tt- ((3-7 -a' ~7-(t44r rfor - fo.r 2 Similarly, it is easy to verify that (ct F'-y:r [)(bt/'rFr)-':b i+(tr k + ")Fr (5.76) CA t 61, I/ (.t

64 where a" and b" are arbitrary constants. The upper bounds of the local rounding errors introduced during the n+l n+l n+l numerical computations of Ui n i, and Vi may now be estimated. Let n+l n+l n+l Ui,j, Gi,j and Vi,j be the exact solutions of FDAEs (3.35), (3.36) and;n+l (3-37), respectively, and the corresponding computer solutions be Ui j, ~n+l ~n+l 8i,j and Vij. Therefore, according to FDAEs (3.35), (3.36) and (3.37), 1J- j 44 U tj /\\ (-bT j: I tt, (Uo\M 4) q J ~I ~. VX-A /L-,j 2; /,.t -:+(-'u,+.. (3,77):.,{ Z.'~l (iT) -... F r;: (3.78) ~it ri -1 (I( I U, (3U*79) In above equations, we have assumed that the numerical values of AT, AXi, AY.Pr, (iUo/6T)i and (Uo0Uo/6X)ni can exactly be read into the computer and no rounding for these numbers is needed. n+l n+l n+l The local rounding errors of Ui e and Vi. are defined as 7_ n+l n+l;n+lll n+l a nn.l n+l the differences (Uij j) l. i.) and (Vi,j - Vi,) respectively provided that Uin n and V,j for all the space grid points (i,j) can exactly be represented on the computer. In other words, assuming L,4 =iJ > L L,J - ]

for all the space grid points (i,j), it is desired to study the upper':z~p~l n+l I I -+l I and V..- Vi ~ l bounds of Uij - Ui|j - 3ij anij i According to equations (3.77), (3.35) and the floating-point arithematic operation rules given by equations (3.73), (3.75) and (3.76),it may be estimated that + -27r~LIJ IfUqn In-frr~Uof'.~ UA-; -,A,- i |+< (2 ) + +(-titi + uj.l t6i ( AY I(AY)2 (3.80) Similarly, Lj- L, L.J I~ uk -" _____\_'__y,4 (az lj 4 &')It. vI +Pr (A~)z (3.81). (AY)t. V64'Vit < (Z ) &-f w e (Y) t (( knWi-4~)D + (5-t)U 7)H - 1 +H hB H: (3-82) ~n~~~~Y where Hi: is defined as the quantity in the brackets [ ] on the righthand side of equation (3.80).

CHAPTER IV NUMERICAL EXAMPLES A. OSCILLATIONS IN BLASIUS FLOW 1. Statement of the Problem Here a fluctuating boundary layer on a stationary semi-infinite flat plate is considered. The fluctuations in the boundary layer are introduced by an unsteady main stream which oscillates periodically about a constant mean value. Let this mean value of the main stream velocity be u, and the corresponding amplitude and the frequency of the oscillation be Auo and w, respectively. Hence, the free stream velocity uo of the problem takes the following form: CL( -Ll t All,, S | n (at) (4.1) Selecting u, as the characteristic velocity of the problem, equation (4.1) may be rearranged in the dimensionless form, u-=- +Au0 S\n(Wi) (4.2) where AUo and W, defined as Auo/U, and w/(u 2/), respectively, are the corresponding dimensionless amplitude and frequency of the free stream oscillation. With UO specified in equation (4.2), the governing FDAEs of the problem can easily be obtained from equations (3.35) and (3 -37). The boundary conditions are given by equations (3.38.a), (3.38.c), (3.38.d) 66

67 and (3.38.f). Since neither blowing nor suction is present in the flow, the value of Vw in equation (3.38.d) is taken to be zero. The initial conditions are specified in equations (3.39.a) and (3.39.b) where Us and Vs are assumed to be the steady Blasius solution. The analytical solutions of the above problem were studied by many authors In general, they can be classified into three groups, namely, high-frequency, intermediate-frequency and low-frequency solutions, according to the characteristic frequency parameter X.W. For high-frequency case (X.W > 10), the problem was investigated by Li.n (31) with a method described in Part (F) of Chapter II. According to Lin's result. the amplitude and the phase lag of the local osc iilations depend only on the combined variable Y\'W/2 and is independent of X. Hence., a single curve is sufficient to represent these quantities for all frequencies. In the intermediate-frequency range (10 > X.W. > 1), Hill and Sterin ng ( 49) studied the probliem analytically and experimentally. The analytical method which they used is actually the method of successive approximaticns described in Part (D) of Chapter II. Lin's (31):results were taken as the solutions of the first order approximation, and the solutions of the second order approximation were computed. They found, in addLtion to variable YNW/2, that the amplitude and the phase lag of the local oscillations also vary with parameter X Ow X

68 For low-frequency case (1 > X.W), the analytical investigations were carried out by Lighthill (35) and Nickerson (50). Lighthill linearized the governing PDEs by using the series solution method which was discussed in Part (C) of Chapter II. Only first two terms in the series expansions of the local velocity were considered, hence, the validity of the results are restricted to the small oscillations. The small oscillation in Blasius flow was also studied by Ghosh (51) for all frequency ranges. Again, only first two terms in the series expansions of the local velocity were calculated. As pointed out by Schlichting (22), the first two terms in the series solution are not sufficient to reveal the net effect of the free stream oscillation on the local velocities in the boundary layer. 2. Numerical Solution for High-Frequency Case In this example, a standard space domain as shown in Figure 2 is adopted. The dimension of the space domain is specified by system no. 2 in Table 1. According to this specification, the maximum value of X is about 8.22 x 104. Again the numerical solutions in the domain that X < 1.96 x 104 are not used, because of the significance of the singular errors introduced in this domain. The dimensionless frequency W and amplitude AUo are chosen to be t x 10-4* and 0.05, respectively. *If the mean value of the main stream velocity, u,is chosen to be 20 fps and air is the medium of the flow,W = n x 104 corresponds to a frequency of 125 cps.

69 Therefore, the characteristic frequency parameter X.W in this example varies from 6.15 to 25.8 which covers most high-frequency range and the upper half of the intermediate-frequency range. The time step AT is selected to be 200, According to the frequency specified in the previous paragraph, 100 time steps are needed to complete each cycle of the oscillation. The flow becomes approximately steady periodic after five cycles of oscillation from the steady Blasius flow, or equivalently, after 500 time steps. In fact, — 400+1 30t _ i6.;X is satisfied for all grid points (i,j) in space domain and for any n between 0 and 100. The above numerical computations are carried out on IBM 7090 according to the program given in Appendix II-A. The total machine execution time is about twenty-five minutes. The numerical solutions shown in Figures 7 and 8 are drawn during the computation of the local velocities at the sixth cycle. In order to draw the necessary information about the local velocity at every degree of the phase angle, the time step AT used in this cycle is changed to 55.555556 and the computer program is modified slightly. Figure 7 shows the numerical solutions of phase angle 8 of the local velocity component U for X.W = 25.8 and 9.0. A comparison is made with Lin's (51) high-frequency solutions. The numerical solutions agree well with Lin's solutions for X.W = 25.8. For X.W = 9.0 which corresponds to a frequency a little lower than the lower limit of the high-frequency

50 High-Freq. Solution by C.C.Lin (31) 40- | XW = 9.0 }Numerical Solution (AUo=0.05) o XW = 25.8 30 20 0 \ -10 0 I 2 3 4 5 6 7 8 9 Figure 7. Phase of high-frequency oscillations in Blasius flow./ Figure 7. Phase of high-frequency oscillations in Blasius flow.

1.2 1.0 0.8 0.6 3 0.4 High -Freq. Solution by C.C. Lin (31) v XW =9.0 Numerical Solution 0.2 0 XW =25.8 (AU =.05) 0 I 2 3 4 5 6 7 8 9 Y, 1/W/2 Figure 8. Amplitude of high-frequency oscillations in Blasius flow.

72 region, a slight deviation appears between the solutions, as expected. The same conclusion can be extended to the local amplitude ratio AU/ AUo according to Figure 8. 3. Numerical Solution for Intermediate-Frequency Case The dimension of the space domain used in this example is specified by system no. 3 in Table 1. The numerical results bounded by the lines X = 2.1 x 10 and X = 10.2 x 104 will later be used for discussion. The frequency W and the amplitude AUo are chosen to be 5.8177641 x 10-5* and 0.05, respectively. The corresponding characteristic frequency parameter X.W varies from 1.22 to 5.94 which covers the lower half of the intermediatefrequency region. The time step AT is chosen to be 300, 360 time steps are needed for one cycle of the oscillation. In the fifth cycle, the flow is almost steady periodic, in fact, 1440+nV L.;-< 4. x10O for all grid points (i,j) in the space domain and for any n between 0 and 360. The numerical results shown in Figures 9 through 16 are drawn during the coimputations of U and V in this cycle. From the first cycle to the fourth cycle, the program given in Appendix II-A is used. This *If up is taken to be 20 fps and air is the medium of the flow, W = 5.817764 x 10-5 corresponds to a frequency of 23.2 cps.

73 program is slightly modified in the fifth cycle in order to print out all the information needed for Figures 9 through 16. The total machine execution time is about ninety minutes. Figures 9 and 10 show the numerical solutions of the phase angle 8 and the amplitude ratio AU/AUo, respectively, for the characteristic frequency parameter X.W = 5.06. A comparison is made with Hill and Stenning's (49) intermediate-frequency solutions and experimental measurements. Both the numerical and the intermediate-frequency solutions slightly deviate from experimental measurements. The numerical solutions together with the intermediate-frequency solutions and the experimental measurements of S and AU/AUo for X.W 2.45 are shown in the Figures 11 and 12, respectively. As far as amplitude ratio AU/AUo is concerned, Figure 12 shows the agreement of the numerical and intermediate-frequency solutions with experimental data. However, for phase angle (, Figure 11 shows numerical solutions agree with experimental data better than the intermediate-frequency solutions. For smaller values of the characteristic frequency parameter, say X.W = 1.40, Figure 13 shows that the deviation of the intermediatefrequency solution from the experimental data becomes more severe,'while the agreement with the numerical solution is still excellent. Figure 14 shows that the numerical solutions of the maximum phase lag have the same trend as the experimental measurements for the complete intermediate-frequency range, while the intermediate-frequency

50 Numerical Solution (A Uo=0.05, XW5.06) 40 -Interrm.-Freq. Solution By Hill and Stenning (49) o o o Measurement (AUo=O.ll XW 4.98) 30 oX I 20 0 0 0 0 -10 I I I! I I 0 2 3 4 5 6 7 8 9 Y 1W/2 Figure 9. Phase of intermediate-frequency oscillations in Blasius flow (xw = 5.6).

1.2 0~~ 1.0 0.8 0 O 0.6 < 0.4 0 - Numerical Solution (Auo=0.05 XW=5.06 Interm.- Freq. Solution By Hill and Stenning 0.2 0 0 0 oMeasurement J (Au O.09XW =4.98) 0 2 3 45 6 7 8 9 Figure 10. Amplitude of intermediate-frequency oscillations in Blasius flow (XW = 5.06).

50.. Numerical Solution (AUo=0.05, XW 2.45) 40 - Interm.-Freq. Solution By Hill and Stenning (49) o o o o Measurement (AUo=0.11, XW=2.49) 30 0 o 20 0 o -10 Lr 1 PI Ie r I I I I I! 0 1 2 3 4 5 6 7 8 9 Y ~~W/2 Figure 11. Phase of intermediate-frequency oscillations in Blasius flow (XW = 2.5).

1.2 1.0 0.8 0.6 0' 0.41,/ -Numerical Solution (AUo =0.05, XW = 2.45) --- --- Interm.- Freq. Solution By Hilland Stenning (49) 0.2 o o o0 Measurement (AuO =0.15,xw=2.65) o I I_ I I I I I I 0 I 2 3 4 5 6 7 8 9 YIW/2 Figure 12. Amplitude of intermediate-frequency oscillations in Blasius flow (XW = 2.45).

50 Numerical Solution (AUo= 0.05, XW= 1.40) 40 Interm.- Freq. Solution } By Hill and Stenning (49) o o o Measurement J(AUo o. 10, XW 1.48) 30 o 20 0 -10 0 1 2 3 4 5 6 7 8 9 Y W /2 Figure 13. Phase of intermediate-frequency oscillations in Blasius flow (XW = 1.40).

0 0 (31) Lin's High-Freq. Solution-' ~~~0 -4 o \ro -8 00 0-=,~-12 - - Interm.-Freq. Solution By i and Snning49) }By Hill and Stenning(49) -16 c) o o Measurement Numerical Solution I I I I I I,! 0 I 2 3 4 5 6 7 8 XW Figure 14. Maximum lagging phase angle of intermediate-frequency oscillations in Blasius flow.

1.2 1.0 0.8 1~ 0.6 0.4 -- Blasius Solution 0 XW 25.8 }Numerical Solution (AUo=0.05) 0.2 a r X W 5.06 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Y/TX Figure 15. Time-averaged velocity profiles of high and intermeidate-frequency oscillations in Blasius flow.

1.2 (Y T_)(-U) U o1.08- =,.- A 3.90, 0.720 0.8 / I 9.00, 1.000 (Free Stream) 0.6 —- 1.94, 0.389 0.4 t? ~- 0.65, 0.131 0.2 sI -~'0.4 0.0 -0.6 50 - 1.5 0 20 40 60 80 1100 120 140 160 180 200 220 260 280 300 32 90 (wr)o Figure 16. Local velocity fluctuations of intermediate-frequency oscillations in Blasius flow (W = 5.818 x 10-5).

82 solutions start to deviate from them at X.W = 3.0. Since Hill and Stenning's intermediate-frequency solution is merely an improved Lin's high-frequency solution by a second order successive approximation, the deviation from the corresponding experimental data at small values of the characteristic frequency parameter is not surprising. It is worth noting that the existing literature (18,31) observes the net streaming effects introduced due to the presence of the oscillation in Blasius flow is of order of (AUo). Since AUo used here is equal to 0.05, therefore, the difference between the time-averaged local velocity U and the steady Blasius velocity profile must be negligibly small. The numerical solution shown in Figure 15 confirms this point. Furthermore, Figure 16 shows that the instantaneous local velocity computed by the numerical method oscillates sinusoidally about its timeaveraged value. B. IMPULSIVE START OF WEDGE FLOW 1. Numerical Solution of the Transient Velocity Profile The boundary-layer flow along a wedge of 450 half-angle is considered. The wedge is initially at rest relative to a stationary fluid and suddenly moves with a constant speed along its center-line. After the motion has started, the free stream velocity uo relative to a moving coordinate system described in Part (A) of Chapter III takes the steady value,

where xc is an arbitrary characteristic length and uoo is the velocity of the fluid particle on the center-line of the wedge at distance Xc from the stagna ion point of the wedge as shown in Figure 17. Figure 17. Wedge Flow Figure 17. Wedge Flow Selecting xc = </uc and arbitrarily setting u. = uc/45, the dimensionless free stream velocity becomes vfuoX 4-5 (4.4) The governing FDAEs and the corresponding boundary and initial conditions of the problem can be obtained in the same way as described in Section (A-1) of this chapter, provided that Us and Vs are replaced by Uo and 0, respectively. A standard space domain is adopted, its dimensions are specified by system no.4 in Table 1. The time step AT is chosen to be 50, 100 and 110 for n = 1 to 100, 101 to 800 and 801 to 1050, respectively, At n = 1050, the flow is considered to be steady, in fact,

84 10u50 1041 LT U <.4 ()4 is satisfied for all the grid points (i,j) in the space domain. The numerical computation of this example is carried out according to the program given in Appendix II-B. The total machine execution time is about forty-five minutes. The numerical solutions of the transient local velocity profile and local skin friction are shown in Figures 18 and 19, respectively. Neither analytical nor experimental solutions are available for these profiles. It is worth noting that the flow becomes fully developed in 4 a very short period, say T = 10.25 x 10.* Figure 20 shows the numerical solutions of the fully developed velocity profile at two different locations. For comparison, Hartree's (26) analytical solution versus the similarity parameter E is also plotted in this figure. The agreement is excellent. 2. Numerical Solution of the Transient Temperature Profile In this example, the velocity field is assumed to be fully developed, the temperature of the fluid and the wedge is assumed to be uniform and equal to zero initially. The transient temperature profiles due to sudden raise of the surface temperature of the wedge from 0 to 1 will be studied. *If air is considered as the medium of the flow and the characteristic velocity u is assumed to be 10 fps, then the dimensionless time T = 10.25 x 10 corresponds to a physical time of 0.164 second.

1.0 0.8 04,-0.6 T= 0.5x10 T- 2.5x 10 4 0 u~ 0.4 ~ — 7.5~~= 5. x 104 ~~OD~~~~~~~~ 0 6 12 18 24 30 36 42 Development prof4id Y/ 32.8 Figure 18. Development of the local velocity profile in wedge flow.

8r — 7 o 50 =0.5x |04 r 2.5x 104:7.5x 104:4 rO10.25x104 (Steady Solution) 3 22 4 X 10-4 5 6 Figure 19. Development of the local skin friction in wedge flow.

1.2 1.0 0.8 / Ueb I 450 (Corresponds to / = 1/3 ) 0 0.6 a4 - Hartree's (26) Solution for P= I/3 o i= 30, X= 58,303.72l Solution 0.2 o i= 15, X = 19,604.05 0 I 2 3 4 5 r =( Ua/Uc)2Il)/2/2] /2x 1/2 (/- ) y Figure 20. Numerical solutions of the fully developed velocity profile in wedge flow.

88 The governing FDAE is given in the equation (3.36) where Prandtl number Pr is chosen to be 0573. The boundary conditions are specified by equations (3.38.b), (3.38.e) and (3.38.g) where 00 = 0 and $w = 1. The initial condition is given by equation (3.39.c) where Gs = O. The velocity components U j and Vi in equation (3.36) are assumed to be the steady numerical solutions of the previous example. The space domain is chosen to be the same as that used in the previous example. The numerical computation of 0 has been carried out for 1400 time steps. For the first 200 time steps, AT is chosen to be 50; while for the last 1200 time steps, AT is equal to 88. At n - 1400, the temperature profile is considered to be fully developed and the instantaneous local temperature satisfies the following relation: L', L < for all the space grid points (i,j). The computer program given in Appendix II-C is used for the numerical computation of this example, Approximately a total of thirty-nine minute machine time is required for the entire calculation. Figures 21 and 22 show the development of the local temperature profile and the local Nusselt number Nu which is defined as ow - eo In this case, no analytical solution is available for comparison.

1.0 t~l.56xl0 4(Steady Solution) r= 8.04 x 104 0.8 r=4.52 x 104,r 1.00 Q 104 C\j - 0.6 50 rr) 00 LO 0.4 Pr =0.73 0.2 0 -0.2 III 0 6 12 18 24 30 36 42 Y/32.8 Figuire 21. Development of the local temperature profile in wedge flow.

7 U00 /~45O UaI 6 N Pr = 0.7 3 5 r -T1.00 x 104 4 O 0 x I I so z 3: 4.52 x 104 1 1.56 x 104 (Steady Solution) I L I l I I 2 3 4 5 6 I x l0 Figure 22. Development of the local Nusselt number in wedge flow.

CHAPTER V CONCLUSION In the preceeding chapters, a finite difference method for computing velocity and temperature profiles of an unsteady, incompressible, laminar boundary layer around a two-dimensional cylinder of arbitrary cross section is developed. Blowing or suction may be present on the wall of the cylinder. The thermal boundary condition at the wall may be either specified wall temperature or specified heat flux at the wall. Two examples, one oscillation in Blasius flow and the other impulsive start of wedge flow are given; the numerical results are compared with the existing analytical and experimental results. From these, the following conclusions are drawn. (1) The method can be used to compute, with high accuracy, the velocity and temperature profiles of an unsteady boundary-layer flow of the foregoing category up to separation point. (2) The method can easily be extended to determine the velocity and temperature profiles at the separation point, where a temporary reverse flow may occur, provided that a forward FDA, instead of a backward FDA, is used for the derivative in the second term on the left-hand side of equation (3.26). (3) At and near a singular point, the numerical solutions given by the method are doubtful because of the singular errors. 91

92 (4) For periodic flow problems, if the period of oscillation is very long or the frequency is very low, a large number of time steps are needed to complete one cycle of oscillation because of the limitation on AT imposed by the stability criterion. Therefore, it may not be advantageous to study oscillating flows of very low frequency by the present method. (5) For starting-ending flow problems, which are difficult to investigate by analytical methods, the present method is indispensable. The method can also be used to study the starting-ending flow problems for large values of time. In this case, the problem becomes pure'ending' type and its solution approaches the steady solution of the corresponding problem.

CHAPTER VI HYDRODYNAMIC STABILITY An attempt is made to simulate Schubauer and Skramstad's experiment by the present numerical method. Disturbances are introduced at a distance Xa away from the leading edge, and they are investigated as they travel downstream. Since the thickness of boundary layer is quite uniform at the downstream of Blasius flow,a rectangular space domain bounded by straight lines X = Xa, Xb and Y = O, Yb as shown in Figure 23 is chosen. Space grids AX and AY are selected to be uniform. For small disturbances, it is assumed that the disturbed velocity at Xa is equal to the velocity of the base flow plus the disturbances introduced by the vibration of a ribbon at Xa. Thus Vwher x (VBa) s F(Y) Sti(Wh) (6.1) where subscript'B' denotes the Blasius solution, function F(Y) the effect of disturbances in the Y-direction and W the dimensionless frequency of the vibration of ribbon. It readily follows from the equation of continuity that Ug;9(=3+ X S ) ( 6.2a) and 95 I ( (6. b) 93

Y Yb AO a__ __ i,j) ~AX~X Figure 23. Rectangular space domain.

95 where C is a constant proportional to the width of the vibrating ribbon. Following Schubauer and Skramstad (52) the distribution function F(Y) and its derivative dF/dy are assumed to have the forms as indicated in Figure 24. This figure is used as boundary condition at X = Xa. The rest of boundary conditions are Y = 0: U = V = 0 (6.3a) Y = Yb: U = 1 (6.3b) where Yb is chosen to be reasonably greater than the maximum thickness of the boundary layer in the region considered. Blasius solution is used as initial condition. According to the boundary and initial conditions mentioned above, the disturbances at downstream of Xa can be calculated by the present numerical method. Two cases are considered here. W is chosen to be 8.0 x 10-5 in both cases, For the first case, X = X1.95 x 1059 Xb = 2 65 x 105, 3 ~~~~~~a b~~~~~~ Yb 750 x 10 AX = 5000, AY = 100, AT = 785.39816; and the distribution of AU at Xa is given in Table 3. According to References (52,53,54)) the foregoing case is in the unstable region. However, the numerical results drawn from the computation of U and V in the eleventh cycle of oscillation show, as indicated by Table 4, that the disturbances are rapidly damped out when they travel downstream. For

96 TABLE 3 DISTRIBUTION OF DISTURBANCES AT Xa AU j AU AU j AU 1.000 20 -.012 39 -.011 58 -. og09 2.010 21 -.013 40 -.011 59 -.oo009 3.020 22 -.014 41 -.011 60 -.009 4.030 23 -.015 42 -.011 61 -.009 5.040 24 -.014 43 -.011 62 -.009 6.050 25 -.013 44 -.011 63 -.009 7.045 26 -.012 45 -.011 64 -.009 8.040 27 -.012 46 -.010 65 -.009 9.035 28 -.012 47 -.010 66 -.008 10.030 29 -.012 48 -.010 67 -.007 11.025 30 -.012 49 -.010 68 -.oo6 12.020 31 -.012 50 -.010 69 -.005 13.32 -.012 51 -.010 70 -.004 14.010 33 -.012 52 -.010 71 -.003 15.005 34 -.012 53 -.010 72 -.002 16.000 35 -.012 54 -.010 73 -.001 17 -.005 36 -.011 55 -.010 74.000 18 -. 008 37 -.011 56 -.009 75.000 19 -.010 38 -.011 57 -.009 76.000 = (j-1)AY. Y A = Distance from ribbon to the surface of the flat plate F(Y) - -- -- ~dF l I4 v~ dY 0~~~~~ ~~F, dF/dY Figure 24. Function F(Y) and dF/dY.

9'( TABLE 4 DAMP OF DISTURBANCES INTRODUCED AT X = 1.95 x 105 Y AU AU I;AU I'AU' =1.95xlO X=2.1540x105 X=2.40xlo5 X=2.65x1 05 100 1.00 x 102 1.109 x 10-4 1.15 x 10-5 4.50 x 10-6 700 4.00 x 10-2 1.05 x 10 -3 4.55 x 10-4 1.60 x 10-4 the second case, Xa 4 10 X 00 x 105, Yb = 9.00 x 103 AX = 5000, AY = 120, AT = 785.39816; and the same distribution of AU at X used in the first case is adopted here. According to References (52,55,54) the flow should be unstable. However the numerical results corresponding to the fourteenth cycle of oscillations (Table 5) indicate the stability of the flow. TABLE 5 DAMP OF DISTURBANCES INTRODUCED AT X = 4.00 x 10 Y Au I AUJ JAUII x=4.OOxl05 X=4.30x105 x=4.65xo105 120 1.00 x 10-2 0.30 x 10-5 0.23 x 10-5 600 5.00 x 10-2 1.20 x 10-4 4.45 x 10-5 Since Prandtl's approximations are certainly valid in describing the boundary layer flows, why the boundary layer equations yield stability

98 for apparently unstable cases is an interesting question to answer. To explore the point further, the classical small perturbation theory of stability is applied to the boundary layer formulation Au +y~ > o, /aa=D CL,(6.4) a no ___/L~~~ z (6.5) Thus introducing U z u Z.AU VL,ntJ (6.6) into equations (6.4) and (6.5), and neglecting the higher order terms gives da= 1, ba'e O) (6.7) -u" (6-7) -Iat Q - pkV peAilIaron/ (6.8) clearly, a base flow that satisfies equation (6.7) and the boundary conditions (6.9) is not possible. Therefore equation (6.8) cannot be made a characteristicvalue problem whose first characteristic-value is the critical Reynolds number. Thus the comparison of the inviscous Orr-Sommerfeld equations based on the Navier-Stokes and boundary layer equations and expressed in terms of the usual nomenclature

99 it(u9 -a ( a- ),(6.10) reveals the importance of the missing terms At x c3 f46 (6.12) in the latter equation. These terms are related to the terms / - 9v'/ (6.13) of the y-momentum. It may then be concluded that, although the vertical velocity component is a small quantity of order 1/Re, it produces a momentum transfer across the boundary layer where both the disturbance quantities and the mean flow properties vary rapidly. Thus the net effect on the transport process turns out to be much larger than the magnitude of the vertical velocity which produces the process.

APPE&C'LX I DERIVATION OF INEQUALITY (3.66) Rewriting inequality (3 65) in the fcrm, CJ \ t -F; &ij:t t ( l 4i Ai t-k iItrK and rearranging it by the help of equation (3.62) gives ~c —e ( 4~ f -f r _ti i + + I ) for any k < n, Applying inequality (I.1) successively for k = n,n-l,...,O, yields t. he n.' A t. rtAK i[ i t * -,! _ t - 4 r -:; K" "? tf i.~ i'. or 1 -1\ t i, ~: \ y - (T -D.nM.e the time doma n i is fi~nite, fc'r a~n value of n It follcws then (, t -- +, 0I 0,v ( I.4) 100

101 Finally, by the use of inequalities (I.3) and (I.4), inequality (1I.2) may be rearranged to give hich is th e statement of inequality (tK.66).) which is the statement of inequality (3.66)

APPENDIX II COMPUTER PROGRAMS 102

103 A. CALCULATION OF VELOCITY OF OSCILLATION IN BLASIUS FLOW $ COMPILE MAD, PUNCH OBJECT READ FORMAT VA,MA,MB,NlN2,T,T2T3,T4,Y1,Y2DTE,WDX(2)...UX(38) VECTOR VALUES VA=$8I5/5F12.1/(10F6.2)*$ PAUSE NO. 1 BACK SET HIGH DENSITY TAPE 9 REWIND TAPE 9 READ BCD TAPE 9,VR,TN(l)...TN(7) VECTOR VALUES VR=$7C6*S WHENEVER TN(3).NE.$RCDEFG$ PRINT COMMENT $WRONG TAPE$ PRINT ON LINE FORMAT WRONG VECTOR VALUES WRONG=$H*WRONG TAPE MOUNTED PLEASE CORRECT**$ SKIP6. UNLOAD TAPE 9 P~USE NO. 77777 TRANSFER rO BACK END OF CONDITIONAL EXECUTE SKIP.(TlT2,9) READ HINARY TAPE 9,N,U(l,1)...U(38,41),V(,1)...V(38,41) AA( 1)=0. AB( 1)=O. NA N=N+1 SINA=SIN. (WN*DT ) COSA=COS. ( W* (N- 1 ) *OT) UDU=E*W*COSA oY=Y1 DYP=OY.P.2 DYU=2./DY DYM=2.*DY DYP2=2. /DYP THROUGH S1,FOR J=2,1,J.G.41 S1 AA(J)=1.+E*SINA THROUGH S2,FOR I=21,iI.G.6 THROUGH S21,FOR J=2,lJ.G.20 WHENEVER.ABS.V ( I,J).G.DYD.OR.DT.G.1./(U( I,J)/DX( I )+DYP2) 1.OR.U(I,J).L.O. PRINT RESULTS I J,N,U(tI,J),V(I,J),DT,DX(I),DY TRANSFER TO SE OTHERWISE AB(J)=U(I,J)-DT*(U(I,J)*(U(I,J)-U(I-1,J))/DX(I)+V(I,J)*( 1 U( I,J+I)-U( IJ-1) )/DYM-( U( I,J+1)-2.*U( I,J)+U( I,J-1 ))/DYP) 2 +DT*UDU V(I,J)=V(I, J-1)-(AB(J-1)- AA(J-1)+AB(J)-AA(J) )/DYD/DX(I) S21 END OF CONDITIONAL THROUGH S2,FOR J=2,1,J.6.20 U(I-1,J)=AA(J) S2 AA(J)=AB(J) THROUGH S3,FOR I=7,1,I.G.15 THROUGH S31,FOR J=2,1,J.G.40 WHENEVER.ABS.V(IJ).G.DYD.OR.DT.G.1./U(I,J)/DX( I )+DYP2) 1.OR.U(I,J).L.O. PRINT RESULTS I,J,N,U(I, J),V( I,J),DTDXI),DY TRANSFER TO SE OTHERWISE

AB(J)=U(I,J)-DT*(U(IJ)*(U(I,J)-U(I-,J))/DX(I +V(I,J)*( 1 U(I,J+1)-U(I,J-1) )/DYM-(U(I,J+1)-2.*U(I,J)+U(I,J-))/DYP) 2 +DT*UDU V(I,J)=V(I,J-L)-(AB(J-1)-AA(J-1)+AB(J)-AA(J))}/YO/DX(I) S31 END OF CONDITIONAL THROUGH S3,FOR J=2,1.J.G.40 l)(I-1,J)=AA(J) S3 AA(J)=AB(J) DY=Y2 DYP=DY.P. 2 DYD=Z./DY DYM=2.*DY DYP2=2./OYP THROUGH S4,FOR J=:1,1J.G.21 AC(J)=AA(2*J-1) S4 U(15,J)=U(15,2*J-1j THROUGH S5,FOR J=2291,J.G.40 AC(J)=AA(41) S5 U(15,J)=U(15,21) I=16 THROUGH S6, FOR J=2,1,J.G.40 WHENEVER.ABS.V( I J).G.DYD.OR.DT.G.1./(U( I )/DX( I )+DYP2) I.OR.U(I,J).L.O. PRINT RESULTS ItJ,N,U(I.J),V(ItJ),DTtDX(I.),DY TRANSFER TO SE OTHERWISE AB(J)=U(I,J)-DT*(U(I,J)*(U(I,J)-U(I-1,J))/DX(I)+V(IJ)*( 1 U I,J+I)-U(I,J-1))/DYM-(U(IJ+I)-2.*U(I,J)+U(I,J-1))/DYP) 2 +DT*UDU V(I,J)=-V(IJ-1)-A( AJ-1 )-AC(J-1)+AB ( J)-AC(J))/YD/DX I) S6 END OF CONDITIONAL THROUGH S7,FOR J=2,1,J.G.40 U(I-1,J)=AA(J) S7 AA(J)=AB(J) THROUGH S8,FOR I=1i71,I.G.38 THROUGH S81,FOR J=2,1,J.G.4O WHENEVER.ABS.V(I,J).G.DYD.OR.DT.G.1./(U(I,J)/DX(I)+DYP2) I.OR.U(I,J).L.O. PRINT RESULTS IJ,NU(I,J),V(IJ),DTOX(I),DY TRANSFER TO SE OTHERWISE AB(J)=U(I,J)-DT*(U(I,J)*(U(I,J)-U(I-l,J))/DX(I)+V(I,J)*( 1 U( I,J+1)-Ut(I,J-1))/DYM-(U(I,J+I)-2.*U(IJ)+U(IJ-1))/DYP) 2 fDT*UDU V(I,J)=V(I,J-1)-(AB(J-1)-AA(J-l)+AB(J)-AA(J))/DYD/DX(I} S81 END OF CONDITIONAL THROUGH S8,FOR J=2,I,J.G.40 U(I-1,J)=AA(J) S8 AA(J)=AB(J) THROUGH S9,FOR J=2,1,J.G.40 S9 U(38,J)=AA(J) THROUGH SIO,FOR I=2,1,I.G.6 S10 U(I,21)=AA(41) THROUGH S11,FOR I-7,1,I.G.38 Sll U(I,41)=AA(41) WHENEVER N.E.N1 PRINT RESULTS N,U(1,1)...U (38,41 ) WHENEVER MA.E.O EXECUTE SKIP.(T3,T4,9) OTHERWISE

105 HACKSPACE RECORD OF TAPE 9 END OF CONDITIONAL WRITE HINARY TAPE 9,N,U(1,1)...U(38,41),V(1,1)...V(38,41) END OF FILE TAPE 9 TRANSFER'ro SE OR WHENEVER N/NZ.E.MB TRANSFER TO NA OR WHENEVER N/N2.G.MB MB=N/N2 PRINT RESULTS N,U25,1)...U(27,41),U(36,1)...U(38,41) TRANSFER TO NA END OF CONDITIONAL DIMENSION U(1558,VB),V(1558,VB) AA41 ),AB(41)AC(41),DX(38)TN7) VECTOR VALUES VB=2,1,41 INTEGER I,J,N,N1,N2,TN,MA,MB,T1,T2, 3,T4 SE UNLOAD TAPE 9 END OF PROGRAM

106 B. CALCULATION OF VELOCITY OF IMPULSIVE START OF WEDGE FLOW $ CCMPILE MAC, PUNCH OBJECT READ FORMAT VA,MA,MB,NI,N2,TT,T2,T3,T4,YltY2,DTDX(2).. DX(30) VECTCR VALUES VA=$8I5,3F12.1/(10F6.2)*$ PAUSE NO. 1 BACK SET kIGH DENSITY TAPE 9 REWINC TAPE g READ BCD TAPE 9,VR,TN(1)...TN(7) VECTCR VALUES VR=$7C6*$ WHENEVER TN(3).NE.$BCDEFG$ PRINT COMMENT $WRONG TAPE$ PRINT ON LINE FORMAT WRONG VECTOR VALUES WRONG=$H*WRONG TAPE MOUNTEC PLEASE CORRECT**s SKIP6. UNLOAC TAPE 9 PAUSE NO. 77777 TRANSFER TO BACK END CF CONDITIONAL EXECUTE SKIP.(TI,T2,9) READ BINARY TAPE 9,N,U(1,1)...U(30,45),V(1,1)...V130,45), 1 UDU(2)...UDU(30) N A h=N+1 ZERO.( AA(' 1)..;'AA(22)tAB( 1i)) DY=Y1 CYP=CY.P.2 CYD=2./DY DYM=2.*DY CYP2=2./DYP THROUGH S11,FOR J=2,1,J.G.22 WHENEVER.ABS.V(I,J).G.DYD.OR.DT.G.1./(U(I,J)/DX(I)+DYP2) 1.CR.U(I,J).L.O. PRINT RESULTS I,JtN,U(I,J),V(I,J),DT,DX(IDOY TRANSFER TO SE ThE R W I S E AB(J)=U(I,J)-DT*(U(ItJ)*(U(I-,J)-UtI-i,J))/DX(I)+VII,J)*( 1 U(I,J+1)-U(I,J-1))/DYM-(U(I,J+1)-2.*U(I,J)+U(I,J-1)f/DYP) 2 +DT*UDU(I) V(I,J)=V(I,J-1)-(AB(J-1)-AA(J-1)+AB(J)-AA(J))/DYD/DX(I) Sil END CF CCNDIT-IONAL TH R OUGH' S 1, F"CR'-' j 2,,JG'.-J.G 22 U T-1,J)=AA(J) S1 AAIJ)=AB(J) SPRAY.(U(2,23), AA(23)..AA(45)) THRCUGH S2,FOR 1=3,1,[.G.? THROUGH S21,FOR J=2,1,J.G.44 WHENEVER. ARBS.V(I, ji)"'G. DY-O.R.DT.G../ (U ( I, J )/OX I )'+DYP2 ) 1.CR.U(I. J).L.O. PRINT RESULTS I,J,N,U(I,J,V (I,J),OT,DX(I),DY TRANSFER TO SE OTHERWISE AB(J)=U( I,J)-OT*(U( I,J)*(U(I,J)-U(I-I,J))/DX(I)*V( I,J)*( 1 ULI,J+I)-U(IJ-1))/DYM-(U(I,J+I)-2.*U(I,J)+U(I,J-1))/DYP) 2 +CT*UDU(I) V I,J)=V(IJ-1)-(AB(J-I)-AA(J-1)+AB(J)-AA(J))/DYD/DX(I)

107 S21 thC Cf CCNOITIONAL TH1ROUGH S2, FOR J=2,1,J.G.44 U( I-J1 )=AA(J) S2 AA(J )=A 3J) CY=Y2 CYP=CY.t.2 LYD —2./OY CYM=2 * CY CYP22. /DYP AA(45):U( 7,45) THRCUGP S3,FCR J=I,1,J.G.23 AC(J )=AA(2*J-1) S3 U{ 7, J )U ( 7, 2*J-1 ) SPRAY.(U(7,23),AC(24)...AC(45),U(7,24)...U( 7,45 ) THROUGH S5,FCR J:2,1,J.G.44 WHtENEVER AS.V( I,J).G.DYD.OR.DT.G.1./(U( IJ ) /DX(I)+DYP2) I.OR.U( 1,J).L.O. PRINT RESULTS I,JNU(IqJ (I,),V(I,JOTtDX(l) TRANSFER TO SE CTHERW I SE AB(J)=U( IJ)-DT(( (IJ)*(Ut ItJ)-U( -1,J) )/OX(I) +V( 1,J)*( 1 U( I,J+1)-U( I,J-1) )/DYM-(IU IJ+i)-2.*U( I,J)+U(IJ-1) )/DYP) 2 +CTUCU ( I) V( I,J)=V(I,J-I)-(AB(J-)-AC(J-1)+ AB(J )-A C{JI))/DYD/Dxl I ) S5 ENO CF CONDITIONAL THRCUGH S6,FCR J=2,1,J.G.44 UC I-1,J)=-A(J) S6 AA(J)=AH(J) THRCUGH S7,FOR 1=9, l,I.G.30 THROUGH S71,FOR J=21,tJ.G.44 WHENEVER.ABS.V( I,J).G.DYD.OR.DT.G.l./(U( I,J)/DX(I +DYP2) 1.OR.U(I,J)-.L-O. PRINT RESULTS I,J,N,U IJ),V( I,J),OT,OX[)tDY TRANSFER TO SE CTHERWISE AB (J)U( 1,J)-DT*IU{ IJ)*{Ut IJ -U(I-1,J ))/OX(I I+VI IJ)*( I U( IJ+1)-U( IJ-1) )/DYM-(U I,J+l)-2.*U( IJ)+U(IJ-1) )/DYP) 2 +DT*UDU(I) V( I,J)=V(I,J-1)-(AB(J- 1 ) -AA( J-1)+AB(J)-AA(J ) ) /DYD/OX ( I ) S71 END CF tCNCITIONAL THROUCH S7,FCR J=2,l,J.G.44 ('I-1,J)=AA J ) S7 AA(J)=AB(J) THROUGH S8,FCR J =2, 1,J.G. 44 S8 30, J ) =AA J ) WHENEVER N.E.N1 PRINT RESULTS N,U(2, 1)...U(30,45) WHENEVER MA.E.O EXECUTE SKIP.(T3,T4,9) CTHERWIoE BACKSPACE RECORD OF TAPE 9 END' CF'ONDITIONAL WRITE BINARY TAPE 9,NU(l,l)...U(30,45),V(ltl)...V(30,45), 1 UCU()...UCIJ(30) END CF FILe TAPE 9 CR WHENEVER N/N2.E.MB TRANSFER TO NA CR WHENEVER N/N2.G.MB PB=N/N2

108 PRINT RESULTS N,I( 251 )...U(25,45),UI 30,)...U (30,45) TRANSFER TC NA ENC CF CUNCITIUNAL CI[ENSIGN U 1350,VB) V(l 350,VB),DX(30),AA{ 45),AB(45),T 1) >JI( 30) VLtCTru VALUES VB=2, 1,45 INTEGER I [J,N,NIN2, TNtMAMB, r T2, T3tT4 SE LNLGAC TAPE 9 tNC CF PROGRAM

lo9 C. CALCULATION OF TEMPERATURE OF WEDGE FLOW AFTER SUDDEN DROP OF WALL TEMIERATURE S COMPILE MAD, PUNCH OBJECT FTRAP. READ FORMAT VA,MA,MB,Nl,N2 T2,TTT3,T4,Q1,PR,Y1,Y2,DU,DX(2) 1...DX(30) VECTOR VALUES VA=S815, F6.1/ ( F6.2) *$ PAUSE A;O. 1 BACK SET HIGH DENSITY TAPE 9 REWIND TAPE 9 READ BCD TAPE 9,VR,TN(1)...TN(7) VECTOR VALUES VR=$7C6*$ WHENEVER TN(3).NE.$BCDEFG$ PRINT COMMENT $WRONG TAPE$ PRINr ON LINE FORMAT WRONG VECTOR VALUES WRONG=$H*WRONG TAPE MOUNTED PLEASE CORRECT**$ SKIP6. UNLOAD TAPE 9 PAUSE NO. 77777 TRANSFER TO BACK END OF CONDITIONAL EXECU[E SKIP.(T1,T2,9) READ BINARY TAPE 9,N,U(1,1)...U(30,45),V(1,1)...V(30,45), 1 Q(1,i)....-(30,45) NA N=N+. SPRAY.(Q1, AA(2)...AA(22)) DY=Y 1 DYP=DY. P.2*PR DYM=2.*DY I=2 THROUGH S11,FOR J=2,1,J.G.22 SlL AB(J)=Q(I,J)-DT*(U(IJ)*(Q(I, J)-[(I-1,J))/DX{ I)+V( IJ)*( 1 (( I,J+i)-Q( I,J-1) )/DYM-(Q( I,J+1)-2.*Q( I J)+Q( I,J-1) )/DYP) THROUGH S1,FOR J=2,1,J.G.22 Q (I-1,J)=AA ( ) S1 AA(J)=AB(J) SPRAY. [QI,AA (23)...AA(45)) THROUGH S2,FOR I=3,1,I.G.7 THROUGH S21,FOR J=2,1,J.G.44 S21 AB(J)=Q(I,J)-DT*(U(IJ)*(Q(I,J)-Q(I-IJ))/I)X(I)+V(I,J)*( I Q( I J+1)-Q( IJ-1l))/DYM- ( Q ( I,J+1 )-2.*Q( I,J ) +Q( I,J-1 ) )/DYP) THROUGH S2, FOR J=2,1 J.G.44 Q(I-1,J)=AA(J) S2 AA(J)=AB(J) DtY=YZ DYP=DYP. 2*PR DYM=2.*DY THROUGH S3,FOR J=l,l,J.G.23 S3 Q(7,J)=Q(7,2*J-1) SPRAY. (( 7,23),Q(7,24)...Q(7,45)) THROUGH S4,FOR I=8,1,I.G.30 THROUGH S41,FOR J=2,1,J.G.44 S41 AB(J)=Q(I, J )-DT*(U(I,J)*(Q(I, J)-Q(I-l,J) )/OX(I)+V( IJ)* 1 ( (I,J+L)-Q(I,J-l) )/DYM-(Q( I,J+1)-2.*Q(I,J)+Q(I,J-l))/DYP) THROUGH S4,FOR J=2,1,J.G.44 Q(I-1,J)=AA(J)

110 S4 AA(J)=AB(J) THROUGH S5,FOR J=2,1,J.G.44 S5 Q(30,J)=AA(J) WHENEVER N.E.N1 PRINT RESULTS N,Q(2,1)...Q(30,45) WHENEVER MA.E.O EXECUTE SKIP.(T3,T4,9) OTHERWISE BACKSPACE RECORD OF TAPE 9 END OF CONDITIONAL WRITE BINARY TAPE 9,N,U(l,1)...U(30,45)tV(l,l)...V30O,45) 1,Q(,ll)...Q(30,45) END OF FILE TAPE 9 TRANSFER TO SE OR WHENEVER N/N2.E.MB TRANSFER To 4NA OR WHENEVER N/N2.G.MB MB=N/N2 PRINT RESULTS N,Q(25,l)...Q(25,45),Q(30,1)...Q(30t45) TRANSFER TO NA END OF CONDITIONAL DIMENSION U(1350,VB),V(135 BQ(1350,VB)DX(30)AA(45),AB45 1 TN(7) VECTOR VALUES Vb=2,1,45 INTEGER I, JN,,N1,N2,TN,MA,MBTT2,T3,T4 SE UNLOAD TAPE 9 END OF PROGRAM

REFE RENCES 1. Prandtl, L. "iUber Flussigkietsbewegung bie sehr kleiner Reibung." Verh. III. Int. Math. Kongr., Heidelberg (1904), 484. 2. De Sante, D. F. and Keller, H. B. "Numerical Studies of Transition from Laminar to Turbulent Flow Over a Flat Plate." AFOSR 723, 1961. 3. Der, J. Jr. and Raetz, G. S. "Solution of General Three-Dimensional Laminar Boundary-Layer Problems by an Exact Numerical Method." 30th IAS Annual Meeting, New York, 1962. 4. Fromm, J. "A Method for Computing Nonsteady, Incompressible, Viscous Fluid Flows." University of California Los Alamos Scientific Laboratory, Report No. LA-2910, 1963. 5. Yang, K. T. and Kelleher, M. D. "On Hydrodynamic Stability of TwoDimensional Unsteady Incompressible Laminar Boundary Layers." Tech. Note 64-22, Notre Dame University, 1964. 6. Schuh, H. "Uber die'ahnlichen' Losungen der instationaren laminaren Grenzschichtgleichung in inkompressibler Str6mung." 50 Jahre Grenzschichtforschung (Editors: Gbrtler, H. and Tollmien, W.), Vieweg, Braunschweig, 1955. 7. Yang, K. T. "On Certain Similar Solutions to Unsteady Laminar Boundary-Layer Equations in Low-Speed Flow." J. Aero. Space Sciences, 25 (1958), 471. 8. Hayasi, N. "On Similar Solutions of the Unsteady Quasi-TwoDimensional Incompressible Laminar Boundary-Layer Equations." J. Phys. Soc. Japan, 16 (1961), 2316. 9. Hayasi, N. "On the Quasi-Two-Dimensional Laminar Boundary-Layer." J. Phys. Soc. Japan, 15 (1960), 208. 10. Hansen, A. G. Similarity Analysis of Boundary Value Problems in Engineering. Prentice-Hall, Englewood Cliffs, N. J., 1964. 11. Cheng, S. I. "Some Aspects of Unsteady Laminar Boundary-Layer on a Flat Plate." Quart. Appl. Math., 14 (1957), 337. 12. Hassan, H. A. "On Unsteady Laminar Boundary-Layer." J. Fluid Mech., (1960o), 300. 111

112 13. Hayasi, N. "On Semi-Similar Solutions of the Unsteady Quasi-TwoDimensional Incompressible Laminar Boundary-Layer Equations." J. Phys. Soc. Japan, 17 (1962), 194. 14. Blasius, H. "Grenzschichten in Fliissigkeiten mit kleiner Reibung." Z. Math. Phys., 56 (1908), 1. 15. Goldstein, S. and Rosenhead, L. "Boundary-Layer Growth." Proc. Camb. Phil. Soc., 32, (1936), 392. 16. Watson, E. J. "Boundary-Layer Growth." Proc. Roy. Soc., A 231 (1955), 104. 17. Cheng, S. I. and Elliott, D. "The Unsteady Boundary Layer on a Flat Plate." Trans. ASME, 79 (1957), 725. 18. Kestin, J., Maeder, P. F. and Wang, H. E. "On Boundary-Layers Associated with Oscillating Streams." Appl. Sci. Res., A 10 (1961), 1. 19. Hori, E. "Unsteady Boundary Layers (4th Report, Calculation of Boundary-Layer Around a Body with Arbitrary Shape by Use of a Series Expansion Method)." Bulletin of JSME, 5 (1962), 461. 20. Rosenhead, L. Laminar Boundary Layers. Clarendon Press, Oxford, 1963. 21. Van Dyke, M. Perturbation Methods in Fluid Mechanics. Academic Press, New York, 1964. 22. Schlichting, H. Boundary-Layer Theory. McGraw-Hill, New York, 1960. 23. Schlichting, H. "Berechnung ebener periodischer Grenzshichtstr6mungen." Phys. Z., 33 (1932), 327. 24. Schuh, H. "Calculation of Unsteady Boundary-Layers in TwoDimensional Laminar Flow." Z. Flugwiss, 1 (1953), 122. 25. Pohlhausen, K. "Zur naherungsweisen Integration der Differentialgleichung der laminaren Grenzschicht." Z. angew. Math. Mech., 1 (1921), 252. 26. Hartree, D. R. "On an Equation Occuring in Falkner and Skan's Approximate Treatment of the Equations of the Boundary-Layer." Proc. Camb. Phil. Soc., 33 (1937), 223.

1.1% 27. Yang, K. T. "Unsteady Laminar Boundary-Layers Over an Arbitrary Cylinder with Heat Transfer in an Incompressible Flow." Trans. ASME Series E, 81 (1959), 171. 28. Rozin, Lo A. "An Approximation Method for the Integration of the Equations of a Nonstationary Laminar Boundary-Layer in an Incompressible Fluid." NASA TT F-22, 1960. 29. Yang, K. T. "On the Calculation of Unsteady Incompressible Laminar Boundary-Layers Over Arbitrary Cylinders." Tech. Note 62-11, Notre Dame University, 1962. 30. Hayasi, N. "On the Approximate Solution of the Unsteady Quasi-TwoDimensional Incompressible Laminar Boundary-Layer Equations." J. Phys. Soc. Japan, 17 (1962), 203. 31. Lin, C. C. "Motion in the Boundary-Layer with a Rapidly Oscillating External Flow." Proc. 9th Int. Congr. Appl. Mech., 4 (1956), 155. 32. Gibson, W. Unsteady Laminar Boundary-Layers. PhD Thesis, Massachusetts Institute of Technology, 1957. 33. Lighthill, M. J. "The Response of Laminar Skin Friction and Heat Transfer to Fluctuations in the Stream Velocity." Proc. Roy. Soc., A 224 (1954), 1. 34. Wu, J, C. "On the Finite Difference Solution of Laminar BoundaryLayer Problems)" Proc. Heat Transfer and Fluid Mechanics Institute, 1961, 955 35. Raetz, G. S. "A Method of Calculating Three-Dimensional Laminar Boundary-Layers of Steady Compressible Flows." Report No. NAI-58-73, Northrop Aircraft Inc., 1958. 36. Hellumrs, J. D. Finite Difference Computation of Natural Convection Heat Transfer. PhD Thesis, University of Michigan, 1960. 37. Wilkes, J. 0, The Finite Difference Computation of Natural Convection in an Enclosed Rectangular Cavity. PhD Thesis, University of Michigan, 1963. 38. O'Brien, G., Hyman, M. and Kaplan, S. "A Study of the Numerical Solution of Partial Differential Equations)" J. Math. Phys., 29 (195), 2253

114 39. Richtmyer, R. D. Difference Methods for Initial-Value Problems. Interscience Publishers, New York, 1957. 40. Douglas, J. Jr. "On the Relation Between Stability and Convergence in the Numerical Solution of Linear Parabolic and Hyperbolic Differential Equations." J. Soc. Indust. Appl. Math., 4 (1956), 20. 41. Lax, P. D. and Richtmyer, R. D. "Survey of the Stability of Linear Finite Difference Equations." Communs. Pure and Appl. Math., 9 (1956), 267. 42. Forsythe, G. E. and Wason, W. R. Finite Difference Methods for Partial Differential Equations. John Wiley and Sons, New York, 1960. 43. John, F. "On the Integration of Parabolic Equations by Difference Methods." Communs. Pure and Appl. Math., 5 (1952), 155. 44. Douglas, J. Jr. "A Survey of Numerical Methods for Parabolic Differential Equations." Advances in Computers, 2 (1961), 1. 45. Douglas, J. "On the Numerical Integration of Quasi-Linear Parabolic Equations." Pacific J. Math., 6 (1956), 35. 46. Lees, M. "Approximate Solution of Parabolic Equations." J. Soc. Indust. Appl. Math., 7 (1959), 167. 47. Hildebrand, F. B. Introduction to Numerical Analysis. McGrawHill, New York. 19%6. 48. Wilkinson, J. H. Rounding Errors in Algebraic Processes. Her Majesty'S Stationery Office, London, 1963. 49. Hill, P. G. and Stenning, A. H. "Laminar Boundary-Layers in Oscillatory Flow." Trans. ASME Series D, 82 (1960), 593. 50. Nickerson, R. J. The Effect of Free Stream Oscillation on the Laminar Boundary-Layers on a Flat Plate. ScD Thesis, Massachusetts Institute of Technology, 1957. 51. Ghosh, A. "Contribution a L'etude de la Couche Limite Laminaire Instationnaire." Publications Scientifiques et Techniques, No. 381, 1961.

115 52o Schubauer, G. B. and Skramstad, H. K. "Laminar Boundary Layer Oscillations and Transition on a Flat Plate." NASA Report No. 909, 1947. 53- Schlichting, H. "Zur Entstehung der Turbulenz bei der Plattenstr6mung." Nach. Gesell. d. Wiss. z. Gottingen, Math. Phys. Klasse, 1933, 181. 54. Schlichting, H. "Amplitudenverteilung und Energiebilanz der kleinen St6rungen bei der Plattenstromung." Nach. Gesell. d. Wiss. z. Gottingen, Math. Phys. Klasse, 1 (1935). 55. Barakat, H. Z. Transient Natural Convection Flows in Closed Containers. Ph.D. Thesis, University of Michigan, 1965.

UNIVERSITY OF MICHIGAN 3 9015 0282o 7782 THE UNIVERSITY OF MICHIGAN DATE DUE ~lzo //. Y