ABSTRACT THE STABILITY OF PLANE POISEUILLE FLOW SUBJECT TO A TRANSVERSE MAGNETIC FIELD By James Anthony Kutchey The stability of an electrically conducting fluid flow- ing between parallel planes subject to a transverse magnetic field is investigated for infinitesimal three-dimensional dis- turbances. Primary interest is in the effect of the magnetic field and magnetic fluid parameters on the critical point of the neutral stability curves. The governing stability equa- tions include perturbations for both the velocity and magnetic fields and result in a sixth order coupled set of linear ordinary differential equations. This set represents an eigen- value problem that is transformed to an initial value problem and solved numerically using a fourth order Runge-Kutta integra- tion scheme with a Special filtering technique. The results indicate a strong dependence of the critical eigenvalues on both the magnetic field strength and a fluid property, the magnetic Prandtl number. The effect of both is to greatly stabilize the flow. Inclusion of the Span-wise wave number does not affect the eigenvalues other than in a manner predicted by Squire's Theorem. The numerical results are also compared to previous data obtained by asymptotic expansion techniques. THE STABILITY OF PLANE POISEUILLE FLOW SUBJECT TO A TRANSVERSE MAGNETIC FIELD BY James Anthony Kutchey A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR 0F PHIIOSOHIY Department of Mechanical Engineering 1970 ACKNOWLEDGMENTS The author wishes to express his sincere appreciation to his Major Professor, Dr. Merle C. Potter, Associate Professor of Mechanical Engineering for his valuable guidance and encouragement throughout this study. Without his personal interest and concern this thesis representing the consumation of the author's graduate work would not have been possible. Special thanks are conveyed to Dr. Charles R. St. Clair, Jr., Chairman, Mechanical Engineering Department for his advice and support in this graduate program. Thankful acknowledgment is extended to other members of the Guidance Committee: Dr. James V. Beck, Dr. J.S. Frame, and Dr. K.M. Chen. The author is grateful for financial support received from 8 NSF Traineeship and for supporting funds from the Mechanical Engineering Department and Division of Engineering Research. To his wife, Arlene, the author dedicates this dissertation for her understanding, cooperation, and patience during the entire graduate program. ii CHAPTER II. III. IV. TABLE OF CONTENTS List Of Tables ......O........... ........... ... ..... List of Figures ............ ......... ......... ...... NomnCIature ..... .0 ...... ...........O 0.. ........... INTRODUCTION ........... ............... ...... ....... 1.1 Review of Literature .......................... 1.2 Description of the Problem ....... ...... ....... STABILITY EQUATIONS ................................ 2.1 Fundamental Equations ............. ..... . ...... 2.2 Non-Dimensional Variables ......... ...... ...... 2.3 The Primary Flow and Magnetic Field Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . O . . . . 2.4 The Linearized Equations for Small Disturbances . . . . . . . . . . . . . . . . . . . . . . . O . . . O O O . O . O 2.5 Eigenvalue Pr0blem .....................O...... NUMERICAL SOLUTION OF THE STABILITY EQUATIONS ...... Introduction .................................. Starting Conditions ............ ............... Growth of the Eigenfunctions .................. Purification Scheme ........................... Iteration Scheme for the Eigenvalues ... ..... .. Criterion for Guessing the Eigenvalues ........ Range of Parameters Considered ................ Special Case of Zero Hartmann Number .. ..... ... wwwuwwww O G>\JO\UIF‘U)h3P‘ RESULTS, CONCLUSIONS, AND RECOMMENDATIONS FOR FURTHER STUDY ..0.l......................... ..... .0. 4.1 Numerical Results . . . . . . . . . . . . . . . . . . . ...... . . . I 4.2 comelus ions . . . . . . . . . . . . . . . . . . . . O ............ . O 4.3 Recommendations for Further Study ...... ....... iii 13 15 18 19 19 19 20 22 25 27 27 28 30 3O 33 34 BIBLIOGRAPHY ILJUSTRATIONS (Tables and Figures) ................. APPENDICES A NUMERICAL TECHNIQUE TO SOLVE THE GOVERNING STABILITY EQUATIONS ............. B mMmTER mOGRAM ......O.... ............. .0 3-1 Description of the Computer Program ......................... ..... B-2 Listing of the Computer Program ..................... ....... .. iv 35 38 58 62 62 LIST OF TABLES Table Page 1. Comparison of Critical Eigenvalues ................. 38 2. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M = 0. 001, Pm =10- -6 cocoa-0000000000 ooooooooooooo oo 39 3. Eigenvalues for the Stability Equations for the Case of Neutral Stability (ci = 0.0) M=1.,O Pm =10-6 ....... ............... 39 4. Eigenvalues for the Stability Equations for the Case of Neutral Stability (ci = 0.0) M = 10 0, Pm =10-1 0.00.00.00.00... ooooooooooooooooo 40 S. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1" = 0.0) M32.0’ Pm -10-6 ................... .............. [+0 6. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M=2.0, Pm B10-1 ..................... ............ 41 7. Eigenvalues for the Stability Equations for the Case of Neutral 10Sgability (ci =0.0) M=3.’O Pm = ................. ................ 41 8. Eigenvalues for the Stability Equations for the Case of Neutral 10Sliability (ci =0.0) =3. 0, P6 = oooooo oooooooooo ooooooooooooooooo 42 9. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M=4.0, Pm =10.6 ................... ........ ...... 42 10. Eigenvalues for the Stability Equations for the Case of Neutral Siability (ci = 0. 0) “84.0, Pm=10- ................ ................. 43 11. Eigenvalues for the Stability Equations for the Case of Neutral Stability (ci =0.0) M = 5.0, Pm = 10'6 ...... .......... ................. 43 12. 13. 14. 15. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M=5.0’Pm=10-1........................O........ Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M=6.0, Pm =10-6................................. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1.. = 0.0) M=6.0’ Pm =10-2 .................. ............ .. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M=6.0, Pm ~-=10"1 44 44 45 45 Figure 10. ll. 12. LIST OF FIGURES Primary Flow Configuration..... ................... .. Dimensionless Primary Velocity Profiles for various Hartmann Numers ................ .......... 0 Wave Number 0 Versus Neutral Stability for Wave Number 0 Versus Neutral Stability for Wave Number a Versus Neutral Stability for Wave Number 0 Versus Neutral Stability for Wave Number 0 Versus Neutral Stability for Wave Number 0 Versus Neutral Stability for Wave Number a 'Versus Neutral Stability for Plot of Eigenfunctions Critical Point for R for Reynolds Number M = 0 Reynolds Number R for M=1.0 ..................... Reynolds Number R for M I 2.0 Reynolds Number R for M'3.0 ..................... Reynolds Number R for M3400 ......... ........ .... Reynolds Number R for M=5.0 ........ ...... ....... Reynolds Number R for M=6.0 ...............o V‘V +i¢i at the M = 0.001”...................... Plot of Eigenfunctions V - V: + iwi at the Critical Point for P5 = 10'6 and 10"1 for M=S.O ............................................ Plot of Eigenfunctions ¢ = ¢ + 1¢i atlthe Critical Point for Pm = 10"6 and 10 for “85.0 ...... ..... ................... ....... . ..... . vii Page 46 47 48 49 50 51 52 53 54 55 56 57 c = c + ic, r 1 0 2:1 ’111 Ell UL 3‘1 hi NOMENCLATURE Channel half width Coefficients used to obtain combined functions Magnetic flux density Complex propagation Speed of wave Wave speed Amplification rate Electric flux density Electric field Body force Applied magnetic field with com- ponents (Hx’ H , 0) Y Magnetic perturbations with com- h h ponents (hx, y’ z) Conduction current Square of combined wave number _ o 5 Hartmann number - uHa (E; Magnetization vector Polarization vector Magnetic Prandtl number = Vac Pressure U a Reynolds number = —%— Change in R viii <1 X:Y:z :n ‘O m W1 Magnetic Reynolds number = UmaHU Time Test function Average velocity Primary velocity with components (U, 0, 0) Velocity perturbations with com- ponents (vx, vy, Vz) Coordinate axes Wave number in flow direction Change in a Spanwise wave number Electrical conductivity Magnetic permeability Permittivity (dielectric constant) Space charge density Vorticity of primary flow = V X V Vorticity of perturbed flow = v X 3 Mass density of fluid Kinematic viscosity of fluid Complex eigenfunction for magnetic perturbations Complex eigenfunction for velocity perturbations Del operator = i §;+j B._+ Ra'- BY 52 2 2 2 Laplacian = l—z’ + 3—2 + 5—5- ax BY 32 ix A Step size = Ay Superscripts ( )' Differentiation with respect to y (') Differentiation with reSpect to time ( )* Dimensional quantities (A) Magnitudes of perturbation quantities Subscripts ( )w Quantities at wall ( )p Primary flow and magnetic quantities CHAPTER I INTRODUCTION 1.1 Review of Literature The term "transition" as generally used in the field of fluid mechanics applies to the observable change in flow pattern when well-ordered laminar motion becomes turbulent. The theory of stability using perturbation techniques with an assumed form of infinitesimal disturbances attempts to predict the value of the critical Reynolds number or point of insta- bility for a prescribed main flow. The earliest work in the area of transition was per- formed by Reynolds (1883) when he used dye injection for flow in pipes. Theoretical studies to predict the transition from laminar to turbulent flow were begun by Lord Rayleigh (1880, 1887) and Lord Kelvin (1887). Based upon this work, independent studies by Orr (1907) and Sommerfeld (1908) led to the now well known Orr-Sommerfeld stability equation for two-dimensional flow, which is given by iv 1 (U-c)(¢" - 02¢) - U"¢ = -aR (m - 2 02¢" + 04¢) where U Velocity c B Complex wave prOpagation Speed Wave number Q ll ¢ = Eigenfunction R = Reynolds number This equation is based upon the assumption of infini- tesimal periodic distrubances. Neglecting the effects of viscosity, Lord Rayleigh (1914) was able to show that any velocity profile that possesses an inflection point is unstable. Much later, Tollmein (1935) proved that this was not only a necessary but sufficient condi- tion for the amplification of small disturbances. Prandtl (1914) postulated the existence of a viscous boundary layer and was able to define transition, separation and drag coefficients on bodies. Incorporating the viscous boundary layer into stability theory, Prandtl (1921) considered flow over a flat plate and included the effects of the largest viscous terms near the wall. This work along with calculations performed by Tietjens (1925) gave the startling result that the introduction of viscosity into the equations did not produce damping as was presumed but amplification for sufficiently large Reynolds numbers for particular wavelengths of the dis- turbances. Tollmein (1929) demonstrated that the effect of viscosity must be taken into account not only near the wall but also in the critical layer. The critical layer is a narrow region surrounding the critical point at which the main flow velocity and the wave prepagation velocity are equal, that is, U = c. In addition, he also showed that the influence of viscosity leads to instability only if the main flow velocity profile is other than a straight line. The method developed by Tollmein, based on Asymptotic Theory, provided the mathematical basis for later progress in the stability area. Lin (1945, 1946, 1955) was able to provide a firm mathematical basis for the asymptotic expansion theory and was able to explain the nature of the functions near the critical point. He discussed what he called the inner viscous layer which includes the critical point, and the outer viscous layer, a wall viscous layer. The asymptotic expansion method was used almost exclu- sively until the advent of modern high-Speed digital computers which permit the use of more accurate numerical techniques. Even now, however, the asymptotic method can be effectively used to predict the type of functional behavior to be expected, prior to obtaining numerical solutions. The stability of plane Poiseuille flow was investigated by Thomas (1953) with a numerical scheme. He obtained a value for the minimum Reynolds number, Rcr for neutral stability it’ of 5780, which is based on maximum channel velocity and the half-width. This value has been shown to be more accurate than Lin's (1945) value of 5300 or Stuart's (1954) value of 5100 based on asymptotic techniques. Potter (1965) studied the stability of plane Couette- Poiseuille flow by asymptotic expansions and later (1967) per- formed numerical calculations for symmetrical parabolic flows. The values obtained for RC were in close agreement with those rit of Thomas. The point of instability as determined theoretically and the physically observable transition point from laminar to turbulent flow often differ considerably. An explanation for these differences was thought by some to be due to the fact that the derivation of the Orr-Sommerfeld equation is based on the assumption of two-dimensional disturbances only. Squire (1933) showed that if three-dimensional disturbances are considered the flow is more stable, that is, a higher value of Rcrit is predicted,than for two-dimensional disturbances. The distance between the point of instability and the actual transition point depends upon the degree of amplifica- tion present and the intensity of fluctuations present in the primary flow. Schlichting (1933) performed calculations for boundary layer flow over a flat plate and investigated the parameters in the interior of the neutral stability curve (Ci > 0) to help explain the actual mechanism of disturbance amplification. More recently, Shen (1954) repeated Schlichting's calculations, and Stuart (1956) investigated the amplification of unstable disturbances by accounting for the effect of the non-linear terms in the equations. Reynolds and Potter (1967) considered the instabilities of channel flow for disturbances of finite amplitude. Except for early pipe flow measurements by Barnes and Coker (1905) and Ekman (1910), who succeeded in maintaining laminar flow for fairly high Reynolds numbers (40,000), experimental verification of the results of stability theory was slow in coming. Rosenbrook (1937) found agreement with the inflexion point theorem due to Rayleigh and Tollmein. Some of the most significant experimental work.was per- formed by Schubauer and Skramstad (1947) and Dryden (1947) who performed very precise measurements for boundary layer flow over a flat plate (with very low free stream fluctuations), and showed the influence of free stream disturbance intensity on the critical Reynolds number. Emmons (1951) observed that any disturbance which trig- gers transition may be "local in time" and once initiated, the turbulent Spot moves downstream growing steadily in all di- rections. This phenomenom was studied by Schubauer-Klebanoff (1956) . Transition from laminar to turbulent flow in a boundary layer is now believed to take place within 4 stages. At the first stage, infinitesimal two dimensional waves called Tollmein- Schlichting waves, begin to amplify and become unstable. The two dimensional waves become three dimensional and result in hairpin eddies at the second stage. In the third Stage low Speed turbulent streaks or bursts (Emmons' spots) originate near the wall, and finally in the fourth stage the burst rate becomes constant and the transition to fully turbulent motion is completed. Morkovin (1958) reviews some of the recent advances in the study of transition and discusses the mechanisms involved in the above mentioned stages. Stability theory yields a critical Reynolds number that corresponds to stage one. Since the third stage is the first point at which large scale variations take place, this is often considered to be transition by many engineers. These dif- ferences along with the slower reSponse times of earlier instru- mentation serve to explain some of the discrepancies between theory and experiment. Stability predictions in channel flow yield critical Reynolds numbers that also correspond to infinitesimal dis- turbances but the stages of transition are not as apparent as in boundary layer flow. Free steam disturbances or distur- bances which result from wall roughness amplify and lead to the transition described above but the effect is now propagated throughout the flow and the entire channel becomes turbulent. In the area of magnetohydrodynamics (MHD), the velocity and magnetic field equations for an incompressible, viscous and electrically conducting fluid moving in the presence of a magnetic field have been derived by Batchelor (1950). The effect of a magnetic field on thermal instabilities was investigated independently by Thompson (1951) and Chandrasekhar (1952). For a complete discussion of this problem and others in the field of hydrodynamic and hydromagnetic stability the reader is re- ferred to Chandrasekhar (1961). Stuart (1954) considered the stability of viscous flow between parallel planes in the presence of a co-planar magnetic field. Lock (1956) investigated a similar problem but con- sidered the effect of a magnetic field perpendicular to both the confining parallel planes and the flow direction. Both Stuart and Lock simplified the resulting sixth order set of equations and solved the fourth order Orr-Sommerfeld equation by asymptotic expansions. To simplify the analysis, Lock utilized Squire's theorem, as detailed by Michael (1953) for the case of MHD flows. 1.2 Description of the Problem The purpose of the present study is to investigate the stability of an electrically conducting fluid flowing between parallel planes Subject to the influence of a transverse mag- netic field. The coordinate system is oriented with the origin at the centerline, the x-axis in the direction of the flow and the y-axis perpendicular to the bounding plates as shown in figure 1. The planes are assumed to be non-conducting and located at y = _-l_- a. The governing equations, developed in the next chapter result in a coupled set of ordinary linear differential equa- tions, a fourth order equation on y, the velocity perturba- tion, and a second order equation on ¢ the magnetic field perturbation. The above problem was considered by Lock in 1956. He had to make several simplifying assumptions in order to obtain a solution by asymptotic expansions. The reduced equation that he solved was the following modified Orr-Sommerfeld equation: (U-c) (1" - azw) - [W = - £1;in (1.2.1) where the only effect of the magnetic field is to modify the primary velocity profile U(y). The eigenvalues that appear in most stability equations normally occur in the combinations of (a2 + 82), aR, and cr’ and thus, Squire's Theorem is applicable. In the govern- ing coupled equations for this problem 0 appears in the magnetic equation as a separate coefficient indicating that the Solutions may be subject to an influence of B: the Span- wise (z-direction) component of the disturbance wave. Chawla (1969) in studying the effect of rotation on the stability of flow over a flat plate, found that the rotational effects were directly coupled to a and Squire's Theorem was not applicable. Although Lock initially assumed three-dimensional disturbances, his simplifying assumptions delete the magnetic equation and he is then justified in applying Squire's Theorem. This study provides a solution of the complete set of stability equations for various magnetic quantities and fluid prOperties and also provides bounds for which Lock's assump- tions are justified. CHAPTER II STABILITY EQUATIONS 2.1 Fundamental Equations Consider first the interaction of the electric and magnetic fields which are given by Maxwell's equations. These equations written for a non-relativistic reference frame are v - D = pe (2.1.1) V -B’ = 0 (2.1.2) v x E = - "13’ (2.1.3) v xfi = 3+5 (2.1.4) (.71 II where Electric flux density ml ll Electric field on II Magnetic flux density H = Magnetic field p8 = Space charge density 3 = Conduction current density and the notation (°) implies if , t being time. In addition to these equations, the current conservation equation v-If +5, = 0 (2.1.5) is often used; although, it is not independent of Maxwell's equations and follows directly from (2.1.1) and (2.1.4). 9 10 The conduction current is given by Ohm's Law 3’ = 065 +\7 x1?) - peV (2.1.6) where g is the electrical conductivity, and V is the velocity vector. The constitutive equations for a linear medium are —0 as + '13 (2.1.7a) 61 ll W1 II 1161' + fie) (2.1.7b) where the polarization P and magnetization. M; vectors can be neglected for conducting fluids; n is the magnetic permea- bility and e is the permittivity, also known as the dielec- tric constant. For most materials, and in particular for all conduct- ing liquids and gases u is that of free Space no. The assump- tion that e = so, however, can be used only for plasmas, but by using the current conservation equation (2.1.5) this assump- tion need not be made. Following the procedure outlined by Chandrasekhar (1961) some simplifications can now be made to arrive at the equations generally used to solve MHD flow problems. The assumptions are: 1. The non-relativistic approximation has already been stated and is consistent with the Newtonian form of the equa- tions of motion that will be used. 2. Without an externally applied electric field, the electric fields originate only from the induced effects and -o are of the same order of magnitude as ‘V X B appearing in ll Ohm's Law. 3. High frequency phenomena are not considered so that the displacement current D is neglected in equation (2.1.4) compared to 3, the conduction current. In fact for metals, the displacement current is meaningless and need not be mentioned. 4. In Ohm's Law, which determines the conduction current, the Space charge pe may be neglected. For liquid conductors and in dense, collision-dominated plasmas, which may be treated by a continuum model as an ordinary conducting gas, the Space charge effects become unimportant. Ohm's Law is then written as -o 3’ = U(E + vx'fs’) (2.1.8) where it is further assumed that the conductivity is constant with frequency and independent of the magnetic field. The small electric field and the negligible diSplacement current imply the main interaction is between the magnetic field and the fluid, hence the magnetohydrodynamics. The Navier-Stokes equations governing the motion of an incompressible fluid are L a a 1 2w 1 a V + (V - v)V = E'Vp + v V V + E'F (2.1.9) -—~O where the body force term is F = J X B, p is the pressure, p is the mass density and v is the kinematic viscosity. The equation of continuity is v - V = o . (2.1.10) 12 The equations now relevant to the problem are: Maxwell's equations v - fi' = 0 (2.1.11) v xfi = 3 (2.1.12) v xE = nil (2.1.13) Ohm's law: 3' = 063 +1} x L{13) (2.1.14) Continuity: V . V = 0 (2.1.15) Navier-Stokes: 2 V+(V-V)V=‘£(3xfi)- Vp+wx7 (2.1 l P P Eliminating the electric field E between equations (2.1.13) and (2.1.14), together with (2.1.12) leads to the magnetic equation . —o --0 2—0 V (V x H) #0 This equation along with (2.1.15), (2.1.16) and (2.1.12) are sufficient to determine all the variables, V, H, and p. 2.2 Non-Dimensional Variables To write the equations in non-dimensional form, choose the average velocity Um, the channel half width (a) and the applied magnetic field strength Ho as reference quantities. The dimensionless variables will then be .16) V H (2.1.17) l3 at 17* X1 v =6... Xi =2._ m 7t * Umt __R_ (2.2.1) t = p - a U2 p m ...-k «at if = E. 3 = Ls H H O 0 where xi = (x, y, z) the cartesian coordinates, and the asterisks denote the previously used dimensional quantities. Introducing these in the governing equations yields . 2 V+N-v)U=-vp+%VZV+'EE-(3XH) (2.2.2) m H-vx (Vxfi)=%—v2fi (2.2.3) In 3 = v x ii (2.2.4) v ~11 = 0 (2.2.5) v - \7’ = 0 (2.2.6) Uma where R = Reynolds number = —;- M = Hartmann number = u Ha(SZ~-)15 pv Rm = Magnetic Reynolds no. = Uma pg It should be noted here that the above dimensionless groups are not unique and that others may be formed by multiply- ing togehter various combinations of the above. One Such para- meter that will be used later is the magnetic Prandtl number R -.JE Pm vuc R . 2.3 The Primary Flow and Magnetic Field Distributions The solution for steady, two-dimensional motion of a conducting fluid between parallel planes subjected to a transverse magnetic field is well known and has been provided 14 by Hartmann and Lazarus (1937). For parallel flows, the velocity components are (U, 0, 0) and for this problem the magnetic field has the components of y and only. Equation (2.2.3) yields 2 l d Hx d — = - —-— . .1 '-“ dy (UHy) (2 3 ) Rm dyz (H , H , 0), all functions X y ———l=0 (2.3.2) Since the normal component of pH must be continuous and have the same value at both walls and since the applied magnetic field induction is in the y direction, equation (2.3.2) also be tinuity reSpect shows that Hy is a constant H. This result can obtained from equation (2.2.5) which represents con- of the magnetic field. Differentiating the x -component of (2.2.2) with to y yields 2 2 d H 3 'E‘ "235- - 1% (2.3.3) m dy dy Eliminating Hx between (2.3.1) and (2.3.3) results in d3U 2 dU "j; = M a; (2.3.4) dy The solution to the above equation is M(cosh M - cosh My) U = (2.3.5) M cosh M - sinh_M 15 For zero magnetic fields, this equation should reduce to the standard parabolic profile, U = 1.5 (l-yz). This is found to be the case when the first two terms of the series expansions are substituted for the hyperbolic functions. To determine the induced magnetic field in the direc- tion of flow, substitute equation (2.3.5) into (2.3.3). The resulting equation and the boundary conditions that Hx must be zero at both walls since there is no applied magnetic field in this direction yields Rm (Sinh My - y sinh M) Hx = M(cosh M-l) 2.4 The Linearized Equations for Small Disturbances The primary velocity and magnetic field quantities are now considered to have superimposed on them three-dimensional infinitesimal disturbances. They are then written V=V +3 i 4p i (2.4.1) H = H + h p where tha (U, 0, 0) <1 :1 U U 2? B x x <: :r: <: o N V V :1 u 9 :J‘ :3“ Substituting (2.4.1) into equation (2.2.3), Subtracting off the original equation for the primary flow and neglecting the squares of all the disturbance quantities yields 3' = dip-v)? + (fi-vfip - (VP-ml? - ('v’-v)ifp + 11— 23 (2.4.2) In 16 The pressure gradient term is first eliminated from the momentum equation by taking the curl of both sides. The dis- turbance equation is then found as above, and the vorticity \ equation for the perturbations results. :0 -o -b —+ a _ v-o . -—b -0. l 2-0 5 + (VP-WE + (v-vmp - (op v)v + (g va + R v s 2 +L RRm where 6p = vorticity of main flow = v X § = vorticity of perturbed flow = v X v 3=(VXH) P P '5 =(vx'fi) V p -—1 ~*.'-‘ “2"-”.‘°-“-‘."' 2.4.3 {(Hp v)J+(h v)Jp (Jp v)h (J <1)le ( ) In addition, the following continuity equations result v o 3 = O :71 II C v 0 (2.4.4) (2.4.5) Consider further that the assumed three-dimensional disturbances take the separated form vX = Vx(y) exp[i(ax + 82) vy =10) exp[i(m< + 82) v2 °z(y) exp[i(ax + 32) hx = hx(y) exp[i(ax + 82) h " (My) exp[i(orx + 82) 3‘ II 3"? A "C v exp[i(ax + 52) iact] iact] iact] iact] iact] iact] (2.4.6) The assumed form of the disturbances implies a Spatially periodic wave with complex amplitudes where a and a are dimensionless wave numbers and are real quantities. The complex “I!“ .' 9. 77“]... 1 A , . l7 wave speed c is given by C = Cr 4. “1 (2.4.7) where c1 is the amplification rate. A positive or negative ci implies growth or decay respectively of the perturbations. This study is concerned with neutral stability, that is, ci 8 0. Introducing the assumed form of the disturbances into the component form of equations (2.4.2) and (2.4.3) and then eliminating vx, vz, hx’ and b2 with a procedure similar to that outlined by Chawla or Stuart (1954), yields the following 1 v 1 ,, 2 11¢ - 34’- - (U-c)¢ WT“: (¢ - K ¢) (2.4.8) (U-c) (1);" - K211) - UN + i: “iv - 2K2¢" + KW) (2. .9) 2 4 ' P{h(¢" ' K20) - “i (0"'- K 0') - h"¢} Rm (sinh My - y sinh M) M (cosh M - l) where h - H B x M (cosh M - cosh My) U a M coshM - sinh Mk K2gmz+az M2 P-— RR m and primes indicate differentiation with respect to y. Equations (2.4.8) and (2.4.9) can be solved simulta- neously or can be combined to form a single sixth order equa- tion. Combining the equations requires tedious algebra and the resulting coefficients are extremely long and unwieldy. No insight or simplification is gained by such a combination 18 so that it is better to solve the two equations simultaneously. The necessary boundary conditions result from the no- slip velocity conditions at the wall and from continuity, namely, €- 0 €- n c> (a '< n l+ H and the nature of the magnetic field perturbations (2.4.10) The problem is now completely Specified. 2.5 Eigenvalue Problem The system of equations and boundary conditions derived above, represent an eigenvalue problem with the characteristic values a, B, c, R and M. The wave propagation Speed, c as stated earlier is complex, with the imaginary part being the exponential growth or decay rate of the assumed disturbances. Since neutral stability curves are desired ci is set to zero. To solve for the characteristic values it is necessary to Specify some of them, say M, B and cr and solve for the remaining two, namely a and R from equations (2.4.8) and (2.4.9). The neutral stability curve (0 vs. R for constant M) will have a minimum R, called the critical Reynolds number crit’ for which a disturbance is neutrally stable. Reynolds numbers greater than Rc result in growth for that parti- rit cular disturbance and Reynolds numbers smaller than Rcrit result in decay. Associated with Rcrit there is also a ) critical wave Speed (cr crit’ and a critical wave number ’1’ . ' CI'LC JI CHAPTER III NUMERICAL SOLUTION OF THE STABILITY EQUATIONS 3.1 Introduction The numerical solution to equations (2.4.8) and (2.4.9) will generate three independent solutions because the solutions are started at one wall with three boundary conditions already satisfied. These solutions cannot each satisfy all of the boundary conditions at both walls. However, a prOper linear combination of these functions will yield the total eigen- function which must then satisfy the three boundary conditions at the opposite wall. The integration of the equations is begun at the lower wall (y = -l) and proceeds Step by step across the channel to the upper wall. A fourth order Runge-Kutta technique, de- tailed in appendix A is used to solve the equations. The three independent solutions are each initialized at the lower wall and integrated simultaneously at each step across the channel. 3.2 Starting Conditions To use the Runge-Kutta integration scheme to solve a differential equation of order n the problem must be trans- formed to an initial value problem where the function and its n-l derivatives are initially specified. For equations 19 20 (2.4.8) and (2.4.9), the boundary conditions (2.4.10) provide starting values for 'i’ w; and $1 (i = l, 2, 3 and repre- sent the separate solutions) with values for the remaining derivatives $2, ¢;" and ¢i somewhat arbitrary. The highest order derivatives, namely 11v and ¢;, do not require ini- tialization since they are determined in terms of the lower order derivatives. The "arbitrary" starting conditions mentioned above must be chosen so as to insure independent functions, at least at the start of the integration. A purification scheme maintains independence as the integration proceeds. Assigning a non-zero value to one of the three unSpecified derivatives for each of the three solutions should help keep the solutions linearly independent. Specifically $3, 45", and $1 are given non- zero values. The purification scheme requires identification of the fastest growing'solution. This was accomplished by checking the ratio of -l at the first integration step for each of the independentisolutions. The values obtained for i = 1, 2, 3 were approximately 400, 300, and 200 respectively and, as Should be expected, were unaffected by the relative magnitudes of the starting values assigned to $3, V3 and ¢i. 3.3 Growth of the Eigenfunctions It is well known that during the numerical integration one of the two independent functions for the fourth order prob- lem and at least one of the three independent functions for the 21 sixth order problem grows very rapidly. Kaplan (1964) referred to these functions as the "growing solutions" and called the others the "well behaved solutions". The growing solutions stem from the viscous portion of the stability equations and the well behaved functions originate from the inviscid portion. The governing stability equations can be written as [hi - i’i' - (U-c)¢] = —£- [¢" - K2 3 (3 3 1) a aRm ¢ ° ° 2 ' - 2 4 [(U-c)(¢" - K W) - u"¢] = if [11V - 2K 1" + K ¢ 2 (3.3.2) - 3:31" 01(4)" - K20) - 5; (0m ' K205 - Mel] m The terms in the square brackets on the left hand side represent the inviscid part and those in the brackets on the right hand side the viscous part. In this particular problem, for the larger Hartmann numbers and magnetic Prandtl numbers the three solutions exhibited three distinct growth rates, which could be termed the largest growing solution, intermediate growing solution and the well behaved solution. As an example, for the case of M = 4.0 and Ptn = 10‘.1 the three functions exhibited growths across the channel on the order of 10280, 10140 and 106 respec- tively. This extreme growth limited the range of Hartmann numbers that could be considered in this Study as it is apparent that the limitations of the computer are Soon exceeded. For M = 0 (.001 for the computer program) the equa- tions effectively reduce to fourth order and result in only one growing solution which exhibited a growth on the order of 22 1036 across the channel, which is in agreement with the observations reported by Reynolds and Potter (1967). The final combined solution, of course, does not exhibit this rapid growth which indicates that only a very small portion of the growing solutions are required to form the eigenfunctions. 3.4 Purification Scheme The very rapid growth rate exhibited by some of the solutions causes some difficulties other than machine over- flow. Initially all three solutions are linearly independent but as the integration proceeds, this independence is observed to disappear rapidly. Kaplan (1964) States that this loss of independence, which is impossible for an exact solution, arises because of the approximate nature of the numerical integration. Errors are introduced because any numerical method being applied at small but finite steps has associated with it a truncation error and in addition a digital computer carries only a fixed number of digits, resulting in round off error. Kaplan further concludes that the arbitrary initial conditions for the well-behaved solution contain a small portion that is also an initial condition for the growing solutions. As the integra- tion proceeds it grows much more rapidly than the behaved solu- tion and soon dominates it. Even if this small portion is not present initially, it is effectively introduced by truncation errors . 23 Kaplan's conclusions appear justified although the actual mechanism that causes one solution to "pollute" the other may be open to question. The goal is, therefore, to keep the solutions independent over as large a range as possible. To reduce the truncation error associated with any numerical scheme, an obvious answer is to choose a step size that is as small as possible con- sistent with the particular limitations of machine storage and speed. Performing all arithmetic operations in double pre- cision should greatly reduce the round-off error. This is in fact, found to be true. However, since all the functions are complex, adding double precision to the program significantly increases machine computation time and storage requirements. Kaplan suggests an alternate approach that does not require double precision, namely, a filtering technique. His method consists of subtracting from the well-behaved solution a portion of the growing solution at every step of the integra- tion. This procedure called "filtering" prevents the growing solution from ever dominating the well-behaved solution and thus maintains the needed functional independence. Kaplan's shceme was implemented in the computer program for this problem. The CDC 6500 computer nominally carries about 15 significant digits in single precision, which is equivalent to double precision on IBM equipment and in parti- cular the IBM 7090 used by Kaplan. These computations were in effect then performed in "double precision". In fact, for the cases with Hartmann number equal to zero or one no suppression 24 scheme was required. Test cases of M = 2 or 3 revealed that Kaplan's scheme was required. For M > 3 however, the Scheme was no longer sufficient to maintain independence. In this range the growth becomes very large, and the second solution becomes an intermediate growing solution which now pollutes the well- behaved solution. A double Suppression scheme Suggested by W. Reynolds of Stanford in a private communication to M. Potter was used to find the two remaining solutions. The filter used consisted of a ratio of the inviscid solutions (left hand side of equation (3.3.2)), namely, inviscid part of the well-behaved solution Filter = , , _ , , 1nV1sc1d part of the growing solution (3.4.1) where, (inviscid) m icyR[ (1mm); - K2111?" u'wm] (3.4.2) m 1, 2, 3 The above ratio determines the fraction of the growing solu- tion to be "extracted" from the behaved solution, so the amount subtracted off is the product of this ratio and the value of the growing solution at the particular integration step. This product was equal to about 20% of the behaved solution. The final Scheme consisted of extracting from solutions two and three a portion of the fastest growing solution (one) and then extracting from solution three a portion of the inter- mediate growing solution (two). Once incorporated this method assured complete independence of the three separate solutions for all growth rates even to the point of exceeding the over- flow limits of the machine. 25 3.5 Iteration Scheme for the Eigenvalues Integrating across the channel, three independent solu- tions are generated at each step. Upon reaching the opposite wall, the functions are linearly combined to form the total eigenfunctions that must satisfy the boundary conditions. The three conditions that must be satisfied are W = W. = ¢ = 0. Consider the following set of combined functions at the wall: alwl +- a2¢2 + a3¢3 = (V (3.5.1) all; '+ 8215 + 8315 = t; (3.5.2) alt; +- azyg +-a3¢g - 13 (3.5.3) a1¢1 + a2¢2 +a3¢3 = Qw (3.5.4) Equations (3.5.2-4) are used to solve for the coeffi- cients a1. Note that there is no Specific boundary condition for '3’ hence the choice here is arbitrary. As long as 'w is non-zero, the value assumed merely changes the normalization factor and still represents a valid solution. This system can now be written as " . '1 r 1 F 1 *1 '2 *3 a1 0 *1 '2 '3 a2 ' 1 H1 ('2 '31 183. 1°. For functions that are linearly independent the determinant of the matrix containing known values of v', V" and ¢ will be non-zero. When the double suppression scheme, discussed in the previous section, was not used this determinant 26 was effectively zero and indicated the functions to be linearly dependent, that is, one column is a multiple or combination of the other two. For independent functions, the Si can be determined by finding the inverse of this matrix or as was done for the computer program by simply writing out the solution. The ai, once determined can now be substituted into equation (3.5.1). With the correct eigenvalues Ww will be zero; hence ww will serve as the test function (T). If T)<(conjugate T) is less than 10-12 the convergence criteria is satisfied and the eigenvalues used to generate the functions are assumed to be the correct ones. If convergence is not attained let T1 = T. Increase a by T1, recalculate T and let T = T. After setting a to its original value, 2 increase R by 1%, calculate T and let T3 = T. The finite difference approximations for the change in T with respect to y and R are a-T— .. 2 i 3 5 r Ad Ad ( ' '3) T - T AT- ~ 3 1 (3.5.6) 6R AR These are substituted into the complex equation 5361 + SEAR + T = 0 (3.5.7) as 5R 1 from which Ag and AR can be calculated. The new values for the eigenvalues are 27 0[new - Gold + AU (3.5.8) R new old + AR R It is apparent that for every iteration or new "guess" of eigenvalues the equations must be integrated across the channel three times. It was found that if the initial guesses are relatively good, convergence is obtained in about three iterations. 3.6 Criterion for Guessing the Eigenvalues The initial estimates for some of the eigenvalues were obtained from the data presented by Lock (1954). Where his data was not applicable (M > 3, Ph I 10.1), guesses had to be made with extreme care. Hopefully, the iteration process would iterate to the correct eigenvalues, although this did not always occur because of the small radius of convergence. After a few points were obtained and plotted, the next estimates could be obtained by extrapolation until the complete stability curve was generated. 3.7 Range of Parameters Considered Hartmann numbers in the range of zero to six were run 6 and 10-1. Test for magnetic Prandtl numbers (Pm) of 10- cases for magnetic Prandtl numbers less than 10-6 were run with no change in eigenvalues. Physically, Prandtl numbers in the range of 10-6 to 10“4 are characteristic of liquid metals and it was hoped that the study could be extended to Pm = 100, which is characteristic of ionized gases. At these values, 28 however, the functional growth rates became excessive, result- ing in machine overflow. As an example of the large variations associated with a change in Pm refer to the case of M = 6 as plotted in Figure 9. Included on this figure is a neutral stability curve for Pm = 10..2 in addition to 10-6 and 10-1. The change -6 for P = 10 to 10-2 is only about one-ninth the in Rcrit m change from 10"2 to 10-1. It is apparent then that a further increase of only one order of magnitude with its corresponding increase in functional growth rates soon puts the values out of manchine range. Several values of B were tried at different Hartmann numbers and magnetic Prandtl numbers but, as will be discussed in the next chapter, no B-effect was found to exist except that predicted by Squire's Theorem. 3.8 Special Case of Zero Hartmann Number The general equations are sixth order, but for the case of zero magnetic field the equations reduce to the standard fourth order Orr-Sommerfeld equation. Rather than write a separate computer program to solve this case it was felt that if the general program could be used it would provide a better check on proper program operation. Hopefully, then, as M approaches zero the program Should yield the known results for plane Poiseuille flow. This was found to be the case when M was set equal to 10-2 or 10.3, with the same eigenvalues resulting for either value. A small 29 value was necessary since setting M equal to zero exactly yields an indeterminate result (0/0) in the calculation of the primary velocity profile. It also results in a reduction in the order of the equations thereby causing an indeterminate result when satisfying the boundary conditions. CHAPTER IV RESULTS, CONCLUSIONS, AND RECOMMENDATIONS FOR FURTHER STUDY 4.1 Numerical Results Prior to obtaining any points on the stability curves, the subroutine which calculates the primary velocity profile and induced magnetic field quantities was run as a separate program. The non-magnetic case, as was discussed in the previous chapter was run at M = 10-3. The results of this run were compared to the parabolic profile for plane Poiseuille flow and found to agree to within at least five significant numbers. This profile along with others for M > 0 are well known, see for example, the original work by Hartmann and Lazarus (1937) or any current textbook on MHD. Dimensionless plots of these curves are presented in Figure 2. The data generated for the neutral stability curves are presented in Tables 2 through 15 and are also plotted along with Lock's data in Figures 3 through 9. A summary of the critical wave numbers and Reynolds numbers from this study and from Lock's data are presented in Table 1. Typical eigen- functions are plotted in Figures 10, 11 and 12. To reduce the equations to a form that could be solved by asymptotic expansions, Lock, based upon Stuart's paper (1954), used the assumption that for most conducting liquids, 3O 31 the magnetic Prandtl number Pm ( Rm/R R Vuo) is small. Provided therefore, that R is not too large, Rm will be small compared with unity. He then neglected the magnetic effects and reduced the sixth order set of equations to the following fourth order, modified Orr-Sommerfeld equation which he then solved: 2 . . (II-cm)" - a1) - u") = - :1; 11" (4.1.1) In reducing the problem to this form, the only effect of the magnetic field is to modify the primary velocity profile. Comparison of Lock's data as presented along with the present data is surprisingly good. For Pm = 10-6 his critical points are a little high for M = 0 and 1, agree almost exactly at M = 2 and for M 2 3 are low compared to this study. As M increases the deviations increase, which indicates a greater dependence of the solutions on the magnetic terms of equations (2.4.8) and (2.4.9) that were neglected by Lock. As discussed in the problem description, an examina- tion of the equations indicates that B, the spanwise component of the disturbance wave.should affect the final eigenvalues such that Squire's Theorem is not applicable. This was not found to be the case, however. The term 4' in equation Q lr-I- (2.4.8) that leads to this conclusion does not Significantly affect the final solution. This was verified by several runs for various values for M and Pm with a non-zero (0.5) value of 3, all of which produCe the same eigenvalues as the B = 0 ORDER NO. DATE ISSUE W.O. N0. YEAR 3 . ,8 [ ($4 4.1 Co ri ht I ’ N . . SCHOOL V/n/L (-211 1 :3, U , ° P°S [ \J/ 4‘2 I I I Not C/R i I M. No. #153977; )éflflH—ow Total 4 ,3 , Po,” to cm W)& 1 W, W PHYS. DESC. 2 1 4 14 I Pos. to School LC C/R SCH. 1 2 25 B PUBL. N0. N0.PAGES PRICE 1 PCS. AUTHOR NAME (LAST, FIRST, MI.) 1—2 25 31 33 36 37 = P 39 57 Format A 16 22 32 35 50 54 59 77 FcTr'Efét B 2.1 7/101631/01/ l 135 l L I 1 jg Dal/ACl/fgl- K sJfl/Am 611‘s: ‘14: l I l’ Cust. No. Adv. Pav No. Adv. Pav Amt. —2 2-8 10 l6 17 22 45 50 51 57 0 A I l I I l 1 O l 4 2 0 0 a a l A n 1 I 1 l a - - CHG. QTY. Service or Product U anmp'nfc AUTH- 1—2 9—11 13 15 16 22 23 27 N0. AMI, FROM Micro and Abst. 0 9 X, . PD00421 00001 @ $20.00 Ea. Micro and Abst.. 0 9 XI I . . PD00421 00002 @ $25.00 Ea. Publ. Abst. only 0 9 XI . I I PD00421 00003 (a $15.00 Ea. Publish Diss. & Store 0 9 X. . . . P1300429 00000 Npgariivg (an12 oo Ea t 0 9 x I 9000422 00000 fizzygaid meta. Copyright 0 I I I I PD00422 00001 Paid in ... Extra Pages in Abst. 0 9 x. . . P1300433 00000 (a $2.40 per page. Extra Listings in DA 0 9 XI .. I I PD00434 00000 @ SLM. — 0 9 A l I 0 9 4 2 6 Format Positive Copies to . . "B' C/R Office and/or School 0 9 I I I I ADDRESS Mfix @ WW 546% 770 ‘ . ) f/ 1 I“ X c/r p. as rec'd 04' I‘ 54/ made by ed. / ._____ stripped by ed. W "a" pp. by author by ed. (7 __ ~ nnbc CS 5 0C Lt. & dk.type W I ., . NOCbg Pap. prnt. 5g 1 7 1 I , Q / NOC _;_Mntd. illus. _ 1 r 0 Ditto Glos. photo H [m 1 7 I Carbon Color photo / -4 CONDITION: 3/5 / I H u M; (I) (11$ . Foreign Text all____ some P—fix\; Moi/n.5,”- 32 case. The influence of increasing Hartmann number is, of course, to stabilize the flow. For a given Hartmann number, a further stabilization of the flow occurs with an increase in the magnetic Prandtl number. These effects are presented graphically in Figures 3 through 9. Consider, for example, the case of M = 5.0 which is depicted in Figure 8. The critical Reynolds number increases from 132300 to 169400 for -6 -l P = 10 and 10 respectively. This increase due to the m change in Pm stems from the fact that the eigenfunctions for the velocity perturbations W and the magnetic field perturbations ¢ at PIII = 10-1 are now of the same order of magnitude. Hence, for this range the magnetic terms in the stability equations significantly affect the final eigenvalues. The magnitudes and behavior of the eigenfunctions for this case are shown in Figures 11 and 12 where V and ¢ are plotted reSpectively. The 4 curves have been normalized to 1 + 0i at the centerline as is commonly done in the current literature. Figure 9, which presents data for M = 6.0 includes in addition to the neutral stability curves for Pm = 10”6 and -1 - 10 a curve for P = 10 2. The shift in RC m , for a change 1711: in Pm from 10-6 to 10.2 is only about 6000, from 169500 to 175800, but a value of Pm = 10-1 Significantly increases Rcrit to a value of 230700, indicating the large influence of the magnetic terms in the stability equations. 33 As discussed in Chapter III the growth rates of the independent eigenfunctions increase rapidly with increasing Hartmann number and magnetic Prandtl number. The growth rates with this numerical scheme become so large for M > 6 or Pm > 10.1 that machine overflow occurs and no data could be obtained. This phenomenon is evidently also due to the ¢ eigenfunctions. 4.2 Conclusions Based on the results obtained and the preceeding dis- cussion it can be concluded that for small magnetic Prandtl numbers (10-6) characteristic of liquid metals the standard Orr-Sommerfeld equation along with the modified velocity pro- file gives satisfactory results in determining neutral sta- bility. For other conducting liquids or gases where the magnetic Prandtl is larger, the effects of the magnetic field must be fully accounted for. The effect of an increasing magnetic field strength is to greatly stabilize the flow and also to increase the critical wave number causing the in- stability. Squire's Theorem is applicable to this problem since 6 does not affect the solutions, and K2 which equals 2 2 a + B can be replaced by q2 in equations (2.4.8) and (2.4.9). 34 4.3 Recommendations for Further Study This study could be further generalized to include conducting non-newtonian fluids which would include liquid metal amalgams, uranium slurries, seeded polymers, and others. The effects of geometry changes such as a finite channel with conducting side walls that would lead to different boundary conditions on the magnetic field could also be studied. Also of interest would be an investigation to compare the results of this study with stability curves obtained for ci # O. i1 J ....3 .-J! B 18 LIOGRAPHY BIBLIOGRAPHY Barnes, H. T. and Coker, E. C., "The Flow of Water Through Pipes". Proc. Roy. Soc. London, A 74, 341 (1905). Batchelor, G. K., "On the Spontaneous Magnetic Field in a Conducting Liquid in Turbulent Motion". Proc. Roy. Soc. London, A 201 (1950). Chandrasekhar, 8., "On the Inhibition of Convection by a Magnetic Field". Phil. Mag. Ser. 7 43, (1952). Chandrasekhar, 8., "Hydrodynamic and Hydromagnetic Stability". Clarendon Press (1961). Chawla, M. D., "The Stability of Boundary Layer Flow Subject to Rotation". Doctoral Dissertation, Michigan State University (1969). Collatz, L., "Numerische Behandlung von Differentialglei- chungen". Springer-Verlag, Berlin (1951). 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. Appl. Mech. I, Academic Press (1948). Ekman, V. W., "Archiv fur Math". Astr. och Fys. VI, No. 12 (1910). Emmons, H. W., "The Laminar-Turbulent Transition in a Boundary Layer". Part 1, Journal of Aeronautical Sciences, 18 (1951). Hartmann, J. and Lazarus, F., Math.-fys. Medd. 15 (1937). Kelvin, Lord, Philosophical Magazine. 5th Series t. xxiv (1887). Lin, C. C., "On the Stability of Two-Dimensional Parallel Flows". Part I, II, III, Quart. Appl. Math. 3, (1945, 1946). Lin, C. C., "The Theory of Hydrodynamic Stability". Cambridge University Press, Tendon (1955). 35 36 Lock, R. C., "The Stability of the Flow of an Electrically Conducting Fluid Between Parallel Planes under a Transverse Magnetic Field". Proc. Royal Soc. A 233 (1956). Michael, D. H., "Stability of Plane Parallel Flows of Electrically Conducting Fluids". Proc. Camb. Phil. Soc. 49, (1953). Morkovin, M. V., "Transition from Laminar to Turbulent Shear Flow - A Review of Some Recent Advances in its Under- standing". Trans. ASME, (1958). 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). Potter, M. C., "The Stability of Plane Couette-Poiseuille Flow". Doctoral Dissertation, University of Michigan (1965). Later published in J. Fluid Mech. (3), 24 (1966). 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 Phy . 2;, 23 (1922). Rayleigh, Lord, "On the Stability or Instability of Certain Fluid Motions". Proc. Lond. Math. Soc., 11, 51 (1880) and 12, 61 (1887). Rayleigh, Lord, "Further Remarks on the Stability of Viscous Fluid Motion". Phil. Mag3 and Journal of Science, 6 Series,_2§ (1914). Reynolds, 0., "An Experimental Investigation of the Circum- stances which determine whether the Motion of Water shall be Direct or Sinous, and of the Law of Resistance in Parallel Channels". Phil. Trans. Roy. Soc. (1883). Reynolds, W. C. and Potter, M. C., "Finite-amplitude In- stability or Parallel Shear Flows". J. Fluid Mech. 13;, _2_7_ (1967) . Rosenbrook, G., "Instabilitat der Gleitschichten im schwach divergenten Kanal". ZAMM 11, 8 (1937). 37 Schlichting, H., "Zur Entstehung der Turbulenz bei der Plattenstromung". Nachr. Ges. Wiss. Gottg, Math. Phys. Klasse, (1933). Also, ZAMM, 12_(l933). Schubauer, G. B. and Klebanoff, P. 8., "Contributions on the Mechanics of Boundary Layer Transition". NACA TN 3489 (1955) and NACA Rep. No. 1289 (1956). 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 Layer Oscillations and Stability of Laminar Flow". J. Aero. Sci., 14 (1947). Shen, S. F., "Calculated Amplified Oscillations in the Plane Poiseuille and Blasius Flows". J. Aero. Sci., Zl.(1954)° Sommerfeld, A., "Ein Beitrag Zur Hydrodynamischen Erklaerung Der Turbulenten Fluessigkeitsbewegungen". Atti del Congr. Internat. dei Mat., Rome (1908). Squire, H. B., "On the Stability for Three-Dimensional Disturbance of Viscous Fluid Flow between Parallel Walls". Proc. Roy. Soc. London, A 142, (1933). Stuart, J. T., "On the Stability of Viscous Flow Between Parallel Planes in the Presence of a Co-Planar Magnetic Field". Proc. Roy, Soc. A, 221 (1954). Stuart, J. T., "On the Effects of the Reynolds Stress on Hydrodynamic Stability". ZAMM Sonderheft (Special issue) 32-38 (1956). See also, J. Fluid Mech., 4 (1958) and J. Fluid Mech., 2 (1960). Thomas, L. H., "The Stability of Plane Poiseuille Flow". Physical Review (4), 21 (1953). Thompson, N. B., "Thermal Convection in a Magnetic Field". Phil. Mag; Ser. 7, 42 (1951). Tietjens, 0., "Beitrage zur Entstenhung der Turbulenz". ZAMM, 2 (1925) . Tollmein, W., "The Production of Turbulence". NACA Tech. Memo 609 (1931), originally published in Nachr. Ges. Wiss. Gott., Math. Phys. Klasse (1929). Tollmein, W., "General Instability Criterion of Laminar Velocity Distributions". NACA Tech. Memo 792 (1936), originally published in German in (1935). I LLUSTRATIONS (Tables and Figures) 38 Table 1. Comparison of Critical Eigenvalues M? ngk Present Study 3.2.212 22.2.1.2 Pm = 10-6 :13 1 a R a R 0 1.03 4000. 1.02 3847. - - l 0.98 6960. 0.97 6782. 1.00 6926. 2 0.93 20000. 0.925 20160. 0.90 21180. 3 0.96 46000. 0.95 48230. 0.95 54600. 4 1.04 79400. 1.04 86340. 1.00 105100. 5 1.15 116700. 1.20 132300. 1.10 169400. 6 - - 1.30 169500. 1.20 230700. Note: Lock's original data has been converted to corre5pond to the same non-dimensional variables as used in this study, that is, R based on average velocity rather than centerline velocity. 39 Table 2. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c. = 0.0) M = 0.001, p = 10-6 1 m c a R r Remarks 0.9780 26050. 0.2750 Upper Branch 1.0710 10050. 0.3400 1.0980 5772. 0.3800 1.0500 3924. 0.3997 1.0200 3847. 0.3960 Critical Point 0.9500 4138. 0.3783 0.8500 5426. 0.3417 0.7500 8305. 0.2975 0.6500 14620. 0.2500 0.5756 26190. 0.2100 Lower Branch Table 3. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c = 0.0) M = 1.0, P - 10'6 i m c a R r Remarks 1.0190 19580. 0.2900 Upper Branch 1.0480 10530. 0.3300 1.0340 7866. 0.3470 1.0000 6920. 0.3506 0.9700 6782. 0.3474 Critical Point 0.9500 6853. 0.3436 0.9000 7365. 0.3310 0.8500 8388. 0.3148 0.7500 12320. 0.2757 0.6700 19250. 0.2400 Lower Branch 40 Table 4. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M = 1.0, Pm = 10'1 c a R r Remarks 1.0330 15910. 0.3040 Upper Branch 1.0300 7641. 0.3490 1.0000 6926. 0.3510 Critical Point 0.9000 7369. 0.3310 0.8500 8388. 0.3152 0.7500 12410. 0.2760 0.6500 21920. 0.2311 Lower Branch Table 5. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c1 = 0.0) M = 2.0, PM - 10-6 c a R __ r Remarks 1.0070 40070. 0.2440 Upper Branch 0.9900 22830 0.2719 0.9500 20380. 0.2741 0.9250 20160. 0.2720 Critical Point 0.9000 20440. 0.2685 0.8000 25630. 0.2446 0.6500 53040. 0.1931 Lower Branch 41 Table 6. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c. = 0.0) M = 2.0, P = 10'1 1 m c a R r Remarks 0.9726 66050. 0.2200 Upper Branch 0.9986 40240. 0.2450 0.9900 25510. 0.2676 0.9500 21430. 0.2730 0.9000 21180. 0.2679 Critical Point 0.8500 22870. 0.2579 0.7500 31716. 0.2291 0.6500 53850. 0.1932 Lower Branch Table 7. Eigenvalues for the Stability Equations for the Case of Neutral Stability (ci = 0.0) M = 3.0, Pm -= 10-6 c a R r Remarks 1.0380 134200. 0.1880 Upper Branch 1.0560 74400. 0.2140 1.0300 54770. 0.2267 1.0000 49650. 0.2294 0.9500 48230. 0.2271 Critical Point 0.8600 54680. 0.2138 0.7500 79570. 0.1878 0.6500 136500. 0.1584 Lower Branch CI. 1t 1C8]. P011111: .....s: Table 8. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c. = 0.0) M = 3.0, p = 10-1 1 m _ c a R r Remarks 1.0140 143400. 0.1870 Upper Branch 1.0280 105600. 0.2000 1.0170 66980. 0.2190 0.9900 57850. 0.2238 0.9500 54600. 0.2235 0.9000 55550. 0.2185 0.8500 60540. 0.2102 0.7500 84930. 0.1867 0.6500 143900. 0.1577 Lower Branch Table 9. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c, = 0.0) M = 4.0, p = 10-6 1 m c a R r Remarks 1.1570 178500. 0.1800 Upper Branch 1.1490 117900. 0.1970 1.0900 88910. 0.2069 1.0400 86340. 0.2057 Critical Point 0.9500 93620. 0.1969 0.9000 104100. 0.1892 0.8000 143800. 0.1697 Lower Branch 43 Table 10. Eigenvalues for the Stability Equations for the Case of Neutral Stability (Ci = 0.0) M = 4.0, Pi = 10-1 C a R r Remarks 1.1120 186000. 0.1800 Upper Branch 1.0900 132200. 0.1930 1.0400 108000. 0.1986 1.0000 105100. 0.1975 Critical Point 0.9500 108900. 0.1930 0.9000 119000. 0.1865 0.8000 159700. 0.1676 Lower Branch Table 11. Eigenvalues for the Stability Equations for the Case of Neutral Stability (ai = 0.0) M a 5.0, PM - 10-6 C a R r Remarks 1.2980 188100. 0.1820 Upper Branch 1.2540 148900. 0.1910 1.2400 141400. 0.1928 1.2300 140800. 0.1926 1.2000 132300. 0.1945 1.1500 127400 0.1944 Critical Point 1.1000 131000. 0.1912 1.0000 148900. 0.1813 0.9000 190200. 0.1666 Lower Branch 44 Table 12. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c, = 0.0) M = 5.0, P = 10‘1 1 m c a R r Remarks 1.1900 195700. 0.1810 Upper Branch 1.1700 184700. 0.1827 1.1500 175400. 0.1840 1.1000 169400. 0.1834 Critical Point 1.0500 175000. 0.1799 1.0300 174600. 0.1789 0.9820 186600. 0.1740 Lower Branch Table 13. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c, - 0.0) M = 6.0, p - 10"6 1 m c a R r Remarks 1.466 199030. 0.1840 Upper Branch 1.378 176800. 0.1880 1.350 171700. 0.1885 1.300 169500. 0.1881 Critical Point 1.250 171300. 0.1865 1.200 172800. 0.1844 1.150 179300. 0.1813 1.100 196400. 0.1758 1.050 216000. 0.1701 Lower Branch 45 Table 14. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c. = 0.0) M = 6.0, P. = 10"2 1. m c a R r Remarks 1.450 211000. 0.1810 Upper Branch 1.380 185700. 0.1855 1.300 175800. 0.1863 Critical Point 1.250 179200. 0.1843 1.150 188400. 0.1791 1.100 199700. 0.1752 1.050 220600. 0.1692 Lower Branch Table 15. Eigenvalues for the Stability Equations for the Case of Neutral Stability (c. = 0.0) M = 6.0, p = 10-1 1 m c a R r Remarks 1.280 237400. 0.1765 Upper Branch 1.250 230800. 0.1771 1.200 230700. 0.1756 Critical Point 1.150 234500. 0.1734 1.111 245000. 0.1700 Lower Branch “hangar 1 , E _ I 46 .4 cowumuawwmcoo 30am xumewum 3: ll /% W /z .H opswfim MN //, 1.6 0.8 0.6 0.4 M = 0 47 (Parabolic) L P 1 l l 0.0 0.2 0.4 0.6 0.8 1.0 Y Figure 2. Dimensionless Primary Velocity Profiles for Various Hartmann Numbers 48 o n 2 you huwawnmum Hmuusmz now m uwnevz mUHoczwm msmum> a monasz m>m3 muoH x m mH HM o n _ _ _ _ musum ucmmwum Moog .m snowflm 49 la 1. o.H u z MOM muHHwQMum Hmuusmz now a penssz mvaoc%mm mawum> o sagasz m>m3 .q madman n-0H x m mm om ma 0H «a NH OH m _ _ _ _ _ _ _ _ /’ / OH Em I Cl / #004. I. I/ /’ / a /' III- 'I": I III‘ n.o w.o ¢.o H.a 50 m-oH x m om ow On 00 om oq om ON fl _ _ — a 7 _ — OH u a .II. H- n m L 0H Em III 0| I 1004 i I L . 1|I\\\\ 0.0 m.o w.o m.o o.H 51 o.m u z ...8 533.3 2.382 08 a oaH OOH om _ m .onm cm _ oo 4 umnaaz wcHoczmm msmu0> o nonaaz u>m3 On — .o ouawwm no.0 mm.o mw.o mm.o mo.H 52 0.3 u 2 you zuwawnmum Hmuuswz u0m m hoessz mwaocxmm msmpo> o amassz m>w3 .m muswflm m-oa x a mmH mqa mma mNH mHH moH ma mm a _ _ _ _ a a a 5 73 u m -.I /' E / I 0.. - xuoa m.o m.o o.H H.H N.H 53 o.m u 2 you huHHMQMUm Hmuusoz now m Monasz mwaocmom msmuw> a nmnEDZ m>m3 .w wpswwm m-oH x a was was mas mos mma was mma mNH _ — _ _ a _ _ _ /' OH u Em Ill \L _ - I.\II m.o o.H H.H N.H m.H 54 0.0 u 2 you zuHHHnwum Hmuusmz you m ponesz mnHocmmm msmuo> o nonssz m>m3 .m muswwm nuoa x m mam mmw mNN mam mom mod mma mna q _ d _ _ a _ _ /’ I a 7.07 00/ HI m I S n 9m I: / NI I /o ./ \ l mo.a ma.H mN.H mm.H m¢.H 55 soo.o u z 001 “seem Hauauaso «:0 0m 391 + as u a mcofluoczucwwflm mo uoam .0H musmwm 0.H 0.0 0.0 0.0 ~.0 0.0 N.0n 0.0- 0.0. 0.0- 0.Hu _ _ _ _ _ _ . a _ 0.0 00.0.11 .1 N.0 H9 Illnllull «0.0-1. l.d.0 us .Illlllll \I a 0 as a 0 a 0 as 0.0 fulllllllnl IIIIIIIIIIIIIIIIIIIIIIIIIII. 1'111.\ m6 N0 0:. .l 0.H 56 No.0: H0.0- 0.0 H0.0 No.0 .5 0.m n 2 you HuoH van ouoa u m a u now ucwom Hmowuwuu msu um .afl.+ a u a mcofluocsmcowflm «0 uon .HH muswam % O 0.H 0.0 «.0 ~.o 0.0 N.0n «.0- 0.0: 0.0- o H- _ _ _ a _ _ . 0.0 57 E H 0.m m 2 you Hu0H new onoH I m pow ucHom HmoHuHuo one on .9H + 0 u e mcoHuoczwcmem 0o uon .NH oustm / .93 x m.0 mo u0uomm >0 coomvou use m>onm Ou uwHHEHm 0 0H u m now mo>p30 "00oz [I m 0 H 0.0 0.0 «.0 N.0 0.0 ~.0n 0.0- 0.0: 0.0- p h b H p . p F H 1 . H . H a . H I 0H 5m HI m.Hu 0.Hu n.01 0.0 m.0 0.H 0.N APPENDICES APPENDIX A Numerical Technique to Solve the Governing Stability_Equations A standard fourth order Runge-Kutta integration scheme was chosen in preference to a predictor-corrector scheme to simultaneously integrate the coupled second order and fourth order equations. The predictor-corrector scheme calculates the value of the highest order derivative based on a curve passed through the values at the previous three steps. All lower order derivatives are then calculated from this. Because of the large growth found to exist for some of the eigenfunctions and the suppression scheme used, variations in the highest order derivative (in particular, the fourth derivative) caused significant changes in the final eigenvalues. The Runge-Kutta method is not subject to this since each derivative is essentially calculated separately at each step. In fact, the fourth derivative at each step need not be calculated, so to save time and machine storage this derivative was eliminated from the computer program. The accuracy of the results were verified by running the non-magnetic plane Poiseuille case and comparing the results against Thomas (1958) and Reynolds and Potter (1967). Rather than write a separate program for the single fourth order equa- tion the present program was run for a Hartmann number of 10'3. This was effectively zero and gave excellent results, which agreed to three significant figures to the presently accepted 58 59 values. Several runs were made for various step sizes, with 0.01 chosen as the optiuunn Doubling this value resulted in eigenvalues that were in error by about four percent, while halving this value resulted in a change in eigenvalues of only 0.4%. Based upon these results, a step size of 0.01 provided accurate values without excessive computation time or machine storage requirements. For this step size, the truncation error 10 associated with the Runge-Kutta scheme is of the order of 10- The scheme as illustrated by Collatz (1951) is outlined below: At Y VOO=W(Y) v10 = A ¢'(y) AZ .. V20=T¢ (y) 3 v30=%-¢"'(y) U00 = 0(y) U10 = A ¢'(y) Using equations (2.4.8) and (2.4.9) 2 V _ A_. .19 2 2 . g— {[i a Rm(U-c) + K 195' - i 0’ Rmhxy - me'} n: l 60 4 2 V - L f (——V20 —10 V ——U10 U ) 1 24 A2 3 A 3 00, A 9 00’ y C) I 4 _ g: {[1 a R(U-c) + 2K2 +M2M" + 2 i QMZH y' 2 - [1 an K (U-c) + 1 ozRu"+1<4 +02M2h2 - i aMZh'N - i a M2(U-c)¢' +-[q2M2h(U-C) - i a M2(U'-h"/Rm)]¢} where A = AY = step size Values at Y + g! v =v +lv +lv +lv +£- 01 00 2 10 4 20 8 30 16 Cl 11 10 20 4 30 2 Cl 3 3 = +— + -- V21 V20 2 V30 2 G1 V31 = V3o + 2 Cl 1 1 a . — + — U01 U00 + 2 U10 4 F1 U11 g U10 + F1 As before 2 v - A. _1_ F2 2 “A ’V01’”01’Y) __ 2“ 2"21 V11 U11 6) l 2 2;qu ’—A—’V01’ A ’U01’Y) Values at Y +-AY V02 =V00 +V10 +V20 +V30 +62 v12=v10+2v20+3v30+402 v22=v20+3v30+6c2 v32=v30+4c2 U02 =U00 +U10 +F2 0 =0 +2F 12 10 2 61 Again as before F3 - 3: “1:2. ’ V02’ ”02’ Y) (; IRA4f(2V22 :1—2—,V ,E-l—g,U ,Y) 3 24 A2 ' A 02 A 02 and F = %-(F1 + 2 F2) F' = 3L0?1 +14 F2 + F3) c=%§(801+802-03) G' ='% (9 G1 + 12 02 - G3) G"=2(Gl+202) G"'= §'(Gl +-4 02 + G3) Thus the functions and their derivatives at Y + AY are W(y + Ay) V00 + V10 + V20 + V3O + G . =.l_ . 4 (y + Ay) Ay (v10 + 2 v20 +-3 v3o + G ) 2 ¢"(y + Ay) = 2 (V20 +-3 v30 + G") (Ay) w"'(y +Ay) = 6 3 030 +0'“) (by) and ¢(y + by) = ”00 +~U10 + F . =.1_ . ¢ (y + Ay) Ay (010 + F ) Now, if desired, 0""(y + AV) and ¢"(y + Ay) can be calculated from the above. Implementing the formulas in the sequence given, the integration is started at one wall and proceeds across the channel to the other wall where the combined functions are formed. APPENDIX B COMPUTER PROGRAM B-l Description of the Computer Program The program initially reads the required values for -, 100p sizes, internal program counters, step size, convergence criteria, the initial guesses for the eigenvalues, and also ~ -4 sets boundary conditions at the wall. Subroutine RUNGE is now called which in turn calls Subroutine VCA1 to calculate the value of the mean velocity and induced magnetic field and their derivatives which are known inputs to Runge. The equations are now integrated to determine the value of the function and its derivatives at the next step. Returning to the main program, the number one solution which is the growing solution is filtered from solutions two and three. To further insure independency, solution number two is filtered from number three. The above procedure is repeated until the filtered solu- tions are obtained across the entire channel. At the next stage the filtered solutions are corrected to account for the portion filtered out and then printed out if this option has been selected. Next the growing solution is modified to compress the magnitude range of the functions. 62 63 The three independent solutions are now combined linearly to find the total function and the derivatives at the opposite wall. A test function based on 1w is checked to see if the boundary condition namely 0w = 0, is satisfied. If not, the eigenvalues are each incremented by a small amount and a new guess is made for the eigenvalues. The process is repeated until convergence is achieved or the iteration counter exceeds its present value. The program incorporates a printout monitor which allows the selection of 3 levels of printout. The value of the monitor is read in and can be assigned values from one to three. A value of one provides express runs and gives the eigenvalues and values of the test function at the end of each iteration; two prints out the above plus the final combined eigenfunctions; and a value of three prints out all the eigenfunctions. When using the two or three option an additional parameter is read fin that Specifies the number of increments between printout points. Explanation of the Input Cards for the Main Program Card Program Number Column lESE Format Designation 1 1-10 Printout monitor I 10 MON 11-20 No. of data cards I 10 KP 21-30 Max. No. of iterations I 10 IT 31-40 Steps between printout pts. I 10 JJ 41-50 Step size F 10.5 DEL! 2 1-10 Percent change in cr F 10.5 PAL 11-20 Percent change in R F 10.5 PCR 21-30 Maximum change in cr or R F 10.5 PCH 31-40 Convergence criteria E 10.5 TIP Card Program Number 991393_ .lEEE Format Designation 3 1-10 Alpha F 10.5 AL 11-20 Reynolds number F 10.5 R 21-30 cr F 10.5 CR 31-40 ci F 10.5 CI 41-50 Beta F 10.5 BE 51-60 Hartmann No. F 10.5 HM 61-70 Magnetic Prandtl No. E 10.3 PM Any number of data cards as shown by card number 3 above may be added. B-2 Listing of the Computer Program The program developed to compute the eigenvalues for any given case is listed in the next section. Actually two versions of this program were used to obtain the data points used in plotting the curves presented. They differed only in regard to which of the eigenvalues was held fixed while iterat- ing on the other two. The program listed here is called MHDA which holds a constant and iterates on cr and R. This program was used on the lower leg of the stability curve up to and slightly beyond the critical point. Since the upper leg of the curve is relatively flat the second version was used which fixed cr and iterated on a and R. This procedure minimized the number of "bad guesses" which result in no convergence and also to conserve computer time by more evenly distributing the points along the curve. 65 The programs were run on the CDC 6500 and took a little over 4 secondsper pass (one third of an iteration) with con- vergence obtained in normally three to four iterations. In the program that follows, the dollar sign statement separator and multiple replacement statements are used to reduce the size of the source deck. 100 (‘00 101 598 59‘) 601 ()0? 179 178 17.8 173 66 PROGRAM HHDAI INPUT,OUTPUT,TAPE2=INPUT,TAPE3=OUTPUT) THIS PROGRAM IIERAIFS 0N CR AND RFY N0 WITH ALPHA CUNSIANT DIMFNSIUN RRIZUIQ3192I701190I3I95131911513)9AI393IoBI39319NI3). IXIBIqII3loIESII3I CHMMUN P(20103)00PI20193’907PI20193),!)39I20193)9VI20193IgnVI20193) l9C1,C29C39C49UOI3)oALofifgCflgRofllegcoKoXM CIIMPLFX P.0P.DZP,03P.V.DV.UI.UZ,UO.RR,RS.A,B,BI.BZ.P.3.T.O,X. H.011.012.CW.A1.A2.A3,A4.A5.A6.A7.A8.Cl.C2.C3.C4 RFAI I" READ IN PR'IGRAM PARAMETERS AND INITIAL GUESSES READI794001 MnNnglegJJoIH'LY RFADIZoLOII PAL.PCR.PCH,TIP N=2.0001/DFLY+I 3 NM=N-I 3 KK=0 $ C=DELY KK=KK01 $ IP=MD=KL=0 S ITF=I IFIKK-HH 600.600.250 LUNTINUE RFAF(?,40?) AL.R.CN.BE.HM.PM CF=RI§ALICHI S CI=AINAGICHI S XM=HN3€1HN WRITFI398IQ) HMQPMQCIQBEQDFLY MD=ML+1 IFIIP-O) 598,598,599 hRIIF—I3QHOI} III; II'IMIVI‘II) (50396029601 HPITFI398031 3 C“ 1n 100 HRIIFI3y807I R'ALOCR K=AI-"AI+BF*HF 5 Y="I.0 5 lIII=-I.0 1 RM=PN$R CH=CR 9.0 INUARY CIINI‘ITIUNS ----- P=DP=V=U AI Y=+"1 V=PHI=VACNFIIC pFRIHRBAIIIINSQ P=P§I=VFLIICIIY PI'RIUKHAIIHNS SFI QI'l-RIINC, CONDIIIHNQ I’II'II=PI1.71:1"1931=DPIIoII=I)PIIy?I=IIpII,3I=IO.U.U.U) Vll.l)=V(l.2)=v¢1,3)=(0.0.0.0) I‘VIIQII‘llcl"13nylot-I3IH 3‘ DVIIp7I=IIVIIq'5I=I().OgU.OI 1179(197'I=II.(1910I)) 3 “79¢Io1"“)?leo7.)‘-‘I”.UoUoOI ‘1‘5I’(112)=II.091.0) S I'3PII9I)=lJE-P(lg3)=(“.UoUoOI CI=IU.09|.0) $ C7=CI*AL”‘RM 5 C3=CI*AI’=R $ L4=ILI*AL*XM HS .300 J=IoNM $ I.I.=J CFNHIAIF 3 SOLUIIUNS AT EACH SIFP CAIL QI’NGFILlonHMgRMI L=|+I 5 ZILI=ZIJI+IIELY FXT'UICI THE GRUHINI‘, SHIUTHJN (900.1) FRHM 511le 2 ANN 4. HU=AHSH¢EALHNHIIII IFIHH-I.FI7‘)I 178,179,179 linll)=U(HlH‘IloL-IOOI H0 20” I=703 IFIUU-IoEIZ'DI 1719172917? ”RILQII=UOII)/Ufl(l) CU Tn 173 RRILoII=IUOII)Il'OII))*II.F-100) PILoII= PILgII-RPILQII’3 plLvI) IJPILvII': anLoII-RRILQII” (IPIL'II “291101I=D?PILQI)-RRILQII‘DZPILOI) 299 300 632 303 431 507 311 344 319 503 174 178 307 505 310 345 506 67 O3PIL9II=O3PIL9II-RRIchI*O3PIL91I VILoII3 VIL!II‘RRIL0II* VILOII OVILoII8 OVILoII-RRILQII‘ DVILOII CONTINUE FXIRACI SOLUTION 2 FROM 3 RRI191I=U0I3IIU0I2I PIL93I=PILv3I‘RRIL9II’PIchI 5 DPIL93I=DPIL93I-RR(L9II.DPIL92I 07PIL93I=DZPIL03I’RRILQII*OZPIL92I D3PIL93I=O3PILo3I‘RRILQII‘O3PILIZ) VIL93I=VIL93I'RRILoII*VIL92I 5 OVIL93I=DVIL93I-RRIL9II‘OVILvZ) CONTINUE ° REPAIR IHF EXTRACTIONS IFIHUN-II 50895089632 J=N 5 NN=N~2 5 RSI2)=RRIJ92I 5 RSI3I=RRIJ93I DO 303 HH=IoNN 5 J'J-I 5 DO 303 13293 PIJQII’PIJoII-RSIII’PIJOII 5 O3PIJ0II=O3PIJQII‘RSIII‘U3PIJ911 DPIJoII=OpIJoII-RSIII‘OPIJolI 5 OZPIJQII=DZPIJ9II-RSII)‘DZPIJQII VIJoII=VIJvII'RSIII‘VIJQII 5 OVIJoII=OVIJ0II‘RSIII*OVIJQII PSIII‘RSIII+RFIJvII J=N 5 RSIII=RRIJ91I 5 I=3 DU 631 HM=19NN 5 J=J‘I PIJoII=PIJ9II-RSI11*pIJoZI 5 O3PIJ911:03PIJQII‘RSIII’U3PIJ02) OPIJoII=UPIJle-RSIII*OPIJ9?I 5 DZPIJyII=DZPIJ91I‘R$III‘QZPIJoZI VIJoII=VIJ0I)-RSIII*VIJ92) 5 OVIJ9II=I)VIJ'II’KSI1I‘OVIJ021 RSIII=RSIII+RRIJ911 PRINT OUT THE INDEPFNDENI EIGENFUNCIIONS IF DESIRED GO TO (50395089507I9N0N HPITFI39807I 5 DO 319 I=Io3 5 HRIIEI30805) I 5 DO 311 J=IQN9JJ NRIIFI3,806I ZIJIgpIJvlIvOPIJoIIQOZPIJoIIoUBPIJ,II HRIIFI39816I 5 00 344 J=IO~9JJ HRIIFI30806I ZIJI'VIJQIIODVIJOII C(INIINHE CONTINUE 5 NCSN/2+I 5 J=N MODIFY IHF GROWING SOLUIION IFIHU-IoEIZSI 17501749174 A8=9IN9II‘II.E'100I 5 A5=PIN,7I/A8*II.F-IOUI 5 GU TU [76 A5=PIN92IIPINOII DO 307 KJ=19NM PIJgIIS PIJoZI-A5* PIJplI 5 091.191): IH’IJoZI"A")lflK DPIJ,“ DZPIJ,II=DZPIJ97I’A5*0?PIJ91I 5 OBPIngI=O3PIJ92I‘A5#U3PIJ9II VIJoII=VIJ92I-A5‘VIJ91I 5 OVIJoII=DVIJ92I-A5‘UVIJ9II J=J~I GO TO (5069506950510HON HRIIEI3980“) 5 I31 5 HRIIEI39805I I 5 DO 310 J=19N7JJ HRIIEI39806I ZIJIoPIngIpanJolI,O?PIJ0II903PIJ9II HRIIEI39816) 5 DO 345 J=19N9JJ HRIIEI39806I ZIJIQVIJQII90VIJ911 CONTINUE UFIERHINE THE TOTAL SOLUTION AI UPPUSIIE HALL AIIQII=VINQII 5 AI1921=VIN92I 5 AII93I=VIN93I 442.1): ”PIN'II 5 AIZoZI= OPIszI 5 AI293I= UPIN93I AI3qII=OZPIN91I 5 AI3121=OZPIN,ZI 5 AI3,3I=OZPIN.3I 692 691 694 693 305 609 610 612 613 614 525 502 503 504 615 616 611 509 511 315 68 HI=PIN911 5 82=P(N.2) 5 O3=pIN93I CALCULATE IHE COEFFICIENTS NEEDED TO COMBINE IHE FUNCTIONS Cl=A(3.l)*(Al1,213Al2.3)-A(2.2)*A(1.3)i-A(392)#(A(1,1)#A(293) I'AI29II3AI193II‘AI393I‘IAII92ItAI29II‘AII9II*AI292)I IEIUU‘IcEIZSI 69196929692 CI=CI¢IIoE-100I XIII='IAI29ZI'AII93I’AII9ZI‘AI293II/CI XI213-IAII91I‘AI293I’AI29II*AI193II/CI X131=QIAII92I*A129II‘AII9I)*AI292II/CI IEIUU'IoE1251 69396949694 XIII=XIII9I10E'IOOI 5 XI7I=XI2I*IIoE-IOOI 5 XI3I=XI3I‘II.E-IUUI CHECK IOIAL SOLUTIONS AT THE HALL HIII=XIII*AII9II+XI2I*AI1921*XI3I‘AII93’ NIZI=XIII*AI?9II*XI2I*AI292I+XI3I‘AIZ93I HI3I=XI1I*AI39II+XIZI9AI392I+XI3I’AI393I-I1.090.0I HRT=0o0 5 DO 305 I=I93 HRI=HRI+ABSIREALIHIIIII+ABSIAIHAGIHIIIII HRIIEI39HIII HRI IFIHRT'IoE-IOI 61096109609 HRIIEI39HOR) HIII9NIZI9HI3) 5 HRITEI39813I 5 GO TO 509 IP‘IP+1 COMPUIE TEST FUNCTION IIIPI=XIII‘BI+XIZI*H?+XI3I#O3 IFSIIIpI=REALIIIIpI*CONJGITIIPI)I IFITESIIIPI-IIPI 61196119612 IEIKL‘I) 61496139613 HRIIEI398I7I IIIPI9IESIIIPI WRIIEI39814I 5 KL‘O IIFRAIION OE EIGENVALUES GO TO I5029503950419IP OCR=pALECR 5 CR=CR+OCR 5 GO TO 101 CR=CR-OCR 5 RDngR*R 5 R=R+RO 5 GO I“ 10] OII=ITI2)-III)I/OCR 5 R=R-RO l)I2=ITI3I-TI III/RP OFN=AIHAGIICONJGIOIII1*OTZI OCR=AIHAGITIII‘CUNJGIOI?II/OEN RO=AIMAGIICONJGIIIIIII*OIII/OFN IEIAHSIOCRI.GE.IPCH*CRII OCR=pCHVCR*LCR/AUSIUCRI IFIABSIRUI'IPCH*RII 616,616,615 RU=PCH*R*Rl/AHSIROI R=R+RO 5 CR=CR+OCR 5 19:0 5 IIE=IIF*I 5 KL=I 5 CU I“ 101 NRIIEI39UI5) 5 HRIILI3981?) IIIPI9IFSIIIPI PRINT OUT THE FINAL COMBINED FUNCTIONS IE UESIREO IFIMON-Z) 10095119511 NRIIEI39809I 5 HRIIEI39810I 5 DO 315 J=I9N9JJ AI=XIII*PIJ9II+XIZI*PIJ92I+XI3I#pIJ93I A7=XI1199PIJ9II‘XIPI‘OPIJ9?I+XI3I‘OPIJ93) ggzxgl)sD?P(J.l)+xt2)vn2P(J.2)+X(31¢OPPIJ.3) A4=XIII*O3PIJ9II+XI21*O3PIJ9?)+XI3I*O3PIJ93) HRIIEI39820) ZIJI9AI9A79A39A4 HRITEI39BZII 5 OD 316 J=I9N9JJ A6=XIII‘VIJ9II*XI2I‘VIJ92I+XI3IEVIJ93I 69 A78X111*DV1J911+X121‘DV1J921+X(31‘DV1J93) 316 HR11E1398201 Z1J19‘69AT 510 CONYINUE 8 GD 10 100 400 FORMAT161109 F10.5| 601 F0RHAT13F10.59EIO.51 402 FORMAT16F10.59 £10.31 801 FORMATIIQHOITERATION Nn..12) 802 FORMA1115HOREYNDtDS ND. =9F20.13910X97HALPHA =.F8.4,10X, 14HCR =9F20.161 ‘ 803 FOHMAT120HOEXCFSSIVE [TERAT10N1 806 FURMAT(27HOCORRECTED GRDHING SOLUTION) 80> FORMAT (IQHOSDLUTIUN N0. 9129/91H093X91HY.10X96HP1J911918x97HDth, 111916X98H02P1J91)916X,RH03P1J911916X98H04P1J9111 806 anvattlu .F5.2.5t613.3.612.3)) 807 FORMAT124HOREPA1RED ElGENFUNCTIDNS) 808 FORMAT‘1H093HVW=9515959E15059I'7HODP51H=9E15059515059/9 18HODZPSIH=9E15o59615.51 809 FORMAT125H1F1NAL COMBINED FUNCTIONSD 810 FURNAT‘5X91HY912X94H91J1919X,5HDP¢J1918X,6HDZP(J1918X96H03P1J19 1 18X96HD¢PIJ11 811 FORNAT(27HOCUHULATIVE ERROR AT HALL 89515.51 612 FORMA111H0911HP AT HALL 892515.5910X915HTEST FUNCTION =,E15.5) 813 FORMAI(30HOFUNCTIDNS 00 NOT SATISFY B.C.1 814 FORNAT1/l1 815 FORMAT127HOCDNVERGENCE TESI SATISFIEDD 816 FORMAT11H093X91HY910X96HV1J911918X,7HDV1J911916X98HDZVIJ9111 819 FORMATQZleTHE FOLLOWING CASF IS F0R9/916H HARYHANN N0.=,F7.3. 15X918HVAGNFTIC PRANDTt 395129395X94HC1 =9F6o39SX96HBE1A 3,Fb.3, 2 I912H STEP SIZE =9F6.39///1 620 FORMAT!1X,F5.2951£12.39E12.311 821 FORMAT!1H094X91HY912X94HV1J1919X95HDV1J1918X96HDZV1J11 750 CONTINUE END SUBROUTlNE RUNGE 1J9Y9HM9RM1 THIS SUBRUU11NF USES A FOURTH ORDER RUNGt-KUTIA SCHEME D1MENS1DN G131 COMMON 9(20193190Pl2019319UZP120193).!)39120193).V(201.3).!)V1201.3) 19C19C29C39C99U01319AL9BF9CN9R9DELY9C9K9XM COMPLEX V009V109V209V309'100911109FK19GK19V019V119V219V319l101911119 1 FK29GK29V079V17.9V229V329U029U129FK39GK39FK9GK9FKP9GKP96KZpobK399 2 991199112P911399V9DV9U19U29U09CV19C19C29C39C’09619029A39A‘99A‘19A79A9 REAL K CALCULATE THE REOU1RED VELUC1TY AND HAGNET1C QUANTIT1FS CALL VCA1 (Y9HH9RM9H9HP9HPP9T9TP91pp1 Y=Y+DFLYIZOO CALL VCAL (Y9HM9RH.Z.ZP,ZPP9$9SP,SPP1 603 599 302 605 606 607 70 Y-Y+DELY/2.o CALL VCAL IY,HH,RM,U,UP,UPP,H,HP,HPPI DO 301 I=l,3 FUNCIIONS A1 Y voo=PIJ,II 9 VIOsCtoPIJ,II s V20-C¢C#DZP¢J,II/2. v30=ctc*c*DBP(J,I)/6. s UOO=VIJ,II s UlOsC‘DVIJ,II FKlsI(£20!H-CHI9KItuoo-cztttvoo-RH9VIOICItc*C/2.O GK1=((C391H-CHI9Z.9K+XHI#2.9V20/(etcI+catTtZ.tVIOIC-IC3*K91H-CH) l9c3thP+KtK+ALtlLt19ItXH~C4tTP)tvoo-CatiH-CH)9U10/C+(ALtlttxntT 2*(H-CNI-ChtIHP-TPP/RHI)9U0019C994/24. FUNCTIONS A1 VooeLYIZ. v01=voo+o.59v10+o.zstvzo+o.1259v30+o.0625*cxl v11-v10+v20+o.759v30+o.5*6K1 s v21=v20+l.59v30+l.5*GKl V31=V30+2.0#GK1 s UOlsUOO+0.5#Ulo+O.2S#FKI s UllsUlo+FKl FK2=(1C2?!z-CH)+K)#UOl-CztstVOI-Rntv1llC)*c*C/2.o GKZ=¢IC3¥Iz-CHI+2.*K+XM)*2.#VZII(th)9C6*S*2.#VII/C-(c3*K*(Z-CH) 1+CB92PP+K0K+AL9AL959stXM-catspItVOI-C491Z-CHItUlI/c+(ALtnttxnts Ztiz-CHI-cattzp—SPOIRM))*UOII-C#*4/24. FUNCTIONS AT v+nELv ‘ vnz=vooov10+v209v30+cn2 s v12-v10+2.o*v20+3.09v30+4.0*GK2 V22=V20+3.0*V30+6.096K2 s V32=V30+k.0tGKZ U02=UOO+UIO+FK2 s UIZ=UIO+2.0*FK2 FK3=t«czacu-cu)+x)*uoz-C2*H*voz-Rntv12/C)tC9C/2.0 GK3=IIC3*IU-CH)+2.*K+XMI#2.#V22/(CVCI+Ch#H#2,‘V12/C-IC3*K*IU-Cw) 1+C3*UPP+K#K+AL‘AL#H*H*XM-C4#HP1*V02-C6*IU-CHI*U12/C+IAL*AL*XH#H 2uIU-CH)-C4¢(UP-HPPIRMII*UOZ)*C**6/2¢, FK=IFK1+2.0*FK2)/3.0 s FKP=(FKI+6.0#FK2+FK3)/3.o GK:€8.*GK1+8,'GK2-GK3)/15, s GKP=I9.*GK1+12.*GK2-GK31/5. GK2P=2.*GK1+4.*GK2 s GK3P=IGK1+4.*GK2+GK31#2./3. S L=J+1 VALUES 0F FUNCTION AND ITS DERIVATIVES A1 Y+C PtL,l)=voo+v10+v20+v30+GK s DPIL,I)=(v10+2.*v20+3.tv30+GKP)lc DZPIL,II:IV20+3,'V30+GK2PI#2./IC*C)303PIL,I)IIV3O+GK3PI#6.IIC*C*CI VIL9II=UOO+UIO+FK s DVIL,I)=IUIO+FKP)/C CALCULATE SOLUTION To INVISID EON To BE USED FOR FILTERING U01!1=((U-CHI*(02PIJ+191I-K*PIJ+191)I-UPP#P(J+1911I‘AL*C1=R DFIERHINE GROWING SOLUTION IFIJ-l) 604,603,604 IFtI-B) 604,599,604 CONTINUE DO 302 H8193 AA=REALIDPIL9M1*ICDNJGIDPIL9HIII) BB=REALIPCL,M)9ICONJG(PtLonitII GIMI=SORTIAAIBBI GI=AMAXIIGIII961219GI3II IF!GI.E0.GI1)I GO TO 604 s IFIGl-GIZII 60696059606 M=2 9 GO TO 607 M23 00 608 N=J9L Al=P(N,M) s PIN,H)=PIN,II s PIN,1)=A1 AZ=DPIN,MI s DP!N.H)=DPIN,II s DP!N.1)=A2 A3=OZPIN,H1 s DZP(N.M)=DZP(N.I) s DZPIN,ID=A3 AbenspIN,MI s D3PIN,MI=03PIN,II S 03°IN,1)=A4 non 600 604 301 0 1 71 A6IVIN,NI S VIN9N18VIN911 5 VIN9118A6 A7-DVIN,HI 8 DVIN9NI=UVIN911 S DVIN9118A7 IFIN-JI 60094009608 A9=UOIHI S UOINIIUOIII S UOIIIIA9 CONTINUE ' CONTINUE CONTINUE RETURN END SUBROUTINE VCAL TV.HN,RN.U.UR,UPP,N,NR,NRPI THIS.SUBROUTINE CALCULATES THE VELOCITY AND INDUCED MAGNETIC FIELD QUANTITIES AND THEIR DERIVATIVES AT THE REQUIRED v STATIUNS COSHIXI=IEXPIXI+EXPI~XII/Z.O SINHIX)=IEXPIXI-EXPI-XI)/2.0 C=CDSHIHMI s S=SINHIHMI IFIYI 32.31.31 =-v s CV=COSHIHM8YI s SY=-SINHIHM*YI s v=-v s 60 TO 30 CY=COSHIHM*YI s SY=SINHIHN9YI CONTINUE 08HM*C-S s U=HM*IC-CYI/D s UPz-HMtHMtsvln s UPP=-HM#*3*CYIU H:RMIHM*ISY-Y*S)/(C-l.1 s HP=RMIHH*IHM#CY-SIIIC-1.1 HPP=RM$IHM*SYI/IC-l,1 RETURN END DATA CARDS 1 1 13 10 0,01 001 0.01. 002 1.5-12 015 190000. 0178 0.0 090 6.0 lee-2 "'il’iififlfifilfliflhjufififlfljflfiflifilfifliflfiflfilflmS