'E‘EEE; SEEEEE 5"”! GEE E‘GEEEEEIIERY LEE’EZE’ En ‘LGEE’ idol SUB.) ECE T0 ROTATEEEN or: 0 . P E .masts E'QP E’E’ie EDGQE'BfiS 0E 2’53. 0. [‘11:1 N {‘10: \ vvvvvvvv ”area-p" EEEE (15.73 951% .i ’tz'i-JEEY 11.1 1 u L .4 Mangat flags ChawEa 3%6? , THESIS llllflllll“HIE“NEEHHIIEHIUIEIIHHEHIIEEll . 310385 9629 This is to certify that the thesis entitled THE STABILITY OF BOUNDARY LAYER FLOW SUBJECT TO ROTATION presented by Manga] Dass ChawIa has been accepted towards fulfillment of the requirements for DoctoraE degree in Mech. Engr. JWI/W; Major professor November 14, 1969 Date 0-169 LIBRARY Michigan State ' University “—o—_._- .u . ,.- <. 73%49 2.. ABSTRACT THE STABILITY OF BOUNDARY LAYER FIOW SUBJECT TO ROTATION By Mangal Dass Chawla The stability of three dimensional Tollmien-Schlichting waves moving in a boundary layer flow on a flat plate which is rotating about its leading edge is investigated. Primary interest is in the effect of rotation on the critical point of the neutral stability curve. The problem is formulated in terms of a sixth-order differ- ential equation, analogous to the Orr-Sommerfeld equation, with homogeneous boundary conditions. The eigen-value problem thus formed is then transformed into an initial value problem.which is solved numerically. The numerical technique develOped first solves for the primary velocity profile and its derivatives. Then using a Runge- Kutta fourth order method and a Special filtering technique the eigen-values are determined by iteration. The results indicate that as the flat plate rotates into the flow it is destabilizing and as the plate rotates away from the flow it is stabilizing. There is also a Span-wise wave number' influence which amplifies this effect. The results compare favour- ably with known experimental results. THE STABILITY OF BOUNDARY LAYER FLOW SUBJECT TO ROTATION By Mangal Dass Chawla A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1969 4.47-70 A DEDICATION This thesis is dedicated to the author's loving wife Santosh who supplied encouragement and displayed great patience and tolerance during the entire course of the graduate study without which the author would not have successfully completed this work. This is also dedicated to the author's elder son Rajeev and the newly born twins Rakesh and Renu. ii ACKNOWLEDGEMENTS The author wishes to express his sincere thanks and deep gratitude to Dr. Merle C. Potter, Associate Professor of Mechanical Engineering for his valuable guidance, advice and assistance through- out the course of this study and in the preparation of the thesis. Thanks are also due to Dr. M.C. Smith, Dr. P.K. Wong, Dr. J.L. Lubkin, Professor L.V. Nothstine and Dr. R.F. McCaulay for serving as members of the doctoral guidance committee. The author wishes to express his gratitude to the Department of Civil Engineering, Michigan State University for providing a graduate teaching assistantship from Fall 1965 to Spring 1966, employing the author as Assistant Instructor from Fall 1966 to Spring 1967 and for allocating the computer time. Thanks are also due to the National Science Foundation for providing the author witha Special Graduate Research Assistantship from July 1967 to June 1968 and for making funds available for computer analysis. Thanks are also due to the Division of Engineering Research and the Department of Mechanical Engineering for providing Research and Teaching Assistant- ships from July 1968 to March 1969 and from March 1969 to June 1969 reSpectively. Thanks are also due to Mrs. N. Barnes and the author's wife Santosh, who typed and did the proof reading of the manuscript reSpectively. iii The author is also indebted to the Canadian Commonwealth Scholarship Committee which got the author started on his graduate work by granting a scholarship for the M.Sc. Engineering degree. iv TABLE OF CONTENTS List of Tables ...................................... List of Figures . ............. ....... ................ Nomenclature ........................................ CHAPTER I. INTRODUCTION ........................................ 1.1 Review of Literature ........................... 1.2 Purpose of the Present Investigations .......... 1.3 Conditions of this Study ......... . ............ . II. INFINITESIMAL DISTURBANCES ON A BOUNDARY LAYER FLOW . III. IV. 2.1 Governing Equations ......... . ...... ...... ...... 2.2 Non-dimensional Variables ........ . ............. 2.3 The Primary Flow ......... . ..................... 2.4 Three Dimensional Disturbances ................. 2.5 Mathematical Formulation of the Three Dimensional Linear Stability Problem ........... 2.6 Boundary Conditions .......... . ................. 2.7 Eigenvalue Problem ............................. 2.8 Closure ........................................ NUMERICAL SOLUTION OF THE PROBLEM ..... .............. 3.1 Solution for the Primary Velocity Profile ...... 3.2 Approximate Solutions of the Orr-Sommerfeld Equation for Zero Angular Rotation ............. 3.3 Approximate Solutions of the Orr-Sommerfeld Equation for Non-Zero Angular Rotation ......... 3.4 Numerical Integration of the Orr-Sommerfeld Equations ...................................... .(a) Initialization ............................ (b) Growth of Eigen-functions ................. (c) Purification Scheme - Filtering the Contributions of the Growing Solution ..... (d) Iteration Scheme for the Eigenvalues ...... (e) Criterion for Guessing the Eigenvalues .... CONCLUSIONS AND SUGGESTIONS FOR FURTHER RESEARCH .... 4.1 The Primary Velocity Profile ................... 4.2 Eigenvalues for the Case of Neutral Stability for Zero Angular Rotation .... .................. V oooou—s I-' 10 12 13 15 19 20 22 23 23 25 26 27 27 29 31 33 35 36 36 36 4.3 Effect of the Angular Rotation Q and Wave Number 3 .............. . ....................... 4. 4 Summary of Conclusions ....... . ................ 4.5 Recommendations for Further Research ........... BIBLIOGRAPHY ........................................ ILLUSTRATIONS (Tables and Figures) .................. APPENDICES A NUMERICAL TECHNIQUES ...................... A-l Numerical Technique to solve for the Primary Velocity Profile ............. A-2 Numerical Technique to solve the Governing Stability Equations ........ B FLOW DIAGRAM .............................. C COMPUTER PROGRAMS ........... .............. C-l Use of the Computer Programs ......... C-2 Description of the Computer Programs . C-3 Listing of the Computer Programs ..... vi 38 39 4O 41 45 79 79 81 85 89 89 91 92 LIST OF TABLES Table Page 1. Values for the Primary Velocity Profile .............. 45 2. Definition of Eigenvalues used by Different Authors .. 48 3. Comparison of Critical Eigenvalues for Zero Rotation . 48 4. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci = 0. 0), for Q = 0.00 and B = 0.00 ............... ....... . ............... 49 5. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C, = 0.0) for Q = -0.01 and a = 0.20 ...... ..... ....}........... ...... .....,. 50 6. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) for 0 = -0.10 and B = 0.00 ...... . ...... .. ....... ....... ........ ... 50 7. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C. = 0.0) for Q = -0.10 and a = 0.10 ........ . ..... .} ............... . ........ 51 8. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) for 0 = -0.10 and B = 0.20 ................... .... ........... . ..... 51 9. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci =0.0) for 0 = -0.10 and 680.30 OOOOOOOOOOOOOO ......OOOOOOOOOO00.0.0...0 52 10. Eigenvalues for the Orr- Sommerfeld Equation for the Case of Neutral Stability (Ci =0. 0) for Q= 0.01 and B = O. 20 ........ . ............................... 52 ll. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci =0.0) for 0 = 0.10 and B= 0.10 ................... . ........... . ........ 53 12. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) for Q = 0.10 and B = 0.20 ............ ...... .... ........... . ...... 53 13. Eigenvalues for the Orr- Sommerfeld Equation for the Case of Neutral Stability (Ci =0.0) for 0 = 0.10 and B= 0.30 ......... . ........... . .................. 54 14. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci =0.0) for Q = 0.10 and B= 0.40 .......... . .................. . .......... 54 15. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci =0. 0) fOr' 0 = 0.10 and B = 0.50 ..... ..... ....... . ....... ... ........ 55 vii 16. 17. 18. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) for O = 0.10 and B = 0.60 ........................................ 55 Eigenvalues for the Orr-Sommerfeld Equation.£or the Case of Neutral Stability (C1 = 0.0) for O = 0.50 and B = 0.20 ........................................ 56 Eigen-functions near a Critical Point on the Lower Branch of the Neutral Stability Curve (C1 = 0.0) for 0 = -0.10 and B = 0.20 ........................ 57 viii LIST OF FIGURES Figure 10. 11. 12. 13. Primary Flow ......................................... Numerical Integration of the Orr-Sommerfeld Equation . Comparison of the Neutral Stability Curves for the Blasius Boundary Layer over a Flat Plate for 0 = 0.00 and B = 0.00 ........................................ Comparison of Eigen-functions ¢ = ¢r + i ¢i with those of Radbill and Van Driest, CorreSponding to the Data of Schubauer and Skramstad for the Upper and Lower Branches of the Neutral Stability Curve (C1 = 0.0) for' 0 = 0.00 .......... ... ..................... . ..... Comparison of the First Derivative of the Eigen- function ¢' = ¢; + i ¢£ with those of Radbill and Van Driest, Corresponding to the Data of Schubauer and Skramstad for the Upper and Lower Branches of the Neutral Stability Curve (C1 = 0.0) for 0 = 0.00 ... Plot of Wave Number 0 against Reynolds Number R for the Neutral Stability Curves for different values of Angular Rotation. Q for B = 0.20 ..... ........... Plot of Wave Velocity Cr against Reynolds Number R for the Neutral Stability Curves for different values of Angular Rotation Q for B = 0.20 ......... Effect of Angular Rotation O on Critical Eigen- values Rerit’ (Cr) and acrit for AB = 0.20 ... Effect of varying Angular Rotation. O 'on the Neutral Stability Curves for B = 0.10 ....................... crit’ Effect of Varying Angular Rotation. O on the plot of Wave Velocity c against Reynolds Number R for the Neutral Stability for B = 0.10 ...... ............ .... Effect of Varying Angular Rotation O on the Neutral Stability Curves for B = 0.30 .......... ............. Effect of Varying Angular Rotation. 0 on the plot of Wave Velocity C against Reynolds Number R for the Neutral Stability for B = 0.30 .......... ........ Plot of Eigen-functions ¢ = ¢r +qi 91 near a Critical Point on the Lower Branch of the Neutral Stability Curve (C1 = 0.0) for n = -0.10 and B = 0.20 ..... . ............... . ....................... ix Page 60 60 61 62 63 64 65 66 67 69 70 71‘ 14. 15. 16. 17. 18. 19. 20. Plot of the First Derivative of the Eigen-function d' = ¢; + i ¢£ near a Critical Point on the Lower Branch of the Neutral Stability Curve (C1 = 0.0) for 0 = -0.10 and B = 0.20 ........ ..... ... ........ Plot of Eigen-functions Y = Yr +-i Y1 near a Critical Point on the Lower Branch of the Neutral Stability Curve (C1 = 0.0) for Q = -0.10 and a=0.20..... ................................ . ....... Effect of varying Wave Number B on the Neutral Stability Curves (C1 = 0.0) for Angular Rotation 0 = -0.10 .. ........ . .......... ... ........ ... ......... Effect of Varying Wave Number B on the plot of Wave Velocity Cr against Reynolds Number R for the Neutral Stability for O = -0.10 ....... . ......... Effect of varying Wave Number B on the Neutral Stability Curves (C1 = 0.0) for Angular Rotation O = 0.10 ............................................. Effect of Varying Wave Number B on the plot of Wave Velocity Cr against Reynolds Number R for the Neutral Stability for Q = 0.10 ...... . ........... Plot of Wave Number B against Critical Reynolds Number Rcrit for Angular Rotation Q = 0.10 and -0.10 72 73 74 75 76 77 78 r 1 C r C. 1 _ d D dy E = exp[i(ax + Bz) - 1 act] r, r1, r2, r3 1’ S2 NOMENCLATURE Complex propagation speed of wave Wave speed Amplification factor Differential operator Components of body forces Dimensionless stream function Indices Differential operators (i = 1 to 9) Pressure of undisturbed flow (Non-dimensional) Pressure of undisturbed flow (Dimensional) Pressure after perturbation Perturbation pressure Reynolds number based on boundary layer thickness 6 Reynolds number based on the distance from the leading edge Change in R Radial distance from the origin Constants Constants xi C) Cl v = V +-v' 4» Test function based on combined D0 at the plate Values of T during iteration Time (Non-dimensional) x - Component of velocity of primary flow Free stream velocity x - component of velocity after perturbation Perturbation in X - component of velocity Amplitude of u' perturbation of velocity Dimensional components of velocity along coordinate axes (i = 1,2,3) Dimensionless components of velocity along coordinate axes (i 1,2,3) y - component of velocity of primary flow y - component of velocity after perturbation Perturbation in y - component of velocity Amplitude of v' perturbation of velocity 2 - component of velocity of primary flow xii w = W +-w' g) — Q. m - U2 a w = g +-w' w! 0 dz W1 :11] X X (1) ,X(2) E. 1 x. 1 Y Z 21,22,23’24’25’26 xiii z - component of velocity after perturbation Perturbation in z - component of velocity Amplitude of w' perturbation of velocity Pressure term including the effect of centrifugal forces (Dimensional) Pressure term including the effect of centrifugal forces (Non- dimensional) Pressure term w after perturbation Perturbation in pressure term w Amplitude of w' perturbation in pressure term Coordinate parallel to plate in the direction of free stream Multiplying constants Dimensional coordinates X, Y and Z respectively (i = 1,2,3) Non-dimensional coordinates x, y and z reapectively (i = 1,2,3) Coordinate prependicular to plate Coordinate perpendicular to the planes of X and Y New variables equal to U, 0, a, DB, RE and 0% respectively 4| Wave number Change in a Spanwise wave number Boundary layer thickness Displacement thickness Transformation variable Laplace operator Vorticity Density of fluid Kinematic viscosity of fluid Time (Dimensional) Square of the combined wave number Angular rotation (Dimensional) Angular rotation (Non-dimensional) Complex eigenfunction - y dependent factor in perturbation stream function VI Individual eigenfunctions at the plate Combined eigenfunction at the plate Complex stream function - y dependent factor in perturbation stream I function u Individual stream-functions at the plate xiv Combined stream function at the plate XV CHAPTER I INTRODUCTION 1.1 Review of Literature The study of transition from laminar to turbulent flow for parallel shear flow started at the end of the ninteenth century. This was first experimentally studied by Reynolds (1883) by inject- ing dye in flow through pipes. 0n the basis of these experiments Reynolds (1895) stated that the motion of the water was direct or sinuous according as the mean velocity Umfg kl;- Do where D = Diameter of the pipe p = Density of water 0 = Absolute viscosity of water k = Dimensionless number or 9——ng § k. The number k was found to lie between 1900 and 2000 for a circular tube. Thus Reynolds found that steady direct motion in round tubes was stable if LEM < 1900 and unstable if E—Um—D— > 2000. 0. Later his theoretical analysis gave the value of k = 517 for flow 1 between parallel plates. Then with the result that the critical velocity was inversely preportional to the hydraulic mean depth, he inferred that for a circular pipe k = 1034. This was below the experimental value. Reynolds stated that this discrepancy was due to not meeting all the kinematical conditions. The corre3ponding values for round and flat tubes were later found to be 470 and 167 respectively by Sharp (1905). In the meantime Lord Rayleigh (1880, 1887) studied the problem of stability in perfect fluids. Lord Kelvin (1887), objecting to the neglect of viscosity by Rayleigh, found that the mathematical difficulties were enormous when viscosity was considered. Later Orr (1907) and Sommerfeld (1908) each independently made a complete theoretical analysis of the entire stability problem. The stability equation thus derived is now known as the Orr-Sommerfeld equation. Orr discussed the stability of both perfect and viscous fluids. Sommerfeld studied couette flow and found no instability. Both Orr and Sommerfeld assumed sinusoidal, periodic, infinitesimal dis- turbances to study stability. Rayleigh (1914) took the Orr-Sommerfeld equation for the two-dimensional frictionless case, given by (U - c)(¢" - a2¢) - U"¢ = 0 (1-1-1) C II where mean velocity c = complex wave velocity 01 = wave number ¢ = eigen-function, written in the form as" = m + 030 and found that at a distance from the boundary where U = c, U" must be zero, that is there must be an inflexion point if the flow is to be unstable.He called the layer at this point the inner friction layer. Rayleigh thus found that the steady motion of an inviscid fluid in two dimensions between fixed parallel planes was stable pro- vided the primary velocity U, parallel to the walls and a function of y only, was such that U" does not change sign anywhere. In other words, Rayleigh proved that the existence of a point of inflexion constituted a necessary condition for the occurrence of instability. Later Tollmien (1935) showed that this was also a sufficient condition for the disturbances to grow and create instability. A similar result, found by Tietjens (1925), shows that if the basic velocity profile is approximated by multiple straight line segments stability results if the corners are convex and instability results if the corners are concave. Prandtl (1914) defined transition, separation and drag co- efficients on bodies. He first tripped a boundary layer with a wire. The discovery of a boundary layer opened new fields for the study of transition to turbulence. Prandtl (1914, 1921) started the analysis of the stability equation retaining only the larger terms involving viscous effects close to the wall and found that with this hypothesis all profiles made up of line segments were unstable for all wave lengths. This result was later confirmed by Tietjens (1925). Heisenberg (1924) applied this assumption to the basic profile, taking it as parabolic, and concluded instability. Tollmien (1929) arrived at the conclusion that viscous effects must be taken into account near the wall and at the level where U = c. He later found the critical Reynolds number to be 420 for flow over a flat plate at zero incidence. Tollmien's method, also known as the Asymptotic Theory, was later applied successfully by Schlichting (1933, 1935) to the velocity profiles in the boundary layer, when approximated by a straight line or a parabola. Lin (1945, 1946, 1966) gave the asymptotic theory a firm mathematical footing. He explained the behaviour of the functions near the singularities. Later, Shen (1954) working with the asymptotic method improved the accuracy of the stability calculations. Synge (1938) surveyed the whole subject of Hydrodynamic Stability and explained in detail the physical problem and its mathematical formulation. He described the application of the method of the exponential time-factor to the couette flow with general and plane disturbances and to the Poiseuille motion in a tube of general section or in a circular tube with disturbances having rotational or axial symmetry. The asymptotic method of solution has been successfully applied to many stability problems. Recently, the advent of large digital computers, have made numerical schemes for solving the Orr-Sommerfeld equations very popular. Numerical methods are usually preferred to the asymptotic methods as they are more accurate and relatively less time consuming, however, asymptotic methods are still employed, especially in the zones of large Reynolds number where numerical solutions show poor convergence and when general properties of the solutions are to be studied. Gill (1950) discussed the integration of the differential equations on the computers and also developed a.form that gave the highest accuracy with minimum machine instructions. His method later became very popular and has been successfully used by many authors. Thomas (1953) investigated the problem of the stability of plane poiseuille flow and verified the work of Pekris (1948) who had used the method of asymptotic expansions. Kaplan (1964) employed a numerical scheme to study the stability of laminar incompressible boundary layers with compliant boundaries. He presented the numerical results for a variety of compliant boundaries. He forms an initial value problem from the boundary value problem at some distance from the compliant boundary before starting numerical integration as did Nachtsheim (1964). Kaplan studied the problem in detail when some roots grow very large and thus due to the limitations of the machine computations, start affecting the other roots which should otherwise stay independent. In Kaplan's terminology, the roots were called "the growing roots" and "the well behaved roots" respectively. He outlined a method of filtering the contributions of the growing roots from the well behaved roots as the integration proceeds and later repaired the well behaved roots. This method as advocated by Kaplan was later advantageously used by Lee and Reynolds (1964) and Reynolds and Potter (1967). Lee and Reynolds discussed other numerical methods and compared them to the variational solution of the Orr-Sommerfeld equations. They, of course, discussed other numerical integration methods where one may or may not use the filtering scheme used by Kaplan. Potter (1965) studied the stability of plane Couette-Poiseuille flow by asymptotic expansions and later in 1967 of symmetrical parabolic flows numerically. Radbill and Van Driest (1966) investigated the stability of laminar boundary layers by quasilinearization and numerically. 0n the experimental side work was done by Nikuradse (1933). He excited the boundary layer with sinusoidal oscillations near the leading edge of a flat plate in water. Dryden (1936) measured the velocity fluctuations in a laminar boundary layer, but these fluctuations were said to be imposed by the mean flow because they were irregular. Schubauer and Skramstad (1947) experimentally proved that induced oscillations exist in the boundary layer. Dryden (1947) and Schubauer and Skramstad (1947) experimented with very low free stream turbulence obtaining very large transition Reynolds numbers. It is generally believed that there are three stages in the transition process before the turbulence sets in. At the first stage infinitesimal Tollmien-Schlichting waves become unstable. At the second stage two-dimensional Tollmien-Schlichting waves become three- dimensional. At the third stage burstsfrom the boundary start. The flow will be termed laminar up to the start of the third stage. The fourth stage is characterized by the constant burst rate from the wall. The critical Reynolds number Rcrit given by the stability theory is the lowest Reynolds number at which a shear layer can bear and amplify an infinitesimal Tollmien-Schlichting wave. This is the start of the first stage as mentioned above. Morkovin (1958) writes that transition to turbulence is generally found to occur at Reynolds numbers much greater than RC which indicates that rit’ certain investigators may refer to the third stage of the transition process as the beginning of transition. Emmons (1951) believes that sometimes disturbances that trigger transition may be "local in time". But once initiated, a small turbulent spot moves downstream and grows steadily in all directions, sweeping out a tongue like wedge with time. Morkovin (1958) writes that in the laminar layer the un- steadiness and vorticity of a lump of fluid decays with time during its travel, that is, the unsteadiness does not feed on the energy of the mean shear motion. Morkovin further states that many investigators feel that three-dimensionality, however small, must lead to small loops in originally straight vortex lines of the steady motion and/ or the unsteady regular two-dimensional velocity distributions. These loops would inevitably stretch into longer and longer loops by the differences in mean velocity at different heights of shear layer, thus increasing w' fluctuations. In the area of present study Conard of Princeton University (1962) investigated the effects of rotation on the stability of laminar boundary layers on curved walls. He assumed these walls to be rotating about an axis perpendicular to both the flow direction and the radii of curvature. He considered Taylor-G3rt1er waves and arrived at the conclusion that positive rotations (rotations for which the absolute velocity of the fluid exceeds its relative velocity) were destabilizing and negative rotations were stabilizing. Regarding Tollmien-Schlichting waves Conard stated that when he considered two- dimensional Tollmien-Schlichting waves, the curvature and rotation terms dropped out of the equations for the order of his approximations. Had he considered three-dimensional disturbances the rotation in- fluence would have been retained. Experimentally, the only work available is due to that of Halleen and Johnston of Stanford (1967). Their results are similar to those obtained by Conard and now numerically verified by this study. 1.2 Purpose of the Present Investigations The purpose of the present study is to investigate the stability of boundary layer flow subject to rotation with the axis of rotation perpendicular to the flow direction and parallel to the flat plate (see figure 1). The flow is subjected to three-dimensional Tollmien- Schlichting waves including wave number a and Spanwise wave number B. Primary emphasis has.been placed on neutral stability, that is, the conditions for which these infinitesimal disturbances neither decay nor grow with time. The result establishes the beginning of the first stage of transition. New information which this study furnishes is the influence of angular rotation 0 on this neutral stability of the boundary layer. Various spanwise wave numbers 5 are included for each n. 1.3 Conditions of this Study In the numerical computations for the present investigation, the boundary value problem is transformed into an initial value problem with calculations starting at Y = 1.5 5. This is the point which will be used as "infinity". It is dictated by the step size and the storage limitations of the machine. Another limiting condition for the present study is that convergence to the true eigenvalues is assumed to be complete when the test function T based on ¢' at the plate divided by ¢' at one step behind is less than 10-3. Under these conditions, this study determines the critical Reynolds number Rcrit for different values of 0 and B at which instability due to three-dimensional, infinitesimal disturbances first sets in. The values of Q and B along with the results are discussed in Chapter IV. CHAPTER II INFINITESIMAL DISTURBANCES ON A BOUNDARY LAYER FLOW 2.1 Governing Equations Consider laminar flow of an incompressible fluid over a flat plate which is rotating with an angular velocity 5 about its leading edge in a counterclockwise sense. The coordinate system, shown in figure 1, is attached to the rotating plate so that the flow can be considered steady when measured in this rotating system. This requires the fluid at infinity to be rotating with the same angular velocity. The Navier-Stokes equations, governing the motion, can be written as ._ ._ z. Bl1i - aui - 1 5.2 B ui l -—-—- —.." = - — L - d 2 51' +uj ax F1 pax +vai ax " (Z‘OX ‘) j i J j + 2 3.. Eh_l (2-1-1) where the bars denote dimensional quantities and hi = (0,0,5). (2.1.2) The last two terms in equation (2.1.1) account for the non-inertial coordinate system chosen. We may combine the pressure term and the centrifugal acceleration term to obtain ‘O'p—I 3% +a_%‘a X 32) = - 3:— (2.1.3) 1 i ax, 9 10 .. - =3. _ where 0) = B - %‘0 X r‘z. (2.1.4) Neglecting body forces and introducing (2.1.3) in (2.1.1), there results 3:1 - 5-1 a; 5231 - - —.— + —.- = - .. 01- L11 axj 3x1 + " aijaij + 2 eijk j 0k (2'15) 50. 1 . (2.1.6) Equations (2.1.5) and (2.1.6) are the equations which govern the flow. 2.2 Non-Dimensional Variables To non-dimensionalize choose u as the reference velocity 0 and 6 as a characteristic length, then the dimensionless variables will be 1 ... IF“:- “7"fo I I) ll CO’ 18 L)! 8:: P (2.2.1) U '1' t .__J°_ 6 =0. ‘” 2 U 0° J Introducing (2.2.1) in (2.1.5) and (2.1.6) we obtain the Navier Stokes and continuity equations in dimensionless form as 2 3:. Elihu. laui e no (222) at 3 3x 3x. R 5x 5x ijk j k ' ' j 1 13 and an. f=0 (2.2.3) 81 where the Reynolds number is Uaé R =-7;- . (2.2.4) The Reynolds number can be changed by varying either 09, 0, v, or Umé, keeping the other variables constant. For flow over a flat plate, as in the present case, both Uup and v will be kept fixed. The boundary layer thickness 6 varies as the square root of the distance from the leading edge, that is 6 ¢ g§_ (2.2.5) U w or a a L (2.2.6) V/R x UQX where R.x = -;-. Taking 5 as a point at which the velocity is 99.9% of the free stream velocity U”, the constant of proportionality in (2.2.6) is approximately equal to 5.0. Using this relationship for 6 in our definition of R we obtain R = 5 fiiflx . (2.2.7) T Knowing the Reynolds number R at which the flow becomes unstable one can then determine the distance from the leading edge along the plate at which the infinitesimal waves become unstable. 12 2.3 The Primary Flow In order to observe the influence of Q on the primary flow, let us write equations (2.2.2) in cartesian coordinates as L... La]- L=-a—+— + 2.3.]. at u BX v BY Bx AU 2V 0 ( ) and Ali. 31+ u=_aQ+.1. _ at u Bx v BY BY R Av 2u 0 (2.3.2) The equation of continuity is 53-+-53 = 0. 2.3.3 Bx BY ( ) Differentiating (2.3.1) with respect to y, (2.3.2) with respect to x and subtracting gives 2 2 2 2 3.2t.- 3.2.) +.u a_2_ _ u a;!.+.a2.aa.- 52.32. (atay atax axay 6x2 BY ax ax ax +.v D__ - v a___ +.a_.a_.- 52.03.: 5,2 ayax ay ay ax By -A(§;- g? + 20g+§$ (2.3.4) The vorticity is g = 53-- 51. (2.3.5) Inserting this in (2.3.4) gives —£ +2 +53. =-l + 53-+-5!' 2.3.6 C(ax BY) R Ag 2max 0y) ( ) Continuity then allows us to write 2Q =.l R Ag . (2.3.7) 13 Equation (2.3.7) shows that the coriolis term does not enter the vorticity equation. By an examination of equation (2.3.1) the effect of 0 is small and thus neglected because of the presence of the small velocity component v. The effect of Q in equation (2.3.2) would be to influence the BD/BY term. Hence 0 does not affect the primary velocity profile and the usual Blasius boundary layer results. 2.4 Three Dimensional Disturbances Consider infinitesimal perturbations in the variables u, v, w and w, so that after perturbing, these variables change to u = U + u' v = V +'v' w = w' (2.4.1) w = w +-w' where U, V and g. pertain to the primary flow. The component of primary velocity parallel.tC)z—axis is zero and the one parallel to y-axis is comparatively small. It is not unreasonable to assume V to be of order 6. Insert equations (2.4.1) in (2.2.2) and (2.2.3), neglect non-linear terms and write the resulting equations in cartesian coordinates to yield E+Ua£+u'afl+v'afl+vag_v+w'ag=_aw_ at ax BX BY BY BZ BX +%Au' + 2v'n (2.4.2) av_'+uay_'_+vax_i+u:al+vv51+w:al=-aal Bt Bx BY Bx BY B2 BY +fiAv' - 2u'n (2.4.3) l4 £+Uafl+vfll =-Bw_' Bt BX BY BZ +%Aw' (2.4.4) where the Laplacian A is 2 2 2 6237+57+57 (2.4.5) ax BY 62 and a£+axl+a2l=o (2.4.6) ax BY 62 Lin (1966) has stated that, "It can be shown, by careful examination of all the terms in the first order disturbance equations of continuity, motion and energy, that the basic flow may be treated as essentially parallel. The analysis is fairly complicated, but the conclusion has been justified”. The resulting equations are then I I I EL+Uau_+Vea_Q=-aw_+lAu'+2vin (2.4.8) Bt Bx BY BX R .I I I a_v_+Uay_ --aw_+lAv' - 2u'0 (2.4.9) at BX 837 R 1 I ' ' aw' + U aw = - aw_ +1 Aw: (2.4.10) at ax 62 R and I I I al“—-+51’—+5W— = 0. (2.4.11) Consider further that the perturbations can be expanded in the form of three-dimensional disturbances, given by 15 u' = 0(y) . E B v' = 0(y) . E w. = My) . E P (2.4.12) w'=0®) E J where E = exp[i(ax + Bz) - 1 act]. (2.4.13) The above implies a spatially periodic wave with complex amplitudes, B, 0 and- 0 in the "x, y and 2 directions respectively. a and B are wave numbers which will be real. C is a complex wave velocity, given by C = C +'i C, (2.4.14) r 1 where C1 is the amplification rate. Ci positive means growth and C1 negative implies decay. Ci equal to zero yields neutral stability, that is, neither growth nor decay of the perturbations. Actually, only the real parts of equations (2.4.12) are used to give the physical velocity and pressure perturbations. 2.5 Mathematical Formulation of the Three-Dimensional Linear Stability Problem Introduce equations (2.4.12)in (2.4.8) to (2.4.11), cancel E throughout and rearrange. The following system of second order linear differential equations results: 16 -20)V+101R&> 3 [D2 - (d2 + 32) - iaR(U - 6)]fi = R(DU [02 - (62 + 92) - mam - em; = 21200 + RD6) (2.5.1 a to d) [D2 - 022 + 52) - iaR(U - c)]& = iBRfi 1(afi + 56:) +00 = 0 where D is a differential operator given by (2.5.2) ale- ‘4 In carrying out the above operations the set of partial differential equations, (2.4.8) to (2.4.11), resulted in a set of ordinary differential equations, (2.5.1). The above equations can be expressed as six, first order, linear differential equations by introducing (21,22,23,24,z5,z6) in place of (0,0,0,00,R&,00). These six equations are = \ DZ1 Z4 1322 = -1(o[z1 + ez3) DZ = z 3 6 2 2 (2.5.3) 024 = [(0 + a ) + 161110149321 + Ram-mm2 + iaZS 2 2 . . 025 = -[Qy 1+ e ) + nyR(U-c)]22 -1(az4 + 826) - 2an1 _ 2 2 . DZ - [05 + a ) + L1R(U-c)]Z + 182 6 3 5 J These can be conveniently written in the matrix form by D[z] = [A][Z] (2.5.4) where matrix {- 'fi Z1 z = z2 23 (2.5.5) Z4 Z5 2 L 6d and matrix r- . fl 0 0 0 1 0 0 A = -ia 0 -1a 0 0 0 0 0 0 0 0 1 b R(DU-20) 0 0 i0: 0 (2.5.6) -212n -b 0 -ioz 0 ~15 0 0 b 0 1B 0 B .J Here 2 2 , b = (a + B ) + 1aR(U - c) (2.5.7) Matrix A is not a matrix of constants but it contains U and DU, which are functions of y. The system of equations (2.5.1 a to d) can also be expressed as two differential equations, one fourth order and one second order. Multiply (2.5.1a) by a and (2.5.1b) by B and add; we get [D2 - (a2 + 62) - iaR(U-c)](afi +-50) = aR(DU-Zfl)v + 1110.2 + 32m, (2.5.8) Differentiate (2.5.8) with respect to y, multiply (2.5.1b) by -iQ12 + 52), add the two and use (2.5.1d). This yields 18 [{Dz-Kz-iaR(U-c)}(D2-K2) + iaR(D2U-20D)]0 + 2K2Rflfi = 0 (2.5.9) where K2 =62 +32. (2.5.10) Again eliminating the pressure term & from (2.5.1a and c) by cross- differentiation, we obtain 2 2 , . .. .. [D - K - 1ozR(U-c)](Bu - aw) = BR(DU - 2mv (2.5.11) Eliminate 0 using (2.5.1d) in (2.5.11) to obtain [02 - k2 - iaR(U-c)](K2'u - 160%) = 521mm - 20%) (2.5.12) Writing (2.5.9) and (2.5.12) in terms of the new variables Z1 to Z6 and rearranging [(Dz-KZ)2 - iaR(U-c)(D2-K2) + iaR(D2U - 201))322 + zxzmzl = 0 P and 2 2 (2513.. 2 K2 . §_ . K: 2 2 [{(D - ) - iozR(U-c)}D - 1 a R(DU - 20)]22 + 1 a—[D -K - and b) iaR(U-c)]zl = 0 J We can also express the equations as a sixth order differential equation by eliminating Z between (2.5.13 a and b) resulting in 1 6 5 4 3 2 _ [F7D + F6D + F50 + F4D + F3D + F2D + F1]Z2 — 0 (2.5.14) where P7 = 1.0 7 F6 = 0.0 F5 = -[3K2 + 2iaR(U-C)] (2.5.15a to f) F4 = -21aR DU F3 = [3K2 + iaR(U-c)][1(2 + mam-6)] . . 3 F2 210R(IZDU + D U) J 19 F1 = MR 040 -K-6 - 216R x“(U-c) + 02R2K2(U-c)2 (2.5.15g) + a2R2(U-c)D2U + 2B2R20(DU - 20) The systems of equations given by (2.5.1), (2.5.3) or (2.5.4), (2.5.13) and (2.5.14) are equivalent. All of these represent the three-dimensional linear stability problem.with the boundary con- ditions discussed in the next section. 2.6 Boundary Conditions The problem outlined in the previous section will be fully defined if the boundary conditions are Specified. The.boundary conditions are such that part of them are known at one boundary and part at infinity. The problem, therefore, falls in the category of a two point boundary value problem. In the study under consideration, velocity perturbations, u', v', w' should vanish at the flat plate, because the particles stick to the plate. Also u', v', w' should tend to zero at large distances from the plate. Equations (2.4.12) then imply that B, 0, 0 I 0 at y = 0 a and (2.6.1) t1, {7, 6:7 =9 0 as y=9co. J Thus the problem of three-dimensional, linear stability is completely defined by equations (2.5.1 a to d) with boundary con- ditions (2.6.1). Or, in terms of the 2's the problem is defined by equations (2.5.3) or (2.5.4) with boundary conditions and (2.6.2) 20 The boundary conditions to go with the form given by (2.5.13 a and b) are 21, 22, DZ2 and (2.6.3) 21, 22, D22 = 0 as y =»m .J = 0 at y = 0 where boundary condition on DZ results from continuity. 2 Finally the boundary conditions that go with equation (2.5.14) are Z2, D22 = 0 at y = 0 and Z2, DZ2 =10 as y = m along with two remaining boundary conditions, one at the plate and one at infinity, which are quite complicated. They can be found by using Z1, Z and D2 = 0 (or 21, Z 2 and DZ2 =»0) in equation 2 2 (2.5.13a). 2 . 7 Eigenvalue Problem The linear stability problem involving rotation is thus given by: l.) the set of four equations (2.5.1) with boundary conditions (2.6.1), in the four unknowns B, 0, 0 and &, or by 2.) a system of six ordinary linear differential equations (2.5.3) with six boundary conditions (2.6.2), in six unknowns Zl’ Z2, 23, 24, 25’ Z6, or by 3.) two higher order but linear differential equations (2.5.13) with boundary conditions (2.6.3), in two unknowns Z and 22, l or by 21 4.) a sixth order, ordinary linear differential equation (2.5.14), in one unknown 22 with four well defined and two involved boundary conditions. Every one of the above mentioned systems with the corresponding boundary conditions, represent an eigenvalue problem with the char- acteristic values, c, a, B, R and .0. The relation among these parameters would result from the Eacular equation. Symbolically, we will have F(C.a.B.R,0) = 0 (2-7-1) Specifying R and .0, c would depend on the real dimension- less wave numbers a and B. Defining c as in (2.4.14) if C1 is positive, one will observe from (2.4.12) that velocity perturba- tions u', v', w' would grow exponentially with time and the motion would become unstable. 0n the other hand for negative Ci’ u', v', u w would decay exponentially—with time and the motion would be stable. With C equal to zero neutral stability would result. 1 From (2.7.1), for a given 0 c = f(a!,B,R) (2.7.2) Since this is a complex equation, 0 ll fi(a:BaR) (2'7°3) and 0 ll fr(01,B,R) (2.7.4) For no growth or decay, fi(a!,B,R) = 0 (2.7.5) 22 For a specified B one could then determine neutral stability curves. These are found by picking one of the eigenvalues, say Cr’ and solving for the other two from equations (2.7.4) and (2.7.5). The neutral stability curve will have a minimum R, called the critical Reynolds number Rc for which a disturbance is neutrally stable. rit’ result in growth for that Reynolds numbers greater than R . cr1t particular disturbance and Reynolds numbers smaller than Rcrit result in decay. Hence one can find an xcrit which separates a region of growth from a region of decay. Associated with Rcrit we also have a (C ) .for each B and each 0. . and an . r cr1t O[c r1t 2.8 Closure If we had considered a two-dimensional disturbance (B = 0), the equations (2.5.14 & 2.5.15) governing the stability of the infinitesimal disturbances (Tollmien-Schlichting waves) would reduce to those governing the stability of flow over a flat plate for which 0 = 0. Thus we see that rotation influences the linear stability problem by interacting with the spanwise component of the three- dimensional disturbances. CHAPTER III NUMERICAL SOLUTION OF THE PROBLEM 3.1 Solution for the Primary Velocity Profile In order to solve the stability problem, presented in Chapter 2, and given by equations (2.5.1), (2.5.3), (2.5.13) or (2.5.14), it is necessary to know the primary velocity profile which is, as shown in Chapter 2, independent of 0. Consider the dimensional boundary layer equations, which are '- '- 2 E§§+v§=05§ (3.1.1) BY and ._ .- 3+ ‘1 = 0 (3.1.2) Bx BY with boundary conditions ._ ._ .. n u = v = 0 at y = 0 and '2 (3.1.3) u.= Uno as y =.m. _ Introducing a transformed stream function and coordinates as f(T\) =' ‘1 (3.1.4) vaon and n = 7/ 2:: U... (3.1.5) with velocity components 23 _=ai_=3iafl= f' u BY B“ BY Um (n) (3.1.6) and -=-ai=_1. I V B52 2 2130f -f) (3.1.7) :2 we obtain f f"‘+ 2f"' = 0 (3.1.8) where a prime denotes differentiation with respect to n. The new boundary conditions are and (3.1.9) f'=1 as 11:00 .J The numerical technique used to solve equation (3.1.8) with boundary conditions (3.1.9) is outlined below. We can express equation (3.1.8) as three first order differ- ential equations as follows: f' = 2 (3.1.10) 2' =‘W1 (3.1.11) 1 I a - _. w1 2 f w1 (3.1.12) The starting conditions are Also some guess must be made about the value of W at n = 0. 1 With this information, and starting from the plate, the computations are carried out in four steps and then the predictions about the 25 values of the functions f, z and "W1 are made at the next step. The scheme used was the fourth order Runge-Kutta. It is outlined in a reference by Weeg and Read. That value of W was accepted, which gave U = 1.0 at the farthest distance from the plate. The value of W1 came out to be 0.332 which gave U = 1.0 at y = 3.29. This corresponds to n 2 16.5. Further details about the computer program are given in appendix C with a list of the program in appendix C-3- In the same program the derivatives of U which appear in the governing equations were computed. 3.2 Approximate Solutions of the Orr-Sommerfeld Equation for Zero Angular Rotation The systems of equations (2.5.1), (2.5.3), (2.5.13) or (2.5. 14) are the different forms of the Orr-Sommerfeld equation. For 0 = 0, consider the set (2.5.13). The first equation yields 4 0 z2 - [2x2 + iaR(U-c)]D222 + [x4 + iQIR K2(U-c) + 101R 020122 = 0 (3.2.1) This is the same equation that Radbill and Van Driest (1966) de- veloped, the only difference being that they had B = 0, so that K reduced to a. At large distances from the boundary (y > 5) we can let U = 1.0 and DZU = 0.0, thereby giving an equation with constant coefficients as 0422 - [2&2 + iaR(1-c)]D222 + [6“ + ia3R(1-c)]22 = 0 (3.2.2) This equation, valid only for large y, permits solutions of the form ery. These approximate solutions then enables us to start 26 our forward integration procedure as we move from the free-stream to the wall. To find these approximate solutions substitute ery into equation (3.2.2), and simplifying, we are left with a poly- nomial of the fourth degree which is r“ - [26.2 + iaR(l-c)]r2 + [a4 + ia3R(1-c)] = 0 (3.2.3) Letting equation (3.2.3) reduces to a quadratic equation in S with roots S1 and 82. The four roots of equation (3.2.3) are then ’iH/EI and 'i,/§;' which would be complex. However, those roots which are admissible in the present case are bounded at infinity. There are only two such roots, r1 = - ./S_1- and r2 = - /-S—; The approximate solution of (3.2.1), therefore, is the linear combination 22 = exp(r1y) +'B exp(r2y) (3.2.4) where B is a constant. To solve equation (3.2.1) with variable U(y) and D2U(y) a numerical scheme, which uses the solution (3.2.4) at large y, outlined in section (3.4) is used. The solution yields the eigen- values a, R and C for Specified B and 0. 3.3 Approximate Solutions of the Orr-Sommerfeld Equation for Non- Zero Angular Rotation For '0 other than zero consider the single equation given by equation (2.5.14) 27 6 5 4 3 2 _ [F7D +F6D +an +1740 +F D +F D +F1]22 -0 (3.3.1) 3 2 where F1, i = 1 to 7, are given by equation (2.5.15). This is a linear, ordinary differential equation of the sixth order with coefficients depending on the known primary velocity U(y) and its derivatives. Equation (3.3.1) again admits solutions of the form ery for large y. Substituting this in (3.3.1) and simplifying yields the polynomial 6 5 4 3 2 F7r +-F6r +-F5r + Far + F3r + Fzr +F1 - 0 (3.3.2) where U = 1.0 and all its derivatives are zero. This gives six complex roots, three of which have negative real parts. These are selected because they are the only roots which would be bounded at infinity. Let those roots be represented by r1, r2 and r3. Then the approximate solution of (3.3.1) becomes 22 = exp(r1y) +'B exp(r2y) + Clexp(r3y) (3.3.3) where B and C1 are constants. Eigenvalues a, R and C with B and 0 specified are again determined from the numerical solution of equation (3.3.1) where the coefficients are functions of y. The solution (3.3.3) is used in the numerical scheme for large y. 3.4 Numerical Integration of the Orr-Sommerfeld Equations (a) Initialization The numerical integration is started from the free stream and carried out step by step until the wall is reached. The initial 28 values of the independent solutions generated, which are combined when the wall is reached, are found from the approximate solutions of sections (3.2) and (3.3). These solutions are $1 = exp(riy) (3.4.1) where i = 1,2 if n = 0 and i = 1,2,3 if n 7‘ 0. The higher order derivatives are found simply by differentiating equation (3.4.1). A distance of y = 1.5 was used as the starting point for the forward stepping integration scheme. Thus, in the case of zero angular rotation, the boundary value problem given by the Orr-Sommerfeld equation (3.2.1) and boundary conditions ¢, D¢ = 0 at y = 0 and (3.4.2) Q5, D0 = 0 for large y .1 changes to an initial value problem where we integrate from the free stream to the wall. For 0 ¥ 0, having found ¢i and all its derivatives, Y1 is computed by (2.5.13a), DYi is calculated by differentiating (2.5.13a) and then DzYi is obtained from (2.5.l3b). These equations are given by 4 2 , 2 . Y = [D ¢ - {2K +-1aR(U-c)}D ¢ - ZLaR.fl D¢ + {K4 + idR K2(U-c) + 101R D2U]¢]/(2K2R (2) (3.4.3) 29 5 2 . 3 . 2 DY = [0 ¢ - {2x 1+ kyR(U-c)}D ¢ - u1R(DU +-zn)D ¢ + {K4 + iaR K2(U-c) + iaR D2U}D¢ + nyR{K200 + D3U}¢]/(2K2R n) (3.4.4) DZY = i g§£03¢ - {K2 + n1R(U-c)}D¢ - iB2R(DU - 20)¢/a] K + {K2 + h1R(U-c)}Y. (3.4.5) Using the derivatives of ¢i from (3.4.1) in (3.4.3) to (3.4.5) we can initialize Vi, DYi and DZYi. Hence, in the case of O f 0, the boundary value problem given by Orr-Sommerfeld equations (2.5.1), (2.5.3), (2.5.13), or (2.5.14) with boundary conditions as discussed in section (2.6) is transformed to an initial value problem. The numerical integration of the equations, started from y = 1.5 towards the plate, utilizes a Runge-Kutta fourth-order method. This is explained in detail in appendix C. The integration proceeds to the wall where the independent solutions are combined to form the eigen-functions. When the correct eigenvalues are used the boundary conditions can all be satisfied. (b) Growth of Eigen-functions During the numerical integration it is observed that for n = 0, one of the two independent solutions grows more rapidly than the other. In the case of rotation, two of the three solutions grow faster than the third one. The functions with the fastest growth are called the "Growing solutions” and the other solutions are called the "Well-behaved solutions”. The growing solution stems from the inviscid part of the Orr-Sommerfeld equation and the 30 well-behaved solution results from the contribution of viscosity. Mathematically, if ¢1 and ¢2 are the growing and the well behaved solutions respectively, then D¢ 2D ...; >12 91 P2 The Orr-Sommerfeld equation can be written as, for Q = 0, (from equation (2.5.13a)) [(U-c)D2 - 22(0-6) - 1320]., = - $02 - x2)2]¢ (3.4.6) and, for {2 f 0 (from equations (2.5.l3a and b)), we have 2 . ‘- [{(U-c) (DZ-K2) - (020-211 mm + 2 1 E- m] = :71; [(DZ-K2)2]6 and _ ' (3.4.7) 2 2 2 U (u-c)1) + Emu-2mm + i(U-c) 1;— Y] = - 3033421));5 + 1 E-(D2-K2)Y]_.J For flows at large Reynolds numbers, which are of interest in this study, in equations (3.4.6) and (3.4.7), terms in the square brackets on the left hand side represent the inviscid part and those in the square brackets on the right hand side, the viscous effects. At the start of the integration i.e. at y = 1.5, all the eigen-functions are linearly independent, but as the integration proceeds, the inviscid solutions grow faster than the well behaved solutions and since the computations are carried out up to a certain number of significant places on the computer, the growing solutions start affecting the viscous solutions. The functions then do not remain independent. To reduce this effect, all the calculations in the present study were made in double precision. Additional 31 steps to make the functions linearly independent all the way to the plate, are explained in the next section. (c) Purification Scheme-Filtering the Contribution of the Growing Solution The problem in the numerical solution is in the numerical integration of the governing equations. For 0 = 0 equation (3.4.6) can be written as Similarly for Q # 0 equations (3.4.7) can be written as 1 _ (L3¢ +-LhY) +-R L5¢ - 0 i. and (3.4.9) 1 _ (L6¢ + L7?) + R(L8¢ + LgY) 0 J where L1,...,L9 are differential operators which can be determined by comparison with equation (3.4.6) or (3.4.7). When R is of the order of 104, the equations (3.4.8) and (3.4.9) become very singular. In any particular case, and in the absence of any rotation, there will be two linearly independent solutions of equation (3.4.8) satisfying the boundary conditions and an appropriate linear combina- tion of these solutions. This general solution must satisfy the remaining two boundary conditions at the plate. Similarly in the presence of a rotation the general solution is an appropriate linear combination of three linearly independent solutions of equation (3.4.9) satisfying the three boundary conditions in the free stream. Here this general solution must satisfy the remaining three boundary conditions at the plate. It is well known that for O = 0 one of 32 the two solutions grows very rapidly as the integration proceeds from the free stream to the plate. Similarly for n # 0 the growth of two out of the three solutions is very rapid. It is difficult to generate the remaining solution which would be inde- pendent because any small round off throws in some small multiple of the growing solution. This grows as the integration proceeds and may later dominate and thus pollute the well behaved solution. Reynolds and Potter (1967) did find this to be so in a solution of Poiseuille flow where they started at the centre with inviscid conditions and proceeded towards the walls. 0n reaching the walls they observed that the well behaved solution was, to eight digits, a multiple of the growing solution and was of the order of 1010 whereas initially it was present to the order of one part in 108. Since the final form of the solution is of the same order as the well behaved solution meaning thereby that a very small amount of the growing solutions is required, the problem is to numerically generate a well behaved solution which does not contain any multiple of the growing solution. In order to generate the linearly independent solutions the scheme employed is to purify the well behaved roots of the contribu- tion of the growing roots. Such a scheme was first used by Kaplan (1964) and later by Lee and Reynolds (1964). In the main program given in appendix C, solutions are arranged in order of their growth. For example the solutions associated with root r1 has the largest growth. After initializing the independent functions and calculating 4 their derivatives, the values of D ¢ and D2? are computed by 33 solving simultaneously equations (3.4.7) for each root for 0 # 0. But for O = 0 D4¢ was computed from equation (3.4.6). Then proceeding forward by the Runge-Kutta fourth order method, these values are found at the next step. Then the effects of the growing solutions are filtered from the behaved solution which are later modified to repair the extractions. The above is schematically given in the flow chart in appendix B. The computer program that does this forms a part of the main program listed in appendix c-3. (d) Iteration Scheme for the Eigenvalues After generating the individual independent functions they should be suitably combined so that for the correct eigenvalues, the boundary conditions at the plate are satisfied. If the guesses on the eigenvalues are off, the correct values are determined by iteration. This is explained below for Q # 0. The combined 0, D6 and Y at the wall are ¢w1 + RUM)”2 + 11am”3 = 1w (3.4.10) 11¢,”1 + X(1)D¢w2 + X(2)D¢W3 = new (3.4.11) Y + X(1)‘Yw2 + 1((2)1IWB = 1’ (3.4.12) W1 W For the right choice of eigenvalues a, R and Cr with fixed 0, B and Ci = 0, the functions QW’ DQW and Yw should be zeros, in order that the no-slip condition and continuity be satisfied. To set up the iteration scheme set 6 and YW equal to zeros. Solve W for X(1) and X(2) from (3.4.10) and (3.4.12). Using these values in (3.4.11) the test function T results, which is = [111W]r T = [Déwji r [(DQW) 1 [(DQW) T 3 one step from wallJr one step from.wall]i 34 2 . giving T = /1:1. + T: ' (3.4.13) -3 - If T s 10 , the eigenvalues are assumed correct. If T > 10 3, let T = T. Increase a by 1% and let T = T. Then increase R l 2 by 1% and let T3 = T. Then T - T 31 e _______Z 1 (3.4.14) on An and T - T T 3 1 SE = T (3.4.15) Au and AR are calculated from the complex equation a: +3: + = 30M 6R AR 11 0. (3.4.16) The new values of a and R are new old and (3.4.17) = + Rnew Rold AR J. with these values of a and R, a new T is computed. If T < T1, then this forms a converging sequence, and the values of a and R would hopefully converge towards the true values. If the initial guess is close, the true values would be reached in about three iterations. Sometimes this might oscillate and then converge. But if the first guess is bad, T may fail to converge. This would be observed by a very slow convergence of T with comparatively large changeSin the values of a or R. For 0 = 0, combined 0 and D0 at the wall are 35 W1 ... X(1) . ¢W2 = 6W (3.4.18) and D¢W1 + x(1) . D¢W2 = Dow (3.4.19) To satisfy one boundary condition at the plate, set QW = 0, solve for X(1) from equation (3.4.18) and using this in (3.4.19) gives the test T. The iteration procedure is the same as discussed for the 0 # 0 case. (e) Criterion for Guessing the Eigenvalues To start the solution, the eigenvalues are first found for the O = 0 case. As first guesses, some values of Radbill and Van Driest (1966) were used. Once some eigenvalues are determined, other guesses are obtained by extrapolation. This is done until a complete neutral-stability curve results. Actually, we are only interested in the critical eigenvalues so when these are found for 0 =‘0, B is assumed equal to 0.2, then taking 0 = 0.01 and -0.01 and taking some of the guesses equal to the zeroed in values corresponding to 0 = 0, new values are obtained. Then by careful guessing the complete stability curves for these values of 0 are completed. Similarly for the other values of’ Q, the true values of a, R, and Cr are found. For B 0.1, 0.3, 0.4, etc., values are guessed from the corresponding B 0.2 and so on. With these, the various curves are completed. CHAPTER IV CONCLUSIONS AND SUGGESTIONS FOR FURTHER RESEARCH 4.1 The Primary Velocity Profile In order to solve the equations governing the growth of the infinitesimal disturbances, it was necessary to have the values of the primary velocity U and its derivatives at all the steps of integration. For this the third order differential equation (3.1.8) was solved. The results for U, U' and U" are given in table 1. The third and fourth derivatives were also computed at the first few points. The values of U checked very closely with those of Howarth (1938). 4.2 Eigenvalues for the Case of Neutral Stability for Zero Angular Rotation The numerical technique used here is the same as with angular rotation. Zero rotation is studied to allow a check on the accuracy of the present technique and to provide information necessary for iterations for other neutral stability curves. The eigenvalues are given in table 4 and plotted in figure 3. In the same figure mean curves due to Schlichting (1935), Lin (1945) and Kaplan (1964) are also shown. There is good agreement.in the lower limit of the neutral stability curves except that Schlichting's curve is a little lower and Kaplan's a little higher than the present calculations show. 36 37 For the upper branch the present curve falls midway between Schlichting's and Lin's curves and slightly below Kaplan's. The difference with Kaplan's curve may be due to not taking U = 1.0 and D20 ' 0 at the start of integration but using the actual values computed from the primary velocity profile described in section 4.1, and also starting at a larger distance from the plate than Kaplan did. Another difference was a larger step size used in the present case. The present calculations were also checked by superimposing the curves of Shen (1954), Kurtz (1962), and the data points of Radbill and Van Driest (1966) (see figure 3). The experimental points due to Schubauer and Skramstad (1947) are also plotted in figure 3. In this connection one should remember the different characteristic length used by each author to non-dimen- sionlize various quantities. Table 2 gives the relation between the eigenvalues used by different authors. The critical eigenvalues for zero rotation as computed by different authors by various techniques are given in table 3 for immediate comparison. In order to compare the eigen-functions ¢ = ¢r + i ¢i and ¢ 0.33156 for the Lower Branch) as taken by Radbill and Van Driest the same values of Cr (0.33052 for the Upper Branch and (1966) correSponding to the experimental data of Schubauer and Skramstad (1947) were taken to compute the values of a and R. Curves of ¢ and ¢' with the corresponding curves of Radbill and Van Driest are shown in figures 4 and 5 for comparison. At this point, for ~O and C1 = 0.0, Squire's theorem (1933) was verified; it requires that: 38 (GR)Two Dim. = (aR)Three Dim. 2 _ 2 2 (a )Two Dim. - (0’ + 6 )Three Dim. (C ) = (C ) r Two Dim. r Three Dim. The verification of this theorem provided a partial check on the numerical calculations. 4.3 Effect of the Angular Rotation Q and Wave Number B The numerical scheme used was employed to investigate the effect on the neutral stability curve of angular rotation 0 and disturbance angle B. Eigenvalues for various 0 and B are given in tables 4 through 17. Neutral stability curves are shown in figures 6, 9, 11 for a verses R and in figures 7, 10, and 12 for Cr verses R. Critical eigenvalues of the neutral stability curves, R (C ) and a . , are lotted a ainst crit’ r cr1t p g crit’ O for constant B in figure 8. The results show (see figure 8) that increasing 0 has a destabilizing effect on the boundary layer and as 0 decreases it stabilizes the boundary layer. The same result was found by Conard (1962) in a study of the effect of rotation on Gbrtler waves. The results also agree with the exper- imental data presented by Halleen and Johnston (1967). Although they did not specifically state it, it is clear from their mean velocity profiles for small and large Reynolds numbers that positive (Idestabilizes whereas negative 0 stabilizes the flow. It is also observed that (C ) . increases raduall with r cr1t g y the increase in 0 at first and later decreases and acrit initially decreases and then increases. These effects are less pronounced for B = 0.1, and much larger for B = 0.3; so it 39 appears that increasing B increases these changes although it may be that there is a critical B for which the smallest critical Reynolds number occurs. The various spanwise wave numbers chosen were B = 0.0, 0.1, 0.2 and 0.3 for O = -0.10, B = 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 for O = 0.10 and B = 0.2 for O = 0.50. It can be observed from a verses R and Cr verses R curves (see figures 16 and 17) that for negative 0 increasing B stabilizes the boundary layer, and from Cr verses R curve (see figure 17) one can also notice that as B increases, Cr decreases. For positive 0 the increase in B destabilizes the boundary layer (see figures 18 and 19). From the correSponding Cr verses R curve (see figure 19) increasing B increases Cr. Figure 20 shows the effect of span- wise wave number B on the critical Reynolds number for angular rotation O = 0.10 and -0.10. The eigen-functions near the critical point on the lower branch of the neutral stability curve for 0 = -O.10 and B = 0.2 are shown in figures 13, 14, and 15, and are also presented in table 18. 4.4 Summary of Conclusions The results obtained in this study can be summarized as follows: 1. In the presence of angular rotation 0, spanwise wave number B must be non zero, if the effect of 0 is to be observed. 2. For a constant B, increasing 0 is destabilizing, whereas decreasing Q has a stabilizing influence on the boundary layer. In particular, for B = 0.2, with O = -0.1 R , = 1897.91, for O = O R , = 1506.8Z for cr1t cr1t 40 n = 0.1 R , = 1347.76 and with 0 = 0.5 R . = 1166.69. cr1t cr1t 3. For a constant 0, increasing B has a destabiLizing effect for positive values of 0 and stabilizing effect for negative values of 0. 4.5 Recommendations for Further Research The present study was conducted to investigate the influence of rotation on the stability of an incompressible fluid undergoing the parallel shear flow, the laminar boundary layer. The angular rotation rate 0 was varied between -0.1 to 0.5 for different values of spanwise wave number B. It is prOposed that future research in this area can be extended along the following lines: 1. To continue the present work to include a greater range of 0 say from -1 to 5 for B varying from 0 to l. 2. To extend the present work for C1 # 0 in order to study the growth and the decay of the three-dimensional disturbances. 3. To restudy the present work and also as suggested in paragraphs 1 and 2 above for a compressible fluid. 4. To conduct research as done in this study and suggested in paragraphs 1 and 2 above for a non-Newtonian fluid. 5. To take into account the temperature variations and thus to study the flow stability when there is a temperature boundary layer superimposed on a velocity boundary layer. 6. To verify the results obtained in the present analytical work by the use of experiments on a flat plate. BIBLIOGRAPHY BIBLIOGRAPHY Collatz, L., "Numerische Behandlung von Differentialgleichungen". Springer-Verlag, Berlin (1951). Conrad, PLW., "The Effects of Rotation on the Stability of Laminar Boundary Layers on Curved Walls". Report FLD9, Mechanical Engineering Department, Princeton University, (1962). Dryden, H.L., "Air Flow in the Boundary Layer Near a Plate". NACA Tech. Rep. 562 (1936). Dryden, H.L., "Some Recent Contributions to the Study of Transition and Turbulent Boundary Layers". NACA TN 1168 (1947). Also, "Recent Advances in the Mechanics of Boundary Layer Flow". Adv. hppl. Mech. I, Academic Press, (1948). Emmons, H.W., "The Laminar-Turbulent Transition in 8 Boundary Layer". Part I, Journal of Aeronautical Sciences, 18 (1951). Gill, 8., "A Process for the Step-by-Step Integration of Differential Equations in an Automatic Digital Computing Machine". Proc. Camb. Phil. Soc., 41 (1951). Halleen, R.M. and Johnston, J.P., "The Influence of Rotation on Flow in a Long Rectangular Channel - An Experimental Study". Report MD-18, Thermosciences Division, Department of Mechanical Engineering, Stanford University, Stanford, California, (1967). Heisenberg, W., "Uber Stabilitfite und Turbulenz von Flfissigkeitsstrfimen'. Ann. d. Phys., Lpz., (4), 14, (1924). Hilderbrand, F.B., "Introduction to Numerical Analysis". McGraw-Hill Book Co., Inc., New York (1956). Howarth, L., "On the Solution of the Laminar Boundary Layer Equations". Proc. Roy. Soc., London, A 164, (1938). Kaplan, R.E., "Solution of Orr-Sommerfeld Equation for Laminar Boundary Layer Flow over Compliant Boundaries". ASRL-TR-ll6-1, Cambridge Aeroelastic and Structures Research Lab. Report, M.I.T., (1964). Kelvin, Lord, Philosophical Magazine. 5th Series, t. xxiv (1887). Kurtz, E.F., Jr. and Crandall, S.H., "Computer-Aided Analysis of Hydrodynamic Stability". J. of Math. and Physics, XLI, (1962). 41 42 Lee, L.H. and Reynolds, W.C., "A Variational Method for Investigating the Stability of Parallel Flows". Tech. Rept. No. FM-l, Mech. Engineering Dept., Stanford University, Stanford, Calif., (1964). Lee, L.H. and Reynolds, W.C., "0n the Approximate and Numerical Solution of Orr-Sommerfeld Problems". Lin, C.C., "0n the Stability of Two-Dimensional Parallel Flows". Parts I, II, III. Quart. Appl. Math. 3, (1945, 1946). Lin, C.C., "The Theory of Hydrodynamic Stability". Cambridge University Press, London (1966). Milne, W.E., "Numerical Solution of Differential Equations". John Wiley & Sons, Inc., New York (1953). Morkovin, MLV., "Transition from Laminar to Turbulent Shear Flow - A Review of Some Recent Advances in its Understanding". Trans. ASME, (1958). Nachtsheim, P.R., "An Initial Value Method for the Numerical Treat- ment of the Orr-Sommerfeld Equation for the Case of Plane Poiseuille Flow". NASA TN D-2414, (1964). Nikuradse, J., "Experimentalle Untersuchungen zur Turbulenzentstehung". 22.1114, 12 (1933). Orr, W.M'F., "The Stability or Instability of the Steady Motions of a Perfect Liquid and of a Viscous Liquid". Proc. Roy. Irish Academy, (A) 27 (1907-09). Pekeris, C.L., "Stability of the Laminar Parabolic Flow of a Viscous Fluid between Parallel Fixed‘Walls". Physical Review (2), 74, (1948). Potter, M.C., "The Stability of Plane Couette-Poiseuille Flow". Doctoral Dissertation, University of Michigan (1965). Later published in J. Fluid Mech. (3), 24_(l966). Potter, M.C., "Linear Stability of Symmetrical Parabolic Flows". The Phy. of Fluids (3), 19_(1967). Prandtl, L.P., "Uber den Luftwiderstand von Kugeln". Gottinger Nachrichter (1914). . Prandtl, L.P., "Bemerkungen Uber die Entstehung der Turbulenz". ZAMM; l (1921), and Phys. 2., 2§_(1922). Radbill, J.R. and Van Driest, E.R., "A New Method for Prediction of Stability of Laminar Boundary Layers". North American Aviation, Inc., Anaheim, California (1966). Rayleigh, Lord, "0n the Stability or Instability of Certain Fluid Motions". Proc. Lond. Math. Soc., ll,_§l (1880) and_lg,_gl_ (1887). 43 Rayleigh, Lord, "Further Remarks on the Stability of Viscous Fluid Motion". Phil. Mag. and Journal of Science, 6 Series, 2§_(19l4). Reynolds, 0., "An Experimental Investigation of the Circumstances which determine‘whether the Motion of Water shall be Direct or Sinuous, and of the Law of Resistance in Parallel Channels". Phil. Trans. Roy. Soc. (1883). Reynolds, 0., "On the Dynamical Theory of Incompressible Viscous Fluids, and the Determination of the Criterion”. Phil. Trans. Roy. Soc., 186 (1895)A. Reynolds, W.C. and Potter, M.C., "Finite-amplitude Instability of Parallel Shear Flows". J. Fluid Mech. (3), 21 (1967). Schlichting, H., "Zur Entstehung der Turbulenz bei der Platten- sterung". Nachr. GeS. Wiss. G6tt.,Math. Phys. Klasse, (1933). Also, ZAMM, 13 (1933). Schlichting, H., "Amplitudenverteilung und Energiebilanz der Kleiner Storungen bei der Plattenstromung". Nachr. Ges. Wiss. Gottingen, Math. Phys. KlasseBd I (1935). Schlichting, H., "Boundary Layer Theory". McGraw-Hill Book Co., Inc., New York (1962). Schubauer, G.B. and Skramstad, H.K., "Laminar Boundary Layer Oscillations and Transition on a Flat Plate". National Bureau of Standards Research Paper 1772. Also, "Laminar Boundary 1 Layer Oscillations and Stability of Laminar Flow". J. Aero. Sgi,, 14 (1947). Shen, S.F., "Calculated Amplified Oscillations in the Plane Poiseuille and Blasius Flows". J. Aero. Sci., 21_(l954). Sommerfeld, A., "Ein Beitrag Zur Hydrodynamischen Erklaerung Der Turbulenten Fluessigkeitsbewegungen". Atti del Congg. Internat. dei Mat., Rome (1908). Squire, H.B., "0n the Stability for Three-Dimensional Disturbance of Viscous Fluid Flow between Parallel Walls". Proc. Roy. Soc. London, A 142, (1933). Streeter, V.L., "Handbook of Fluid Dynamics". McGraw-Hill Book Co., Inc., New York (1961). Synge, J.L., "Hydrodynamic Stability". Semi-centennial publications of the Amer. Math. Soc. 2 (Addresses)(1938). Thomas, L.H., "The Stability of Plane Poiseuille Flow”. Physical Review (4),‘21 (1953). Tietjens, 0., "Beitrfige zur Entstehung der Turbulenz". ZAMM, 2 (1925). 44 Tollmien, W., "The Production of Turbulence". NACA Tech. Memo. 609 (1931), originally published in Nachr. Ges. Wiss. GBtt., Math. Phys. Klasse (1929). Tollmien, W., "General Instability Criterion of Laminar Velocity Distributions". NACA Tech. Memo. 792 (1936), originally published in German in (1935). Weeg, G.P. and Reed, G.B., "Introduction to Numerical Analysis". Computer Laboratory, Michigan State University, East Lansing, Michigan. Zurmfihl, V.R., "Runge-Kutta-Verfahren Zur Numerischen Integration von Differentialgleichungen n-ter Ordnung". Z. FUr A» Math. und Mech., 2§_(1948). I LLUSTRATIONS (Tables and Figures) 45 Table 1. Values for the Primary Velocity Profile y U U' U" 1.500 0 999981858 0 000276314 - -0.003992197 1.490 0.999978886 0.000319064 -0.004569974 1.480 0 999975456 0.000367969 -0 005224437 1.470 0.999971503 0.000423839 -0 005964703 1.460 0.999966953 0 000487582 -0.006800817 1.450 0 999961722 0.000560211 -0.007743826 1.440 0 999955716 0 000642854 -0.008805861 1.430 0.999948828 0.000736768 -0.010000208 1.420 0 999940938 0.000843347 -0 011341399 1.410 0 999931913 0.000964138 -0.012845291 1.400 0 999921602 0.001100852 -0.014529150 1.390 0 999909837 0.001255382 -0 016411745 1.380 0.999896429 0 001429816 -0.018513429 1.370 0.999881167 0.001626453 -0.020856227 1.360 0 999863817 0 001847822 -0.023463924 1.350 0.999844119 0.002096699 -0.026362145 1.340 0 999821781 0 002376124 -0.029578440 1.330 0 999796483 0.002689425 -0.033142352 1.320 0.999767868 0.003040234 -0 037085497 1.310 0.999735541 0.003432510 -0 041441616 1.300 0.999699065 0.003870560 -0.046246638 1.290 0.999657961 0.004359064 -0 051538719 1.280 0 999611699 0.004903091 —0.057358276 1.270 0 999559697 0.005508128 -0.063748003 1.260 0.999501314 0.006180099 -0.070752875 1.250 0 999435850 0.006925391 -0.078420134 1.240 0.999362539 0.007750872 -0.086799250 1.230 0.999280541 0.008663918 -0 095941871 1.220 0.999188943 0.009672431 -0.105901736 1.210 0.999086747 0 010784860 -0.116734578 1.200 0.998972869 0.012010221 -0.128497983 1.190 0.998846134 0.013358116 -0.141251231 1.180 0 998705265 0.014838744 -0.155055105 1.170 0 998548881 0.016462923 -0.l69971660 1.160 0.998375490 0.018242093 -0.186063966 1.150 0.998183483 0.020188331 -0 203395809 1.140 0.997971125 0.022314351 -O.22203136l 1.130 0 997736552 0.024633514 -0.242034810 1.120 0.997477764 0 027159816 -0 263469950 1.110 0 997192617 0.029907891 -0.286399741 1.10:) 0 996878817 0.032892994 -0.310885823 1.0 9o 0. 996533915 0 . 036130990 -0 . 336988001 1.080 0. 996155300 0.039638328 -0.3647 63695 1.0713 0.995740194 0 043432019 -0 394267350 1-060 0. 995285647 0.047529600 -0.425549825 1 050 0. 994788530 0.051949096 -O.4586S7749 1-0 0.01() 0.0cu3 0.717446605 0.705652885 0.693605706 0.681310375 0.668772561 0.655998261 0.642993787 0.629765735 0.616320958 0.602666544 0.588809782 0.574758143 0.560519242 0.546100817 0.531510700 0.516756784 0.501847001 0.486789294 0.471591587 0.456261764 0.440807640 0.425236938 0.409557269 0.393776104 0.377900760 0.361938374 0.345895888 0.329780031 0.313597303 0.297353958 0.281055995 0.264709139 0.248318835 0.231890236 0.215428194 0.198937252 0.182421637 0.165885254 0.149331679 0.132764161 0.116185609 0.099598598 0.083005362 0.066407792 0.049807438 0.033205504 '0.016602852 0.000000000 47 1.166531332 1.192130354 1.217217027 1.241754569 1.265708403 1.289046298 1.311738478 1.333757731 1.355079488 1.375681894 1.395545859 1.414655088 1.432996106 1.450558254 1.467333679 1.483317309 1.498506804 1.512902513 1.526507397 1.539326960 1.551369161 1-562644319 1.573165009 1.582945956 1.592003921 1.600357584 1.608027420 1.615035584 1.621405780 1.627163145 1.632334124 1.636946351 1.641028533 1.644610339 1.647722285 1.650395638 1.652662311 1.654554774 1.656105968 1.657349221 1.658318180 1.659046742 1.659568999 1.659919185 1.660131633 1.660240739 1.660280938 1.660286681 .584241153 .534915513 .481805906 .425126001 .365103575 .301978964 .236003444 .167437585 .096549572 .023613526 .948907829 .872713480 .795312487 .716986311 .638014382 .558672681 .479232417 .399958797 .321109894 .242935624 .165676837 .089564518 .014819103 .941649912 .870254700 .800819314 .733517461 .668510583 .605947835 .545966151 .488690409 .434233671 .382697510 .334172393 .288738149 .246464470 .207411485 .171630363 .139163964 .110047512 .084309313 .061971471 .043050646 .027558809 .015504013 .006891171 .001722839 .000000000 48 Table 2. Definition of Eigenvalues used by Different Authors Author Characteristic Length Eigenvalues equal to those used in the present study a R C r * Present 6 = 3 6 a R C Calculations r _ * 3 2_ Kaplan 6 - 3;5 6 3.5 ak 3.5 Rk (C )k . _ 6 3 3 L1“ 5 ‘ 0.28673 3.4876 “1. 3 4876 R1. (Cr)L Radbill and * Van Driest 6 3 ar 3 Rr (Cr)r * Schlichting 6 3 a 3 R (C ) s s r s Schubauer and * Skramstad 6 3 ass 3 Rss (Cr)ss Here 6 = Boundary Layer Thickness * 6 = Displacement Thickness The subscripts k,];, r, 8, ss refer to the corresponding author's value. Table 3. Comparison of Critical Eigenvalues for Zero Rotation Author Critical Eigenvalue equivalent to the one derived in present study for neutral stability for Q = 0.0 R a C r Present study 1506.87 0.87553 0.396 Kaplan 1546.28 0.94285 0.40 Lin 1264.48 1.11549 0.411 Radbill and Van Driest 1548.75 0.89307 0.395 Schlichting 1725 0.834 0.42 Tollmien 1260 1.101 0.425 49 Table 4. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) {l = 0.00, B = 0.00 Cr R a Remarks 0.31000 8822.88 0.83457 Upper Branch 0.33000 6099.27 0.90570 1 * 0.33052 6042.47 0.90743 0.35000 4269.23 0.96918 0.37000 3004.11 1.01801 0.39000 2082.86 1.03359 0.39500 1876.92 1.02365 0.39800 1745.12 1.00871 0.39900 1694.89 0.99753 0.40000 1632.00 0.98411 0 39900 1516.99 0.91674 0.39600 1506.87 0.87553 Critical Point 0.39300 1517.60 0.84520 0 39000 1537.65 0.81961 0.38000 1636.50 0.74971 0.36000 1936.01 0.64408 0.34000 2368.45 0.56093 * 0.33156 2599.29 0.53015 0.33000 2645.57 0.52468 0.32000 2973.92 0.49131 0.30000 3826.21 0.43135 V 0.27000 5848.15 0.35514 Lower Branch * These Cr values correspond to the data of Schubauer and Skramstad and also discussed by Radbill and Van Driest. 50 Table 5. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C. = 0.0) n - -0.01, B = 0.20 1 Cr R a Remarks 0.380 2547.46 1.00998 Upper Branch 0.395 1894.79 0.99681 “ 0.399 1602.89 0.91967 0.397 1573.87 0.87877 0.395 1573.80 0.85301 Critical Point 0.392 1588.38 0.82187 0.390 1603.66 0.80397 0.380 1714.79 0.73055 0.360 2056.52 0.62051 0.340 2561.35 0.53372 0.320 3290.80 0.46055 0.300 4360.18 0.39721 lower‘granch Table 6. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) O = -0.10, B = 0.00 Cr R a Remarks 0.360 3584.44 0.99753 Upper Branch 0.380 2517.68 1.03414 “ 0.385 2297.35 1.03720 0.399 1509.12 0.91266 0.396 1502.38 0.87287 Critical Point 0.394 1509.57 0.85276 0.390 1534.40 0.81763 0 385 1579.23 0.78084 0.382 I 1611.08 0.76093 0.380 1634.22 0.74843 0 360 1933.97 0.64310 0.340 2366.12 0.56011 LoweruBranch 51 Table 7. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) a = -0010, B = 0010 Cr R a Remarks 0.360 3546.84 0.98813 Upper Branch 0.380 2487.38 1.02115 0 0.385 2265.57 1.02229 0.397 1605.28 0.92642 0.396 1589.35 0.90659 0.395 1583.74 0.89174 Critical Point 0.393 . 1584.21 0.86767 0.390 1598.88 0.83834 0.385 1640.87 0.79870 0.382 1673.20 0.77765 0.380 1696.52 0.76465 0.370 1837.18 0.70702 0.360 2012.64 0.65812 9 0.340 2475.87 0.57647 Lower Branch Table 8. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) C) = -0.10, B = 0.20 Cr R a Remarks 0.360 3414.36 0.95491 Upper“Branch 0.380 2355.89 0.96814 0.385 2091.83 0.94831 0.386 1917.42 0.89783 0.384 1897.91 0.86574 Critical Point 0.383 1899.82 0.85411 0.380 1921.51 0.82578 0.375 1983.35 0.78857 0.360 2263;17 0.70519 0.340 2814.62 0.62342 LowerwBranch 52 Table 9. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) 0 = -0.10, B = 0.3 Cr R a Remarks 0.340 4321.70 0.83378 Upper Branch 0.350 3585.03 0.84277 A 0.353 3352.32 0.83768 0.354 3258.25 0.83201 0.354 3136.17 0.80835 0.353 3132.25 0.79574 Critical Point 0.350 3189.73 0.77222 0.340 3529.45 0.72037 Lower Branch Table 10. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) O = 0.01, B = 0.20 C R a Remarks 0.380 2575.67 1.01568 Upper Branch 0.395 1931.37 1.01037 “ 0.400 1528.68 0.89431 0.397 1522.26 0.85194 Critical Point 0.395 1529.97 0.83026 0.392 1549.64 0.80212 0.390 1566.49 0.78532 0.380 1678.94 0.71390 0.360 2012.36 0.60334 0.340 2501.23 0.51361 “ 0.320 3208.42 0.43576 Lower Branch 53 Table 11. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C = 0.0) n = 0.10, a = 0.10 1 Cr R a Remarks 0.360 3637.51 0.99358 Upper Branch 0.380 2550.08 1.03241 “ 0.395 1915.97 1.02993 0.400 1470.11 0.89652 0.398 1467.13 0.87058 Critical Point 0.390 1509.61 0.79669 0.385 1555.36 0.76031 0.380 1610.19 0.72806 0.360 1904.46 0.62120 0.340 2325.14 0.53542 LowerUBranch Table 12. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) O = 0.10, B = 0.20 Cr R a Remarks 0.360 3815.29 0.98934 Upper Branch 0.380 2669.23 1.03518 “ 0.395 2027.86 1.04605 0.408 1514.59 0.99625 0.406 1349.66 0.87664 0.405 1347.76 0.86270 Critical Point 0.400 1359.86 0.80694 0.398 1370.38 0.78813 0.390 1429.58 0.72321 0.380 1529.96 0.65490 - 0.360 1799.36 0.53817 0.340 2157.41 0.43102 0.320 2552.18 0.30861 V 0 .310 2610.45 0 . 24343 Lower Branch 54 Table 13. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) O = 0.10, B = 0.30 Cr R a Remarks 0.360 4130.97 0.97924 Upper Branch 0.380 2863.98 1.03502 “ 0.395 2185.14 1.05851 0.415 1174.49 0.82865 0.413 1172.76 0.79993 ' Critical Point 0.410 1176.75 0.76225 0.400 1215.53 0.65587 0.398 1225.59 0.63611 0.390 1267.62 0.55607 0.380 1278.87 0.41023 Lower “Branch Table 14. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) O = 0.10, B = 0.4 Cr R a Remarks 0.360 4651.27 0.95944 Upper Branch 0.380 3147.21 1.02830 P 0.395 2388.88 1.06347 0.405 1992.76 1.07599 0.410 1818.64 1.07766 0.420 1505.08 1.06666 0.430 1204.89 1.01011 0.432 998.49 0.84961 0.430 970.33 0.78864 0.428 953.66 0.73934 0.427 946.79 0.71588 0.425 932.32 0.66718 0.423 911.27 0.60709 0.4225 900.88 0.58400 0.422 877.11 0.54589 55 Table 15. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C1 = 0.0) 0 = 0.10, B = 0.5 Cr R a Remarks 0.430 1397.06 1.08694 Upper Branch 0.440 1145.52 1.04910 0 0.445 1019.47 1.00449 0.450 840.74 0.86321 0.453 680.46 0.62868 0.455 649.15 0.58823 0.457 625.03 0.55958 0.460 595.32 0.52685 0.462 578.07 0.50902 0.463 569.99 0.50091 0.465 554.66 0.48590 0.470 519.77 0.45346 0.475 487.52 0.42516 Table 16. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (Ci 3 0.0) 0 = 0.10, B = 0.6 Cr R 0 Remarks 0.440 1304.46 1.10778 Upper Branch 0.450 1083.42 1.08000 “' 0.460 882.96 1.01217 0.470 680.09 0.84493 0.473 626.43 0.77933 0.475 597.04 0.74155 0.477 572.17 0.70953 0.480 541.15 0.67013 0.483 515.33 0.63799 0.485 500.19 0.61955 0.490 467.40 0.58049 0.495 439.59 0.54843 0.497 429.45 0.53700 0.500 415.06 0.52093 0.510 371.87 0.47394 56 Table 17. Eigenvalues for the Orr-Sommerfeld Equation for the Case of Neutral Stability (C. = 0.0) n = 0.50, a = 0.20 1 Cr R a Remarks 0.360 3789.20 1.03303 Upper Branch 0.380 2640.32 1.08096 0 0.400 1819.80 1.09706 0.410 1184.75 0.94655 0.400 1166.69 0.85660 Critical Point 0.395 1177.47 0.82488 0.390 1193.60 0.79741 0.380 1236.30 0.75090 0.370 1289.69 0.71261 4 0.360 1351.71 0.68037 lower Branch 57 666+66H.H- ~66--H.~ ~66-66~.H- 666+N66.H- 666-666.N 666+~em.a 6666.6 666+RRH.H- ~66-H66.- N66-666.H- 666+6~6.H- 666-mmm.~ 666+Hmm.a 6666.6 666+~6~.H- ~66-~64.H ~66-mme.-- 666+666.H- 666-66~.~ 666+6~m.a 6666.6 666+6NN.H- ~66-666.H ~66-Ham.a- 666+~6H.H- 666-mam.a 666+e6m.a 6666.H 666+mm~.H- ~66-MN6.H ~66-mea.a- 666+6HH.H- 666-6~e.a 666+666.H 6666.H 666+66~.H- N66-~6N.H ~66-666.H- 666+6~H.H- 666-mom.a 666+~66.H 6666.5 666+66N.H- ~66-H-.H 666-~6~.6- 666+H~H.H- 666-com.a 666+666.- 6666.H 666+66N.H- 666.666.6 666-464.6- 666+6HH.-- m66-~ma.a 666+6H6.H 6666.H 666+6m~.H- 666-666.6 666.66H.e- 666+6aa.a- 666-464.6 666+666.H 666H.H 666+mm~.H- 666-Hmm.~ 666-66~.6- 666+e6a.a- 666-466.6 666+6em.- 66~H.H 666+6~N.H- 666-666.6 666-666.m- 666+666.H- 666-66N.e 666+Hmm.a 666H.H 666+mH~.H- m66-666.m 666-e66.6- 666+666.H- 666-64~.6 666+6~m.- 666H.H 666+H6~.H- 666-666.6 666-H66.6- 666+646.H- 666-~66.m 666+e66.a 666H.H 666+R6H.H- 666--N.6 666-666.m- 666+666.a- 666-666.6 666+66N.H 666~.H 666+mkH.-- 666-666.m 666-466.m- 666+666.H- 666-666.m 666+m6N.H 66-.H 666+46H.H- 666-66H.m 666.6H6.N- 666+666.H- 666-H66.m 666+66N.H 666~.H 666+H6H.-- 666-66e.~ 666-66~.~- 666+-6.H- 666-666.~ 666+mNN.H 666~.H 666+m~H.H- 666-6~6.N 666-5e6.-- 666+e66.a- 666-H~m.~ 666+66~.H 666~.H 666+66-.H- 666--66.- 666-m~e.a- H66-RH6.6- 666-Nma.m 666+66H.H 666m.H 666+466.H- 666-e66.a 666-6~6.H- H66-66e.6- 666-6N6.H 666+66H.H 66~m.- 666+646.H- 666-6H6.- 666-~mm.a- ~66-NH6.6- 666-~em.a 666+665.- 6666.4 666+666.-- 666--6~.H m66-6a~.-- H66-666.6- 666-66N.H 666+6~H.a 6666.4 666+H66.H- 666-~66.H 666-H6H.H- H66-466.6- 666-666.H 666+66-.H 6666.H 666+6~6.H- m66-H6H.H 666-666.H- H66-mma.6- 666-666.6 666+666.H 6666.5 666+666.H- 666-666.H 666-66N.6- H66-666.6- 666-666.6 666+646.- 66N6.H H66-6H6.6- 666-6-.N 666-666.6- H66-566.6- 666-mme.6 666+~mo.H 6666.4 H66-66~.6- 666-66m.m 666-6N6.e- H66-6He.6- m66-~6e.m 666+666.H 6666.H H66-Hem.6- 666--H6.6 466-46H.e- H66-666.6- 666-666.H 666+e~6.- 6666.H H66-m~6.6- 666-666.e 666-666.~- H66-666.6- ~66-me-.H 666+666.H 6666.H :- 63 we we .6 we .... 646N6.6 u 6 . ~6.H~6H u a . 6~.6 u 6 . 6m.6 u 66 6H.6- u c 666 .A6.6 u woo 6>666 acaawcrum Hmuusoz BLu mo noamum “630A Baa co ucHom Hmowuwno n Home maofiuoanwusowwm .wH mHan 58 000+00H.m 000+50N.m 000+0mm.m 000+m¢~.m 000+00H.m 000+H0m.0 000+m00.¢ 000+00m.0 000+0N0.¢ 000+000.m 000+00m.m 000+0mm.~ 000+Nm0zw 000+m~m.m 000+500.N 000+005.H 000+000.H 000+00N.H 000+N00.H H001000.0 H00105m.m H00100m.m H001500.N N00:0H0.N H001000.Hn H001000.Nu H00:H5N.¢u Hooummm.mn H00:000.0u H00:050.5u H00100m.0n H001000.mu 000+000.Hu 000+000.Hu 000+50H.Hu 000+0N0.Hn 000+000.Hu 000+mmm.Hn 000+00H.H1 H001000.01 H00:000.0u H0015mm.0n H0010H0.Nu H001H0N.Hu «001005.01 ~00:0N0.H N0010NH.¢ ~0010H5.¢ N0010N~.¢ ~0010Hm.m N00:¢~¢.N N00:¢55.H N001H00.H N0015mm.H 5001000.H N00:000.H N00150H.~ N00:000.N «00:005.~ N00|m0m.~ 5001000.m N001~0H.m ~00100H.m ~0010NH.0 N00:0¢0.0 N0015mm.~ ~001a00.~ N0015¢0.N ~00105¢.N N00:000.~ H0010H¢.m H001000.m H001H0¢.0 H001Nm5.0 H001050.N H00150N.N H001000.H H00-000.H ~00105m.0 ~001H00.m ~001H00.~ N001000.H 0001500.¢ 0001000.H 000:00n.~ 0001H00.Hu 000:00¢.Mu 000:000.0u 000100H.mu N00100N.Hn N00-000.Hn N00:000.Hu ~001000.~u N00:00~.~u N00:00¢.~u N00:000.Nu 500-050.01 N00100m.~u «001000.51 ~001000.Nu N00:000.Nu N00:000.Nu N00:05H.~u N001000.~1 N001000.Hu 000+00H.¢ 00¢+0NH.¢ 000+0N0.¢ 000+~m0.0 000+0~5.m 000+0mm.m 000+NON.m 000+000.m 000+000.~ 000+HNO.~ 000+N00.N 000+0¢H.~ 000+0N0.H 00¢+0m0.a 000+~00.H 000+¢5N.H 00¢+N50.H H00:H05.0 H0010Hm.0 H0010HH.0 H00100¢.m H00:005.H N0010¢0.N H00:HOH.Hu H00100¢.~u H001005.01 H0010a0.0a H0010N0.mu H001005.0a H00:000.5n H00:¢¢N.0n H001H00.0n a00:amm.0u Hooumm0.mu 000+HN0.H1 ~001000.0u N001~00.Nn N001000.Hu N00100H.Hn 000-000.01 000-050.N 0001N00.0 000-050.0 000100m.0 000-00m.m 0001000.0 N00:0N0.H N00:000.H N00:0¢0.H ~00:5¢0.H N001000.H ~001H00.H ~00|N00.~ N0015H0.H 0001Hmm.m 0001~50.0 0001H00.m 0001500.0 0001000.0 H0010N0.0 000:0Nm.5 000-0H0.5 000-000.0 0001000.m 000105¢.m 000-00m.¢ 000:5Hn.0 0001000.0 0001000.m 000100N.0 H001-0.5 H00:000.0 H00100¢.m 000+0N0.H 000+NOH.H 000+05H.H 000+¢¢N.H 000+000.a 000+500.H 000+~N0.H 000+N50.H 000+0H0.H 000+000.H 00¢+00m.a 000+0~0.H 000+0m0.H 000+550.H 000+500.H 000+NH5.H 000+0~5.H 000+005.H 00¢+005.H 00¢+0¢5.H 000+005.H 00¢+mm5.a 000+m~5.H 000+HN5.H 000+0H5.H 000+0m0.H 000+000.H 000+500.H 000+000.H 000+~m0.a 000+0H0.H 000+000.H 000N.0 0000.0 000N.0 0000.0 00mm.0 0000.0 0000.0 0000.0 0000.0 00N0.0 0000.0 0000.0 0000.0 0000.0 00mm.0 0000.0 0000.0 0000.0 0000.0 00N0.0 0000.0 0000.0 0000.0 0005.0 00~5.0 0005.0 0005.0 0005.0 0000.0 0000.0 0000.0 0000.0 0000.0 0000.0 00N0.0 .6.0666 59 0H010-.Nn 000+0NH.H 000+550.N 000+0H0.N 000+mmm.m 000+005.m 000+mH0.0 000+0m~.0 000+~00.0 000+N00.0 000+H00.0 000+000.m 000t000.0 H00:~00.0 000+000.H H00100m.m H00-000.0 H001000.~ H00105m.au H0010m~.0n 000+000.Hn 000+0H0.Hu 000+0Nm.aa 000+0N0.Hn 000-000.01 H00:00H.5u H001500.mn H0010NH.01 H00:00H.5n H001~N0.0u H001000.Ha N001005.m H0015mm.m H00:0N0.0 H001000.0 H001000.m 0001005.~n 000+N00.H 000+000.H 000+N00.N 666+66H.m. 00¢+00m.m 000+0H0.m 000+~mm.0 000+00H.0 000+05H.0 00¢+aa~.0 000+0HN.0 5H0:H00.H 000:05H.0u N0010mm.~u N00:000.0n N00100H.0n N00:05N.5u N00:0H0.5n «001000.01 ~0015H5.5n ~001000.5n «00:50Hn0n N00-mmo.mn 000+000.0 N001000.H N00:050.0 ~0010H5.0 H001000.H H00150H.~ H00:550.~ H001000.m H001050.0 H00100N.m H00100H.0 H001H0m.0 0000.0 0000.0 0000.0 0000.0 0000.0 000H.0 00NH.0 000H.0 000H.0 000H.0 000N.0 00~N.0 .6-6660 Axis of Rotation 60 IN 6 Boundary Layer Thickness Flat Plate Figure 1. Primary Flow ¢, D¢, Y = 0.0 as y a 0 y H ———- U(y)=>1.0 —————— as y=0 (Numerically at y = 3.29) *——-—‘ Initial Conditions 1.5 . Numerical ._____.. d 1.0 integrations 0 /Edge of the Boun ary Layer 91’92’93 1 U(y)' Boundary Conditions ¢ = D¢ = Y = 0.0 0 I'II/f/ll/f/lllff/IIIIflli Flat Plate Figure 2. Numerical Integration of the Orr-Sommerfeld Equation. . 051N411... ... 0.0 u G How oumam umam m Hw>o uo5mg 5um0c3o0 630mmam map 000 mm>uso zufiaflsmom Hmwuooz 050 no acmfiumaaoo .0 005000 0000 0000 0005 0000 0000 m 0000 0000 000m 000a .r 1 a I - m I _ N... 0mmH swam IIIII I. . Aawuamawuomxmv a...llllI.JrnlMW.HHhHhH R . I 50mH vaumEmuxm 0 umsmnsnom on. ..... II Ilu.urlnllu. . -I . 68H wean-€266 l..l.l .. - 6 6 663 662.6 66> .2362- a 663 63 oIIII iiii II r N0mH Nuuax 66$ 6263-6IEI uneaumanoamo 3035 III .l).0 c unawmq a 300.5 wanna—60x0 05m non-602000 .. / . anewuwanoamo unmmmum 0 III. 0 .\.., 1| 0 .qolo \ p / $4.0 I/ 7:16,: $660 wan-r6366... Ii. |\. / Ilowl o./o. u I [6, o /.. V L 3.00.3 nuusx IIIIIII./OIOI.IICI . /0../ o ... / IU’I/ .L/./ / ../QI a/ II, I. o o to llloofi CU III 330 .5 III /./... /-.. 66. I.F#.ocou ‘6 .o _\ 3660 6266M . \ 1 / . $660 63 Ike/II .\N 0.0 n 0 Infill a [4,4, \ AIN.H IIIHHMU\ n0 -0 62 00.0" 6 HOW 50.0" H00 0>Hso 5uHHHnmum Hmuusmz 0:0 00 mmnoomum H0304 0:0 H0000 0:0 How wmumamuxm 0cm Hmsmnnzomu Ho moan 000 cu wcchomm0HHoo .um0HH0 am> 0cm HHHnumm mo 0monu suH3 He H + us" 0 mcoHuocsmI Ic0me mo oomHHmnEoo .0 0H30Hm % H N.CI H 6 6.H 6.H 6.6 6 \.II. 6 H w H H . H m H H 1 O . . \\ \\ 0. I \\ \\ \\ L .mcoHumHsoH-wo uc0m0ua 05 Ho 0050 mm \\ I 06-00 0:0 3 um0HH0 96> .0 HHHnwmm .How \\ H6 m0soomun H09? 000 H033 noon .Hom H9 "00oz \\ I x... \\ .6 If O.H II, . ”I” 66666.6 I a .66.~666 I a 747, . . a . H // 66 6 I 6 66666 6 I 6 77 526000 .0000 III /// um0Hun 65> 0:0 HHHn—cmx \\ 6H6mm.6 I 6 .66.666~ I 6 own, 6 .\ 66.6 I 6 .66Hmm.6 I we xnnnx/ 6.IIJV\.\\\\\ soc-6H0 .033 III /// \ / / l||\\ \ / l I‘ \\ \ 666666 66> 666 HHHceee VI, .“ 2) Iv‘.‘§‘ r 6.6 I c 666 H6.6 I H66 6>666 huHHHnmum Hmuunmz may mo mwnucmum umsoq cam Hmam: mnu uom wmumsmuxm tam umsmnsnam Mo 6669 6:6 06 waHucoammuuou aummHun cm> a HHHLvmm mo mmonu nuH3 m9 H + we u .6 aoHuucDMIcmem 036 we m>Hum>Humn umHHm 6:6 mo somHummaou .m mustm // I . .mcoHumHHUHmo / ucommua msu mo mmosu mm mem wfiu / 666 A666H6 666666 66> 666 HH66666 “M6 mmnucwun ummn: van um3OH suon How .6 "muoz 9 ‘ 63 / - H 6; 666666 66> 666 HH66666 66 . H 66.6 I 6 .66666.6 I 6 \H - nocmum “mam: liIl. /A/ \\ 66 I\7//\\\ //6\ 66666 6 I 66 6666 I 6 /oAWI..6 \ 6H666.6 I a .6N.666~ I 6 6 666666 66> 666 HH66666 66 66.6 I 6 .66H66.6 I 6 nosmum um3QH .IIII I 1 o.~ C :oHumuom HNszc< mo mmsHm> ucwuwwwHw pow mm>uso muHHHnmum Hmuuzmz 6:6 you m Hmnssz mvHoczwm umchmm b umn552 m>63 no uon .o mustm ooom ooow 0006 0006 m 000m 0006 ooom ooou OOOH v _ _ _ b _ _ _ A H q _ g _ o.o Ho.o H.o Il ~.6I6 65 030m C :oHumuom umstc< mo mmsHm> ucmummch now mm>hso quHHnmum Hmuusmz msu now m umgasz mvHochmm umchwm o huHuon> o>mz «o uon .m mustm 6666 6666 6666 6666 6666 6666 6666 666H . . _ .6 m _ m 6.6 o.o HO.OI H.O o.o u no u no u ~.o u a now H B van H A UV m mmaHm>meHm HmoHuHuo :o C aoHumuom umsta< wo uomwwm .m mustm uHHu . B m.o w.o 6.0 0.0 n.o H 6 _ _ _ q . q . . HMO .H 6. 66 6.6 6.6 6.6 h /. H m \ . WV, 6| H.0I 6H66 6 66m6l 666.66 6666 82 . _ _ H _ 0.0 I1.H.o uHuo u 6Huo A 66 . uHuom. H... Al N.o "W C AI 6.6 -.. 6.6 4| m.o 67 0000 wm>uso 66HHHnmum kuusmz 656 so 0000 0005 _ 7 0000 c :oHumuom umHsmc< wchnm? mo uumwmm 000m 0006 H.6I6 000m pom .m mustm 000w 000H «.0 m.0 6.0 m.0 0.0 6.0 0.0 m.0 0.H H.H OH.O H m HOW uhuwHHDmum Human—DMZ 05H HON M— .HOQEDZ mUHOCFAwM umcfimwfi .HU huHuon> 6:53 mo uon .23 no G HHOHuwuom ustwc< muffin? mo 3th .0H mustm 6666 6666 6666 .6 83 6 _I _. 6 _ 66.6 6 6 I c ...666 8 6 Hoe! 166.6 H.$ . 6 .I666 H6 I 6 I666 69 66.6 I 6 666 mm>650 muHHHnmum Hmuuswz 656 00 C :oHumuom umHsmc< wcqum> mo uomwwm .HH mustm m 000m 0000 0006 0000 000m 0006 000m 000m 000H 6 6 _ _ _ 6 _ _I H.0I 0.0 u c H.0 N.0 0.0 0.0 0.0 0.H I. 6.16:. .i .. i 0m.0 u m 600 muHHHnmum Hmuuswz mnu How M 660552 mvHoczmm umchwm LH0 >uHuon> mpmz no uon 656 no G :oHumuom answc< wchumP mo uomwwm .NH muzmHm 0000 000m 0000 m 000m 000m 000H 0 _ _ p _ _ _ Nm.0 H H H H H H 6.6 . H6- -I 66 6 Il.0m.0 70 II 0m.0 .I 66.6 6.6 I 6 g! N¢.0 71 0N.0 u 0 0:6 0H.0I u c How A0.0 u H00 m>6=0 zuHHHnmum Hmuusmz 656 «0 560660 H p 66304 may do ucHom HMUHuHHU m Ham: .6 H + 6 u 6 6:0Huoasmuamem mo uon .MH mustm 6.H 0.H 6.0 T _ #— HS 0 N.0I 0.0 N.0 6.0 0.0 0.0 a 0.H N.H ¢.H 0.H m.H 72 H 66.6 I 6 666 36. I c 666 8.6 I .66 m>u=0 muHHHnmum Hmuusmz mnu mo nucmum um3og wnu co ucHom HmoHuHuo m H 6 666: “6 H + .6 u .6 :oHuuc=MIcmem mnu mo m>Hum>Humn umuHm can we uon .6H mustm 11 OTHI 6 73 6 66.6 I 6 666 OH.OI u G How A0.0 fl .Uv m>HDQ huHHHDmum Hmhuamz 2.3 NO SonHm H033 636 no usHom HmoHuHuo a mum: H? H + M? u 9 6:0Huua:MIcmem no uon .mH mustm \\ Ill ooH! m.H 0.H mwo 0 -I-I Q5 0H.0I u C .SOHumuom umstc¢ How A0.0 n H00 mw>u:0 muHHHnmum Hmuusmz mnu co 0 660852 m>m3_wcH%um>.mo uomwwm .0H musmHm m 74 00mm 0000 - 4 00mm 0000 p 00mm 000m 000a [02H [H.H 75 0H.0I u c .Hom 6uHHHnmum kuuswz msu How m umnasz mvHoa6om udemww .H 0 6uHuon> 6263 «o 6.on ma... :0 0 “6.0852 65.63 wcH6um> mo uommmm .6H onstm 0006 0000 000m 0006 m 000m 000w 000H 0 .6 _. _ _ _ a _ 6.... L Ildm.0 I10m.0 4.6090 H.oI u c I Il06.0 76 6H.6 I 0 66666666 6666666 666 H6.6 I 660 696.650 zuHHHnmum Hmnusmz on“. so a 606—832 0263 06:66.6? mo 666.000 6666 6666 6666 6666 6 6666 6666 6666 r _ P _ _ _ _ u q u d d a 1 N0 66.6 I 66 6.6 I 6 6.6 H.0uC .6H 666666 003 _ d m.0 H0 L1N.O 66.0 [0.0 10.0 l0.H FHA LI6.H 0.48 0.44 0.42 0.40 77 O = 0.1 I J l I I I I I I 1000 2000 3000 R 4000 5000 Figure 19. Effect of varying Wave Number 5 on the plot of Wave Velocity C against Reynolds Number R for the Neutral Stagility for 0 = 0.10 0'6'—I 0.5..- 0.4"' 0.3..- 0.2-— 78 O = 0.1 l l J l J 0.0 1 l l ’1 "7 1000 2000 R 3000 4000 5000 crit Figure 20. Plot of Wave Number 5 against Critical Reynolds Number R i for Angular Rotation 0 = 0.10 and -o.1o cr t d uL-u-m..u _ . l APPENDI CES .235" ‘ "6’1": I..——___-__ .... APPENDIX A NUMERICAL TECHNIQUES A-l Numerical Technique to Solve for the Primary Velocity Profile This numerical scheme solves the equation (3.1.8). It is a fourth order Runge-Kutta method as presented by Weeg and Reed. The accuracy of this scheme was checked by comparing the values of the primary velocity U with those given in Schlichting (1962). The agreement between the two is very good. The technique used is given in appendix:C-3 and explained below: The third order differential equation (3-1.8) is written as three first order differential equations 6 F' -2 n n I 6 - -1 Zn Wn (A 1 ) W'=-lZW n 2 n n .J with the initial (boundary) conditions at the plate being Fn 8 0.0 -I Z ' 0.0 (A-1-2) n W .003 n ..J The value of W“. was guessed and was later found by iteration to be 0.332. The functions F, 2, and W at one step away from the plate are 79 and 80 + 2 a + a 1 Fn +6(a1 + 2 a2 3 4) + 2 b +-b 1 zn+6(b1+2b2 3 4) + 2 c l Wn + 6(c1 + 2 c2 3 +-c4) c for i = l to 4 i are computed as follows: a a =h(Zn +—) b "hmn+—) a C E. .l. ..l 2(zn+2)(wn+2) 0 II 82 1 h(Zn + '2—-) b2 3 h(Wu +.2_9 0" II P. 32 .°_2 '2(Zn+2)mn+2) a3 h(Zn + T) b h(W +—3) 4 n 2 0‘ II a C h .1 .61 "2(zn+2)mn+2 ). (A-l-3) (A-1-4) (A-l-S) (A—1-6) where h is the step size. Using this technique the solution is generated to free stream. That value of Wh is adOpted which gives U 8 1.0 at the farthest distance from the plate. 81 A-2 Numerical Technique to solve the Governing Stability Equations The numerical integration of the Orr-Sommerfeld equation (3.2.1) for O - O and equations (2.5.13) for Q 7‘ O is successfully completed by using a standard Runge-Kutta integration method for second and fourth order equations as illustrated by Collatz (1951). The .0 = 0 case was solved separately because the order of the equation reduces by two and thus the numerical technique changes for the two cases. Equations (2.5.13) are solved separately but simultaneously so that the order of the equations remains less than fourth. Alternatively we can solve one sixth order ordinary linear differential equation in one variable as explained by Zurmfihl (1948). The accuracy of the results is checked by comparing with the results of Kaplan (1964), Lin (1945), Schlichting (1935), Radbill and Van Driest (1966) and also by comparing with the experimental data of Schubauer and Skramstad (1947). This is further verified by varying the step size. The scheme as used in the final form for the simultaneous solution of equations (2.5.13) is outlined below: At y V00 ‘ ¢(y.I) V10 = -h ¢'(y,I) v20 = 22¢"(y.1> h3 v30 :3 - .6— ¢"'(y,1) U 00 Y(y.1) U “h Y'(Y,I) 10 Using equations (2.5.13), Now at Then as Now at 82 2 h 6 1 2 f(’ ha Vao’ h V10’ V00’ "00’ Y) h2 1 933 5—[1 i- ¢"' - fix + iR(U-c)}¢' - 1 2 (U'-20)¢' oz K + {K + 1R(U-c)}‘!] 4 II 2 1 24 f(h2 v20’ ' h V10’ V00’ Uoo’ V) II4 2 Euiaw-c) + 2x}¢" + 2 1R 0 ¢' - {K + iR(K(U-c) + u")}¢ - 2K R n ‘1’] 'h/2 01 V11 v21 V31 = U -I-\II-‘ 01 before 2 F = ‘— 2 2 b4 Gz'fl §_ 1 f(' 3 V31’ ' EV11’V01’ U01’ y ’ 2) D“ ....I 2 h. f‘ 2 V21’ ' ' 2) h V11’ V01’ U01’ y Y'h: +'G = V +'V 30 2 00 10 +-V +-V V 20 02 +I3V = V + 2V 30 10 20 + “G V12 83 V 3 V +-3V +'6G 22 20 30 2 v32 = v30 +-402 U02 3 Uoo +U10 +F2 Again 2 -L .6— 1 F3 ' 2 f(' ha V32’ h V12’ V02’ ”02’ y ' h) 4 h 2 1 Ga ’ 24 f(h2 V22’ ' h v12, V02’ "02’ y ' h) Then F = l(F +I2F ) 3 1 2 F'='1'(F +41? +1?) 3 1 2 3 G ='l- (86 +Iac - G ) 15 1 2 3 1 'a— - G 5(9c;1+1262 G3) II= G 2(G1 + 262) c"'=£(c +16 +6) 3 1 2 3 Thus, the functions and their derivatives at y - h are ¢(y'h) = voo +V10 + v20 +'Vao +’G 1 I =_._ I ¢ (y-h) h (V10 + 2v20 +-3v30 +IG ) II = .2— II ¢ (y-h) hz (V20 +r3v30 +'G ) III = _ L II ¢ (y-h) h3 (V3o +-G ') and - = + Y(y h) U00 U10 + F 84 Y's-h) = - hi (u10 + F') This is done until the plate is reached and the combined functions are formed. Flow Diagram APPENDIX B Solve Equation (3.1.8) for primary velocity profile and store results Polynomial from.0rr- Sommerfeld Equation for Inviscid Flow 01.: 0 _1_‘_0. Read -(— New Data Read OIBICr,a,R To (1 ) : (1) 0 = 0 YES II _ I NO Polynomial Polynomial gives gives two three bounded bounded roots A roots I I 2 I = 3 Initialize Individual functions at y = 1.5 E=J m II Start Numerical Integration using Runge-Kutta Fourth Order Method. Step size D = 0.02. No of Steps N - 75 6T of o ‘XES ‘41NO ileum f 6 Calculate D4¢ from Calculate D4¢ and DZY from.O-S O-S equation (3.2.1) equations (2.5.13a) and (2.5.13b). and use this equa- Add these two equations for later I 85 I 86 :: 11 = 2 (A ¢II" + B) 8 Calculate rn+1 = n+1 I 11 (A ¢IIII + B) = n+1 I 1 Compute for 11 P = (¢)I=1 ¢n+1 . ¢n+1 ' rn+1Pn+1 I g I _ I ¢n+l ¢n+1 rn+an+l I"; . . ;. ° IIII _ "II ¢n+1 ¢n+1 rn'I'an'I']. I 0 3 0 ' IE I For I1 Compute YES I W .. Yn+l “n+1 l::I+1Vn-I-1 I a I _ I Yn+l Yn+l rn+lvn-I-l II = II _ II YD‘I’I n+1 rII+1VII+1 . _ l y = + II I1 1F 11 = J§_____I ILNO f Ln = n+l']** n < N' n = N I = 2 _ n = N, i = I1 ' S = (r ) | I #6 11 ! For i = 11 %-1 = 056-1 ' SiPn-l I . I _ I ¢n-l ¢n-l 8iPn-l 61:. = «>32. "3'1""? n-l I 87 Cont'd. II I O 0 A NO For i = 11 ' Compute YES Yn-l = Yn-l " SiVn-l I = I _ I Yn-l Yn-l siVn-l =‘s +1: II = II _ II .. Yn-l VII-'1 siVn-l i n 1 81 = Sifin-l I 11 = I +1 V I =53 II 1 150 ; _ YES YES r - 1 n_= N Set 11J - 2 0 II For 1 g 11 I In n-l, In > 1 = - s P t a 1 ¢n-1 ¢n-l i n-l LNOL n > 1 91’6-1 = ¢II-1 ' SIPI'I-l I I s a +1 n = l . . . . i i -l ‘ IIII = IIII _ 8 PM” u I ¢n-l ¢n-l i n-l V Yn-l g Yn-l - siVn-l YI'l-l . VII-1 ' 336-1 0 = 0 II = II ._ II Yn-l Yn-l 8ivn-l .6 __ < N0 Compute A and B from YES ¢11+A¢12+B¢13 = 0 I Y11+AYth --l-B‘i’13 = 0 Compute Agfgom ¢11+“? Compute Compute Q1 . (¢I1+A¢I2+B¢1'3)n=n Q2 = (¢11+A¢12+B¢13)n=N-1 I Q 1 = (¢I1+A¢I2)n=u Q2 = (¢11+A¢12)n=N-1 Cont'd. To (1) Test T (Q1): ((21% Tr (Q2): ’ T. (Q2)1 T 3"T: +'T1 I T < 10'3 I I 4N0 ' YES To Read new data Iterate on a and R for fixed C until conver- gence 15 satisfied I Convergence satisfied I Print O’Bacraa3R APPENDIX C COMPUTER PROGRAMS Cel Use of Computer Programs Eigenvalues for the Orr-Sommerfeld equations (3.2.1) for -0 = 0 or (2.5.13 a and b) fbr O # 0 were obtained by computer programs written in fortran IV with the use of the CDC 3600 computer. A description of these programs is given in section C-2 and listings of the main program and subroutines in section C-3. Input to the main program is by two common cards and the rest by data cards. These are detailed later in this section. In these cards we specify the step sizes outside and inside the boundary layer, total number of steps in the boundary layer, number of data cards, the value of the monitor which controls output, tolerance for the roots of the polynomial, percent changes to be.made in a and R, the maximum changes in these variables to be made during iteration, the desired tolerance in the convergence of the eigenvalues, specifications for the angular rotation. 0 and spanwise wave number 3, initial estimates for the eigenvalues Cr’ a and R and the maximum no. of iterations to be performed. Output from the program is monitored. Monitor = 1 gives express runs and only prints out the eigenvalues along with the convergence obtained after every change of eigenvalue during iteration.l Monitor = 2 prints the above plus the eigen-functions. Monitor = 3 prints the same as monitor = 2 plus additional intermediate results, and monitor = 4 89 —~.- - 90 gives the maximum print out. Monitor 3 and 4 find more use during the debugging process. Once the program is debugged, monitor = 2 shows the smoothness of the eigen-functions. Then monitor = l is used for finding the true eigenvalues. Explanation of the Input Cards for the Main Program Card Number Column Item Format 1 1-10 Step size D1 outside the D 10.2 boundary layer 11-20 7 Step size D2 inside the D 10.2 boundary layer 21-30 No. of steps in the I 10 boundary layer 31-40 No. of data cards I 10 41-50 Monitor I 10 51-60 Tolerance for the roots D 10.2 of the polynomial 2 1-10 Percent change in a D 10.3 11-20 Percent change in R D 10.3 21-30 Maximum change in a or R D 10,3 31-40 Tolerance in the conver- D 10.3 gence of eigenvalues 3 l-10 Angular rotation O D 10.3 ll-20 Wave number 3 D 10~3 21-30 Real part of wave D 10.3 velocity Cr 31-45 Wave number a D 15.4 46-60 Reynolds number R D 15.4 61-65 Number of iterations I S Any number of data cards as shown in input card number one may be added. .1: t' 1“- "H A .1133 91 C-2 Description of the Computer Programs Main Program: This program is labelled CHANLA. It reads the input data, solves and stores the values of the primary velocity and its derivatives, then starting with the first estimate of the eigen- values computes the polynomial, finds out its roots by calling the subroutine polyrt, initializes the individual eigen-functions in the free stream, starts the numerical integration by the Runge- Kutta fourth order method using the filtering technique for extracting the contributions of the growing solutions from the well behaved solutions and then repairs the extractions. After this the combined functions ‘¢, Y and ¢' are formed. At this stage the subroutine Dsole is called to solve the linear algebraic equations formed with the boundary conditions. Then the test function T is developed to check the convergence of the eigen- values. If the test fails to converge within the specified tolerance, iteration of the eigenvalues is started, first increasing a then R by a very small amount and then the values of a and R are predicted for the next step. The whole process is repeated until the values converge to the true values. If after about four iterations the convergence is found to be slow or absent, the initial guess is revised. Subroutines: POLYRT. I This program was taken from the Michigan State University, Computer Library. It was modified for relative accuracy and for use in double precision arithmetic. This solves a complex poly- nomial for its roots. It makes use of other subroutines namely 92 Improve, Reduce, Root and Basic which are also stored in the MBU library. V-CAL Using Runge-Kutta fourth order method it solves for the primary velocity and its derivatives and stores these results for use in the main program at the different steps of numerical integration. D-Sole This subroutine solves a set of n linear algebraic equations with complex coefficients, all in double precision. CMULT and CEXP Subroutine CMULT multiplies two complex numbers whereas CEXP gives the exponential of a complex number, all in double precision. C-3 Listing of the Computer Programs The programs developed to compute the eigenvalues for any given case are listed in this section. One run on the CDC 3600 with one set of estimates of the eigenvalues took approximately 10 seconds. In the programs that follow dollar sign ($) is often used. This is a legal statement separator for the CDC 3600 computer. -.-_“I J u C C C C C 425 C C 22 C 460 461 ll 93 PROGRAM CHAULA THIS PROGRAM SOLVES THE FOLLOWING DIFFERENTIAL EQUATIONS (1) ONE EQUATION OF FOURTH ORDER. (2) TWO EQUATIONS IN 2 VAR. OR ONE IN ONE VAR. OF 6TH ORDER THREE DIMENSIONAL LINEAR STABILITY PROBLEM. DIMENSION PR(15393)0DPR(ISSOBIQDZPRIISBOBIO IDBPR(15303)ID4PR(15393)005PR(193)9VR(15393)IDVR(15393)0 ZDZVR(15303)IFR(7)ORR(4)IORI6)9TR(6)OUOR(3)IRRII76)ORR2(76)O 3TBR(6)0TCR(6)0AR(6I6)OBRI606)0XR(6)'WR(6)0UB(153)IDUB‘IBS)‘ 4DZUB(ISSIODSUBISIOD4UBISIIZI153)ITI6)IY(76I DIMENSION PI(15303)9DPI(15393)CDZPIII5303)Q IDBPII15303)0D4PI(15303)OD5PIIIQBIOVI‘153C3IODVIIIS393I' 202VI‘15303)0FI(7)0RI‘4IOOII6)OTII6)9UOI(3’CRII‘76)ORI2(76I9 3TBII6)ITCI(6)IAI(696)IBI(696)9XI(6)9WI(6) DOUBLE PRECISION PRIDPRIDZPROD3PR904PRIDSPRIVRIDVRIDZVRI IFRQRROQROTROUORQRRIORRZ9TBRQTCR'AROBROXRQWRORSRIORSRZOBIRO ZBZROBSRICIRICZRICSRIPURIFIRIERODERIDZERIDSERID4ERIURIDURO BDTRIIDTRZOVOORIVIORIVZORIVSORIVOIRIVIIRIVZIRIVBIRIVOZRO 4V1ZRIVZZRIVBZRIGKIRIGKZRIGKSRIGKRIGKPRIGKZPRIGKSPRIUOOR DOUBLE PRECISION PIIDPIIDZPIIDSPIOD4PIID5PIIVIODVIODZVII IFIIRIIQIITIIUOIIRIIIRIZITBIITCIIAIIBIIXIIUIIRSIIORSIZIBIII 282198310CIIOC2I9C310p01OPIIOEIODEIQDZEIODBEIID4EIOUIODUIO 3DTI190TI2IVOOIIVIOIIVZOIIVBOICVOIIIVIIIIVZIIOVBIIIVOZIO 4V12I0V22IIV321IGKIIOGKZIIGKSIIGKIIGKPIIGKZPIOGKSPIOUOOI DOUBLE PRECISION UIOR’UOIROUI1ROU02ROU12R9U03ROUI3ROFKIRO IFKZRQFKanFKRQFKpROAIROFKZIOFKBIOFKIOFKPIOAII‘UIOIOUOIIO ZUlIIIUOZIIUIZIOUOBIIUISIIFKIIIUBIDUBIDZUBIDSUBOD4UBIUPODUPI 302UPID3UPQZIT0YIDIoDZoDoCIPPIIRALIRBEIRCEIKIRIALIOMOCROCII 43EIPALIPCRIDALIRDIPCHITIPIDELTAZIWRTIDENIALN COMMON /2/ UBODUBODZUBOZQDSUBOD4UB STEP SIZES AND NO. OF STEPS. READ 425001ODZIMBIKPIMONIDELTAZ NBOoS/DI+IoO/DZ+I S IN‘N READ 4200PALOPCROPCHITIP FORMAT (2010.2o3110.2010.2) CALCULATING UBIDUBIDZUBIDSUBOD4UB DISDl/ZoODO 5 DZIDZ/ZoODO 5 N=2*N-I VCAL SOLVES FOR THE PRIMARY VELOCITY DISTRIBUTION. CALL VCALIDIIDZINI GO TO (46594650465022) MON PRINT 460 PRIMARY VEL. U 0F TEXT 8 U8 OF COMPUTER PROGRAM AND 50 ON. FORMAT I10X0*UB(J)*I15XI*DUB(J)*I15XI*D2UBIJ)*99XI*Z(J)*I IBXQ*J*I PRINT 46IOIUBIJ)ODUBIJ’ODZUBIJ)OZIJ)OJOJ‘ION’ FORMAT (3(5X1F1509IOSX0F60395X9I5) PRINT ll ' FORMAT (IOX9*03UB(J)*032XO*D4UB(J)*927XI*J*) PRINT 120((D3UBIJ)ID4UBIJIIJ)0J8193) 12 465 500 21 490 400 405 454 420 see 661 662 421 887 422 65 66 67 471 510 94 FORMAT (2(5XID32025)05X015) PRINT 4250010020M30KP0MON00ELTA2 PRINT 4200PAL0PCR0PCH0TIP CONTINUE Dl=20ODO*DI $ 02‘20ODO*02 S N8005/Dl+100/02+1 FORMAT (60100302I5I €180.00 6 PPI=301415926535897932384626430 KKIO PRINT 210IN FORMAT (*I*0*EIGEN VALUE ITERATION NO OF STEPS N=* I3) KK IS A COUNTER FOR DATA CARDSo IF (MON-2) 49004000400 PRINT 135 KKtKK+l - IP80 S M030 3 ITE=1 IF (KK‘KR) 40504050410 CONTINUE INITIAL GUESS FOR EIGEN-VALU550 READ 4640OM0BE0CR0AL0RIIT FORMAT (30I00302015040I5) FORMAT (40100303I5) PRINT 886 FORMAT (30X0*R*035X0*AL*035X0*CR*) IF (OM 0E0. 00000) 6610662 IR=2 S ILtIR*2 $ GO TO 421 IR=3 $ ILBIR*2 IT IS AN ITERATION COUNTER. MDzMD+I PRINT BBT0R0AL0CR FORMAT (24X03035025) IF (MD-IT) 42204220400 CONTINUE KiK OF TEXT IS EQUAL TO K OF COMPUTER PROGRAM. KIAL*AL+BE*BE IF (MON-3) 471065065 PRINT 66 FORMAT (3X0*0MEGA*03X0*ALPHA*06X0*BE*011X0*CR*010X0*REY NO*) PRINT 670OM0AL0BE0CR0R FORMAT (1X0090201X0Dl20502(1X0D1205)01X001205) M=1 5 UPaUB(M) 3 DUP=DUB¢MI $ DZUP=D2UB(M) $ DBUP=DBUB(M) D4UPBD4UBIM) IF (OM 0E0. 00000) 5050510 CONTINUE FUNCTIONS -COEFFICIENTS OF THE 6TH DEGREE POLYNOMIAL- FRII)820OD*BE*BE*R*R*OM*(DUP‘2000*OM)+AL*AL*R*R*(UP-CR)*02UP 1+K*AL*AL*R*R*((UP-CR)**2-CI*CI)“20OD*AL*R*K*K*CI-K**3 FIIl)=“AL*AL*R*R*CI*02UP-20OD*(K*AL*AL*R*R)*(UP-CR)*CI-20OD 1*AL*R*K*K*(UP-CR)-AL*R*D4UP FR(2)¢OoOD $ FI(2)820OD*AL*R*(K*DUP+DBUP) OOOOO O 95 FR(3)=3.0D*K*K+4.0D*K*AL*R*CI-(AL*R*(UP-CR))**2+(AL*R*CI)**2 F1(3)=4.0D*K*AL*R*(UP-CR)+2.0D*AL*AL*R*R*C1*(UP-CR) FR(4)=0.0D S FI(4)=‘2.0D*AL*R*DUP FR(5)=-3.0D*K-2.0D*AL*R*CI S FI(5)=‘2.0D*AL*R*(UP-CP) FR(6)=0.0D $ FII6):0.0D FR(7)=1.0D $ FIC7)=0.0D GO TO (100100100201) MON 201 PRINT 202 202 FORMAT (9X0*FR(I)*025X0*FI(I)*) PRINT 2030(FR(I)0FI(I)018107) 203 FORMAT (2(4XOD25.15)) POLYRT GIVES THE ROOTS OF THE POLYNOMIAL. 10 CALL POLYRT (FROFI060QROQIODELTA2) J=1 5 DO 451 I=I06 3 IF (QR(I) .LT. 0000) 4500451 450 RRIJISORCI) S RICJ1-OIII) $ J=J+1 451 CONTINUE S EQALBO 5 IF (RR(1) .EQ. RR(2)) 4540452 452 IF (RR(1) .EQ. RRC3)’ 4550453 453 IF (RR(2) .EQ. RRIBI) 4559456 454 RR(2)=RR(3) $ RI(2)=RI(3) $ RR(3)=RR(1) $ RI(3)=RI(11 455 EQAL=2 3 GO TO 456 505 ALNaDSQRT(K) S RIAL*R/ALN S ALzALN FUNCTIONS -COEFFICIENTS OF THE 4TH DEGREE POLYNOMIAL. FR(1)=K*K+AL*R*CI*K $ FIIIIIAL*R*(K*(UP-CR)+D2UP) FR(2)=0.0D $ FII2)=0.0D 3 FR(3)=-((2.0D)*K+AL*R*CI) FI(3)=-AL*R*(UP-CR) FR(4)=0.0D 5 FI(4)80.OD FR(5)=1.0D S FI(5)=0.0D CALL POLYRTCFR0FI040QR0QI0DELTA21 J31 5 DO 515 18104 S IF(QR(I) .LT. 0.00) 5200515 520 RR(J)=QR(I) 5 RIIJIBQIII) S J=J+1 515 CONTINUE S EQAL=0 537 IF (MON-3) 20702070205 456 GO TO (207020702070205) MON 205 PRINT 470 470 FORMAT (9X0*QR(I)*022X0*QI(I)*022X0*ROOTS OF THE POLYNOMIAL*) PRINT 2030((QR(I)0QI(I))01810IL) 545 PRINT 206 206 FORMAT (9X0*RR(I)iIZZX0*RI(l)*.22Xo*NEG. ROOTS OF THE 1POLYNOM1AL*) PRINT 2030((RRII10RIIII)0I8101R’ P IS EQUAL TO THE VELOCITY V‘HAT0 AND DP=D(V-HAT). V I5 EQUAL TO THE VELOCITY U-HAT. BOUNDARY CONDITIONS - PaDP8V=0 AT Y=00 AND THESE APPROACH ZERO AS Y APPROACHES INFINITY. CONDITION ON DP RESULTS FROM CONTINUITY. 207 CONTINUE FINDING OUT THE GROWING SOLUTION I.E. WITH LARGEST DP/P. RRI4I=RRIII S RI(4)=RI(1) RALSDSQRT(RR(1)*RR(1)+RI(1)*RI(1)I 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1121 1120 1122 1125 9 96 RBEsDSGRT(RR(2)*RR(2)+RI(21*R11211 IF ‘0" 0E0. 00000) 111901109 RCE=DSQRTCRR(3)*RR(3)+RI(31*R1(3)1 1F (RAL-RBE) 11100111001111 1F (REE-RC5) 11150111501116 1F IRAL-RCE) 11120111201113 RR(1)=RRI3) s RI(1)-R1(3) s RR<3)=RR(2) s R!(3)=RI(2) RRI2):RRI4) s R1(2)=R1(4) s 60 To 1120 IFIRBE-RCE) 11140111401120 RR(4)IRR(2) s R1(4)aRII2) s RR<2)=RR(3) s RI(2)=R1(3) RR(3)=RR(4) s R1(3)8R1(4) 3 GO To 1120 99(1)=RRI3)5RI(1)=R1(3> s RRI3)=RR(4) s RI(3)=RI(4) Go TO 1120 RRI1)=RRI2) S R1II)=RI(2) s [F (RAL-RCE)1118.111801117 RR(2)=RR(4) s R1(21=R1(4) s GO To 1120 RRI2)=RRI3) s R1(2)IR1(3) s RRI3)=RR(4) s RII3)=RII4) GO TO 1120 IF (RAL-RBE) 11210112001120 RR(1)=RR(2) S RIIII-RIIZ) s RR(2)=RR(4) s RII2)=RI(4) CONTINUE IF (MON-2) 11250112201122 PRINT 2030((RR1110R1(1))018101R) CONTINUE ASYMPTOTIC SOLUTION AND INITIALIZING. DO 6 1=IIIR $ J=1 5 M21 A1R=0.0DO s A11=0.000 CALL CXPEIAIRIAIIIPR(J01)0P1(J01)) CALL CMULT(RR(I)0R1(I)0PR(J01)0P1(J01)0DPR(J0I)0DPIIJ0I)) CALL CMULT(RR(1)0RI(I)0DPR(J01)0DP1(J0I)0D2PR(J01)IDZPI(J01)) CALL CMULTIRP(I)0RI(I)0DZPR(J0I)0D2PI(J0I)IDBPR(J01)0 103P1(JT1)) CALL CMULTIRRIII0RII11003PR(J011003PI(J0I)004PR(J0I)0 104P1¢J,111 IF (OM 0E0. 00000, GO TO 6 CALL CMULTIRRIIIORII110D4PRIJOIIOD4PIIJOI)ODSPRIJ0110 IDSPIIJ0111 UpsUBIM) $ DUPsDUBIM) 5 DZUP8D2U8(M1 $ DSUP=D3UB(M) VR‘J0I)=(D4PR(J0I)-(AL*R*CI+2.OD*K)*02PR(J0I1+AL*R*(UP-CR)* 102P11J011+2.OD*AL*R*OM*DP1(J0I)+(K*K+AL*R*K*CI1*PR1J011'AL** 2R*(K*(UP-CR1+D2UP1*PI(J01))/(-2.OD*R*K*OM1 VI(J0I)=(D4PI(J01)-(AL*R*C1+2.0D*K)*D2P1(J0I)-AL*R*(UP-CR)* 102PRIJ0I1-2.OD*AL*R*OM*DPR(J0I)+(K*K+AL*R*K*CI)*PI(J01)+AL* 2R*(K*(UP-CRI+D2UP)*PR(J0I)1/(-2.OD*R*K*OM) DVR(J01)=(D5PR(J01)-(AL*R*C1+2.OD*K)*DBPR(J01)+AL*R*(UP-CR) 1*03PI(J01)+2.OD*AL*R*OM*DZP1(J01)+(K*K+AL*R*K*CI)*DPR(J0I)- 2AL*R*(K*(UP-CR)+DZUP)*DPI(J01)+AL*R*DUP*DZPI(J0II‘AL*R*(K* 30UP+DBUP)*P1(J0I))/(-Z.OD*R*K*OMI DVI(J0I)=(05PI(J0I)“(AL*R*CI+2.OD*K)*03PI(J0I)‘AL*R*(UP-CR) 1*03PRCJ011-2.OD*AL*R*OM*DZPR(J0I)+(K*K+AL*R*K*CI)*DPI(J0I)+ 1. 0.25.1019“. 'a‘_.' I C 97 2AL*R*(K*(UP-CR)+DZUP)*DPR(J01)‘AL*R*DUP*DZPR(J01)+AL*R*(K* 30UP+03UP)*PR(J.I)1/(-2.0D*R*K*OM) D2VR(JII)=-(03PI(J.I)-(K+AL*R*CI)*DPI(JOI1‘AL*R*(UP-CR1*DPR 1(J.I)-BE*BE*R*(DUP-2.0D*0M)*PR(J0I)/AL)*AL/K+(K+AL*R*CI)*VR 2(J0I1‘AL*R*(UP-CR)*VI(J0I1 02V1(J01)=( D3PR(J.I)-(K+AL*R*CI)*DPR(J0I)+AL*R*(UP-CR)*DPI IIJOI1+BE*BE*P*(DUP-2.00*0M)*PI(J0I)/AL)*AL/K+(K+AL*R*CI1*VI 2(J0I1+AL*R*(UP-CR)*VR(J0I) 6 CONTINUE 31 DO 31 J=10IN 3 N=2*J-1 Y(J)=Z(N) 5 C301 3 NN=IN-1 RUNGE KUTTA METHOD FOR THE SIMULTANEOUS SOLN. OF TWO DIFF. EQUATIONS. REFERENCE COLLATZ. DO 55 J=10NN 5 DO 54 I=10IR $ M=2*J-1 UPBUB(M) 5 DUPSDUB(M) S DZUPSDZUB(M) VALUES AT Y. VO0R=PR(J.I) $ VOOIBPI(J01) S V10R=-C*DPR(J.I) V101=-C*DPI(J01) V20R=C*C*D2PR(J0I)/2.0D0 $ V201=C*C*D2P1(J01)/2.000 V30R=-C*C*C*D3PR(J0I)/6.ODO $ V301=-C*C*C*D3PI(J0I1/6.ODO 1F (OM .EQ. 0.000) GO TO 1200 U00R=VR(J0I) $ U0018V1(J0I) $ U10R=-C*DVR(J0I) U101=-C*DVI(J0I) FK1R=((6.0DO*V301/(C**3)-(K+AL*R*C1)*V101/C-AL*R*(UP-CR)* 1V10R/C+BE*BE*R*(DUP-2.0D0*OM)*VOOR/AL)*AL/K+(K+AL*R*CI1* 2UO0R-AL*R*(UP-CR)*U001)*C*C/2.0DO FK11=((-6.000*V30R/(C**3)+(K+AL*R*CI1*V10R/C-AL*R*(UP-CR)* 1V101/C+BE*BE*R*(DUP-Z.ODO*0M)*VO0I/AL)*AL/K+(K+AL*R*CI1* 2UOOI+AL*R*(UP-CR)*UO0R)*C*C/2.0DO 1200 GK1R=(C**4)*((AL*R*CI+2.0D0*K)*V20R*2.0DO/(C*C)-AL*R*(UP-CR 1205 1)*V201*2.0D0/(C*C1-(K*K+AL*R*K*CI)*VO0R+AL*R*(K*(UP-CR)+ 2D2UP1*V0011/24.0D GKII=(C**4)*((AL*R*CI+2.ODO*K)*V201*2.ODO/(C*C)+AL*R*(UP-CR 1)*V20R*2.0D0/(C*C)-(K*K+AL*R*K*CI)*V001-AL*R*(K*(UP-CRI+ 202UP)*VO0R)/24.0D IF (OM .EQ. 0.0001 GO TO 1205 GK!R=GK1R+(2.ODO*AL*R*OM*V10I/C‘2.ODO*K*R*0M*UOOR)*(C**41/ 124.000 GKIISGKII'(2.ODO*AL*R*OM*V10R/C+2.ODO*K*R*OM*UOOI1*(C**4)/ 124.0D0 VALUES AT Y-C/2 M=M+1 UPBUB(M) s DUPsDUB(M) $ DZUP!D2UB(M1 V01R=VO0R+0.500*V10R+0.25D0*V20R+0.12500*V30R+0.0625D0*GK1R V01I=VOOI+0.5D0*V10I+0.25D0*V201+0.125D0*V301+0.0625D0*GK1I V11R8V10R+V20R+0.75D0*V30R+0.5D0*GK1R V11I=V101+V201+0.75D0*V30I+0.5D0*GK1I V21R8V20R+1.50D0*V30R+105D0*GK1R V21I=V2OI+1.500*V301+1.5D0*GK1I "-'¢-:j_~ - 1210 1215 98 v31R=V30R+2.0DO*GK1R s V311=V301+2.0DO*GK11 1F (OM .EQ. 0.000) 60 To 1210 UO1R=UOOR+0.SDO*01OR+0.25DO*FK1R UOI1=U001+0.500*U101+0.2500*FK1I U11R0U10R+FK1R s 0111=U101+FK11 FK2R21(6.ODO*V311/(C**3)-(K+AL*R*CI)*V1ll/C-AL*R*(UP-CR)* 1V1JR/C+BE*BE*R*(DUP-2.0DO*OM)*VOlR/AL)*AL/K+(K+AL*R*CI1* 2U01R-AL*R*(UP-CR)*U011)*C*C/2.0DO FK21==XR<1)*C11+X1(1)*CIR+C21 GO TO 635 630 T8R(IP1=XR(11*81R-XI(11*81I+XR(21*82R-XI(21*82I+B3R TBI(IP1=XR(11*81I+XI(11*81R+XR(21*821+XI(21*82R+831 TCR(IP1=XR(1)*C1R-XI(11*BII+XR(2)*C28“XI(21*C2I+C3R TCICIP)=XR(11*CII+XI(11*C1R+XR(21*C21+XI(21*C2R+C3I 635 TR(IP1=T8R(IP1/TCR(IP1 S TIIIP1=TBIIIP)/TCI(IP) TCIP1=DSQRT(TR(IP1*TR(IP1+TI(IP1*TI(IP11 IF (MON-21 14501260126 126 IF (MD-21 13001400140 130 PRINT 135 135 FORMAT‘1X0*1T*01X0*IP*04X0*OM*013X0*R*015X0*AL*014X0*BE*0 114X0*CR*014X0*TR*015X0*TI*014X0*T*1 GO TO 145 140 PRINT 24 24 FORMAT (4X0*IP*04X0*OM*013X0*R*015X0*AL*014X0*BE*014X0*CR*0 114X0*TR*015X0*TI*014X0*T*) 145 IF (TIIP1‘TIP1 10201020101 101 IF (IP~2) 11001150115 110 PRINT 2801TE0IP0OM0R0AL0BE0CR0TR(IP10TI(IP)0T(1P) 28 FORMAT (1X0I201X0120811X00150811 GO TO 120 115 PRINT 2301P00M0R0AL0BE0CR0TR(IP10TI(IP)0T(1P) 23 FORMATC4X0I208(1X00150811 ITERATION 0F EIGEN-VALUES. 120 GO TO (105019701981 IP 105 DAL=PAL*AL 5 AL=AL+DAL S GO TO 421 197 AL=AL-DAL S RDBPCR*R $ R=R+RD 5 GO TO 421 198 DTR1=CTR(21-TR(111/0AL $ DTII=CT1121-TI(111/DAL 5 R=R'RD DTR2=CTRC31~TR(1)1/RD $ 0T12a(TI(31-T1(111/RD 302 DEN=DTR1*DT12‘DT11*DTR2 5 DAL=(TI(11*0TR2-TR(11*DT121/0EN RD=(TR(11*DTIl-DTR1*TI(111/0EN 000 000 708 109 108 102 211 410 108 IF(DABS(DAL1 .GE. (PCH *AL11 DAL=PCH*AL*DAL/DABS(DAL) IF (DABSIRD1'(PCH*R11 10891080109 RD=PCH*R*RD/DABS(RD1 DTR1=100.0D*RD/R 3 R=R+RD $ DTR2=100.0D*DAL/AL ALBAL+DAL S IP=0 S ITE=ITE+1 PRINT 230IP00M0R0AL0BE0CR PRINT 230IP0OM0DTR10DTR2 GO TO 421 PRINT 211 FORMAT (* CONVERGENCE TEST SATISFIED *1 PRINT 24 PRINT 230IP0OM0R0AL0BE0CR0TR(IP10TI(IP)0T(IP) GO TO 400 END SUBROUTINE DSOLECAR0AI0XR0XI0N0ND1 THIS SUBROUTINE WAS WRITTEN BY DR. W.C.REYNOLDS 0F STANFORD UNIVERSITY. IT WAS MODIFIED FOR USE IN DOUBLE PRECISION. DPRE.COMPLEX SOLUTION OF N SIMULTANEOUS LINEAR ALGEBRAIC E0. DIMENSION AR(606)0AI(60610XR(6)0XI(61 DOUBLE PRECISION AR0AI0XR0X10ER0EI0C0CX0D0Y0Z CONSTRUCT THE TRIANGULAR ARRAY NM1=N-1 DO 50 K=10NM1 NORMALIZE EACH ROW 0N ITS MAXIMUM MODULUS ELEMENT M=N-K+1 DO 8 I=10M C30. 000 DO 4 J=10M CX=AR(I0J)*AR(I0J1+AI(I011*AI(I0J1 IF (CX-C) 40403 CBCX FORMAT (10X0D25.161 CONTINUE C81.000/DSQRT(C1 XRII1=XR(I1*C XI(I1=XI(11*C DO 7 J=10M AR(IOJ)=AR(IOJ)*C AI(IOJ18AI(10J)*C FORMAT (1X01401X0I404I5X0D25.1611 CONTINUE CONTINUE FIND THE OPTIMUM PIVOT ROW. C800ODO DO 12 I=10M CX=AR(I0M)*AR(I0M1+AI(I0M1*AI(I0M1 10 23 12 13 14 16 80 81 20 35 24 30 50 54 25 57 109 IF (CX-C) 12012010 C=CX IP81 FORMAT (5X0I405X0D25.161 CONTINUE STORE THE PIVOT ROW IF (C1 70070013 IF (IP-M) 14020020 00 16 J=10M ER=AR(1P0J) $ EI:AI(IP0J) $ ARIIP0J1=ARIM0J1 AIIIP0J1=AIIM0J1 AR(MOJ1=ER AI(MOJ1=EI ER=XRIIP1 S EI=XIIIP1 XR(IP1=XR( M1 5 XI(IP)=XI( M) S XR(M)=ER S XI‘M)=EI FORMAT (* REARRANGED MATRIX A *1 FORMAT (12011.31 ELIMINATE THE MTH ELEMENTS MM1=M‘1 DO 30 1:10MM1 D=ARIM0M1¥ARIM0M1+AI(M0M1*AI(M0M) ER=IAR(I0M1*AR(M0M1+AI(I0M1§AIIM0M11/D EI=(AI(I0M)*AR(M0M1‘AR(I0M1*AI(M0M11/D Y=XRII)-ER*XR(M1+EI*XI(M1 $ Z=XI(I1‘ER*XI(M1-EI*XR(M) XR(I1=Y S XI(I1=Z FORMAT (*D=*05(1X002501611 DO 30 J=10MM1 Y=AR(I0J1“ER*AR(M0J1+EI*AI(M0J1 Z=AIIIOJ1‘ER*AI(MOJ1-EI*AR(MOJ1 AR(IOJ1=Y S AIII0J132 FORMAT (1X0I401X01406(1X0D20.1211 CONTINUE CONTINUE CARRY OUT THE BACK SOLUTION C=ARI1011*AR(1011+AI(1011*AIII011 IF (C1 70070054 Y:(XR(11*AR(1011+XI(11*AI(10111/C Z=(XI(11*AR(1011-XR(11*AI(10111/C XR(11=Y S XII11=Z FORMAT (5X03(5X0025.1611 DO 60 1820N IM1=1~1 ER=XRII1 S EI=XIII1 DO 57 J=10IM1 ER=ER-AR(I0J1*XR(J1+AI(I0J1*XI(J1 EI=EI-AR(I0J1*XI(J1“AI(I0J1*XR(J1 D=ARII0I)*AR(I0I1+AIII0I1*AI(I011 XRII1=(ER*AR(I0I1+EI*AI(10111/0 XIII):(EI*AR(I0I)-ER*AI(10111/0 '.i“'. . hv' 000 00000000000 60 7O 75 11 55 33 44 0‘ NUO‘ \OQOUI 10 110 CONTINUE RETURN pRINT 75 FORMAT IZZHODSOLE SYSTEM SINGULAR1 PRINT 810I((ARIIOJ)0AI(IOJ110J810N10I=10N1 PRINT 810(1XR(I10XI(I110I=10N) RETURN END SUBROUTINE POLYRTIFR0FI0NORDER0RR0RI0DELTA21 THIS PROGRAM ALONGWITH THE NECESSARY SUBROUTINES WAS TAKEN FROM THE M.SoU. COMPUTING LIBRARY. THIS WAS WRITTEN BY MR0 WILLIAM HOLSTEN OF THE UNIV. OF CALIFORNIA0SAN‘DIEGO. IT WAS MODIFIED FOR USE IN DOUBLE PRECISION AND WITH EXTENDED ACCURACY. THIS SUBROUTINE ALONGWITH SUBROUTINES IMPROVE0REDUCE0ROOT0 AND BASIC DETERMINES A ROOT OF THE POLYNOMIAL WITH COMPLEX COEFFICIENTSOWHICH IS THEN IMPROVED UNTIL THE DESIRED ACCURACY DELTA2 IS REACHED. THEN THE ORDER OF THE POLYNOMIAL IS REDUCED. THIS PROCESS IS REPEATED UNTIL ALL THE ROOTS ARE OBTAINED. DIMENSION RRI9 )0RI(9 )0FR(9 )0FII9 10FFR‘9 )OFFI(9 ) DOUBLE PRECISION RR0RI0FR0FI0FFR0FFI0ZD0T10R KaNORDER+1 DO 11 I=10K FFR(I1=FR(I1 FFI(I1=FI(I1 NORDR=NORDER IFIFFI(K11203302 IFIFFRIK11304403 NORDR=NORDR-1 K=K'1 IF(NORDR~1166066055 STOP 66‘ IF(FFR(K1“1.000120402 ZD=FFR(K)*FFR(K1+FFIIK)*FFI(K1 DO 5 I=10K T18(FFR(I)*FFR(K1+FFI(I1*FFI‘K11/ZD FFI(I1:(FFR(K1*FFI(I1‘FFR(I1*FFI(K11/ZD FFRII18T1 IFIFFR(11170807 IFIFFI(11170907 DO 10 I=10K FFRII1=FFRII+11 FFIII1=FFI(I+11 RR‘NORDR180.ODO RI(NORDR1¢0.0DO 000 111 NORDRBNORDR-I KsK-l GO TO 4 NORINORDR-I R=1.0DO DO 1 I=10NOR JJBNORDR+I-I CALL ROOT(FFR0FFI0JJ0RR(I10RI(I)0R) CALL IMPROVE(FFR0FF10JJ09R(I10RIII100ELTA21 R=DSQRT(RR(I1*RR(1)+RI(I)*RI(1)) CALL REDUCEIFFR0FFI0JJ0RR¢I)0RI(11) CONTINUE RR(NORDR)=-FFR(1) RI(NORDR)8-FFI(I) END SUBROUTINE IMPROVE(TR0T1oNORDoCRoCIoDELTAZ) DIMENSION TR(9 )0TI(9 )0DR(9 10DI(9 ) DOUBLE PRECISION TRoTlooRooloTTRoTTloTl0DDR00010CR0C10 1CRNEW0CINEW COMMON 09.0! DO I I=10NORD ZI=I DR(I)=TR(1+I)*ZI OIII)=TIII+1)*ZI ICNT=~25 TTR=TR