University of Southern California Dissertations and Theses

Convective instability of fluid interfaces

(USC Thesis Other)

Convective instability of fluid interfaces

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content INFORMATION TO USERS This manuscript has been reproduced from the microfihn master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter free, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, b^inning at the upper 1^-hand comer and continuing from to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6” x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. UMI A Bell & Howell Jnfimiiation Company 300 North Zeeb Road, Ann Arbor MI 48106-1346 USA 313/761-4700 800/521-0600 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. CONVECTIVE INSTABILITY OF FLUID INTERFACES By Frank Onice Chandler A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Aerospace Engineering) Aug;ust 1998 Copyright 1998 Frank O. Chandler R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . DMI Number: 9919021 UMI Microform 9919021 Copyright 1999, by UMI Company. All rights resered. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Ann Arbor, MI 48103 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIYERSmr PARK LOS ANGELES. CALIFORNIA 90007 This dissertation, written by Frank Onice Chandler under the direction of k .u . Dissertation Committee, and approved by all its members, has been presented to and accepted by The Graduate School, in partial fulfillment of re quirements for the degree of DOCTOR OF PHILOSOPHY Dean ifoduate Studies D ate.... DISSERTATION CO ____ R ep r o d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . ACKNOWLEDGEMENTS I would like to express my deep gratitude and appreciation to Professor Larry Redekopp for his constant encouragement, guidance and willingness to teach during the entire course of the development of this dissertation. I would also like to thank my family for their patience and endurance during the entire time the work for the research and dissertation was being performed. u R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . TABLE OF CONTENTS Acknowledgements ii Table of Contents iii List of Figures v List of Tables xi Abstract xii Part I. Rayleigh-Taylor Instability 1. Introduction 1 2. Formulation: Analysis of Layered Models 7 3. Instability Results for Layered Models 3.1 An Isolated Unstable Interface 13 3.2 An Unstable Interface Between Two Bounded Layers 17 3.2.1 Problem Definition 17 3.2.2 Two Bounded Layers with Inertial Effects 19 3.2.3 Two Bounded Layers without Inertial Effects 25 3.2.4 Effect of Geometric Constraints 30 3.2.5 Effect o f Surface Tension and Viscosity Ratio 32 3.3 The Effect of a Stabilizing Upper Interface 35 3.4 The Effect of a Stratified Layer Below an Unstable Interface 42 4. Formulation: Analysis of Continuously-Varying, Time-Dependent Density Distributions 52 4.1 Introduction 52 4.2 Model Development 52 5. Results for a Steady, Continuously-Varying Density Distribution 61 6. Results for a Time-Dependent, Continuously-Varying Density Distribution 68 7. Conclusions 72 Part n. Instability o f a Rapidly-Evaporating Liquid Interface 1. Introduction 79 iii R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 2. Formulation 86 3. Computational Results 99 4. Conclusions 110 References 112 Appendix A. Numerical Modeling of Time Dependent Model 115 IV R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . List of Figures Figure No. page Part I. Rayleigh-Taylor Instability Figure 3.1. Variation of the growth rate with respect to wave number for various values of the surface tension parameter Ttilda for an unbounded case (r2=I.0,s2=1.0). Ttilda = 0 is for a case with no surface tension. 15 Figure 3.2. Variation in maximum growth rate with surface tension parameter Ttilda (r2=l, s2=l). 16 Figure 3.3. Variation in the cutoff wave number as a function of the surface tension parameter (r2=l, s2=l ). 16 Figure 3.4. Two layer model with top and bottom boundaries. 17 Figure 3.5. Variation of Growth rate with wave number for a two layer case and the conditions of stress-free on upper and no-slip on the lower boundaries (D=20, s2=l, r2=0.9, Ttilda=0). 22 Figure 3.6. Variation of maximum growth rate for each htilda for a two layer system with a heavier fluid of height h over a lighter fluid of depth d-with stress free conditions at the top boundary (s2=l, Ttilda=0, r2=0.9TD=20). 23 Figure 3.7. Variation of wave number (beta) at max growth rate (omega max) for a two fluid layer system with a heavier on top of thickness htilda and a lighter on the bottom of thickness d-with a stress free on the top and no-slip on the bottom boundaries (s2=l, Ttilda=0, r2=0.9,D=20). 23 Figure 3.8. Variation of growth rate with wave number for two layer case, stress free conditions at the top (D=2.5,s2=l,Ttilda=0, r2=.9). 24 Figure 3.9. Variation of growth rate versus wave number for two layer case, stress free conditions on the top (D=2.5, htilda=40, s2=l, Ttilda=0, r2=0.904762, Chandrasekar check case). 24 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Figure 3.10. Variation of growth rate with wave number for two layers with r2=0.9, s2=l, Ttilda=0, D=20, No slip and stress free at top. 25 Figure 3.11. Comparison of no time derivative model for Stress free vs no slip b.c. at top ( htilda=20, d/h=2.5). 27 Figure 3.12. Comparison of models with no time derivatives and with time derivatives, both stress free on top boundary (htilda=20, dfri=2.5). 29 Figure 3.13. Variation of maximum growth rate with d/h ratio for a constant overall layer thickness for a two layer system, d+h=35. 30 Figure 3.14. Variation of maximum growth rate with htilda for a case where the depth is constant at d+h=35 and for a case where d/h=2.5=constant. The case where d/h=16.5=constant fell on top of the curve for d+h=35. 31 Figure 3.15. Variation of the maximum growth rate with htilda for various d/h ratios. 31 Figure 3.16. Variation in cutoff wave number with surface tension parameter T tilda (r2=0.9,s2=l .0, D=2.5, htilda=20, stress free model). 33 Figure 3.17. Variation of maximum growth rate with surface tension parameter T tilda (r2=0.9, s2=1.0, D=2.5, htilda=20, stress free model). 34 Figure 3.18. Variation of maximum growth rate with surface tension parameter T tilda for two different models (r2=0.9,s2=l .0, D=2.5, htilda=20). 34 Figure 3.19. Two layer model with stress free top (htilda=20, D=2.5, f=0). 35 Figure 3.20. Three layer model in an unbounded domain. 39 Figure 3.21. Variation in omega with wave number beta for a three layer case, unbounded top and bottom, using intrinsic scales (rl=0.8, r2=0.9, delta=1.0, sl=s2=1.0, htilda=20). 39 Figure 3.22. Variation of growth rate vs wave number for a 3 fluid vi R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . configuration as shown in Figure 3.20 (rl=r2=sl=s2=I, delta=0.0). 40 Figure 3.23. Variation of maximum growth rate with density parameter -A for various h . 40 Figure 3.24. Variation in maximum growth rate with htilda for various density parameters. 41 Figure 3.25. Stratified lower layer with boundaries at the top and bottom 43 Figure 3.26. Variation of growth rate with wave number for a configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=l). 46 Figure 3.27. Variation of growth rate with wave number for configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=5). 47 Figure 3.28. Variation of growth rate vs wave number for a fluid configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=10). 47 Figure 3.29. Variation of growth rate with wave number for a configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=100). 48 Figure 3.30. Variation of maximum growth rate with depth ratio D for a fluid configuration as shown in figure 3.25 (s2=l, Ttilda=0). 49 Figure 3.31. Variation of wave number at maximum growth rate with depth ratio D for a fluid configuration shown in figure 3.25 (s2=l, Ttilda=0). 49 Figure 3.32. Variation of the growth rate with beta for a two layer case with a stratified lower layer ( htil=5, m’=5, s2=l, Tbold=0). 51 Figure 3.33. Variation of growth rate with wave number for two layers, bottom layer is stratified. Stress free on top and stratified layer is bounded at bottom (htil= 10, m’=5, s2=l., Tbold=0). 51 Figure 4.1. Continuously varying, time dependent density distribution model. 54 Figure 5.1. Time dependent layer growth rate vs wave number, zero time rates. 62 VII R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . Figure 5.2. Variation of growth rate vs wave number for the static case results (dh=2.5, dl=100.l, htilda=40, r=0.1, tratio=91, 20 coefficients, 40X40). 63 Figure 5.3. Variation in maximum growth rate as a function of htilda, zero time rates. 64 Figure 5.4. Variation in the wave number where maximum growth rate was achieved versus htilda. 65 Figure 5.5. Variation of Vertical velocity and Zeta with tau time units (htilda=5, alph=.43(alph at sigmax), 60X60, dt X n^2/htilda^2<16, n=60). 65 Figure 5.6. Variation of the maximum growth rate and the growth rate at a constant wave number of 0.31 with htilda. 66 Figure 5.7. Variation in the growth rate at 50 tau for various d/h and htildas (alpha=0.31 ). 67 Figure 6.1. Variation of instantaneous sigma- growth rate with time with three different time slopes for htilda variation (dt=.0175 for stability at htilda=2). 69 Figure 6.2. Variation of growth rate with three different wave numbers for a transient case where htilda is ramped from 2 to 10 in 7.024 tau time units. 70 Figure 6.3. Variation of W/Wi with time for three different htilda variation cases. 71 Figure 7.1. Variation of maximum growth rate with upper layer thickness in relation to most unstable wave length. Two layer case as given in figure 3.7 (s2=l, Ttilda=0, r2=0.9, D=20). 73 Figure 7.2. Variation of growth rate with wave number for three different models at an htilda=l. The three layer model has rl=r2=l, sl=s2=l, Ttilda=0, delta=0. The two layer models have r2=.9, s2=l, Ttilda=0, D=20. 74 Figure 7.3. Variation of growth rate with wave number for three different models at an htilda=3. The three layer model has rl=r2=l, sl=s2=l, Ttilda=0, delta=0. The two layer models have r2=.9, s2=l, Ttilda=0, D=20. 74 viii R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Figure 7.4. Variation in the ratio of max growth rate over wave number squared to wave number. 76 Figure 7.5. Variation in maximum growth rate with htilda for three different models. The three layer model has delta ranging from .001 to 1000 with rl=r2=sl=s2=l. The two layer model has stress free conditions at the top boundary and no slip on the bottom and D=20 with r2=.9, s2=l. The bottom curve is for the cosine distribution with density parameter r=. 1. All have Ttilda=0. 78 Part n. Instability of a Rapidly-Evaporating Liquid Interface Figure 1.1. Physical Setup of Propellant Densification Process. 80 Figure 1.2. Physical Setup of Millers model for stability analysis. 81 Figure 1.3. Physical Setup of Palmer’s Base State model for stability 82 Figure 3.1. (a) Neutral curve of Hickman number as a fimction of wave number showing logarithmic behavior at large wave numbers, (b) Neutral curve showing wave number behavior at the critical Hickman number. 100 Figure 3.2. Critical Hickman number variation with Reynolds number for various density ratios (N^=100, Ncr=l0'^, Nb o =L Pr=10, Br=Ma=0). 101 Figure 3.3. Critical Hickman number variation with wave number for various density ratios (N^i=100, Ncr=10'^, Nb o =U Pr=10, Br=Ma=0). 101 Figure 3.4. Variation of the wave number where the critical Re or Hi occurred for various density ratios (Nji=100, Ncr=lO'^, Nb o =1, Pr=10, Br=Ma=0). 102 Figure 3.5. Neutral curves of Hickman number as a fimction of wave number for various density ratios illustrating the vapor recoil phase out (Nn=IOO, Ncr=10‘^, Nb o =U Pr=10, Br=Ma=0). 103 Figure 3.6. Neutral curves for Hickman number as a fimction of wave number for various density ratios illustrating vapor recoil phase out with a different N ^i (N^=20, Ncr=10'^, Nb o =U Pr=IO, Br=Ma=0). 104 IX R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Figure 3.7. Neutral curves with Hickman number versus wave number with various dynamic viscosity ratios (NRe=10'^, Ncr=10'^, Nb o =1, Pt = 10, Br=Ma=0, Np=2.06xl0^. 105 Figure 3.8. Neutral curves with Hickman number versus wave number for various Prandtl numbers (NRe=10'^, Ncr=10‘^, Nb o =U Np= 1 G O , Br=Ma=0, Np=5x 10^). 105 Figure 3.9. Neutral curves of Reynolds number versus wave number near the critical Hiclonan number (Ncr=10‘^, Pr =10, Nb o =U Np=100, Br=Ma=0). 106 Figure 3.10. Variation of Reynolds number with wave number on the neutral curve for the zero Hickman number case (Np=10*, Pr =10, Nb o = 1, Np=100, Br=Ma=0, Ncr=10'^). 107 Figure 3.11. Variation of Reynolds number with wave number for the neutral curve with Hi=0, linear wave number depiction (Np=10*, Pr =10, N bo= 1 , Np=lOO, Br=Ma=0, Ncr=10'^). 107 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . List of Tables Table No. Page Part I. Rayleigh-Taylor Instability Table 3.1 Matrix D for the two layer case with inertial effects and boundaries on the top and bottom. 21 Table 3.2. Matrix for the two layer case without inertial terms with two boundaries and stress free on the top. 28 Table 3.3. Elements of D used in the dispersion relation for the three layer case with unbounded top and bottom layers. 38 Table 3.4. Coefficient matrix of the stratified lower layer without intrinsic scaling for a stress-free condition on the upper boundary. 44 Table 3.5. Coefficient matrix for the stratified lower layer with intrinsic scaling and a stress free condition at the top. 45 Table 4.1. Matrix A for use in the continuously-varying time dependent density distribution model. 59 Table 4.2. Matrix B for use in the continuously-varying time dependent density distribution model. 60 Part n. Instability of a Rapidly-Evaporating Liquid Interface Table 2.1. Elements of the determinant for the stability analysis with heat and mass transfer. 98 XI R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Abstract The linear instability of fluid interfaces in two different physical situations has been analyzed. First, the instability of gravitationally unstable layers has been studied to assess the influence of a variety of environmental and physical property effects on expected growth rates. Secondly, the instability of a rapidly- evaporating interface has been studied to assess the role of the ‘vapor recoil’ mechanism and to identify its regime of applicability. The linear Rayleigh-Taylor instability analysis for several different models has been performed. The effect due to an upper boundary adjacent to and above a heavier upper layer of varying thickness has been determined for either stress free or no slip boundary conditions. The effect of the thickness of the heavier upper layer on the growth rate has been evaluated and the advantage of using intrinsic scales has been clarified. A three-layer system of a stabilizing, lighter-density layer above a heavy core layer, both overlaying an unbounded layer of medium density has assessed the effect of a stable interface on the upper surface of a heavy layer. In addition, the stabilizing effect of a stable, density stratification in the fluid layer underneath a heavier upper layer, interfacial surface tension, and viscosity ratio has been analyzed. The effect of using a creeping flow assumption in the stability formulation reveals that the neglect of the inertial terms can lead to substantial errors in predicting the growth rate. A model was also developed to assess the effect of a continuous, time dependent density distribution. Distributed density profiles yield significantly reduced growth rates compared to sharp interfaces. Linear hydrodynamic stability analysis was performed for a liquid undergoing steady evaporation at very low pressures. Specific models of interest included configurations where cryogenic liquids could undergo an evaporative xii R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . process or a propellant densification process. Results were developed for liquids undergoing the ‘vapor recoil’ effect and applied to the configurations of interest. The results indicated that the configurations of interest were not in the regime for ‘vapor recoil’ to be operative, and that ‘vapor recoil’ phase out had occurred prior to nearing the regime of interest. Xlll R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . EartJ, Rayleigh-Taylor Instability 1. Introduction The study of the motion of superposed, statically unstable liquid layers, specifically Rayleigh-Taylor instability, is important from a number of physical applications in natural processes and industry. Large scale natural processes include those associated with the overturn of the outer portion o f a collapsed core star, or the formation of high luminosity twin-exhaust jets in rotating gas clouds in an external gravitational potential, or the sedimentation of heavy particulate laden liquid intrusions such as in river outflows to the ocean, and sinking slabs of tectonic plates over layers of hot mantle core material. Industrial processes where Rayleigh-Taylor instability is relevant include Laser implosion of deuterium- tritium fusion targets, coating processes such as painting, gelatin emulsions and magnetic suspensions producing photographic films as well as magnetic media, and the internal propellant dynamics o f space vehicles. Clearly the applications for liquid layer research is many and varied. Rayleigh-Taylor instability of liquid layers has received renewed interest for both the natural phenomena and industrial applications mentioned. The understanding of the parameters that effect the instability behavior of natural phenomena or that can control industrial processes so that production efficiency is optimized are important reasons to continue the study of this particular mode of instability. Rayleigh-Taylor instability involves the situation where we have layers of two or more fluids undergoing acceleration perpendicular to the interface with the lighter fluid accelerating the heavier fluid. The fluids can be unbounded or bounded with no-sIip or stress free boundaries. This physical situation sets up Rayleigh-Taylor instability, whereby the heavier fluid ‘falls’ through the lighter 1 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . fluid in columns or fingers. The specific morphology of these columns under special conditions has been described by various authors utilizing nonlinear analysis. These columns can have significant entrainment of ambient fluid after the vorticity has had time to effect the falling column. The earlier works on Rayleigh-Taylor instability involve linear analysis of the interface and start with the work of Rayleigh (1900), who first discovered the phenomenon, and Taylor (1950) who dealt with arbitrary accelerations. Chandrasekhar (1961) extended the analysis by describing the effects due to viscosity and surface tension for two layers of unbounded fluids. Chandrasekhar also provided a mathematical proof, for unstable configurations of the Rayleigh- Taylor problem in the absence of viscosity and smface tension, that there is only one real positive root for the growth rate. Sharp (1984) provides a good overview of the basic phenomenon with a discussion of nonlinear analysis and numerical computations. Menikoff et al. (1977) performed a linear analysis and concluded, in consistency with Chandrasekhar’s proof, that the dispersion relation would admit only one real positive root for the growth rate for unstable fluid configurations with stuface tension and viscosity. They also developed approximations for the variation of growth rate with respect to wave niunber for parametric ranges for viscosities and densities. The effect on the maximum growth rate with large exciu’ sions in surface tension were also investigated and showed a typical s-curve fi"om the upper growth rate without surface tension asymptoting to zero as the surface tension went to infinity. Menikoff et al. (1978) developed the analysis further to lift restrictions on the density or viscosity of either fluid. Jain and Ruckenstein (1976) performed a linear stability analysis of stagnant films under a viscous fluid on a solid surface. They neglected the inertial 2 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . terms in the development of the governing equations and they limited their analysis to film thicknesses smaller than the distance over which the van der Waal intermolecular forces are effective. Batchelor and Nitsche (1991) performed a linear stability analysis of unbounded fluids with a region containing a continuously varying central or core density for several density variations profiles in both stably stratified density distributions and neutral density gradients outside the core layer. These core layer profiles include, for stably stratifîed density distributions, a cubic density profile, a piecewise-linear density profile, and a sinusoidal undisturbed density distribution. For the cases with vanishing density gradients outside the core layer, the profiles considered included a layer of transition from light fluid below to heavy fluid above, an isolated layer of heavy or light fluid, and a heavy and a light sub-layer, one on top of the other. In each case the layer thickness is taken as the length scale in between two semi-infinite regions where the density gradient tends rapidly to zero. They performed asymptotic long wave analysis of these cases (shrinking the central layer to an interface with zero thickness), using similar profiles represented by Gaussian functions and similar profiles represented by piecewise constant density profiles. This allowed a comparison of the effect of layer thickness as well as layer profile smoothness on the growth rate. Their results indicated that an appreciable thickness of the central layer has a stabilizing effect and discontinuities in the density profile have a destabilizing effect. The growth rate for a Gaussian form of a sinusoidal density profile of a thick layer is approximately 1/3 that of a sharp interface and the growth rate for a smooth layer is approximately 70% that of a piecewise continuous layer. Batchelor and Nitsche (1993) reanalyzed the stability of a stationary stratified fluid in a vertical cylinder to various density profiles similar to those of the earlier work and obtained results 3 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . that indicated that the vertical walls of the cylinder dominate some of the instability results. That is, the finite aspect ratio of the domain had a stabilizing effect. In order to develop the morphology of the interface caused by Rayleigh- Taylor instability, a nonlinear analysis must be performed. Several works have developed various estimates of the form of the developing layer intrusions, or fingers, and an initial overview is again described by Sharp (1984). An experimental investigation of the mixing zone between two layers undergoing Rayleigh-Taylor instability is presented in Emmons et al. (1960) and in Read (1984) with photographic depiction of the developing columns. Babchin et al. (1983) considered a weakly nonlinear analysis in a system of a thin film flowing under a heavier film of equal viscosity. They neglected the inertia terms in the model development. The nonlinear results show that, for certain values of the base velocity (shear rates), the instability undergoes a nonlinear saturation that prevents the film from rupturing and has enhanced stability over stagnant films. Yiantsios and Higgins (1989) used creeping flow assumptions and developed a weakly nonlinear analysis producing developing morphology. They used a model of a thin film underneath another fluid with different viscosities but without a base velocity. The analysis shows that the instability leads to a number of steady-state interfacial shapes with accompanying film rupture and pendant shaped drops of the light fluid in the film rising into the upper fluid. Newhouse and Pozrikidis (1990) consider the nonlinear, viscous Rayleigh- Taylor instability of a liquid layer resting on a plane wall below a higher density fluid. Using creeping flow assumptions, the evolution of the interface is developed by means of numerical simulation as a function of surface tension and 4 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . viscosities of the two fluids. The results show that for small or moderate surface tension the rising fluid takes the form of rising plumes. The morphology of the rising fluid, either as a drop with a trailing small column or as a rising column, depends on the viscosity ratio of the two fluids. Small surface tension has the effect of allowing vorticity (generated by the interfacial shear) to cause the rising drop to entrain the heavier fluid underneath the trailing edge of the rising drop. The terminal phase of the column development is not realistic as stated in the results due to the ill-posed nature of neglecting the time derivatives and inertial terms in the derivation. The rate of thinning of the film as well as the speed of the rising drops are also developed. They indicate that the growth rate, as predicted by linear theory, decreases as the viscosity of the upper layer is increased. Also, they obtain the result that the instability development remains essentially linear (same as the linear theory results) up to the point where the perturbation has grown to approximately 20% of the original film thickness or length-scale. Clearly, a number of questions regarding Rayleigh-Taylor instability in various contexts remain unanswered. In many of the studies to date, the focus was on unboimded fluids or on the effect of instability to thin films confined by a single horizontal boundary. The effect of the stabilizing influence of confining boundaries (measured in terms of the length scale) has not been systematically investigated. Also, the stabilizing effect o f multiple boundaries with various boundary conditions (stress free versus no-slip) has not been explored. In a number of the studies to date, the creeping flow assumption was taken to simplify the analysis. The effect o f this assumption on the stability results has not been carefully assessed. The actual effect of surface tension and various layer viscosities on the temporal perturbation growth rate needs to be assessed without restriction o f the creeping flow assumption. The question of whether a central 5 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . core layer of heavy fluid having finite thickness positioned within a stably stratified density profile is also of great interest. Particle-laden upper layer intrusions can, during settling, exhibit an unstable smooth density profile at the interface between the heavier upper layer intrusion and the lower lighter layer as detailed in Maxworthy (1998). The settling particles produce a time dependent density distribution that exhibits a Rayleigh-Taylor instability mechanism as shown in the reference. In order to assist characterization o f this phenomena, it would be very interesting to evaluate the effect of a smoothly varying density profile on top of a constant density fluid that can have temporally developing profiles. The analysis undertaken in this study seeks to provide answers to each of these areas in question. Results are presented in the various sections that follow for different models designed to address some of these questions in isolation. R ep r o d u ced w ith p erm issio n o f th e cop yrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 2. Formulation: Analysis of Layered Models Several multi-layer models are analyzed to elucidate the linear instability characteristics associated with a heavy layer of fluid superposed above a lighter layer o f fluid. For this purpose we consider the planar motion of a viscous, incompressible, non-diflusive fluid. Since we propose to investigate a model incorporating the effect of a stably-stratified fluid positioned below the heavy layer, the linear equations of motion are taken in the form, u^+w. = 0, (2.1) P, + ^P's = 0, (2.2) p,K, = -p^ + (2.3) p^w, = -p^ - ^ (2.4) The coordinate x points in the horizontal direction and z points vertically upward, anti-parallel to the gravitational field. In this form of the equations of motion, the hydrostatic state relating to base-state density Pj(z) and associated pressure distribution p,(z) given by ^ = - p .g (2-5) az has been subtracted finm the vertical momentum equation (2.4). The perturbation density and pressure relative to the hydrostatic state are represented by and p(x,z,t) , respectively. In the present formulation we specifically do not impose the Boussinesq approximation and we assume that the coefficient of kinematic viscosity v = p, I p, is constant in any layer. It is convenient to eliminate dependent variables in favor of the vertical velocity w(x,z,t) from the set of equations (2.1-4). Eliminating the perturbation pressure p and employing the conservation of mass constraint (2.1), we obtain the equation. R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . V^w, + — W j, = +vV^V^w (2.6) P. P, In the special case where the base-state density is a uniform constant throughout any given layer, equation (2.6) reduces to jv'vi; = 0. (2.7) If the layer possesses a base-state density stratification, equation (2.2) can be used to eliminate p from (2.6) and yield the governing equation w, + ^ = 0, (2.8) where the function N*(z) is the square of the Brunt-Vaisala frequency N \ z ) = - g ^ . (2.9) P. The two equations (2.7) and (2.8) define the governing equations for the motion in any given layer depending on the absence or presence of an underlying density stratification. For the purposes of a linear instability analysis, we look for normal mode solutions of the form w{x,z,t) = W{z)e°’' ^ . (2.10) A Fourier mode disturbance with wave number k (of infinitesimal amplitude in this linear analysis) is imposed and the temporal growth rate a and vertical eigenfunction W(z) are sought. For purposes of simplicity and analytical convenience, the base-state density stratification is chosen to have the form /7 /z) N-=g/3. (2.11) For a stably stratified state, P>0 and the Brunt-Vaisala frequency N is real. With this choice for the density stratification, equation (2.8) assumes the form £ > 2 + + = (2.12) R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . where we have introduced the short hand notation D = — . Holding a and v dz fixed and letting ^ X 0 (i.e., considering the limit of vanishing density stratification), equation (2.12) reduces to [D^ -k^)W = 0. (2.13) The parameter K defined as = A ;" + - (2.14) V is introduced here for convenience and later reference. The analysis of specific flow models which follows requires the solution of the fourth-order equations (2.12) and (2.13). With the specific choice (2.11) of the density distribution, these equations have constant coefficients and, therefore, admit the solutions PK(z) % g" (2.15) where, for the limit case represented by (2.13), the four linearly independent solutions are represented by the four values X = ±K, ± k . (2.16) In the more general case represented by (2.12), the four solutions for Â . are given by the quartic equation i X} = (2.17) ' V (JV In the limit where the Boussinesq approximation is invoked for the motion interior to a given layer, the term with |3 appearing explicitly is ignored and characteristic equation reduces to a bi-quadratic with the solutions À = ±' k'+K^ a 2 ~2 V 1 - 4 ^ a (2.18) In general, and certainly for amplified disturbances with a>0, the solutions for X are complex corresponding to propagating internal waves. As is well known fi"om 9 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . internal wave theory, the waves are propagating vertically upward for Im X<0 and vertically downward for Im X>0. A special limit case of the foregoing development deserves comment. Some investigators (see, for example, Newhouse and Pozrikidis (1990)) have ignored inertial effects in their linear analysis of an unstable interface above a shallow layer of fluid. This amounts to ignoring the effect of the time derivative in equation (2.7) and the term involving a in (2.14). As a point of explanation, neglect of the inertial terms in the momentum equations does not automatically imply that the flow is stationary. Unsteady terms still enter through the kinematic matching conditions at a fluid interface. In this limit the solutions for the X involve repeated roots (c.f, equation 2.16). This approximation can be expected to yield reasonably approximate results provided k^»cr/v, which holds when the instability is dominated by moderate wave numbers and small growth rates. However, this approximation will fail when the peak growth rate is shifted toward the long wave end of the spectrum. This approximation will be explored further in a later section in the context of a specific flow model. In the analysis of layered flow models, interface matching conditions are required to connect the solutions applicable in neighboring layers. In the context of the linear instability analysis of a viscous fluid being considered here, the requisite matching conditions are the following. In writing these conditions we assume that the interface is located at the position z=h+r|(x,t), where r|(x,t) denotes the displacement of the interface from its equilibrium position z=h. Continuity of normal velocity: [H ] = 0. (2.19) Continuity of tangential velocity: [H] = 0. (2.20) 10 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Continuity o f tangential stress: [ k “, + '^ J ] = 0- (2 21) Continuity of normal stress: {p. - P* )gri - [ |- p + 2 /av , I] - = 0 . (2.22) The interface distortion ii(x,t) is related to the vertical velocity through the kinematic condition w = 7,. (2.23) In writing the above forms we have used the notation [|(-")i] = ("-)L ^ ♦ (2.24) and p* and p_ denote, respectively, the fluid density immediately above and below the interface. The parameter T in (2.22) is the coefficient of surface tension (in units of force/length). To apply these conditions we need to represent them in terms of normal modes as in (2.10). To this end we define the modal amplitudes ={U(z),P{z),NYe°'*‘ ^. (2.25) Then, using the equation for conservation of mass (2.1 ) and the horizontal momentum equation (2.3), the respective modal amplitudes can be related to the eigenfunction W(z) as follows: U{z) = - D W , (2.26) k P{z) = -K ^)D W . (2.27) The amplitude of the interfacial displacement is related, through the kinematic condition (2.23), to the eigenfunction W by the relation N = -W {h). (2.28) < j The interfacial matching conditions can then be rewritten as follows: 11 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . . _ , (2.29) {(P --P .)g-^ + -2k^)DW^ = 0. The remaining conditions needed to completely specify the solution for a given flow are those at a solid or free surface. In what follows we consider two different cases. First, at an impermeable, no-slip surface we have the conditions W=DW = 0. (2.30) Secondly, at an impermeable, stress-free surface we have the conditions W=D-W = 0. (2.31) The influence of these different boundary conditions on the growth rate of disturbances on a statically unstable interface in proximity to an impermeable surface will be investigated. 12 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3. Instability Results for Different Layer Models. 3.1 An Isolated, Unstable Interface As an initial example we consider the simple case of a single density interface in an unbounded domain. The analysis of this configuration is well- known (see, for example, Chandrasekhar (1961), pp 441-448), but we summarize briefly some results to reveal some special limits for later models and to elucidate some scaling issues. Consider a single plane interface at z=0 separating immiscible fluids of density pc and viscosity pc in the region z>0 from that with density p2 and viscosity pz in the region z<0. The solutions for the vertical velocity eigenfunction W(z) satisfying the boundary conditions at z -> ±oo are: fV^(z) = Ae'*^ + Be~^^, 0<z<oo; (3.1) (z) = C e** + De^'~, - o o < z < 0. (3.2) The designation Kc and Kz denotes the expression for K given in (2.14) evaluated using the kinematic viscosity in the respective layer noted by the subscript. Then, applying the matching conditions (2.29) yields the dispersion or eigenvalue relation 1 Ji Zk' - I k M e - I ^2 ^2 ( ,1 M e =0. (3.3) Tfiç{,k^ + K^) -IfdçOkK^ In this problem there are no extrinsic length or time scales which are naturally available for forming relevant non-dimensional variables. Hence, we form the intrinsic scales 13 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . I s = g (3.4) where g = g{p^ ~ P i)! Pc the reduced gravity for the interface. The reader will note that we have specifically assumed Pc>P2 so that the interface is gravitationally unstable and g >0. Using these intrinsic scales, the following dimensionless variables are defined: = '• 2 = — . Pc Pc T = (3.5) The parameters o f the problem are the density ratio t2, the viscosity ratio si, and the inverse Bond number T based on the intrinsic length scale /c Chandrasekhar (1961) defines the intrinsic scales in terms of g, rather than the reduced gravity g . In the limit of the Boussinesq approximation rz=l and, if we ignore any viscosity stratification, S 2=l also. In this special case p2=Pc and the only parameter of the problem is the Bond number T . Using the intrinsic scaling defined by (3.5), the eigenvalue relation (3.3) takes the form -1 P This relation admits several temporal modes (i.e., several roots for o) at a fixed value of P and at fixed fluid parameters), but we are interested only in the largest positive root which represents the growth rate of the most unstable mode. 14 1 P P" 1 Pc -1 Pi P^- + Û 3 =0. (3.6) R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . Instability occurs by means o f a stationary mode (i.e., c o is real) which, as discussed by Chandrasekhar (1961), is only marginally stabilized as T tends to infinity. That is, the instability is restricted to a band of wave numbers 0<P<|3n , where is the cut-off wavenumber, and Pn decreases toward zero as T increases. These results are exhibited in Figure 3.1 where the temporal growth rate c o is shown as a function of the wave number P for selected values o f the Bond number T . The results shown in Figure 3.1 correspond to the Boussinesq limit rz=I with S 2= l. The peak value of the growth rate and cutoff wave number as a function of T are shown in Figure 3.2 and 3.3, respectively. The results shown in these figures provide limiting values against which results for other flow models studied later can be compared. 0.3 0.2 Q ) da=0 Tti I Ô 0.15 0.1 0.1 0.3 0.5 0.05 2 1.8 1.4 1.6 0.2 0.4 0.6 0.8 1 0 Wave number Figure 3.1 Variation of the growth rate with respect to wave number for various values of the surface tension parameter Ttilda for an unbounded case (r2=l .0,s2=l .0). Ttilda = 0 is for a case with no surface tension. 15 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.3 Maximi m growth rate s t Tlilda=0 is 0.2 )25 0 2 5 8 I 0 ^ o I 0.15 I I 1 5 0.05 10000 0.1 10 1000 100 1 Ttilda, sur^ce tension parameter Figure 3.2 Variation in maximum growth rate with surface tension parameter Ttilda (r2=l, s2=l). 1.8 1.6 I 12 c fc 0.8 I ^ 0.6 0.4 02 10000 1000 10 100 0.1 1 Surface Tension parameter, Ttilda Figure 3.3 Variation in the cutoff wave number as a function of the surface tension parameter (r2=l, s2=l). 16 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3.2 An Unstable Interface Between Two Bounded Layers. 3.2.1. Problem Definition The first physical or environmental effect we investigate is that of finite thickness of the layer of heavy fluid overlying a lighter fluid. We also consider the effect of finite depth of a lighter layer positioned below a heavy layer. To this end we consider the instability of the flow model shown in Figure 3.4. The upper layer of density pc and viscosity pc has depth h. The lower layer of density pz and viscosity pz has depth d. We choose the level of the interface between the two fluids of differing density and viscosity to be at z = 0, and the solutions for the eigenfunction W(z) in the respective layers as ff'.Cz) = + Qe*^ + Re-'^^+Se^‘ ^, 0<z<h, (3.7) (z) = A e'^ +Be^ + , - d < z < 0. (3.8) The eight independent constants appearing in these expressions require the application of eight independent conditions. These include the four matching conditions given in (2.29) and applied at z = 0 and the boundary conditions, either (2.30) or (2.31), applied at the two boundaries located at z = h and z = -d. pc P2 z= h stress free or No slip (Wzz = 0) (U=W = 0) impermeable and no slip (u = w =0) rZ = -d Figure 3.4. Two layer model with top and bottom boundaries 17 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . Since we also use this model to explore the implications of neglecting the inertial terms, an approximation used by Nev&ouse and Pozrikidis (1990), we need to specify the eigenfunctions W(z) ^propriate to this lim it As already noted, the characteristic equation corresponding to (2.13) possesses repeated roots in this case and the eigenfunctions must take the alternate forms lP,(z) = Pe~^ + Qe^ + Rze-*^ + Sze*^, 0 < z < h, (3.9) W-^{z) = Ae~'^ + Be'° + Cze'’ ^ D z e ' ^ , - d < z < 0 . (3.10) Before proceeding to the full dispersion relation for different boundary conditions and different approximations, we consider the normal stress matching condition (see equation 2.29) and use the eigenfunctions (3.7) and (3.8). Substituting and carrying forth some algebraic simplification leads to the equation j M l^ ( k h ) — ^ { k h Y + 2 ^ — (AA)(^,/z)lc (3.11) L Pe^ c Pc ^ c J I C P c C P c ^ c J - 2 ) + 2 (k h \K ,h lR - S)} = 0. The matching condition was written in this form to expose alternate non- dimensionalizations that can be used. Equation (3.11) is written in a form where the thickness h of the heavy layer is used as the reference length scale. The natural time scale in this case is the diffusive time X d = hVvc. With these choices the Rayleigh number = = (3.12) Vc PrVc appears naturally, together with the Crispation number 18 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . Cr = ^ ^ . (3.13) Pc^c Use of this set of length and time scales (h, td ) is certainly reasonable and may, under some circumstances, be preferable. However, one of the principle goals here is to assess the role of the thickness of the heavy layer on Ray leigh-Taylor instability and to clarify when the growth rate of the linear instability asymptotes to the result for an isolated interface in an unbounded domain (c.f., section 3.1). For this purpose we choose to employ rather the intrinsic scales defined in (3.4). This has the advantage of introducing explicity the parameter h = - = R j , (3.14) which immediately yields the layer thickness relative to the intrinsic scale U that characterizes the limit problem when (h,d) -> oo. Of course, as noted in (3.14), these alternate choices for reference scales are not unrelated in terms of the basic parameters characterizing the problem. In fact, the ratio of the viscous diffusion time scale Zd to the intrinsic time scale X c defined in (3.4) is also related to the Rayleigh number: (3.15) ( v , / r ) = Hence, conversion between these two alternate scalings depends directly on the value of the Rayleigh number. Alternately, the Rayleigh number can be viewed in terms of a ratio of length scales or a ratio of time scales. 3.2.2. Two Bounded Lavers with Inertial Effects Considering the case where the boundary conditions on the upper surface at z = h are impermeable and stress free and those on the lower surface at z = -d are impermeable and no-slip, and including inertial effects so that the 19 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . eigenfunctions are (3.7) and (3.8), the dispersion relation is given as det D = |D| = 0, and matrix D is the 8x8 matrix given in table 3.1. In the two layer model, stress free conditions on the top boundary were selected in order to more closely simulate a free surface on the top, and no slip on the bottom boundary was chosen (at a variable distance d from the interface) to be able to show the influence of a bottom boundary. This boundary condition choice should allow the results to be applicable to intrusions where a heavy fluid, or a particle-laden fluid, overlays a lighter layer of fluid. The results for the two-layer model with boundary conditions as stress free on the top and no slip on the bottom is shown in figure 3.5. The general behavior, for a constant h , is that the growth rate, ©, rises from a wave number p of zero up to a maximum and then asymptotically approaches zero as beta goes to infinity. The value of T is zero unless specifically mentioned. The maximum growth rate, ©mw, for each constant h curve, occurs at larger wave number for small h and rapidly decreases to the value P = 0.4 for larger h . The maximum growth rales and the wave numbers where they occur began to coalesce after an h o f approximately 6. The variation of the peak growth rate © m ax as a function of the intrinsic depth h is shown in figure 3.6 for a particular value of the non- dimensional depth D= d /h =20 ( a value sufBciently large that the constraining effect o f the bottom boundary is essentially zero). The wave number corresponding to the maximum growth rate © m ax is given in figure 3.7. Clearly for h >10, the peak (maximum) growth rate becomes independent of h . The variation in the stability behavior at another depth ratio (D= d /h ) is shown in figure 3.8 for a D=2.5. The effect of decreasing depth is evident at lower wave numbers as expected, but has a negligible effect on the peak growth rate which occurs at moderate wave numbers relative to a depth ratio D=2.5. 20 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . '■ 5 . ca. o . o % S. ea. o ' â o ^ o 7 5 a . 5 * '*' j= ») «N ca. I % , o 3 ? o ^ I ê r r s ^ ca. o ^ - 4 - ca. ■3 = f = ' » . - ? ■ ~ t â Ç ' ■ c 3 J. 3 _ % < N i = \ ° â" f t & ' I « - i r J î — " 3 . "Sa. ca. + ^ ® \ ® ca. "• '=- cT ”t£ ! J ü I CN I 'ZZ- g % 21 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . A comparison was made with the results from Chandrasekhar (1961). Chandrasekhar used g (instead of the reduced g ’) in his scaling. A check case was run using the two layer model, with a stress free boundary condition at the top and intrinsic scales with D=2.5, h =40, s%=l, T =0, and r2=0.904762. The check case is shown in figure 3.9. Converting the present results to Chandrasekhar’s scaling we obtain @=.059958 and k=.l 87235 versus Chandrasekhar’s values of tO c= 06086 and kc=.1939. This is in reasonably good agreement, and the differences lie in the resolution capability of the eigenvalue search routine where A@=.0005 and Ak=.02. Slight differences may also arise from the finite depth effect (D=2.5 here) since Chandrasekhar’s results were for a configuration with a single interface in an unbounded domain. 0.3 2 0.15 2 O 0.1 0.05 1 2 1.4 1.6 0 0 2 0.4 0.6 0.8 1 htilda 1 3 5 ....... 7 - C- 9 — o. 11 Wa\/e number Figure 3.5. Variation of Growth rate with wave number for a two layer case and the conditions of stress-firee on upper and no-slip on the lower boundaries (D=20, s2=l, r2=0.9, Ttilda=0). 22 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0. 3 025 ! " & 0.15 I I o - ' I 0.05 12 10 8 6 0 2 4 h tilda-nondimensional liquid height of top layer Figure 3.6. Variation of maximum growth rate for each htilda for a two layer system with a heavier fluid of height h over a lighter fluid of depth d-with stress free conditions at the top boundary (s2=l, Ttilda=0, r2=0.9d>=20). 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 02 0.1 0 12 6 8 10 0 2 4 h tilda-nondimensional height of top layer Figures.?. Variation of wave number (beta) at max growth rate (omega max) for a two fluid layer system with a heavier on top of thickness htilda and a lighter on the bottom of thickness d-with a stress free on the topand no-slip on the bottom boundaries(s2=l, Ttilda=0, r2=0.9J)=20). 23 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0 . 3 02 G B 15 0.15 0.05 1.6 02 0.6 0.8 1.4 0.4 1 0 •nondimensional wave number Figure 3.8. Variation of growth rate with wave number for two layer case, stress free conditions at the top (D=2.5,s2=l,Ttilda=0, r2=.9). 0.3 02 ts 0.15 0.1 0.05 12 1.4 1.6 0.4 0.6 0.8 1 0 02 Wave number Figure 3.9. Variation of growth rate versus wave number for two layer case, stress free conditions on the top (D=2.5, htilda=40, s2=l, Ttilda=0, r2=0.904762, Chandrasekar check case). 24 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . The two layer model o f Figure 3.4 (with inertial terms) was compared with respect to different boundary conditions on the upper surface. Results are shown for no slip as well as stress free boundary conditions on the upper surface for h ’s o f 1,3 and 5. It can be seen that the no slip boundary condition has a slightly stabilizing effect relative to the stress free condition but that the effect is most pronounced at the lower h of 1 and becomes less influential at the higher h ’s. Boundary Condition/ h SF1 SF3 SF5 ....... ^31 ..Æ... nS 3 NS5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Wave number Figure 3.10. Variation of growth rate with wave number for two layers with r2=0.9, s2=l, Ttilda=0,0=20, No slip and stress free at top. 3.2.3. Two Bounded Lavers without Inertial Effects Several authors have neglected the inertial terms in the momentum equations in their derivation of the equations defining the stability of a two layer system (Babchin, et al., (1983), Newhouse and Pozrikidis (1990), Yiantsios and Higgins (1989) and Jain and Ruckenstein (1976)). This has been done by utilization of the creeping flow assumption (Stoke’s flow) and using lubrication 25 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . theory to justify neglecting the inertia terms. This causes significant simplification in the model, but the validity of the approximation has not been assessed. It was desired to determine if this assumption would cause a divergence in the solution and if so where the divergence occurred. When the inertial terms are neglected, we have to use the alternate eigenfunctions as given in (3.9-10). The matrix given in Table 3.1 must now be modified and results in the form given in Table 3.2. The two layer model (with no time derivatives) provides results for the growth rate, c c > , versus wave num ber, P, as given in figure 3.11. Results are shown for a no slip upper boundary as well as stress firee upper boundary. The general behavior of © vs p is consistent with the earlier two-layer model. However, the values of © m ax and the wave number p where © m ax occur are significantly changed. It can be seen that the no slip boundary condition has a slightly stabilizing effect relative to a stress fiee boundary by yielding a somewhat smaller value of © m ax (1.92 vs 2.35). The wave number where the peak growth rate occurs is also slightly different for the two different cases (.085 vs .07). The difference in the behavior of the growth rate versus the wave number for the two layer models with and without time derivatives is shown in Figure 3.12 for a case with stress firee top and no slip bottom boundaries and h =20, D=20. It is clear that the time derivative terms must be included to obtain accurate results at the lower wave numbers. It can be seen that the no time derivative model © diverged firom the time derivative model © by greater than 10% below a wave number of 1. The differences are expected at lower wave numbers by reference to equation (2.14) [compare K to k]. Since the peak growth rates for the case including inertial effects are bounded by a number less than one, as a becomes small the difference between K and k can be substantial. Therefore, 26 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . for the Rayleigh-Taylor instability mechanism, these results bring into question any use of a creeping flow or lubrication theory assumption. I 2 O 2.5 SF NS 2 ■ 'I , 1.5 1 0.5 0 0.25 0.1 0.15 0.05 0 W aw number Figure 3.11. Comparison of no time derivative model for Stress fiee vs no slip b.c. at top ( htilda=20, d/h=2.5). 27 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . f I I ^ = T - U CO. + (S § c n < I ^ ST -o ^ 0 ^ 0 0 _ I O (N I J3 i i \ ® \ ® T 7 7 t ; cn CÛ- 7 ■ c s ë % g o o ins o O CO. C O . I a "3 . 0 "S ^ ^ O O O — “ p * o -2 ? + 1 ■=' I çO o ' ^ o ^ — — I “ “ il" S + ri _ C" 3 ® \ ® ^ 7 ^ « I - g 28 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 2.5 NT/SF « > S 2 O 0.5 0- 0.8 0.9 1 0.7 0.6 0.5 0.2 0.3 0.4 0.1 0 W aw number Figure 3.12. Comparison of models with no time derivatives and with time derivatives, both stress free on top boundary (htilda=20, d/h=2.5). 29 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3.2.4. Effect of Geometric Constraints All of the above results produced eigenvalues for cases with a constant h or a constant D. It is desired to know for a case with a constant depth how the maximum growth rate changes with a changing D ratio. A case is shown in figure 3.13 for a constant depth {d + h ) of 35. These results will be used to compare with the case of a continuously varying, time dependent density distribution. It can be seen that the growth rate decreases slightly as the D ratio increases. The variation of the maximum growth rate with h for a constant D and depth of 35 is shown in figure 3.14. It is seen that the growth rate is mostly dependent upon h , and if a secondary variable such as D is changing for a constant depth, the driver for changing the growth rate is h . The variation in maximum growth rate with h for various ratios of D is given in figure 3.15. It can be seen that for h >4, the variation in © m ax with D (for D>2) has asymptoted and the variation in © m ax is entirely dependent upon h . 0.3 P 02 o > 0.15 0.1 0.05 6 0 8 12 10 14 16 2 4 18 d/h ratio Figure 3.13. Variation o f maximum growth rate with d/h ratio for a constant overall layer thickness for a two layer system, d+h=35. 30 R ep ro d u ced with p erm issio n o f th e cop yright ow ner. Further reproduction prohibited w ithout p erm issio n . 0.3 I m 0.15 E a CO 0.1 sigmax 0.05 6 5 7 8 9 10 3 4 0 1 2 htilda Figure 3.14. Variation of maximum growth rate with htilda for a case where the depth is constant at d+h=35 and for a case where d/h=2.5=constant. The case where d/h=l 6.5=constant fell on top of the curve for d+h=35. 0.3 0.75 ® 02 s I W 0.15 E 0.1 0.05 6 8 10 12 2 4 0 htilda Figure 3.15. Variation of the maximum growth rate with htilda for various d/h ratios. 31 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3.2.5. Effect of Surface Tension and Viscosity Ratio The stabilizing effect of surface tension acting on the statically unstable interface was evaluated using the two-layer model with time derivatives. The variation in the growth rate with surface tension applied is very similar to that shown in figure 3.1( for the unbounded layers), where it is seen that for any surface tension, the region of instability is limited to wave numbers between 0 and the cutoff wave number. The cutoff wave number versus the surface tension parameter T is shown in figure 3.16. The maximum growth rate that resulted for each of the T is shown in figure 3.17. No stabilization occurs at a finite value of T . As can we seen in figure 3.17, the peak growth rate is reduced to approximately half when T is about 50. The variation in the maximum growth rate with surface tension for the two boundary conditions, stress free and no slip, is shown in figure 3.18. This confirms the previous comment that the form of the boundary conditions (no-slip versus stress free) has virtually no effect on the peak growth rate when T is small (or zero) and h is large ( h =20). The effect of viscosity on the growth rate versus wave number curves is shown in figure 3.19. The initial curve used for most cases is shown for r2=0.9, s2=l, h =20, D=2.5. The Boussinesq approximation was made by setting r2=l .0. The change in the baseline (BL) curve is seen to be slight. About the BL curve with Boussinesq approximation (B L ba ) the viscosity ratio was varied by an order of magnitude either side of s2=l. The results are shown for s2=0.1, 1,10. It can be seen, as one would expect, that increasing the viscosity of the lower layer with respect to the upper layer has a damping effect on the growth rate. This result is opposite to that of Newhouse and Pozrikidis (1990). The difference in effect is seen to be explained by the fact that their thin film analysis had the boundary at 32 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . the bottom with the lighter, less viscous layer rising away from the boundary while this analysis has the boundary at the top with the lighter, more viscous layer moving to the boundary. 1.6 1.4 12 I 3 I °» I 0.6 Ü 0.4 100 1000 10000 10 0.1 1 100000 Ttilda Figure 3.16. Variation in cutoff wave number with surface tension parameter T tilda (r2=0.9,s2=1.0, D=2.5, htilda=20, stress free model). 33 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . 0 . 3 Ttida=0 limiti0.28TO) I o 0.15 E 3 I 0.1 0.05 1000 100 10000 10 0.1 1 100000 Ttilda Figure 3.17. Variation of maximum growth rate with surface tension parameter T tilda (r2=0.9, s2=1.0, D=2.5, htilda=20, stress free model). 0.3 1 o > 0.2 E I 0.15 S I E 0.1 a E •i* 0.05 Stress free nc slip 100 10 1000 1 0.1 Ttilda Figure 3.18. Variation o f maximum growth rate with surface tension parameter T tilda for two different models (r2=0.9,s2=1.0J3=2.5, htilda=20). 34 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.4 0.35 0.3 t2= 0.9. s2= s2=1 I 0 - 2 I 0.15 to 0.1 0.05 0 0.4 0.6 0.8 1 1.6 1.4 beta-nondimensional wave number Figure 3.19. Two layer model with stress free top (htilda=20, D=2.5, T =0). 35 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 33 The Effect of a Stabilizmg Upper Interface In the previous section we discovered that the instability of a density interface separating a heavy fluid superposed over a lighter fluid is unaffected by the thickness of the layer o f heavy fluid when h >6. This result raises the question; What might be the effect of a stable density jump across the interface defining the upper level of the heavy fluid layer? The purpose of this section is to explore this question. The configuration we analyze consists of a single layer of heavy fluid having density pc, viscosity pc. and depth h positioned between an unbounded layer of fluid with density pi<pc and viscosity pi positioned above the heavy layer and another unbounded layer having density pi<p2<Pc and viscosity p: below the heavy layer. The base-state structure for this model is shown in Figure 3.20. A new parameter, A say, emerges in this case measuring the stabilizing influence of the density jump across the heavy layer to the destabilizing density jump across the lower interface: that is, A = = (3.16) Pc-P2 1-^2 Pc When the thickness of the heavy layer plays a significant role in determining the growth rate characteristics of the instability developing on the lower interface, we expect that a positive value of the stability parameter A > 0 will exert a stabilizing influence. Based on the results presented earlier, the influence of the stability parameter A can be expected to become negligible when h >6. The vertical velocity eigenfunction for the conGguration shown in figtire 3.20 is given by (3.7) for the heavy layer contained within 0 < z < h and by the following expressions in the unbounded layers contiguous with the heavy layer = z>A, (r) = Ce*" + De*":", z < 0. 36 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . These later expressions satisfy the requirements of boundedness as z -> ± 0 0 . The eigenvalue relation defining the instability characteristics for the configuration under study is obtained by application of the matching conditions (2.29) at the equilibrium levels z = h and z = 0 o f the two interfaces. Carrying through the algebra, and using the intrinsic length and time scales to define non-dimensional variables, the eigenvalue relation is given by det D=0, and the elements of D are shown in Table 3.3. As anticipated, the parameter A appears explicitly in the second-to-last row in the determinant which derives from the normal stress matching condition at z = h. The variation in growth rate with wave number follows the same profiles as previously developed for the two-layer case. Figure 3.21 illustrates the growth rate © versus wave number p for the three layer case with an h =20 and equal viscosity ratios (si=S2=l), and density ratios of ri=0.8, r2=0.9, A=1.0. Setting the density ratios equal to each other and the density stability parameter equal to zero (ri=r2, A=0) consistent with the Boussinesq approximation produces the results shown in figure 3.22. The variation in growth rate with wave number is shown for three parametric h ’s (1, 3 and 5). It can be seen that the maximum growth rate © m ax for large h ’s is very similar to the two layer case and gives a value of © m ax of approximately 0.2875. The three layer model is useful in being able to elucidate the effects of the stability parameter A on the peak growth rate. The variation in © m ax with A for parametric h ’s (0.5 to 10) is shown in figure 3.23. What can be seen is that for h >6, the stability parameter A ceases to have any effect on the peak growth rate. The variation in © m ax with h for various parametric values of A is shown in figure 3.24. As stated for figure 3.23, after h =6, variation by six orders of magnitude in A causes no significant change in © m ax- 37 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . u I i i2 I I 1 ' ^ _ - ^ _ V t ^ 2 T3 •g % 3 \ " ^"^T ' + '< q P 0 » ^ Q r ^ c^ + !< + f * < # ca. # * • ^ ' f a ■ * . T % '&z< '< 3L h> '< 3 5 T ' Ÿ ^ ■S 7 , - ^ 3 M '& S % % _%_ e * ^ + l-i I \ t ^ ' ê S g ■ V « Î •e f f r f O S a. + i< c o - a » « + ® e x “ - j ^ 4 k 4 > g I O - O '§. ® ® 1 'S. Æ 'T '& L _ X C=L O 3 + <N CCL ^ r—^ 1 S ^. = v = - ; = f | I ^ ' < 0 > U cn (< G ) * & % _!aj__ !< I ■ ’ - • i't.'è î ' ii 38 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . — pl-> z = h stress free (W zz= 0) z = 0 pl < p2 < Pc Figure 3.20. Three layer model in an unbounded domain. 0 1 I (0 g E 0.3 02 0.15 0.1 0.05 0 0.8 12 0.4 0.6 1 1.4 0 2 0 beta-nondimensional wave number Figure 3.21. Variation in omega with wave number beta for a three layer case, unbounded top and bottom, using intrinsic scales (rl=0.8, r2=0.9, delta=1.0, sl=s2=1.0, htilda=20). 39 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0 . 3 0.25 I 0 9 I 1 0.15 0 1 o V ' ° 0.05 1.6 0.8 1 1.4 0.6 0.4 0 02 Beta-nondimensional wave number Figure 3.22. Variation of growth rate vs wave number for a 3 fluid configuration as shown in Figure 3.20 (rl=r2=sl=s2=l, delta=0.0). 0.3 htiUa=10 0.25 I 1.5 0.15 2 O 0.1 0.5 0.05 100 10 1000 1 0.1 0.01 0.001 Delta Figure 3.23. Variation of maximum growth rate with density parameter -A for various h . 40 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0.3 0 2 5 .001 S a 0.15 2 O 1000 0.1 0.05 9 10 7 8 6 5 3 4 0 1 2 htilda - nondimensional core layer height Figure 3.24. Variation in maximum growth rate with htilda for various density parameters. 41 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3.4 The Effect o f a Stratified Layer Below an Unstable Interface In this section we explore the influence of a stably stratified fluid positioned below a finite layer o f heavy fluid. Intuitively, one expects that the ‘*stiffiiess” created by stable density stratification inhibiting vertical motions will restrict the growth of disturbances on a statically unstable interface lying above a stratified layer. The analysis and computations presented here are directed toward quantifying this effect. In order to isolate the intended effect and minimize the required algebraic and computational burden, we analyze the two layer configuration shown in Figure 3.25. Also, in order to simplify the analysis we will look for a solution that is valid for moderate to small wave lengths and make the creeping flow assumption so that we can ignore the inertia terms. As noted in section 3.1 the solution for the vertical velocity in the upper layer is given by equation (2.13), and eliminating the inertia terms reduces to a forth order ODE for vertical velocity with repeated roots as in equation (3.9) of section 3.2. The equation defining the vertical velocity in the density stratified layer is given by equation (2.12) of section 2.0. The solution is given by : W 2 (z) = + C g-*- + , (3.18) (3.19) aVjk is related to the variable X appearing in (2.18) by X = ±*p ^ . P ^ is obtained fi"om equation (2.18) o f section 2.0 by neglecting the terms involving a in the expression for K (i.e., neglect of inertial terms). 42 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . The scaling for length and time are taken as I = h, t = (p c -P z k A . This gives the following: a = kh, y = — -, r = ( P c - P z k ^ ' {Pc-P2)gh^ and we d fX y additionally denote D = —, and 5, = Î 1 z = h z = 0 z = -d Figure 3.25. Stratified lower layer with boundaries at the top and bottom The term P ^ can be recast in a convenient form by noting that N is the dç) ôp g — Brunt-Vaisala frequency, N s ---- — = — — , using the parameters noted in P2 P2 figure 3.25. To measure the stabilizing effect of the stratification in the lower layer, we introduce a length scale 5 which measures the vertical distance from the heavy fluid at the unstable interface to a vertical position in the lower layer with equal density. This length scale is shown schematically in figure 3.25. Using this scale 5, two related stratification parameters m s — — — = — and 5p d 43 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . d s m' 5 —m = — can be defined. The square root factor in the definition of P can h h then be rewritten in terms of m ’ m'ya S 2 (3.20) It should be noted that m'as a . stratification parameter, denotes stratification in an inverse fashion (i.e. as m' -> 00, — -> 0). This causes the expression in dz equation (3.19) to be given by. 1±/ (3.21) Development of the equations involving the matching conditions o f equation (2.29) and boundary conditions of equation (2.30) and (2.31) yields a coefficient matrix D. The dispersion relation is given as Z ) = det [D] = 0 (3.22) The coefficient matrix D is given in Table 3.4. Utilizing intrinsic scaling the coefficient matrix of Table 3.4 is modified to the form given in Table 3.5. Table 3.4. Coefficient matrix of the stratified lower layer without intrinsic scaling for a stress-free condition on the upper boundary. 0 0 0 0 e“ de"* t u T gP-oO g-P.oD gP .«0 0 0 0 0 0 0 0 0 e -(2-a)e-* (2+a)«' P _g-p.«o - p y “ « P,e-»-“ " 0 0 0 0 1 1 1 1 -1 - I 0 0 -P . 3- -P. P. 1 - I -1 -1 ^(i+Pi) 4 1 + P!) -2 -2 +2 -2 ^orirP. X -•y»ïP- X w P . X -PjOYP. X (l-fa * )- ( l-fa ')+ A A (3-p!) (3-K ) (3-3Î) (3-P!) 2ay 2ot U U 44 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . t 0 8 J % § oo c s tu J « o o o 7 7 •o c ca. + S ^ I tn I s ‘ â ^ b5 % ) o o o I 7 0 \ o o 7 7 7 ^ ^ 1 i î c i * '$ o a c T - U tu 0> , cQ . 3 rs % a ^ =2. cû. j '-r-^iïî I a î — a o î = > ' “ -<f- + â ‘f •I “ T I ^ X -w " S . « ^ ' I I ô « X " S . 7 7 X ^ —_ _ en < 3 C L I gn *^' ju o * o ' » » | - ^ e û . + § ; c a . R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . The results for the stratified lower layer indicate that, as would be expected, the degree of stratification as depicted by m ’ greatly effects the instability growth rate. Figures 3.26 through 3.29 provide the relationship of growth rate as a function of wave number for m's of 1,5,10 and 100 respectively. They also show the effect ofD=d/h, as a parameter. As can be seen the growth rate as a function of wave number gets larger as the stratification parameter gets larger (smaller stratification). It can be seen that the stabilizing influence of a stabilizing density stratification does enhance the stability of an overlying step or unstable density layer. The effect of D is that as D gets larger, the growth rate for a specific stratification parameter also increases. 0.05 0.045 0.04 0.035 0.03 I '1.5 0.025 S ^ 0.02 0.015 0.01 0.005 2 2.5 3 1.5 0.5 1 0 V W aw number Figure 3.26. Variation of growth rate with wave number for a configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=l). 46 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . 0 . 1 2 0.1 D=3, 0.08 0.06 2 O 0.04 0.02 3 1.5 2 2.5 0 0.5 1 Wave number Figure 3.27. Variation of growth rate with wave number for configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=5). 0.14 0.12 D=. 0.1 « 0.08 2 0.06 0.04 0.02 2 1.5 3 0 2.5 0.5 1 Wave number Figure 3.28. Variation of growth rate vs wave number for a fluid configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=10). 47 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0.14 0.12 2.5, 0.1 1.5 g 0.08 A 0.06 0.04 0.02 3 2.5 0 0.5 1 1.5 2 Waw number Figure 3.29. Variation of growth rate with wave number for a configuration shown in figure 3.25 (s2=l, Ttilda=0, m’=100). The effect of D and stratification parameter m ’ on the maximum growth rate is shown in figure 3.30. As the stratification is decreased in the lower layer, it can be seen that the maximum growth rate increases at each D. Also, for a specific stratification (m ' = 10 say) the maximum growth rate increases as a function of D and reaches an asymptote at a specific D (in this case y = 0.127 at D = 4 ). At higher stratifications the maximum growth rate is reached at a lower value of D (i.e. for m ’ = 1, y ^ = 0.047 at D = 1.5 ). The wave number at which the growth rate reached its maximum is shown in figure 3.31 as a function of D and with the stratification parameter. The non-smoothness of some of the curves is an artifact due to the granularity of the eigenvalue search routine in the wave number. 48 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . 0.14 m’ =10 0.12 0.1 I 0.08 Ô » I 0.06 E ■ § ^ 0.04 0.02 6 4 5 2 3 0 1 D=d/h Figure 3.30. Variation of maximum growth rate with depth ratio D for a fluid configuration as shown in figure 3.25 (s2=l, Ttilda=0). 1.8 m'=10 o * 12 m'=10 0.8 0.6 0.4 0.2 5 2 3 4 6 0 1 D=d/h Figure 3.31. Variation of wave number at maximum growth rate with depth ratio D for a fluid configuration shown in figure 3.25 (s2=l, Ttilda=0). Since intrinsic scales were determined to be the most appropriate scaling for the preceding sections, the stratified eigenvalue relation was recast using 49 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . intrinsic scales. The scaling used was given by equation (3.4) of section 3.1. This introduced the scaled upper height h into the matrix. In the new matrix the following variable replacements are made, a = pA , y = —, T = T h^, D = D = —, and h h P* = l± f w'mAp ^2 > (3.23) With intrinsic scaling and the above parameters incorporated into the eigenvalue relation, the growth rate as a function o f wave number was developed and is shown in figures 3.32 and 3.33 for /x’s of 5 and 10 respectively and a constant stratification value of 5. It can be seen that as the unstable upper layer thickness is increased (h increases) the maximum growth rate also increases from 0.55 to 1.1 for A of 5 and 10. The wave numbers where the maximum growth rates occurred were .27 and .13 for A of 5 and 10 respectively. The intrinsically scaled maximum growth rate value for A % 10 was 1.1 for this stratified case without inertia terms. The two layer case with intrinsic scaling without inertia terms yielded a maximum growth rate of 2.4 for an h =20, at a wave number of approximately 0.07. The result for the stratified case is entirely consistent with the two layered results utilizing the same assumptions in the derivation of the eigenvalue relations. 50 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0 . 6 § 0.5 g, 0.4 c 0.3 c 02 O ) 1 0.8 0 0.6 02 0.4 Beta- nondimensional wave number Figure 3.32. Variation of the growth rate with beta for a two layer case with a stratified lower layer ( htil=5, m’=5, s2=l, Tbold=0). 0 .8 a I 0.6 0 0 .2 0.4 0 .6 0 .8 1 Beta - nondimensional wave number Figure 3.33. Variation of growth rate with wave number for two layers, bottom layer is stratified. Stress free on top and stratified layer is bounded at bottom (htil= 10, m’=5, s2=l., Tbold=0). 51 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 4. Formulation: Analysis of continuously-Varying, Time Dependent Density Distributions 4.1 Introduction The dynamics of sedimenting surface gravity currents and the role instability has on the sedimentation of particles and flow of the interstitial fluid has been studied by Maxworthy (1998) and is of great interest in geophysical flows. For this reason it is of interest to study a temporally developing, heavy intrusive layer flowing over a lighter layer of fluid. To model the effect of a heavy intrusive layer which is developing temporally as a layer over a fluid of finite depth, we will idealize with a model where the density and height of an upper layer are arbitrary functions of time. We will also assume that the fluid is incompressible and that the viscosity is constant throughout the layer. 4.2 Model Development The non-diffusive, linear equations for motion in a vertical plane are given below, with u,w, p and p as perturbed variables and p$ = p;(z/) as the base state density: = 0, (4.1) = -Px + (4.2) Px^, = -Pz - P g + Px^^ (4.3) p , + p > = 0. (4.4) We have subtracted the hydrostatic balance relation -psg for the base state dz in writing equation 4.3. Eliminating the pressure from the momentum equations and assuming constant viscosity ps = Ps v* leads to the vorticity equation which can be written in the form P x -r^ ^ ^ + P > x , = -P x .g + P x^^^^^' (4.5) at 52 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . Using equation 4.4 to eliminate p in favor of w leads to the single governing equation for the vertical velocity. g — +V, Ps .P , df Ps dt dz^ xw. = 1 dp. (4.6) Ps dt The time dependent coefficients do not allow an instability analysis with a constant exponential growth rate. However, since the equation is homogeneous in the horizontal direction we can consider a single Fourier mode, w{x,z,t) = W(z,t)e‘ ^ . (4.7) Our governing equation becomes in terms of the modal amplitude W , dz p , dz ^ 1 dp; 3 'I p, dt dz^ W ,^ Ps „2, - g — k W + v ^ P. W (4.8) The instability problem is now reduced to the solution of a partial differential equation for the modal amplitude with coefficients that vary in both space and time. The selected model for the time varying intrusive upper layer is shown in figure 4.1. The upper layer of depth h will have a smooth density distribution with a maximum increase in density of po r(t) at the upper surface compared to the ambient density po- The density of the intrusive layer varies continuously down to a depth of h where the density reaches po and remains constant. The entire fluid depth is d and both extremities of the layer are assumed to be stress firee (z = 0, d). 53 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . Cosine distribution stress free ÇJVzz-0) at txrih boundaries Figure 4.1. Continuously varying, time dependent density distribution model Since the boundary surfaces are both impermeable = 0 at z = 0,d) and stress free = 0 at z = 0,d) , the vertical structure of the solution for W(z,t) can be represented in terms of a Fourier sine series given as ( 0 s i n ^ ^ j . ( 4 . 9 ) For analytic convenience we choose a base state density distribution of the form p,(z) = po 0< z< d-h Then we arbitrarily specify the coefficients appearing in the governing equation (4.8) in terms of the following Fourier series: _ L ^ = y £ - W c o s f ^ \ Ps ^ o h \ d ) - L ^ = y M ) s i n f u l p , dt f ^ d / £ l = y £ = t ( Ü s i n f ^ ) p , ^ H, \ d J (4.11) (4.12) (4.13) (4.14) 54 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . The variables 7% are scaling factors. The variable is the square o f the magnitude of the Brunt-Vaisala fret uency given by -g d p s Ps ^ to is a time scale over which p changes, and Ho is a characteristic length scale for a time-varying depth h(t). These scales are not all independent by virtue of the specified density distribution (4.10). For example, the total density deficit (por) divided by the maximum value of the density gradient ( given as f Po — — | ) az y 2 h j H 2 taken at the midpoint location (z = d-h/2) provides a value = — = 0.63662. ft 7 C Substitution of equations 4.9, and 4.11 through 4.14 into 4.8 yields the following governing equation for the coefficients a„ in terms of the specified coefficients bm, C m , dm and em arising from the specified density distribution: Z c o s ^ - ^ jZ '^ » )[â„ + ^ „ ')a „ ]s in ^ ^ j (4.15) The factor kn is defined by kn = nrc/d. Using the orthogonality properties of the trigonometric fimctions we obtain a coupled system of ordinary differential equations for the amplitude fimctions a„(t): 55 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . {k^ +kl)(â^ •*'"/--) 1 ^ "^2^0 Ç""* +vX *" + kl„)a^„] (4.16) 2 f f ^ / ^ m \ ^ p - m ^ p - m ^/M -w ^p+w] 1 A # ~ ' J Z J t S ^ " [^ P -" ^ < » -" ~ ^ / H - « ^ i » + " J — •^ “ 0 * 0 I This coupled set of linear ODE’s defines an infinite dimensional system. However, in practice we specify a truncation of the system to obtain an approximate description of the linear instability. As in the earlier sections, we choose to use the intrinsic scaling I I fv^V fv, 1 3 / = — and T 0 = V s’} • to obtain the governing nondimensional equation. p;+—a i \2fo C_ T l2Llift2 + —L i s 2 r„ 2H, t. d„y p *m p * m SL S.Q i 2 /o 2/fo to (4.17) p -jw rp-w p+ « p-w = 0, where +*J)> Y = /(p + m )^ , and a = / 1 The notation (d) now denotes the derivative of a(x) with respect to the non-dimensional time 56 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . variable t. In preparation for numerical solution, the second-order non- dimensional system is rewritten in terms of an extended system of fîrst-order equations using the variables = â, and A, = à,. Equation (4.17) becomes with this substitution. M • I C_ T ^mY p+nf^p+m 2 /o 2H, /o To o2 / p + m — — PL -■ 2 /o 2//„ r P-» 0 ‘0 p*" P-" £e.L2.R^ * • 0'‘^0 - 2 2 /„ p — # » (4.18) = 0, Using the selected forms (4.11-4.14) for the coefficient functions b„, C m , d„ and of equation (4.8) and the orthogonality properties of sine functions, the coefficients appearing in (4.17) are given by: d (4.20) (4.21) 57 R ep ro d u ced with p erm issio n o f th e copyright owner. Further reproduction prohibited w ithout p erm issio n . The above integrals were developed and their solutions used in subroutines for the mechanization o f the solutions of the ODEs. The governing ODEs developed fix>m equation (4.18) by specifying a truncation M and developing all of the equations in p up to M produces two 2Mx2M matrices. The resulting matrix equation is given as A i + Bx = 0, (4.23) where i is the vector of all of the À^'s and à^'s and x is the vector of and A p’s. In order to cast the equation in a form allowing the use of standard integration routines, we rewrite (4.23) as x = Cx, where C = -A B . The solution to this set of equations is developed by using a fourth order Runge-Kutta equation solver (Press, et. al., (1992)) because of its time tested robust capability. The resulting A and B matrices are shown in Tables 4.1 and 4.2. 58 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . s % " Q «% I J I 1 73 I 0 > E f I e / 3 g 8 1 E 2 % 2 5 f J 8 1 % C N + ca. 3 u fN a: < N a; ( N + en. er|<^ + a . — i C 5 es 0 ^ 0 r * O p # o es P i o ► * k * » * • o p * o P * o 5: * a; $ V 3 3 •? C S L :( u a: es ■ l< ^ -Q >c fS es a. G S L < N H es " C " c 5 l T ri a. (T I .s i < 1 I II < •a M I '• 0-|cN > - a: es < !< c I I < I c a . ^ 1 - 0"|«N ca ç T M es * Q - >c a :“ V es < I I I ea. en. en. f f „ -U " - i « 3 •Q I« o e s 3 u I I < I I < 59 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . ■? § I 1 • f I 1 8- •o 1 8 0 1 E S I ea t o $ X . % % i: «S 3 » T 5; « N N 5 Q Z L 0 0 c CO > c a 3 § •S c 8 a *c % 3 U r n c2 ea % 1 2 cs S I 8 1 O II e a " I I e a " •a % T : ï; < N o T n % o T ï: C N $ % a: (N 60 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 5. Results for a Steady, Continuousiy-Varymg Density Distribution The results for a steady, continuously-varying density distribution as shown in figure 4.1 were developed for a twenty (20) coe£5cient truncated series for the vertical velocity W(z). The initial sensitivity of the vertical velocity result to the number of coefficients retained in the series indicated that after 10 coefficients no significant accuracy increase resulted in higher series termination numbers. Coefficient terminations were performed at 5,10,20, 30 and 50. The difference between 20, 30 or 50 was insignificant as compared to 10, however the difference in run times became significantly longer for termination numbers higher than 30. In the derivation of the coefficient fimctions for the bp, C p . dp, and ep in equations (4.18) through (4.21), expressions resulted which had singularities near some of the increasing factors within the sununations. As d/h increased, the value where the singularity might occur happened at a term later in the expansion. Care needed to be taken such that the summations always included at least 10 terms beyond the term which was close to being singular. Therefore, for this reason, series termination was always taken such that there were at least 10 coefficients beyond the coefficient where divergence in a series expression occurred. With a d/h = 2.5, a twenty coefficient summation for the vertical velocity gave around seventeen terms beyond the near divergent term. The numerical model was tested to determine a satisfactory integration ending time in order to calculate the vertical velocity. This was done so that the startup transients would not affect the stability results. The response o f the model was tested as to which coefficient was best to provide an initial kick to minimize the initial transients before the uniform growth rate phase was realized. This was done in order to decrease the total integration time required. Also, as various parametric values 61 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . were used in the different case runs, numerical instabilities were observed. The stability criteria for the numerical model were empirically determined over the ranges of cases studied. For a more detailed description o f each of these numerical modeling developmental issues, see the Appendix A of this report The steady density distribution gave the growth rate as a function of wave number as shown in figure 5.1 for a wide wave number view and figure 5.2 for short wave number behavior. As can be seen firom figure 5.1 and 5.2, these results are similar in form to the results of the two layer case as given in section 3.2 and shown in figure 3.5. The difference lies in the magnitude of the growth rate for the cosine distribution which is approximately 2.5 times lower than the two layer case and the wave number of the maximum growth rate occurring at approximately 0.2 versus a value of about 0.4 for the two layer case. 0.12 o 0.1 s I 0.08 a 1 o g 0.06 I g f 0.04 a E " 0.02 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 Alpha-waw number Figure 5.1. Time dependent layer growth rate vs wave number, zero time rates. 62 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0 . 12 _ 0.08 o > * a c I 0.06 I I 0.04 É ^ 0.02 0.15 02 0.3 0.35 0.4 0.05 0.1 0 alpha- nondimensional wave number Figure 5.2. Variation of growth rate vs wave number for the static case results (dh=2.5, dl=100.1, htilda==40, r=0.1, tratio=91,20 coefficients, 40X40). The results of Batchelor and Nitsche (1991) indicate that for a pseudo sinusoidal shape (i.e., second derivative of an error function profile) separating two regions of equal density, the reduction in growth rate fi’ om a two layer case could be as much as a factor of approximately 3. Their maximum growth rates for these smooth profiles occurred at wave numbers of 5 4 that of the discontinuous profiles. For unbounded fluids with an unstable density on the top and a smooth transition layer between (i.e., error function profile), the differences when compared to a two layer model were not as pronounced. Although the case studied here differs fi'om those considered by Batchelor and Nitsche both in profile shape and in boundary conditions, some similar trends are observed. The maximum growth rate as a function of the layer depth h and the wave number where the maximum growth rate was reached as a function of layer depth are displayed in figures 5.3 and 5.4. These results look very similar in form to the results for the two layer case given in section 3.2, figures 3.6 and 3.7. The maximum growth rate for the smooth profile in figure 5.3 shows a slight decrease 63 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . after a layer height o f approximately 10 versus the constant maximum growth rate of the discontinuous profile. The smooth profile reaches its maximum value after an A > 6, which corresponds to the value of h where the two layer profile has reached its asymptotic value. The variation in vertical velocity and surface displacement as a function of the integration time is shown in figure 5.5. As can be seen the vertical velocity and surface displacement reach constant linear growth rates very quickly in the integration after approximately 15 to 20 time units. Another point to observe is that the surface displacement remains very small, in the linear regime, over a long integration time, for t « 130 time units. 0.16 0.14 0.12 0.1 & 0.08 0.06 0.04 0.02 25 20 10 15 0 5 hb'lda Figure 5.3. Variation in maximum growth rate as a function of htilda, zero time rates. 64 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . 1.6 « 1.4 a» 0.8 0.6 0.4 02 25 20 10 15 0 5 htilda Figure 5.4. Variation in the wave number where maximum growth rate was achieved versus htilda. IE-2 1E-3-. IE-4 IE-7 1E-8-1 IE-9 E-10 5 1E-11nr IE -1 2 -■ OE+0 1E+1 2E+1 3E+1 4E+1 5E+1 6E+1 7E+1 8E+1 9E+1 1E+2 time - tau units Figure 5.5. Variation of Vertical velocity and Zeta with tau time units (htilda=5, alph=.43(alph at sigmax), 60X60, dt x n^2/htilda^2< 16, n=60). The variation in the growth rate with the upper layer thickness is shown in figure 5.6 for a constant wave number of 0.31 and the maximum growth rate (over all wave numbers) for each h . This variation highlights the strong dependence of 65 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . the growth rate on the specified wave number. A wave number of 0.31 was selected because it represents the wave number where the growth rate is maximum for an upper layer thickness o f h =10. This result will be used with a later simulation of a time varying layer thickness. The variation in the growth rate as a fimction o f the d/h ratio (the ratio o f the total depth to the upper layer thickness) with the upper layer thickness A as a parameter is shown in figure 5.7. It can be seen that the d/h ratio has a similar effect on growth rate as h had, however, the growth rate for each d/h still exhibited a slight increase with d/h after the initial rise for a d/h > 8 and for A >4. Viewing the curves at constant d/h shows the dependence of growth rate with h and the slow drop off after A > 6 after the initial steep rise. 0.16 0.14 Signr amax 0.12 Sigma at a= 31 0.1 5,0.08 0.06 0.04 0.02 a 9 10 6 7 5 1 2 0 3 4 htilda Figure 5.6. Variation of the maximum growth rate and the growth rate at a constant wave number of 0.31 with htilda. 66 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0.18 0.18 0.14 S 0.12 § a 0 . 1 8 % 0.08 Ê • f 0.06 h = 2.001 4.001 6.001 0.04 8.001 0.02 18 16 12 8 10 14 2 4 6 0 d/h Figure 5.7. Variation in the growth rate at 50 tau for various d/h and htildas (alpha=0.31). 67 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 6. Results for a Time Dependent, Continnonsly-Varying Density Distribution The results for a time dependent, continuously-varying density distribution as given in Ggure 4.1 were developed for 3 different time dependent density distributions. Each different time dependent distribution ramped the upper layer of higher density from a thickness of h = 2 to a final value 10 over a varying number of time units. The times were chosen to correspond to 1,2 or 3 e-fold times based on the maximum growth rate of a stationary profile for h = 10. The corresponding ramp times were 7.024, 14.048 and 21.072 tau time units respectively. The total combined depth d of the upper and lower layers was kept constant in order to obtain the destabilizing effect of a growing upper layer density distribution. For this reason the instantaneous ^ was changing (from h 17.5 to 3.5) as h changed (2 to 10). The simulation was initialized and set in motion, then it was allowed to run for 50 tau time units to damp out start up transients. At this point the depth of the upper layer was allowed to vary according to the prescribed rate from h =2 to 10. The local (in time) growth rate was calculated by evaluating the derivative of the logarithm of the vertical velocity at a fixed spatial position. The resultant effect on the growth rate for the time varying upper depth is shown in figure 6.1. The results in figure 6.1 indicate for the 7.024 tau case that the growth rate started at a value of 0.14 and had an initial peak at approximately 0.15 and a low of 0.118 . These results are shown for a constant wave number, a = 0.31, corresponding to the maximum growth rate at h =10. The final value of the growth rate of 0.143 was obtained 3 time units after h reached 10. The other h ramp cases illustrated similar initial peaks and secondary imdershoots with 68 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . smaller peaks and larger under shoots being exhibited for the longer h ramp times. Also, the value of the growth rate at the end of the h variation indicated an overshoot in the growth rate with longer times that subsided to the steady state growth rate value of 0.143 after the transients were allowed to settle out. The results of figure 6.1 are surprising in that the growth rate excursions are not extreme but have an oscillatory wave form dip firom the starting value o f 0.14 to the lowest value o f 0.11 for the 21.072 tau case. This would lead one to conclude, at least to first order, that the growth rate might be assumed to be constant over the h variation period. O f course, these results are obtained for a fixed wave number. 0.16 ^t3da va iation ov^_LX t a n ^ 0.14 14. 0.12 21.072 Q ) 0,0.08 I 0.06 g 0.04 0.02 25 15 20 0 5 10 tim e-tau units Figure 6.1. Variation of instantaneous sigma- growth rate with time with three different time slopes for htilda variation (dt=.0175 for stability at htilda=2). The effect of changing the wave number was addressed by repeating the calculations at three different wave numbers. The results are shown in figure 6.2. The variation in the growth rate is shown as a function of time for three separate wave numbers corresponding to the wave numbers of the maximum growth rates 69 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . for h of 2,6 and 10. It can be seen that using the wave number of a = 0.31 for the ending /t = 10 generally yields the maximum growth rate (except for the initial peak) over the total h variation from 2 to 10. Therefore using the wave number of a = 0.31 should be illustrative of the maximum growth rate response. 0.16 Sigr I ax (d/h=3.5) 0.14 0.12 Sigma 0.1 i.0 .0 8 Sigma (a =.77) I 0.06 0.04 0.02 16 0 2 12 14 8 10 4 6 time - tau units Figure 6.2. Variation of growth rate with three different wave numbers for a transient case where htilda is ramped from 2 to 10 in 7.024 tau time units. The variation in the ratio of the instantaneous velocity at a fixed r» level (z, = 0.7619*d) to the initial velocity as a function of time is shown in figure 6.3 for each h ramp profile. This value of z , was chosen to be slightly above the bottom of the interface and below the middle of the interface where - < d - z ^ <h and —— — = 4 2 . The same linear growth rate is achieved by 2 d — z^ each h ramp profile. It can be seen that in general the velocity ratios achieved are slightly displaced from each other by the integrated effect of the dip in the growth rate profile during the h transient (with the deficit in growth rate profile being more pronounced as time of h variation increases). The value of the 70 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . slowest increasing velocity ratio (for the 21.072 tau case) at the end of 25 time units is approximately 29 versus a value of approximately 34 for the constant steady state rate obtained firom the end conditions. This yields a maximum reduction of only 15% in the velocity ratios. The growth rates after the initial transients are the same so this deficit would remain throughout the instability growth. Therefore, it appears that, to first order, a conservative estimate of the growth rate could have been the linear steady state growth rate estimate of 0.143 at the ending h value of 10. 100 I O W/Wi for 7.024 S c c s W/Vi for 14.048 W /W lfjr 21.072 1 o I I 25 10 20 0 5 15 time - tau units Figure 6.3. Variation of WAVi with time for three different htilda variation cases. 71 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 7. Conclusions The effect of the stabilizing influence of confining boundaries has been studied and their effects determined. For a two layer system it has been shown that, utilizing intrinsic scaling, the maximum growth rate of a disturbance increases firom zero at an A =0 up to its maximum value at A = 6 to 8 for the heavy upper layer. The wave numbers corresponding to the maximum growth rates shift firom a value of about unity for A =1 to a constant value of about 0.4 after A ^ 6. The physics of the phenomena is that the effect of the wall is no longer felt after the layer has grown to a distance of A approximately 6 to 8 because the fastest growing wave length is then commensurate with the layer depth. When the allowable cell size reaches the aspect ratio of unity, with respect to the fastest growing wave length , the effect of the wall on the instability mechanism at the interface is no longer dominant. The fastest growing wave length for the two layer model is "Kj s 16. This can be seen in figure 3.7 where the wave number asymptotes to a wave number of approximately 0.39 corresponding to a wave length of 16. We can also see in figure 7.1 that as the ratio reaches one half, the effect on the growth rate due to the upper boundary disappears. 72 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 0.3 0.25 0.1 0.05 0 0.1 0 2 0.3 0.4 0.5 0.6 0.7 0.8 h Figure 7.1. Variation of maximum growth rate with upper layer thickness in relation to most unstable wave length. Two layer case as given in figure 3.7 (s2=l, Ttilda=0, r2=0.9, D=20). The effect of boundary conditions on an upper boundary in proximity to an unstable upper layer was determined for both stress fi'ee and no-slip. The effect of the boundary conditions are most pronounced at the lower A 's of 1 and 3 as shown in figure 7.2 and 7.3 respectively. Equivalent results for the unbounded three layer case are included for comparison with the stress-firee and no-slip conditions. The three layer case introduces another degree of fireedom where the upper interface is allowed to deform and can be seen to be the most unstable of the three boundary conditions. The most stable case and most constrained case is the no-slip boundary condition. 73 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.18 0.16 0.14 0.12 I 0.1 S 0.08 O 0.06 0.04 0.02 0 0 2 0.4 0.6 0.8 1 2 1.4 1.6 1 V\^ve number Figure 7.2. Variation of growth rate with wave number for three different models at an htilda=l. The three layer model has rl=r2=l, sl=s2=l, Ttilda=0, delta=0. The two layer models have r2=.9, s2=l, Ttilda=0, D=20. 0.3 02 I 0.15 2 O 0.1 0.05 0 0 2 0.4 0.6 0.8 1 2 1.4 1.6 1 SF3 3L-3 NS3 Wave number Figure 7.3. Variation of growth rate with wave number for three different models at an htilda=3. The three layer model has rl=r2=l, sl=s2=l, Ttilda=0, delta=0. The two layer models have r2=,9, s2=l, Ttilda=0, D=20. 74 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . The effect of neglecting the inertia terms in the linear stability model was evaluated for the two layer model. It was shown that the growth rates resulting from the derivation diverged from the growth rates of the model which included inertia terms at mid to small wave numbers (in the region where the maximum growth rate occurred). Also the wave number corresponding to the maximum growth rate moved to smaller values for the case with no inertia terms, making it inaccurate as well. The maximum growth rates achieved with the lubrication model were approximately 10 times larger than the model which included the inertia terms. This can be seen by observing the characteristic roots as defined in (2.14)andnon-dimensionallyas p] = p^+o> = p' in (3.5). This formulation of the characteristic roots includes the inertial terms. It can be seen that if the time derivative terms were neglected that it would be in essence neglecting the m term versus the P^ term in (3.5). When co reaches a value of approximately 10% of the P ^ term, then we would expect the formulation without inertial terms to diverge. We see in figure 7.4, using the results for the two layer case with inertial terms for the maximum growth rates and wave 0) P numbers at the maximum growth rates, that exceeds ten percent around an P a s 1. Looking at the results for the two separate formulations in figure 3.12 we see this is indeed the case. This indicates that to adequately determine the maximum growth rate and wave number behavior of linear stability systems, the inertia terms must be included. 75 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 1.8 1.6 0.8 O ) 0.6 0.4 0 0.1 02. 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 w a\e number at max growth rate-beta Figure 7.4. Variation in the ratio of max growth rate over wave number squared to wave number The effect of surface tension on the maximum growth rate and cutoff wave number was determined for a two layer system. It was shown that without surface tension the system is unstable at all wave numbers. However, with the addition of any surface tension, the region of instability shrinks to a specific band of wave numbers from 0 to a cutoff wave number. As the surface tension increased the cutoff wave number asymptotically decreased to zero as the inverse Bond number T tended to infinity. The effect of the upper layer viscosity in relation to the lower layer viscosity was also determined. A case with the upper layer viscosity being greater than the lower layer is less stable than the case with the lower layer viscosity greater than the upper layer. This result is opposite to that of Newhouse and Pozrikidis (1990). The difference in effect is seen to be explained by the fact that their thin film analysis had the boundary at the bottom with the lighter, less viscous layer rising away from the boimdary while this analysis has the boundary 76 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . at the top with the lighter, more viscous layer moving to the boundary. The result o f MenikofFet al. (1977) (in the case o f any interface far removed from any boundaries) indicates that the case of upper layer viscosity greater than lower layer and lower layer viscosity greater than upper layer are equally destabilizing from the equal viscosity case. The effect of a stabilizing upper interface or stable upper layer was explored. It was found that the general instability result was the same as that for a two layer model, however there was a secondary effect of the density parameter A (as defined in equation ( 3.16 ) of section 3.3) in the lower h range of 0 to 6. In general for the lower h ’s, a high A of 10 or 1000 caused a slight stabilizing effect. However after an h> 6, there was no significant effect on the maximum growth rates for a A excursion of 6 orders of magnitude (.001 to 1000). This is consistent with results described earlier where it was shown that the growth rate was independent of the layer depth when h >6. The effects of a continuously varying density distribution were developed for a cosine density variation in an upper layer. It was shown that the general effect is an approximate reduction in the maximum growth rate by a factor of approximately 2.05 (two layer a = .2875 divided by C T jq =.140 for d/h = 17, A = 10). This reduction is in general agreement with Batchelor and Nitsche (1991) for their Gaussian pseudo-sinusoidal profiles where they obtained a reduction factor of 3.2. The effect of a smooth, distributed density profile is shown in figure 7.5 in comparison to the two and three layer cases. It can be seen that the general behavior of the maximum growth rate for the smooth profile with non- dimensional depth is similar to that for the two and three layer cases with the effect on the maximum growth rate asymptoting to a maximum value at an A of 6 to 8. It can also be seen that the smooth profile causes the instability at the 77 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . interface to act as if the layer is only at an A » 1 for an equivalent sharp two layer case of uniform densities. In reality, however, the distributed profile has an A = 8 to 10. It was also observed that the time span over which the linear solution could be expected to be valid occurred over t w 130 time units, well beyond the time of a number of e-folding increases in velocity. 0.3 3L-.001 3L-1 3L-1000 02 2L-D20 cosine 0.1 0.05 6 8 10 12 0 2 4 htilda Figure 7.5. Variation in maximum growth rate with htilda for three different models. The three layer model has delta ranging from .001 to 1000 with rl=r2=sl=s2=l. The two layer model has stress free conditions at the top boundary and no slip on the bottom and D=20 with r2=.9, s2=l. The bottom curve is for the cosine distribution with density parameter r=. 1. All have Ttilda=0. The effects of temporally increasing the depth of the continuously varying density profile in an upper layer was shown to cause a delay in the growth rate, affecting the velocity. However, the growth rate ‘caught up’ with the steady state growth rate and from that point on the velocity grew at the steady state growth rate. This result caused the postulation that, to first order, a temporally deepening, continuously varying density profile effects on velocity, or surface deflection, could be estimated by using the end condition steady state growth rate. 78 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Part II Instability of a Rapidly Evaporating Liquid Interface 1. Introduction The study of rapidly evaporating liquid interfaces is important in a number of physical applications in the distillation industry and in the space industry (propellant processing). Densifîed cryogenic propellants are the propellants of choice in many applications due to the high performance requirements of future hypersonic vehicles. This densification requires special processing of the propellants to achieve the thermophysical properties required of the cryogenic fuels. Several densification techniques have been studied over the years from the early 60’s into the early nineties. These various processes have included freeze- thaw, helium gas/liquid injection, nozzle expansion, auger, and magnetic refrigeration to create densified hydrogen. A review of a number of these processes are given in Sonntag (1992) and Mann (1966). The physics of the propellant densification process need to be well understood in order to maximize the production process for densified propellant. This study seeks to define the effect of evaporation and surface tension of the liquid interface on surface stability which could be representative of conditions occurring in the densification process. The free surface stability while undergoing heat transfer from the upper gas region to the lower liquid region and the resultant mass transfer from the free surface into the gas ullage region needs to be developed to represent the propellant densification process. A model of the configuration is given in figure 1. 1. Problem formulations for this configuration given in Miller (1973) and Palmer et al. (1972,1976) were evaluated for their applicability to the problem 79 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . area indicated above. Generalized formulations for the linear stability of liquid films is given by Spindler (1982) and the jump conditions across the interface by Delhaye (1974). The effect due to surface tension are given by Levich and Krylov (1969). Miller (1973) described a problem o f free surface stability with different temperature gradients than the configuration in figure 1.1. Miller’s example (shown in figure 1.2) was heated from below as might be the case for distillation processes. The liquid at the fiee surface was vaporized into the ullage by heat transferred from the liquid below. P<Pv for a specific cyclic rate w .z Compressible G as Interfece is fixed m = evaporation rate with respect to the coordinate system Tj at interface h(x.t) or Çfx.t) u .x \v e r s u s z Liquid W 1- varfable f T TL = triple point temperature ! Z = -H 1 . Figure 1.1. Physical Setup of Propellant Densification Process Miller analyzed the effect of the vaporization and movement of the interface on stability of a liquid surface. He assumed the liquid and vapor phases were incompressible and included the effects of surface tension, vaporization and interface movement The density differences between the liquid and vapor were also included. The small amplitude stability analysis indicated that surface 80 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . tension is a highly stabilizing force but this could be overcome with vaporization and large density differences between the vapor and liquid states. The initial temperature profiles o f the liquid and vapor states both transferred heat to the interface as shown in figure 1.2. However to simplify the stability analysis Miller assumed no heat transport in the vapor phase and assumed the interfacial temperature remained constant. VAPOR Evaporation rate = const V . v.= o T „v sz LIQUID I , v s 2 w. z Millers Model Figure 1.2. Physical Setup of Millers model for stability analysis Palmer used a mathematical formulation similar to Miller’s and described a phenomenon called “vapor recoil” which caused a significant increase in the amount of heat transferred and a violent eruption o f the liquid surface. He also developed his stability analysis to show the effects of “vapor recoil”, surface tension, fluid inertia, viscous dissipation, and the moving boundary mechanism. Palmer analyzed the hydrodynamic stability of a steady evaporating, pure liquid at reduced pressure using linear stability analysis, shown in figure 1.3. His results revealed that a rapidly evaporating liquid is unstable to local variations in evaporation rate, caused by local surface depressions produced by the force 81 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . exerted on the surface by the rapidly departing vapor (“differential vapor recoil”) and the resulting shear exerted on the liquid surface by the vapor. The coupling of the “differential vapor recoil” mechanism with the Marangoni effect was studied as well as the relative importance on stability o f the moving boundary, fluid inertia and viscous dissipation at the interface. Palmers Base State Interface » fixed with respect to the coordinate system iw .z No heat transfer Is assum ed In vapor n* = evaporation rate = constant iI v.v= o Tg a t Interface Flat interfece u .x linear tem perature profile in the BL constant in the bulk liquid V.VL=0 general up flowing of fluid to m atch evaporation rate g* Figure 1.3. Physical Setup of Palmer’s Base State model for stabili^ Palmer notes that previous work with regards to spontaneous convection in fluid-fluid systems had mostly looked at two mechanisms, surface-tension variations and density stratification, as to their effects as destabilizing mechanisms at the fluid-fluid interfaces. He indicates that mass-transfer rates can be increased by more than threefold with the onset of spontaneous convection. However a third mechanism “differential vapor recoil,” first noted and described by Hickman (1952 ), had shown an increase in evaporation rates as high as 20- 82 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . fold. Hickman performed a number of experiments involving rapid evaporation of liquids exposed to high vacuum and being heated from below. Palmer’s model, figure 1.3, assumed a pure liquid of infinite lateral extent and unbounded from below undergoing steady evaporation at reduced pressure. His model also described a liquid heated from below and ignored the heat transfer away fium the interface similar to Miller’s model. Palmer performed the linear stability analysis by imparting to a base state a small disturbance, exponential in time and periodic in space with arbitrary wavelength. The disturbance is made to statisfy the conservation of mass, momentum and energy subject to the boundary conditions at the fluid interface. His characteristic equation related the Hickman number to the wave number and the other dimensionless groups at the condition of neutral, stationary stability. Palmer’s stability model yielded stability boundaries of critical Hickman number (Nh c ) versus wavenumber. Palmer found three regions of instability: Region 1 - due to the moving boundary mechanism. Region II- combined effects of “vapor recoil” and fluid inertia, and Region III- viscous dissipation. Palmer and Maheshri (1981) performed experiments to quantify the enhanced heat transfer coefficients obtained during the convection caused by the “vapor recoil” mechanism of instability. Palmer and Maheshri found that the rates of heat transfer were enhanced on the order of 20 over the still values. Gillion, et al (1988) performed experiments on determining quantitatively the rate of evaporation under controlled pressures and the rate of microwave energy absorbed of ethyl alcohol. The experiment foimd three regions of evaporation rates. They foimd that the two transition points between the regions correspond to the onset of convection by the buoyancy effect and the instability of the thermal boundary layer by a coupled Marangoni and vapor recoil instability. 83 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . Burelbach, et al. (1988) performed a stability analysis on a physical configuration which differed from Miller and Palmer. Burelbach’s configuration was a thin liquid film on a horizontal heated surface. They were able to develop a partial differential equation (pde), by considering long wave disturbances which described the evolution of the interface h(x,t). This pde developed the actual profiles of the height of the interface while undergoing evaporation up to the point of film rupture as a function of time. A numerical model was developed to solve the strongly nonlinear partial differential equation which resulted. It was solved in conservative form, as part of an initial value problem for spatially periodic solutions on the fixed interval 0<x< 2x/k. Centered differences in space and the Crank-Nicholson routine was used to obtain second-order accuracy in space and time. The difference equations were solved by Newton-Raphson iteration. Numerical solutions of film profiles at various times for an isothermal case were developed. The effect due to sudden depressurization and the phenomena associated with superheated liquids and their resulting turning to froth, vapor with entrained fluid, is discussed in Shepard and Sturtevant (1982), Prosperetti and Plesett (1984) and Higuera (1987). The effect as far as mass surge (with entrained liquid) is different but the initial stages of instability are the same as developed in Palmer (1976). The present study will perform the formulation of the dispersion relation for the stability o f an interface undergoing heat and mass transfer and will endeavor to apply the result to the conditions of a fluid interface in the parameter space describing propellant densification. The formulation will follow Palmer’s (1976) derivation in order to provide a method to calibrate the results with respect to “vapor recoil”. An investigation o f the regime of the conditions germane to propellant densification will be conducted to see if “vapor recoil” instabilities are 84 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . possible. If the regime of propellant densification does not fall within the “vapor recoil” region, then the phase out of the “vapor recoil” mechanism will be investigated. An evaluation of the resulting dispersion relation to oscillatory instabilities will also be made. 85 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 2. Formulation We will use a simplified model to analyze the effects of heat and mass transfer on the stability of an evaporating interface. The model similar to that of Palmer (1976) has a two layer region divided by a horizontal interface, with both layers consisting of the same fluid in different states (shown in figure 1.3). The interface in the base state is flat with a constant evaporation rate at the interface and a linear temperature profile in the liquid layer underneath the interface down to a depth 5 in the liquid. For -6 >z> -oo we have a constant temperature Tuq. At z = -Ô we have a constant up-welling velocity of fluid to provide the mass loss at the interface due to the evaporation rate q * and to maintain the equilibrium interface at the fixed level z = 0. The surface temperature of the interface is T,. The equations of motion of the lower liquid layer for -6 < z < 0 in perturbed variables are: V F = 0, (2.1) DV — P, — = -V p + p , a % + p V 'F . (2.2) We have assumed that p is constant and the relation that the base state density pe satisfies the equation of state (consistent with a Boussinesq approximation) ~ (2-3) where a = — — ) is the coefficient o f thermal expansion and pV dT), P = - P lo^T. (2.4) The energy equation is (2.5) where we have assumed no shaft work, dissipation or compressive heating. We will need to impose a volumetric heating rate to support the base state linear 86 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . temperature profile (or the base state convected part of DT/Dt) in the lower thermal boundary layer. This is not uncommon in stability theory where body forces are added to the equations of motion in order that the model, base state profiles are exact solutions of the governing equations. The perturbations to the chosen base state are then described by free (or homogeneous) solutions to the equations of motion. Taking equation 2.5 and subtracting the resulting base state equation and linearizing we obtain, in perturbed variables, n r + pw = K V ^r, (2.6) Dt where k = is the thermal diffusivity, and p = — — is the linear pCp 6 temperature gradient in the boundary layer just below the interface. The equations in the upper gas phase layer (assuming p = constant, no heat transfer or internal heat generation, and incompressibility) are V F = 0, (2.7) P, ^ = -V p + gp, + . (2.8) Subtracting the hydrostatic base state from equation 2.8 yields n v - P , ^ = -V p + p ,V ^ F . (2.9) Hence the equations describing the perturbed fields in the gas layer are 2.7 and 2.9. With mass flow at the interface, the matching relations to be used at the interface will include continuity of mass flux, normal momentum, and tangential momentum with a no-slip condition (i.e., continuity of tangential velocity). We will also require that the jump in energy across the interface be equal to the latent heat of vaporization. 87 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . The linearized continuity of mass equation across the interface in perturbed variables is p,w,-p,iv, = ( p , - p j ^ , (2.10) where w,, vP , are the perturbed vertical velocities in the liquid and vapor regions immediately above and below the interface and — is the vertical velocity of the dt interface (with ^ denoting the displacement above the base state). The linearized normal momentum jump across the interface in perturbed ... . Pz - A + 2 q 'n (p ;' - p [ ') + 2[p,w„ - P 2W 2J variables is (2.11) +&(pi =0. where q* is the base state mass evaporation rate per reference area and the subscripts 1 and 2 refer to the liquid and vapor state respectively. The second term in (2.11) is the change in normal momentum flux arising from the density change due to evaporation. This term is sometimes referred to as the “vapor recoil” effect. One of the goals of the present effort is to assess the role of this term on the instability of an evaporating interface. The parameter a is the coefficient of interfacial tension at the evaporating interface. The linearized equation describing the continuity of the tangential momentum across the interface, written in terms of perturbed variables, and using the no-slip condition at the interface, is dx = 0. ( 2. 12) The third term, the Marangoni effect, defines the influence of a temperature- dependent coefficient of surface tension on the tangential stress. We assume a linear dependence of c t on temperature so that the parameter is constant. dT The two terms in the bracket come from, respectively, horizontal gradients in the 88 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . perturbed liquid temperature and the distortion of the base state temperature gradient by fiee surface motions. The linearized energy jump relation written in perturbed variables is given as df. + 2 { v iW ,w ,) ^ - k ,^ dz = 0, (2.13) In writing this expression we have used Wb to denote the liquid base state upwelling vertical velocity, and we have used the assumption that kinetic energy is much smaller than the enthalpy change due to evaporation. Hence, energy flux terms are neglected in comparison with the evaporative flux which involves the latent heat of vaporization X ev ap . A relation for the perturbed evaporation rate rf is obtained by using the DC, definition q = p,(w , — w j , and the kinematic condition w, = for the liquid- side interface velocity. Linearizing about the base state yields the relation - dx\ (f,+ P Ç ) = p . [ i P . - ^ j , (2.14) For scaling we will utilize the thermal boundary layer depth 6 for a length scale, 5^/kl for time, kl/5 for velocity, -pS for temperature and plKl/S^ for pressure. This scaling results in the non-dimensional perturbed lower liquid layer equations of motion being: | + f = 0. (2.15) — + — - w = V"e. (2.18) dt dz 89 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . The non-dimensional groups are given by Pr^ = — , Ra Wg = Re^ Pr^ and Re^ = The scaled non-dimensional perturbed equations of motion in the upper vapor region are given as; — + — = 0, (2.19) dx dz 1 du Pr^ dt 1 dw + Re^ A ^ p — = - ^ + 4 - V 'w . (2.20) (2.21) Pr^ ar " " az dz N, The additional non-dimensional groups are the liquid-to-vapor density and kinematic viscosity ratios N = — , and = — . Pv V, Performing the curl of equations 2.20 and 2.21 and eliminating p and u in favor of w results in a single equation for the vertical velocity (2.22) Representing the vertical velocity and pressure in terms of normal modes defined by (w,p) = (»'(z),/>(2 ))e"‘' ^ ’, the fundamental equation (2.22) becomes an ODE for the vertical velocity structure function Wy(z) for the vapor region (D ^-k^) D "-R e ^ - i(ùN^ ^ Pri J W = 0 . (2.23) In the same way, using the z-momentum equation and the normal mode representation we obtain a relation linking the perturbation pressure field to the vertical velocity 90 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . J_ W y-D P y= 0. (2.24) In like manner, performing the curl of equations 2.16 and 2.17, assuming that buoyancy effects in the lower liquid layer are negligible (Ra = 0), and again eliminating p and u in favor of w we obtain f ^ T > 5 1 a Y - Re, + 1 1 dz dt)^ = 0. (2.25) We specifically ignore any effects of buoyancy-induced instability within the liquid layer. The principal focus is on surface-driven instabilities in the presence of evaporation. Using a similar normal mode decomposition for the liquid layer as for the gas layer, we obtain the fimdamental structure equation for the liquid layer D ^ -k ^ - W ,= 0 . (2.26) Taking the z-momentum equation (2.17) with (Ra=0) leads to a relation for the liquid-side pressure perturbation in terms of the liquid-side vertical velocity -R e D -(‘■ -a W^=DP,. (2.27) Now looking at the energy equation 2.18, using normal mode representation for 0 and W, with = RePr and 0 =0(z)e‘ ^ '“'“‘*, we obtain the following equations for the two liquid regions (i.e., in the liquid boundary layer and in the liquid bulk): [£)"-R eP rZ )-(A :^-/G ))]0= -lF i 0 > z > - l (2.28) [l> ^ -R eP r£ > -(;t^ -i© )]0 = O -l>z> -oo (2.29) The scaled interface and matching jump conditions, which are applied at the equilibrium interface level z = 0 for the linearized problem, are as follows: Mass flux: ~ ^ » (2.30) 91 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm ission . Kinematic; Continuity of tangential velocity: Wki - W u - R e P r(iV p - 1 ) ; „ = 0 , (2.31) (2.32) Normal momentum: ^ ( 9 -Ç ) +2Af, Cr f p (2.33) Tangential momentum: N.. (d^ a n r a ' a n Lar: dz^j W j, - .a x ' a z 'j Energy: - 2 + 2 Ar„w,, + Pr9 ^. = 0 Re = (2.34) (2.35) The subscripts 1 and 2 are for the liquid (lower) and vapor (upper) regions respectively. Also at z = -1 , the interface between the thermal boundary layer and the bulk on the liquid side, & l and DQl must be continuous. The additional nondimensional groups appearing in the scaled problem include the Hickman number dT p-0 (2.36) which measures the vapor recoil effect relative to the restoring effect of surface tension; the Crispation number o 5 (2.37) measuring the ratio of normal viscous stresses to interfacial tension; the Bond number S V P i - P x ) (2.38) 92 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . measuring the ratio of gravitational forces (body forces) to that o f interfacial tension; the Marangoni number which measures the lateral variations in interfacial tension to the tangential viscous stresses; and the Brinkman number (2.39) (2.40) which measures the effect of viscous dissipation relative to heat diffusion. Solutions for the vertical velocity in the upper and lower layers take the form, r ( r ) = 2 A « ‘" . (2.41) /•t where /=1,2 for the upper and lower layers respectively. Using the conditions that fV(z) is bounded as z — > ± o o gives the two equations for the vertical velocities as: Wy (r) = ^,2^"*^ + , 0 < z < oo and ^ ( z ) = ^ 21^*^ + for 0 > z > -<». The coefficients X 23 and X 14 are given as, Re+ (2.42) (2.43) X23 — V Re"+ 4 PrJ ReA^,- and =• (2.44) (2.45) In like manner, the homogeneous solutions for 9 take the form in the thermal boundary layer and lower region as, (2.46) where 1= 1 or 2 for the lower and upper boundaries of the thermal boundary layer respectively. For the two regions the homogeneous solutions are given as: 93 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . 0, = , for - 1 > z > -oo, (2.47) and 0 2 = + 5 3 3 6 * ^ * * , for 0 > z > - l . (2.48) The roots o f the characteristic equation are RePr± J(RePr)^ +4(t^ -/co) = ---------- ^ ------ • (2.49) Applying the matching conditions at z = -1 of 0, = 0 2 and D 0, = DQ j , (2.50) and the condition at z = 0 of equation (2.28), we are able to develop the complete solutions for 0 in terms ofthe unknown coefficients o f W(z). The solutions are 0, = for - 1 > z > - 0 0 , (2.51) and ........................ (2.52) for 0 > z > -I. With the representations for W given in (2.42) and (2.43) plus those for 0 given in (2.51) and (2.52), we are able to derive expressions for the matching conditions at the iiquid-vapor interface in normal mode variables. Using the normal mode representation for ^ as, ^ ^ ^ resulting matching condition equations become; Mass flux: W^^-Wy = (-/© q , (2.53) Kinematic: ^ 0 ^ 0 ) = C 0 > (2.54) Continuity of tangential velocity: N^DWy - DW^ + RePr(jVp - l)t%g = 0, (2.55) Normal momentum: 94 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . +2N. Cr 5 N.. DWy-DWj^ (2.56) Tangential momentum: N. (D ^+ e]W , ~ { D ‘ + * ’K - ;o ) = 0. (2.57) Energy: ^ - 2 - ^ N . , D W y +PrD 0i =0, (2.58) where we have used X«vap» K.E. and neglected the kinetic energy term. Solving for the pressure in the vapor and liquid layers using equations (2.24) and (2.27) respectively, we obtain relations for Pl and Pv in terms of the unknown coefiRcients for W. These relations are 7’ .= - [ R e - ^ ] . 4 „ e * = , (2.59) (2.60) Using the solutions for Wl, Wy, 0 ^ (or 8%), Pc and Pv from equations 2.42,2.43, 2.52,2.59 and 2.60, and substituting into the matching conditions of (2.53-58) yields the following set of linear equations for the unknown coefficients. Mass flux: Kinematic: •^p {Ai2 + An ) Kq +(Ko+/(o);o=0, V ^21 ^23 (2.61) (2.62) Continuity of tangential velocity: -N,{kA„ +r,A,,)-{kA„ +r,A^) + RcPl{N, = 0, (2.63) 95 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Tangential momentum; k^rl + ^ n - ; o ) = 0 V «23 ; Normal momentum nentum: —rp(2^^-^i2 + 'k )'^u) + -^ 2/ 2 *^ ^ 1 » V «21 y . . r,.2_2 . .r . 2 / . . „ (2.64) 2Æ I ,o « , . L C r 5 iV .. * 1 4 ^Bo 1 f2 + 2N çi,ri^ W 2 3 ■ * ■ “ Pr ^ 2. (Bj, + Bj2 ) (2.65) Energy: - ^ B r(^ i 2 + n-A* ) + P r ^ 2i + (Re P r- q)B^ ) L + 2JV „i + M + ( ^ ^ + 2Ar„r, + ^ = 0. The various parameters appearing in these algebraic equations are given by (^ 2, = Re Pr A : - /©, = r ^ - Re(Pr-1), „ N h K N , ° RePrA^cr(l-^p~‘) ’ (2.66) 'V =^14 = 2 -ReA^„ V '■ /. = ^23 = ^^R e+ ^Re^+4^k^ g = X, = l^ R e P r+ ^(RePr)^ + 4 (t^ -rcù)j. 96 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . The conditions requiring the continuity of 9 ^ and £70^, at z = -1 yield the set of equations: - B,,e~'^ + ^ 4 , + — =0, and (2.67) «21 «23 B ^,q e- + 522(RePr-g)e-<“'*^-’> - B ,,q e'‘ ' e”* e~'‘ - (2.68) +k — A 2 1 + fi — = 0. « 2 1 «23 The dispersion relation defining the instability characteristics of the evaporating interface is now obtained by forming the determinant of the eight unknown coefBcients given in the eight homogeneous equations (2.61-68). The dispersion relation is given by the determinant 9 = det [D] =0. The elements of the determinant are given in table 2.1. For the temporal instability problem this characteristic equation defines an eigenvalue problem for the complex frequency (0 for a given wave number k and fixed values of the parameters N cr, , N , N , Re, Pr ). Clearly, the parameter space of the problem is quite large and we are practically limited to exploring instability trends for only a limited number of the physical processes represented by this system. 97 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . + I -*e I ::s I % n I -*e + « ^ - y k a o o ? « 2 c n I • o I I C S I S ; 2 2 * 3 3 > % I >» I J T f 5 S ; "* e I ? -; « a : fS a ; $ a: C M a* I 6 e c Q- C Ü £ £ K. fS + a: % a: f N T S o Z 6 c a: £ i X Q * I & § • .§ I -»e Ô S a a : I ■ * « (N + -Q 5 a- < N I o c a : s a: f N ■ Ç J C Ü c a: H h .3 1 -*e £ a C L -AC 6 < N + w oi 0 0 C m o 1 . —• - :a <N V ; a: I c - aT 0 a: « N ■ * e - » e ? 5 a; C M c è a: [a: a; L T < N a |a tN 98 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 3. Computational Results There are several assumptions implicit in the dispersion relation given in table 2.1. These assumptions include the neglect of buoyancy effects (i.e., Ra=0) and that the kinetic energy o f the vapor is much less than the latent heat of vaporization. However, the formulation explicitly includes the possibility of oscillatory modes of instability (i.e., c o is complex and m = c o ^ on the neutral curve where o o , = 0). The initial eigenvalue searches further reduced the matrix elements by searching exclusively for stationary neutral modes (i.e., co = 0). In order to validate the eigenvalue search against some o f the results of Palmer (1976) we set = 0, Nb, = 0 and chose Np = 10'°. The neutral curve in the Hickman number- wave number plane for this condition is shown in figure 3.1a. As can be seen the critical value of Hickman number is = 6x10"*. This neutral curve agrees with Palmer’s (1976) results for the ‘differential vapor recoil’ region of instability. Figure 3.1b is the same neutral curve depicted in the manner that the information will be displayed for the rest of the results. The investigation of the mapping of critical Hickman number as a function of the Re number for various density ratios of 10’ to 10'° is shown in figure 3.2. Although these density ratios are much higher than is realized in many physical applications, this is the only range of density ratios where Palmer presents numerical results. The stability boundaries shown in figure 3.2 agree very well with those given by Palmer (1976). Palmer, however, provided no information regarding the wave numbers at which the occurs. The wave number associated with the vertical asymptotes of the Np = 10’ neutral curve vary approximately from a value of a = 1.4 to 2.4. 99 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ith out p erm issio n . 0 . 0 0 1 0. 0001- I E 0.00001 3 c I I 0. 000001- 0 .0000001-1 0.00000001 100 0.1 1 10 \Nave number 0.001 -3 0.00011 I E 3 C 0.00001 - c 0 1 X 0.0000011 0.0000001 0.00000001 10 20 30 40 50 60 70 80 90 100 0 W iawe number Figure 3.1. (a) Neutral curve of Hickman number as a function of wave number showing logarithmic behavior at large wave numbers, (b) Neutral curve showing wave number behavior at the critical Hickman number. A depiction o f the wave number at which the critical Hickman number occurred is shown in figure 3.3 for N^= 10’ through 10'°. It can be seen that, similar to the Np = 10’ case, the wave numbers for the other Np cases do not vary greatly at the different ends of the neutral curve where the instability becomes independent of Hickman number. A plot of the critical Reynolds number versus 100 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . wave number for the cases shown in figure 3.3 is given in figure 3.4. As can be seen in figure 3.2 and 3.4, the critical Hickman number only occurs over a bounded band of Re for each Np. In the case of = 10*, the Re band is approximately 10'^ to 10"*. 0.001 I 0. 0001- 0.00001 X 0.000001 8 O 0.0000001- 0 .00000001- 1 Np=1(K7 10^ i I lO'O 1010 0.01 0.1 0.000001 0.00001 0.0001 0.001 Reynolds number Figure 3.2. Critical Hickman number variation with Reynolds number for various density ratios (N =100, N c=10'\ Nb o =1, Pr=10, Br=Ma=0). 0.001 -3 I 8 Ü 0.0001 0.00001 ■ 0 .000001- 0.0000001- 10^ 10^ 1010 Np: 10«7 0.00000001 0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9 2.1 2.3 2.5 Wave number Figure 3.3. Critical Hickman number variation with wave number for various density ratios (N =100, Ncr=10'*, Nb ^,= 1, Pr=10, Br=Ma=0). 101 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . I E 3 C 2.5 X C D g 1.5 0*6 I « I 0.5 3 C 0.000001 0.00001 0.0001 0.001 0.01 0.1 Critical Reynolds number Figure 3.4. Variation of the wave number where the critical Re or Hi occurred for various density ratios (N^=100, Ncr=10'*, Ngg=I, Pr=10, Br=Ma=0). The region where vapor recoil becomes inoperative is shown by viewing Palmer’s (1976) work (Palmer’s figure 3) and focusing on the region to the right of the curves for Nr e> 10'‘ and for density ratios of Np <10* (for Nc^=10'*, Ngg=I, Npr=10, Np=10^, Ng =Nm ,=0). The morphology of the neutral curves, as vapor recoil phases out, can be seen by viewing the following sequence of neutral curves near the vapor recoil fade out boundary. A depiction of the neutral curve for Nr e = 10‘ ^ and Np = 5x10* is shown in figure 3.5. Notice that the critical Hickman number is approximately 1.3x10'^ and the span of unstable wave numbers at a Hickman number that is one order of magnitude higher than critical is approximately 16. The region to the right of these boundaries, i.e. Nr e> 10'‘ and Np < 10* is where vapor recoil ceases to be operative as an instability mechanism and this fade out of its effect is shown in the neutral curves for the cases of Np varying from 2.5 to 2.06x10* in figure 3.5. In figure 3.5, the density ratio is lowered to 2.5x10* and it can be seen the critical Hickman number increases to approximately 3.x 10'^ and the span of unstable wave numbers at a decade above the critical Hickman number is still approximately 17 or 18. However, when the 102 R ep ro d u ced with p erm issio n o f th e cop yright ow ner. Further reproduction prohibited w ithout p erm issio n . density ratio is decreased to 2 .2 5 x 1 and then to 2.06x10 \ both neutral curve boundaries of the low and high wave numbers asymptote to infinity away Srom the critical Hickman number. Also, the unstable wave number range, at the Hickman number of a decade higher than critical, decreases to a very small range o f 2 for the density ratio of 2.06x10®, just before the vapor recoil instability fades out. As the density ratio decreases, the neutral curves become nonexistent. An appropriate change of any parameter can cause the neutral curve to reappear as is shown in figure 3.6 for a change in the dynamic viscosity ratio from 10^ to 20 at the critical Np of 2.06x10®. 100-g 2.06X s E 3 C 2.5x c ( 0 I 0.01- 0.001 60 70 40 50 80 90 100 0 20 30 10 Wave number Figure 3.5. Neutral curves of Hickman number as a function of wave number for various density ratios illustrating the vapor recoil phase out (Np=100, Ncr=10*® , Nb o = 1, Pr=10, Br=Ma=0). At the critical Np of 2.06x10®, with the dynamic viscosity ratio equal to 20, the region of vapor recoil can again be made to phase out by varying the density ratio. This is shown in figure 3.6 by varying the density ratio from 2.06x10® to 2.6x10® where phase out again occurs. 103 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 2.6x11'S 3x10 4x10: Ë E 3 C C c o I X 0.01-= 0.0001 10 20 30 40 50 60 70 80 90 100 0 Wave number Figure 3.6. Neutral curves for Hickman number as a function of wave number for various density ratios illustrating vapor recoil phase out with a different (N^=20, Ncr=10 \ Nb „=1, Pr=10, Br=Ma=0). The effect of varying the dynamic viscosity ratio from 10 to 100 is shown in figure 3.7 for the case with Np= 2.06x10®. It can be seen that the phase out of vapor recoil instability occurs at = 100 as we obtained previously. Using a higher density ratio of Np= 5x10® and a viscosity ratio of = 100, the effect of Prandtl number on vapor recoil phase out is shown in figure 3.8, with Prandtl number varying from 100 down to 4.15 where phase out occurred. 104 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 100-3 80 Nti=100 I E i 0.1 0.011 0.001 0.0001 0.00001 30 40 50 60 70 0 10 80 20 Wave number Figure 3.7. Neutral curves with Hickman number versus wave number with various dynamic viscosity ratios (NR,=IO'\ Ncr=10‘* , Ngg=l, Pr=10, Br=Ma=0, N =2.06xl0‘). : 0., I 0.01 4.5 100 0.001 - = 0.0001 10 20 30 40 50 60 70 80 90 loO 0 Wave number Figure 3.8. Neutral curves with Hickman number versus wave number for various Prandtl numbers (NR,=10'\ Ncr=10'\ Nb o = 1, N^=100, Br=Ma=0, Np=5xl0*). The variation of the neutral curve at a constant Hickman number and varying Reynolds number as a function o f wave number is shown in figure 3.9 for two different density ratios (10* and 10’ respectively). The chosen constant Hickman numbers are very near the critical Hickman numbers for the density 105 R ep ro d u ced with p erm issio n o f th e copyrigh t ow ner. Further reproduction prohibited w ithout p erm issio n . ratios selected. The variation o f the Reynolds number as a function of wave number has a smoother variation than the Hickman number versus wave number. Eliminating the effect of the variations in evaporation rate due to surface temperature excursions is illustrated by setting Hickman number to zero. This eliminates vapor recoil as a mode of instability. Figures 3.10 and 3.11 illustrate the neutral curves of Reynolds number as a function of wave number for the no vapor recoil case at a density ratio of 1 O'. It can be compared to figure 3.9 which includes vapor recoil. It can be seen that the wave number of the critical parameter shifts to a wave number of unity (firom a wave number of 1.25) for the case without vapor recoil. Figure 3.11 illustrates that the Reynolds number neutral curve, for Hickman number equal to zero, as a function of wave number appears similar to the neutral curve of Hickman number versus wave number as shown in figure 3.1b and takes on the role of the critical instability parameter. N )=10'«. Hi = 9 x 1 0 ^ 0.01- ë O c = 0.001^ I 0.0001 - 0.00001 2 2.5 1.5 0 1 0.5 Wave number Figure 3.9. Neutral curves of Reynolds number versus wave number near the critical Hickman number (Ncr=10 * , Pr =10, Nb o =1, N =100, Br=Ma=0). 106 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . I E 3 C « 0.1 " 5 I 0.01 1 0.1 0.01 10 Wave number Figure 3.10. Variation of Reynolds number with wave number on the neutral curve for the zero Hickman number case (Np=10*, Pr =10, Nb o = 1 , N^=100, Br=Ma=0, Nc,=10'*). I E 3 C « 0.1 I IT 0.01 70 80 90 100 10 20 30 40 50 60 0 V\6ve number Figure 3.11. Variation of Reynolds number with wave number for the neutral curve with Hi=0, linear wave number depiction (Np=10*, Pr =10, Ngp=I, N^=100, Br=Ma=0, Ncr=10'*). 107 R ep ro d u ced with p erm issio n of th e cop yright ow ner. Further reproduction prohibited w ith out p erm issio n . The regime o f interest that is associated with hydrogen vaporization does not occur at a pressure in sub-torr ranges, therefore the fluid properties and non- dimensional ratios do not occur in the range of values that exhibit vapor recoil instabilities. They can, however, cross over into the regime of the fluid inertia instability. A realistic set of values of the non-dimensional parameters for the hydrogen vaporization case of interest are: 10^ <N^< 10* \< N ^ ^ < 10* «20 10 " < Ncr < 10 * 10-' < A T g , <10*, ^PR = » U ^Br = A ^ A fa = 0 Selected values o f the above parameters were made to see if any common regions of vapor recoil or fluid inertia could be found. They are given as: Np = 10*, = 20, Pr = 1 , Ncr = 10 *, = 1 , = I . These values, which were selected to favor vapor recoil instabilities, did not overlap within the vapor recoil regime. There is some degree of overlap for the fluid inertial region however. In order to achieve an overlap with vapor recoil we would need to have a decrease in the Nr ^ by 10 * and an increase in the density ratio Np by 10^. It should be observed that this could be achieved with a low volatility mineral oil, of the type that is used in the lubrication of vacuum pumps, and decreasing the pressure above the evaporating liquid to sub-torr ranges (10 * to 10'*) as described by Hickman (1952) and Palmer (1976). For liquid hydrogen the bulk vapor pressure is so high that a substantial decrease in presswe produces bulk boiling and an accompanying, as the liquid vapor rises in the bulk, disruption of the surface that would mask any of the instability mechanisms such as vapor recoil or fluid inertia. 108 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . These results have corroborated Palmer’s (1976) results for the ^differential vapor recoil’ region o f stationary instability. The other regions were not specifically evaluated, such as the moving boundary, fluid inertia or viscous dissipation. Since the dispersion relation was developed for a general, complex fiiequency o , it was desired to search for oscillatory modes of instability were Re(a>) # 0. A numerical search for oscillatory modes proved futile, however, in the regime where vapor recoil effects are significant. This search was not exhaustive in regard to the whole of the extensive parameter space and does not rule out oscillatory modes existing in regimes where the vapor recoil effect is unimportant. 109 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . 4. Conclusions The neutral stability curves for an interface were developed for a fluid undergoing heat and mass transfer. The analytic model was developed similar to the model developed by Palmer (1976). The neutral stability curves were developed for the instability region termed “vapor recoil.” Other mechanisms of instability, including the moving boundary, fluid inertia and viscous dissipation, were not examined but were described in Palmer (1976) and Miller (1973). It should be noted that Hickman (1952,1972) was the first to note “vapor recoil” as a way to induce hydrodynamic instability for rapidly evaporating fluids in low pressure, and this can be a significant instability mechanism in certain circumstances. The regime of interest that is associated with hydrogen vaporization does not occur at a pressure in sub-torr ranges, therefore the fluid properties and non- dimensional ratios do not occur in the range of values that exhibit vapor recoil instabilities. They can, however, cross over into the regime of the fluid inertia instability. A realistic set of values for the non-dimensional parameters for the hydrogen vaporization case of interest are given as: = 10^, = 20, Pr = 1, = 1 0*. A /^ * , = 1 , A ^ r c ^ — 1. These values, which were selected to favor vapor recoil instabilities, did not overlap within the vapor recoil regime. There is some degree of overlap for the fluid inertial region however. In order to achieve an overlap with vapor recoil we would need to have a decrease in the Nr ^ by 10 ^ and an increase in the density ratio Np by 10^. Linear instability studies by Prosperetti and Plesset (1984) neglected viscosity and surface tension gradients and found an instability mechanism similar to vapor recoil. They obtained the perturbation growth rate as a function of wave number and found where oscillatory instabilities might occur. For water at 100° C, they predicted oscillatory instabilities might start at wave lengths of 0.014 cm 110 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . and smaller. At this wave length however, viscosity and surface tension should be very dominant players in the damping of instability and these effects were specifically neglected in their derivations. Additional linear instability analysis was performed by Higuera (1987) including the effects o f viscosity and he obtained a growth rate versus wave number relation similar to Prosperetti and Plesset (1984). His results indicate that, to first order, gravity and surface tension must be large enough to get oscillatory instabilities. This leads to the opposite view that if gravity and surface tension are not large enough, then the growth rate should be real for all wave numbers. If this were the case, then it would lend credence to the stationary instability assumption made by Palmer (1976) and this present work. An additional study performed by Burelbach, Bankoff and Davis (1988) found vapor recoil to be an important effect in the rupture o f thin films. They developed a surface height evolution partial differential equation using long wave theory. The numerical solution of this equation allowed them to determine the effect of mass loss, vapor recoil, cooling due to the loss of latent heat, and van der Waals attraction forces on film rupture. This present study developed the dispersion relation with explicit inclusion of complex o o to be able to search for the eigenvalues of o) that could delineate the regions o f oscillatory instability. For the areas of the growth rate as a function of the wave number space searched, no oscillatory eigenvalues were found. This doesn’t prove that they don’t exist, just that none were discovered. I l l R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . References 1. Babchin, A. J., Frenkel, A. L., Levich, B. G., Sivanshinsky, G. I., 1983, Nonlinear saturation of Rayleigh-Taylor instability in thin films, Phys. Fluids, Vol. 26, No. 11, pp. 3159-3161. 2. Batchelor, G. K. and Nitsche, J. M., 1991, Instability o f stationary unbounded stratified fluid, J. Fluid Meek, Vol. 227, pp. 357-391. 3. Batchelor, G. K. and Nitsche, J. M., 1993, Instability o f stratified fluid in a vertical cylinder, J. Fluid M eek, Vol. 252, pp. 419-448. 4. Burelbach, J.P., Bankoff, S.G., and Davis, S.H., 1988, Nonlinear Stability of Evaporating/Condensing Liquid Films, J. Fluid M eek, Vol. 195, pp. 463- 494. 5. Chandrasekhar, S., 1961, Hydrodynamic and Hydromagnetic stability, Oxford University Press. 6. Delhaye, J.M., 1974, Jump Conditions and Entropy Sources in Two-Phase Systems, Local Instant Formulation, Int. Journal M ultiphase Flow, Vol. 1, Pergamon Press, pp. 395-409. 7. Dongarra, J. J., Bunch, J. R., Moler, C. B., Stewart, G. W., 1979, LINPACK Users’ Guide, Philadelphia: S. I. A. M. 8. Emmons, H. W., Chang, C. T. and Watson, B. C., 1960, Taylor instability of finite surface waves, J. Fluid M eek, Vol 7, pp 177-193. 9. Gear, C. W., 1971, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall Inc., New Jersey. 10. Gillon, P., Courville, P., Steinchen-Sanfeld, A., Bertrand, G., and Lallemant, M., 1988, Free Convection and Boundary Layer Instabilities in Evaporating Liquids Under Low Pressure and/or Microwave Irradiation, Pf^sieoChemical Hydrodynamies, Vol. 10, No. 2, pp. 149-163. 11. Hickman, K. C. D., August 1952, Studies in High Vacuum Evaporation. Surface Behavior in the Pot Still, Industrial and Engineering Chemistry, Vol. 44, No. 8, pp. 1892-1902. 12. Hickman, K. C. D., 1972, J. Vacuum Science T eek, Vol. 9, p. 960. 112 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 13. Higuera, F. J . , 1987, The hydrodynamic stability of an evaporating liquid, Phys o f Fluids 30(3), March, pp. 679-686. 14. Jain, R. K., and Ruckenstein, E., 1976, Stability of stagnant viscous films on a solid surface, J. o f Colloid and Interface Science, Vol. 54, No 1, pp. 108- 116. 15. Levich, V.G., Krylov, V S., 1969, Surface Tension Driven Phenomena, Annual Review o f Fluid Mechanics, 1, pp. 293-316. 16. Maa, J. R., 1967, Evaporation coefficient of liquids, In d Eng. Chem. Fundamentals, Vol. 6, No. 4, pp. 504-518. 17. Mann, D. B., Sindt, C. P., Ludtke, P.R., and Chelton, D.B. 1966, Liquid- Solid Mixtures of Hydrogen near the Triple Point, Advances in Cryogenic Engineering. Vol. 11. Plenum Press, N.Y., pp. 207-217. 18. Maxworthy, T., 1998, The Dynamics of Sedimenting Surface Gravity Currents, (submitted to J. Fluid Mech.). 19. Menikoff, R., Mjolsness, R. C., Sharp, D. H., Zemach, C., 1977, Unstable normal mode for Ray leigh-Tay lor instability in viscous fluids. Physics o f Fluids, Vol. 20, No. 12, pp. 2000-2004. 20. Menikoff, R., Mjolsness, R. C., Sharp, D. H., Zemach, C., Doyle, B. J., 1978, Initial Value problem for Rayleigh-Taylor instability of viscous fluids. Physics o f Fluids, Vo\. 12, No 10, pp. 1674-1687. 21. Miller, Clarence A. 1973, Stability of Moving Surfaces in Fluid Systems with Heat and Mass Transport, AIChE, Vol. 19, No. 5, pp. 909-915. 22. Newhouse, L. A. and Pozrikidis, C., 1990, The Rayleigh-Taylor instability of a viscous liquid layer resting on a plane wall, J. F luid Mech., Vol. 217, pp. 615-638. 23. Palmer, H. J. and Berg, J. C., 1972, Hydrodynamic stability of surfactant solutions heated from below, J. Fluid Mech., Vol. 51, Part 2, pp. 385-402. 24. Palmer, Harvey J. 1976,The hydrodynamic stability of rapidly evaporating liquids at reduced pressure, J. Fluid Mech., Vol. 75, part 3 ., pp. 487-511, 25. Palmer, H.J., and Maheshri, J.C., 1981, Enhanced Interfacial Heat Transfer by Differential Vapor Recoil Instabilities, International Journal o f Heat and Mass Transfer, Vol. 4, Pergamon Press Ltd., pp. 117-123. 113 R ep ro d u ced with p erm issio n o f th e copyrigh t ow n er. Further reproduction prohibited w ithout p erm issio n . 26. Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., 1992, Numerical Recipes in FORTRAN, Cambridge University Press. 27. Prosperetti, A. and Plesett, M.S., July 1984, The Stability of an Evaporating Liquid Surface, Phys Fluids, 27(7), pp. 1590-1602. 28. Rayleigh, Lord, 1900,Scientific Papers, Vol H, Cambridge University Press, Cambridge England, p. 200. 29. Read, K. I., 1984, Experimental Investigation of Turbulent mixing by Rayleigh-Taylor instability, Physica 12D, pp. 45-58. 30. Sharp, D. H., 1984, An overview of Rayleigh-Taylor Instability, Physica 12D, pp. 3-18. 31. Shepherd, J. E., and Sturtevant, B., 1982, Rapid evaporation at the superheat limit, J. Fluid Mech., Vol. 121, pp. 379-402. 32. Sonntag, R., Kaviany, M., Squires, B., Park, Y., June 1992, Slush Hydrogen, Gelled-Hydrogen, Gelled-slush Hydrogen, Report o f Dept, o f Mechanical Engineering and applied Mechanics, University o f Michigan. Final Report for NASA. Lewis RC. 33. Spindler, B., 1982, Linear Stability of Liquid Films with Interfacial Phase Change, Intl. Journal o f Heat and M ass Transfer, Vol. 25, No. 2., pp. 171- 173. 34. Taylor, G. I., 1950, The instability of liquid surfaces when accelerated in a direction perpendicular to their planes, I. Proc. R. Soc., London, A201, pp. 192-196. 35. Yiantsios, S. G. and Higgins, B. G., 1989, Rayleigh-Taylor instability in thin viscous films, Phys. Fluids A, Vol. 1, No. 9, pp. 1484-1501. 114 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . Appendix A. Numerical Modeling of Time Dependent Model This appendix provides the results of the investigation of the transient model developed in Part I, Chapter 4.0 during buildup to provide confidence that it was providing results of the actual physical phenomena. The code utilized routines for matrix manipulation from Dongarra et al. (1979) and for the Runge- Kutta routine from Press et al. (1992). The numerical model was tested to determine a satisfactory integration ending time in order to calculate the vertical velocity. This was done so that the startup transients would not affect the stability results. Satisfactory ending times for the integration were determined to be in the range of 50 to 100 tau time units depending on the parameters in the model (i.e., nondimensional depth of the upper layer-htilda, size of the matrix, size of the initial kick). Various responses are shown in figures A.1-A.4. The response of the model was tested as to which coefficient was best to provide an initial kick to minimize the initial transients before the uniform growth rate phase was realized. This was done in order to decrease the total integration time required. The response of the growth rate to a number of coefficients v/hich received the initial kick is shown in figure A.5. The result was that for a 40x40 matrix the best coefBcient to receive the initial kick was y(21). The effect due to series tennination on the growth rate was also investigated. It can be seen in figure A.6 that, for this specific d/h case, a series of 20 or 30 coefficients are required. This number of coefficients in the series gives 10 coefficients past the discontinuity in the denominator of the series in order to have the growth rate approach its full value. Another check involved the response of the calculated growth rate with respect to the magnitude of the initial kick or 115 R ep ro d u ced with p erm issio n o f th e copyright ow ner. Further reproduction prohibited w ithout p erm issio n . perturbation which started the unstable behavior. It can be seen in figure A.7 that for coefficient y(21) very little difference resulted however on later runs it was noticed that an initial kick of Ix allowed the transient to dissipate quickly allowing the actual instability to setup. During the code build-up, double precision was used for every variable. As various parametric values were used in the different case runs, numerical instabilities were observed which could be removed by an appropriate selection of a smaller time increment. The stability criteria for the numerical model was empirically determined over the ranges of cases studied. Gear (1971) indicated that Runge-Kutta routines had instability boundaries but the method was very robust compared to other methods. He also gave an estimate for an R-K instability boundary but this didn’t prove to be adequate. Also, Press et al. (1992) indicated that the instability boundary would be required to be determined empirically. An estimate for the numeric stability boundary was developed and is given as -= T - ^ 0.01. h This criteria was noticed to require modification as the matrix size increased and a modification which worked for larger matrices of size n was <0.01. dt P (40)' An estimate was made o f the accuracy of the code with respect to changes in the time increment away firom the stability boundary. This result is shown in figure A.8. 116 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.08 0.07 0.06 0.05 I 0.04 2 ® 0.03 0.02 0.01 20 40 60 80 100 120 140 160 180 200 0 integration time in tau time units Figure A. 1. Variation in sigma with time for y(21 )=. 1X 10-9 and alph=.05. 0.08- 0.07- 0.06- 0.05- 1 0.04- 2 ® 0.03- 0.02- 0.01- 0- 0 20 40 60 80 100 120 140 160 180 200 integration time in tau time units Figure A.2. Variation of sigma with time for initial kick of y(21)=0.1X10-9 (initial sigma set to 1X10-6 to have good resolution). 117 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.06 0.05 o 0.04 S, 0.03 5* ^ 0.02 0.01 50 20 35 40 45 10 15 25 30 0 5 integration time in tau time units Figure A.3. Variation of growth rate with time for an initial startup of y(31 )=I .E- 10, for dt=.0175 for stability, htilda=2, alph=.31. 0.16 0.14 0.12 0.1 a ( D 0.08 B 0.06 - 0.04 0.02 90 100 40 60 70 80 10 20 30 50 0 time in tau units Figure A.4. Variation of instantaneous growth rate with time for a 60X60 coefRcient matrix and an initial startup of y(31)=1.0E-10 at tau=05. 118 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm ission . 0.08 0.07 0.06 0.05 g I 0.04 C O 0.03 0.02 y(40) 0.01 80 100 120 140 160 180 200 60 0 20 40 t- nondimensional time Figure A.5. Variation of instantaneous sigma with respect to which coefiBcient gets the initial kick of 1.0X10-7. 0.08 y(6) 0.07 0.06 0.05 §.0.04 0.03 0.02 - 0.01 - 0 20 40 60 80 100 120 140 160 180 200 time Figure A.6. Variation of sigma with time for different series truncations, coeffs truncated at 5, 10,20,30 (20 and 30 series results are identical) 119 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . 0.08 0.07 y(21)E-7 0.06 y(21)E-10 0.05 g 1 ,0 .0 4 0} 0.03 0.02 0.01 0 20 40 60 80 100 120 140 160 180 200 time nondimensional Figure A.7. Variation of Sigma with various initial values for y(21 ) % error - 1.4- 2 2.5 3 3.5 4 0.5 1.5 0 1 Detta time increment Figure A.8. Variation in sigma at 100 time units comparison with various delta time selections, htilda=20, alpha=0.2, the sigma at delta t=0.05 was the basis of accuracy calculation, delta t for stability = 4.0. 120 R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n . IMAGE EVALUATION TEST TARGET (Q A -3 ) / / m y.x % 1 .0 l.l 1.25 a s , Lâ lA P 1.4 ■■ M 2.0 1 .8 1 . 6 150mm V O / V , /A P P L IE D ^ ItW IG E . In c 1653 East Main street Rochester, N Y 14609 USA . = % = Phone; 716/482-0300 . = % = Fax: 716/288-5989 O '993. Applied Image. Inc.. A ll Rights Reserved R ep ro d u ced with p erm issio n o f th e copyright ow n er. Further reproduction prohibited w ithout p erm issio n .

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

An experimental and numerical study on the stability and propagation of laminar premixed flames

PDF

A comprehensive guidance and control theory for spacecraft using pulse-modulated propulsion

PDF

Instability and flow vectoring of jets

PDF

The Structure Of Decaying Turbulence In A Stably Stratified Fluid, Using A Novel Dpiv Technique

PDF

Computer simulation of viscous fluid interaction with a dynamic structural system

PDF

Experimental investigations of swirling jets

PDF

Boundary layer transition due to the entry of a small particle.

PDF

Controller validation, identification and learning

PDF

Cluster aggregation and fluid displacement processes in porous media

PDF

An experimental and numerical study of the effects of heat loss and unsteadiness on laminar strained flames

PDF

Multiphoton laser-induced fluorescence for measuring point specific densities of ground state atomic hydrogen in an arcjet plume

PDF

Experimental and analytical investigation of a ring cusp ion thruster: Discharge chamber physics and performance

PDF

Energy distribution in the plume of a helium arcjet

PDF

Large eddy simulation of turbulent convection

PDF

Effects of convection and radiation on flame spread over solid fuel beds

PDF

Adaptive soft-input soft-output algorithms for iterative detection

PDF

A study of geotechnical applications of biopolymer treated soils with an emphasis on silt

PDF

A velocity estimation model for large eddy simulations of high Reynolds number homogeneous isotropic turbulence

PDF

Application of electrokinetic phenomena in: dewatering, consolidation andstabilization of soils and weak rocks in civil and petroleum engineering,and augmenting reservoir energy during petroleum...

PDF

Analysis of the mathematical modeling approach to understanding the renal concentrating mechanism

Asset Metadata

Creator Chandler, Frank Onice (author)

Core Title Convective instability of fluid interfaces

Contributor Digitized by ProQuest (provenance)

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag engineering, aerospace,Engineering, Marine and Ocean,OAI-PMH Harvest,Physics, Fluid and Plasma

Language English

Advisor Redekopp, Larry (committee chair), [illegible] (committee member), Egolfopoulos, Fokian (committee member)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c17-376799

Unique identifier UC11350265

Identifier 9919021.pdf (filename),usctheses-c17-376799 (legacy record id)

Legacy Identifier 9919021.pdf

Dmrecord 376799

Document Type Dissertation

Rights Chandler, Frank Onice

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA