THE UNIVERS I TY OF MI CHI GAN COLLEGE OF ENGINEERING Department of Mechanical Engineering Heat Transfer Laboratory Technical Report No. 1 THE DYNAMICS OF MOVING BUBBLES IN SINGLE- AND BINARY-COMPONENT SYSTEMS Naoyuki Tokuda Wen Jei Yang John A. Clark Herman Merte, Jr. ORA Project 07461 under contract with: NATIONAL AERONAUTICS AND SPACE ADMINISTRATION GEORGE Co MARSHALL SPACE FLIGHT CENTER CONTRACT NO. NAS 8-20228 HUNTSVTLLE, ALABAMA administered through OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR December 1966

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

ABSTRACT The dynamics of a single bubble moving in a quiescent liquid is analyzed for single- and binary component systems. The analysis is made for the region in which the bubble dynamics is controlled by the transport of energy and/or mass subject to thermodynamic phase-equilibrium at the bubble interface. The dynamics of moving bubble in a single-component system is investigated initially. The interfacial temperature remains constant with time for this case. With the application of the boundary layer simplification and approximating the velocity field around the bubble as a uniform flow, two asymptotic solutions of the bubble dynamics for small and large times, are obtained by means of the coordinate perturbation method. The bubble behavior during small times is dominated by diffusion and/or radial convection while at large times it is controlled by diffusion and axial convection. In the analysis, the temperature distribution in the liquid around the bubble is obtained as a function of dimensionless parameters and universal functions. Then the total heat flux over the entire bubble surface is evaluated and related to the interfacial energy balance condition. The resulting equation is integrated to yield the bubble growth or collapse rate. iii

The solutions for small and large times may be joined successfully at an intermediate time. It is disclosed that the dimensionless parameter r= J e governs the bubble Jd dynamics, where K is a constant, Pe is the Peclet number and Ja is the Jakob number. The results agree very well with experiments. For binary component systems, both injection cooling and boiling are treated. The method employed is the extension of that used for the single-component system. Both the interfacial temperature and concentration vary with time. The dynamics of a moving bubble in injection cooling is governed by the parameters V and p= v where Lu is the Lukomskiy number for gaseous phase. The dynamics in boiling binary component is a function of the parameters, = _ and p Pa X/ uXo h e; _ = ^ ^ in which c^ and O/ are the temperature gradients of the phase-equilibrium curves, Lu is the Lukomskiy number for liquid phase, f and P / are the liquid and gases densities respectively, and oc- o is the relative volatility. The analytical results for the injection cooling case agree well with experiments. iv

TABLE OF CONTENTS Page ABSTRACT.............................. iii LIST OF TABLES...................................... vi LIST OF FIGURES.................................... vii LIST OF APPENDICES.......................x NOMENCLATURE..........xi CHAPTERS I. INTRODUCTION 1......... II. LITERATURE SURVEY........4 III. ANALYSIS 13 A. General Formulation for Dynamics of Moving Bubble........................ 13 B. Mechanism of Bubble Dynamics............ 28 IV. SOLUTIONS FOR DYNAMICS OF MOVING BUBBLES.... 36 V. RESULTS AND DISCUSSION...................... 60 A. Dynamics of Moving Bubble in SingleComponent System................ 60 B. Dynamics. of Moving Bubble in BinaryComponent System - Injection Cooling.... 78 C. Dynamics of Moving Bubbles in Boiling Binary-Component System............... 87 VI. CONCLUSION..95 APPENDICES $..9 LAPPENDICES O................................................ 99 B3IBLIOGRAPHY..... 116 v

LIST OF TABLES Table Page 3-1 Mechakisms of Bubble Dynamics Controlled by Thermal Diffusion...................... 31 3-2 Mechanisms of Bubble Dynamics in BinaryComponent Systems........................ 33 5-1 Experimental Results for Bubble Collapse in Single-Component System (21) 73 5-2 K and KelI For Ellipsoidal Bubbles....... 76 vi

LIST OF FIGURES Figure Page 3-1 Uniformly growing spherical bubble moving with velocity Uoo(t).............. 15 3-2 Phase - equilibrium diagram of binarycomponent systems.......... 20 3-3 Phase - equilibrium diagram of noncondensing gas-cryogenic liquid system......... 35 5-1 Universal functions (%, %o, 9 o ~cJO, Xo kX? Xffo and XS) associated with temperature and concentration fields at small and large times.............................. 61 5-2 Universal function (QFO, h0/ tot-o, X?,o and Xt^i) associated with temperature and concenBition fields at small and large times...................... 62 5-3 Universal functions(@,E) 6, S3;. AX,) Vgo0 ) associated with temperature and concentration fields at small and large times.................................... 63 5-4 Universal functions (0go, e81, e0oo, x,o and %W 0) associated with temperature and concentration fields at small and large times.................................... 64 5-5 Functions ( 0o )A,) 0oV /0 7 ) X ) for the case of u = r2 at large times.... 65 5-6 <(GfXaG f I ) for the case of u = Y2 a large times............. 66 5-7 Convergence of bubble-dynamics solution for single-component system at small times... *.....*......... 68 5-8 Convergence of bubble-dynamics solution for single-component at large times...... 70 5-9 Growth of moving bubbles in single-compon+ ent system............................ 72 5-10 Collapse of moving bubble in singlecomponent system....................... 75 5-11 Collapse of ellipsoidal bubbles......... 78 vii

Figure Page 5-12 Growth of moving bubble in injection cooling at small times for -3-*0....... 80 5-13 Effects of 3 on the growth of moving bubble in injection cooling at small times For Xo = 0.01 and X = 0.90............. 81 5-14 Effects of C on the growth of moving bubble in injection cooling at large times. For X = 0.01, Xoo = 0.90 and qg = 0.436............................. 82 5-15 Effects of!' and q:' on the growth of mving bubble in injection cooling at large times. For Xo = 0.01 and X = 0.90.... 83 5-16 Effects of Xo on the growth of moving bubble in injection cooling at large times. For X oo = 0.90 and r = 0.2..... 84 5-17 Effects of X oo on the growth of moving bubble in injection cooling at large times. For Xo = 0.003 and r = 0.2..... 85 5-18 Growth of moving bubble in injection cooling for Xo = 0.003 and Xoo = 0.856.. 86 5-19 Growth of moving bubble in boiling binary-component system at small times.. 88 5-20 Surface temperature of moving bubble in boiling binary-component system at small time............................. 89 5-21 Effects of Po on the surface temperature and growth of moving bubble in boiling binary-component system at small times for =20 90 for I - 2.0............................ 90 5-22 Effects of (3 on the surface temperature and growth of moving bubble in boiling binary-component system at small times for 3o = 2.0......................... 91 5-23 Bubble surface temperature, phase equilibrium diagram, (3 and / for moving bubble in boiling ethylglycolwater solution.......................... 93 viii

Figure Page 5-24 Effects of ~ on the growth of moving bubble in boiling binarycomponent system for = 0.5, and = 0.5 94 A-1 Functions O (, 3 ),'(3)/,/ o0) and Oc0(o)/O(oo)for equation (A-6). 107 A-2 Functions G( 0 0 <) and, (J..0 0'1X* 108 LX ~ ~ S[n 0

LIST OF APPENDICES Appendix Page I FLOW DISTRIBUTION AROUND A BUBBLE........ 99 II DETERMINATION OF K...................... 101 III COMPUTER PROGRAMS.110 A. Programs for Transport Equations..... 110 B. Programs for Bubble Growth Rate...... 113 x

NOMENCLATURE Roman a,b major and minor axis of an ellipsoid, c constant, 1/sec. CD drag coefficient, dimensionless Cp specific heat at constant pressure, BTU/lbm-~F D mass diffusivity, ft2/sec. e surface viscosity coefficient in equation (2-3), lbf/sec f function, dimensionless hfg enthalpy change in phase transition, BTU/lbm Ja Jacob number, defined as,Cp (ITa-T )or / Cp T dimensionless f k thermal conductivity, BTU/hr-ft-~P K coefficient to determine an equivalent uniform flow as defined in equation (3-28) KpQQ K for an ellipsoid as defined in equation (5-8b) dimensionless Lu, Li Lukomskiy number, defined as and - respectively dimensionless M dimensionless parameter, defined as - p7o,0 order symbols "little" and "big" oho EF pressure, lbf/ft2 Ee Peclet number, defined as 2U0o, dimensionless q dimensionless bubble size defined as T- ] for injection cooling Q total heat flux over the entire bubble surface, BTU r dimensionless bubble radius R/Ro, a spherical coordinate in radial direction, fto R bubble radius, ft xi

Re Reynolds number defined as, U Jo, dimensionless t time, sec T temperature, ~F T* dimensionless temperature, defined as(Ts-To)/(TSf-Too) for binary boiling andTF-To)/(T,-) for injection cooling u velocity component in 0 or x-direction, ft/sec or dimensionless velocity U//U U, velocity component in 9 or x direction as defined by equation (3-28), ft/sec UOO Velocity of a moving bubble, or equivalently free stream velocity, ft/sec. v velocity component in r or y- direction, ft/sec V velocity component in r or y-direction as defined by equations (3-26) and (3-27), ft/sec W ratio of major and minor axis a/b, dimensionless We Weber number defined as -, dimcnsionless x distance measured from the front stagnation point of bubble in orthogonal coordinates fixed on the bubble surface, ft. X molefraction of one species, dimensionless'X* dimensionless molefraction defined as (Xs - X) for binary boiling and(X-X X X-X O) for injection cooling y distance measured normal to the bubble surface in orthogonal coordinate fixed on the bubble surface, fto Greek thermal diffusivity, ft2/sec cIO(CIJ slope of boiling and dew curves of phase equilibrium diagram defined as do-(X/T'Jp and O X( Y7 T)p respectively, dimensionless <, ~asymptotic coordinates defined by equation (Ad7), frl dimensionless xii

parameter in injection cooling defined as c /Lt or a dummy variable for X/R, dimensionless dimensionless parameters in binary boiling defined ^P' as (3 = d pX/0 and X 2 =X*X' dimensionless parameter defined as (gap & boundary layer thickness, ft An~ asymptotic coordinates defined by equation (4-28d) for small times and by equation (4-42c) for large times, dimensionless 0^ defined by equation (4-49d) for small times and by equation (4-61c) for large times, dimensionless temperature function defined as(T-TOOy(T,-TO) or angle in spherical coordinates, dimensionless temperature functions, dimensionless /J viscosity, lbm/ft-sec. iy kinematic viscosity, ft2/sec In defined by equations (4-2c), (4-28c) and (4-49c) for small times and by equations (4-13b), (4-42b) and (4-61b) for large times, dimensionless f density, lbm/ft3 C-' surface tension, ft-lbf/ft tangential viscous stress, lbf/ft2 dimensionless time for small times defined as T~ t dimensionless time for large times defined as(k0o/R)t asymptotic coordinates defined by equations (4-2b) (a-28b) and (4-49b) for small times, dimensionless potential function, ft2/sec or a variable in spherical coordinates, ft stream function, ft2/sec similarity variables for energy transport as defined by equations (4-2a) and (4-28a) for small times and by equations (4-13a) and (4-42a) for large times, dimensionless (X I,-) concentration function, dimensionless xiii

similarity variables for mass transport as defined by equation (4-49a) for small times and by equation (4-61a) for large times, dimensionless Subscripts ell ellipsoid eq equivalent M matching point s interfacial sp sphere 0 initial, reference value 00 asymptotic, at a distance from Superscripts 0()' gas-vapor phase ()e, () derivatives with respect to ~ or 0,() time derivatives. xiv

CHAPTER I INTRODUCTION Problems concerning the dynamics of gas or vapor bubbles with or without translatory motion are continuously arising in engineering practice where the knowledge of growth or collapse rate is necessary. The applications of studies in bubble-dynamics include boiling, cavitation, injection cooling and fluidization, among otherso A comprehensive review of the literature is presented in the following chapter. The time-history of the growth of a vapor bubble, either at rest or in translatory motion in an otherwise quiescent liquid may be divided into three stages: the initial, intermediate, and asymptotic stages. This classification is determined mainly by the bubble size. In the initial stage corresponding to very small sizes, the bubble growth is significantly influenced by the inertia of liquid, the surface tension, and the vapor pressure. After the bubble growth is initiated, there is a rapid increase in its size and its growth rate reaches a maximum. The bubble growth rate decreases thereafter as the bubble grows to larger sizes, corresponding to the asymptotic stage of growth. In the asymptotic stage, the effects of viscosity, inertia force (fluid momentum) and surface tension become vanishingly small and the transport of energy and mass appears to be the dominant mechanism. -1

-2For the collapse of a vapor bubble, Florschuetz and Chao ( 21) has identified a dimensionless parameter B to characterize the modes of collapse as (a) liquid-inertia control case, (b) heat transfer control case, and (c) general case in which liquid inertia and heat transfer effects are both retained. A large value of B implies that the collapse process is governed by liquid inertia, a small value of B is associated with a process essentially controlled by the rate of heat transfer at the bubble wallo This work is devoted to analyze the dynamics of a single, moving gas bubble suddenly brought into contact with a volatile liquido Two cases are treated: Singlecomponent system where both phases are of the same component and binary-component systems in which either one or both of the phases may comprise two components. The growth or collapse of a single, stationary bubble in single- and binarycomponent systems were investigated by Yang and Clark (71) by means of the source theory and by Larsen (42) through the use of the thin thermal-layer theoryo This thesis is the extension of their work on the similar bubbles but differing in that the bubbles are in translatory motion resulting from the action of a buoyancy force. The transfer of momentum, heat and mass between the moving bubble and the bulk of the surrounding liquid is considered to take place in a thin boundary layer around the bubbleO Thus, the boundary-layer concept may be applied to simplify the physical problem in which the bubble may either grow or collapse as the result of a simultaneous action of these coupled transport phenomena,

-3The transport equations governing the dynamic behaviors of a moving bubble are solved by means of the generalized coordinate-perturbation method, an extension of the method of Moore (48,49). Two asymptotic solutions, one for small times and the other for large times, are obtained.

CHAPTER II LITERATURE SURVEY A, DYNAMICS OF,STATIONARY BUBBLE If a spherical cavity in a liquid of infinite extent undergoes a step increase in the system pressure, its collapse rate is controlled by liquid inertia. This is wellknown as a classical Rayleigh's problem (56). Plesset and Zwick (53) extended the analysis to include the effect of surface tension. The criterion for the inertia controlled region is identified by Florschuetz and Chaoc (21) as B>1, where the parameter B is defined B = Ja /' In the asymptotic stage of growth the effect of liquid inertia becomes vanishingly small and the bubble growth is controlled by energy and/or mass transport. By means of the similarity transformation technique, Chambre (14) and Birkhoff, et. al.,(7) investigated the bubble growth in onecomponent system and Scriven (61) treated both single- and binary-component systems. Plesset and Zwick (53) and Florschu:etz and Chao (21) employed the "thin thermal layer"e concept to predict the time history of bubble size during the growth and collapse, respectivelyo Forster (23) and Yang and Clark (71) applied the source theory to obtain the approximate solution for the bubble growth in single- and binary-component systems. The results agree well with Scriven s (61)o -4

-5For injection cooling, the growth or collapse of the vapor bubble was predicted by Yang and Clark (71) using the source theory and Larsen (42) by an integral technique. B. DYNAMICS OF MOVING BUBBLE Langston and Eustis (40) correlated the Nusselt number in terms of the circulation number, for Xe gas bubbles in aquaous glycol. A moving bubble is approximately a sphere when its size is less than 1 mm in diameter. As the bubble size increases its shape gradually deforms into an oblate ellipsoid and then assumes a "spherical cap" shape when its diameter exceeds 2 cm. The deformation in bubble shape depends strongly on liquid properties. The principal dimensionless parameters governing the stability of the bubble shape and translatory motion are known as the M number (= z5 ) ) Weber number 2 a C P and Reynolds number (2 CJ RP ) (30). A small spherical bubble or large spherical bubble or large spherical cap bubble moves upward in a very stable, rectilinear manner. However, in the unstable intermediate region, the bubble fluctuates in shape from an ellipsoidal to a spherical one accompanied by a zig-zig or spiralling upward motion. The continuity of the velocity and stress at the bubble interface requires that the gas or vapor inside the bubble be in a circulating motion. This flow pattern provides an apparent slip flow condition. The degree of slip is found to be closely related to the concentration of the active materials

-6at the interface. Consequently, a small spherical bubble in tap water behaves like a rigid sphere since the thin film of impurities at the interface can resist the shearing force and a full internal circulation of the gas inside the bubble is impededo However, in pure water, there is developed a full internal circulation and the drag characteristics approach that of a fluid droplet, For a boundary-layer type flow the flow may separate from the bubble interface leaving a weak circular wake behind The wake is more pronounced in a more viscous and impure liquid. Hadamard (31) and Rybczynski (58) extended the Stokes0(64) creeping flow analysis around a rigid sphere to the droplet and bubble case by stipulating slip-flow boundary conditions along the interface; ie,, 1. Continuous velocity across the interface, and 2o Continuous tangential stress across the interface, The result is then Bond (10) derived the same solution for the gas or vapor bubble. However, not all the bubbles are found to have a full internal circulation. Bond and Newton (11) were among the first who observed the bubble to behave like a rigid sphere particularly when the bubble is small. They used dimensional analysis to derive the critical radius at which the bubble

-7will behave like a rigid sphere RCY/t ~ J(2-2) They were the first to suggest that this is due to surface contamination. Empirical investigation by Garner and Hammerton (26) has resulted in a clear picture of this subject. Observation of internal circulation was done by adding ammonia chloride fog; and the experiment confirmed Bond and Newton's criteria of Eq. (2-2) for a vapor bubble. When impurities were added, the circulation ceased and the bubble motion behaved as predicted by the non-slip Stokes solution, However, when the liquid was pure, circulation started and the bubble behaved as predicted by Hadamardvs solution for the slip condition, The first attempt toward an analytical explanation for such transition was made by Boussinesq (13), who advanced the hypothesis of dynamic surface tension in such a way that a thin layer of high viscosity exists near the surface layer and introduced the following equation. / J2e?ODc2-3) where & = surface viscosity coefficient which expresses the relation between surface tension and the rate at which the surface area of the liquid changes. For small bubbles, Uo' approaches Stokers solution while for large bubbles Uo approaches Hadamardvs solution,

-8However, both BondVs and Bcussinesq's attempted explanations utilizing surface tension and dynamic viscosity layer effects respectively are indeed hypothetical. For example, Bond9s attempt cannot explain the fact that the bubble follows Hadamard's law if the liquid is very pure. Likewise, the existence of high viscous layer in Boussinesq's hypothesis is not observed physically nor experimentally. The analysis of Levich (44) seems to explain the observed phenomena most clearly. He analyzed the diffusion mechanism of the surface active substance along the interface of a bubble with adsorption and forced convection theory and predicted the existence of thin layer of such substance which establish the surface tension gradient, His analysis is limited to the creeping flow range in which such phenomena are profound, At higher Reynolds Numbers, however, the internal circulation always sets in, regardless of impurities, Levich argues that at higher Reynolds Numbers such thin films of contaminated surface are stripped away by inertia effects An interesting study has been made on the interfacial surface phenomena. Young, Goldstein and B1Qck (72) made a successful analysis of an effect of temperature gradient on the bubble motion, It is known that the bubble starts falling instead of rising in the quinscent liquid if the negative temperature gradient exists. They succeeded in making a qualitative analysis, taking into account the inverse surface tension gradient induced by the temperature gradient, The analysis predicts the condition for which the bubble motion stops.

-9Ramshaw and Thornton (55) attributed the observed oscillation of the bubble to interfacial disturbances. Leonard and Houghton (43) and Datta,Napier and Newitt (17) observed that mass transfer retarded the rise of the bubble. Leonard and Houghton attributed this retardation to the electrokinetic effects due to mass transfer StokesQ solution and its extension to bubbles and droplets by Hadamand and Rybczynski are limited to creeping flow range so that Re <,. 1. The analysis for higher Reynolds Number was carried out by Lunnon (46) for a rigid sphere, and by Allen (2), Miyagi (47), Bond (10) and Bond and Newton (11) for air bubbles. Miyagi (47) was apparently the first to introduce the concept of added mass into bubble motion analysis, the concept which takes accout of the portion of the liquid which must be accelerated with a bubble. This made the bubble motion analysis valid for a wide range, including the rapid acceleration period in which a bubble accelerates rapidly from a standstill position to a terminal rising velocity. Miyagi also noted a zig-zag motion of the bubble. The concept of added mass has since been employed by many others for bubble motion analysis. Davidson and Schuler (18) employed the concept for bubble formation analysis at the orifice tip. Han and Griffin (32) for bubble formation on a heated plate and Gaywood and Birkhoff (27) for moving bubbles. The first study which gave a careful attention to a flow condition around a bubble was reported by Davies and Taylor (19).

-10They experimented with large spherical cap bubbles and observed for the first time that the flow field around the front stagnation region of the bubble is irtotational. The irrr:.otational flow distribution around the bubble has proved to be valid unless the internal circulation is impeded by the formation of surface active material layer along the interface. Because of the irrotational flow, the analysis of bubble motion is much simplifiedo Levich (44) argued that the irrotational flow distribution prevails over the entire rising spherical bubble. He obtained the drag coefficient by calculating the viscous dissipation in the irrotational flow. Cha!' (16) and later Moore (48) more rigorously, showed that indeed the actual flow distribution around the rising spherical bubble deviates from that of irroational flow only by 0 (Re) everywhere. Moore (49) corrected the error in the pressure drag term which was committed both by himself and by Chao. He obtained the correct asymptotic solution in terms of Re, thus the 1st term of the corrected result agreed with Levich's result obtained by viscous dissipation analysis. The stability of a deformable bubble is recognized to be due to the interaction of hydrodynamic pressure and surface tension and it was studied by several authors. Saffman (59) analyzed the stability criteria for the oblate shape of the bubble and its zig-zag upward motion by assuming that the irrotational flow field prevails at the front stagnation point as it was observed by Davies and Taylor (19)

-11The essential points of the analysis are (i) The pressure of the gas within the bubble is constant everywhere and the thermodynamic relation of surface tension and pressure will be P+ 4( - ) 4 co nt (2-4) (ii) The bubble surface corresponds to a stream line. Hartonian and Sears (33) modified the Saffman's analysis by considering the following factors. (i) Inviscid and irrotational flow prevails not only at the front stagnation region but essentially over the entire bubble surface, (ii) Changes in the shape of the bubble surface have effects on the unstabilityo The Weber number measures the ratio of the hydrodynamic pressure forces to the surface tension forces, which maintain the bubble shape. If a bubble is rising in pure and relatively inviscid liquid, the Weber number is found to be a critical parameter, If a bubble is rising in impure and more viscous liquids on the other hand, the Reynolds number becomes a critical parameter, Hartonian and Searsv calculation gives We, crit = 1,26 and Re. crit = 202, Moore (48) improved Hartonian and Sears' analysis by including non-linear terms. Tadaki and Maeda (66) empirically correlated the drag coefficient CD and the deformation ratio of major and minor axis in'terms of Re and Mo However, Haberman and

-12Morton (30) conclude that the above sets of dimensionless parameters are insufficient for correlation. The greatest difficulties for analysis of bubble motion are as follows: (i) The shape and motion of the bubble are interdependent. (ii) The exact condition at the interfacial surface is not well established. If the surface active material forms a thin layer, a full slip flow condition is not applicable. (iii) Separated flow introduces complexities beyond the capabilities of present analytical methods.

CHAPTER II ANALYSIS A. GENERAL FORMULATION FOR DYNAMICS OF MOVING BUBBLE This chapter presents a general formulation for the dynamics of a single, moving spherical bubble growing or collapsing in a binary-component system. Consider a vapor or gas bubble consisting of a binarycomponent mixture moving and growing in a quiescent, superheated, binary-component liquid of infinite extent. Initially, the bubble of radius Ro has a uniform composition XQ and temperature TO. On the other hand, the superheated liquid has a uniform composition Xo, and temperature Too throughout. Then the bubble is brought into the liquid. Due to the buoyancy force the bubble ascends in the liquid as it grows or collapses. After the bubble is in contact with the liquid, the interfacial compositions and temperatures become XI and Ts on the gas-phase side and Xs and Ts on the liquid-phase side, respectively. The following restrictions and simplifying assumptions are imposed on the solution, (1) The bubble is spherical in shape; and grows uniformly over the bubble surface. The effect of its deformation to an ellipsoidal shape is presented in Chapter V-A. (2) The liquid is considered infinite in extent. (3) The liquid is incompressible, -13

-14(4) All physical properties such as specific heat, thermal and mass diffusivity and viscosity remain constant. (5) The viscous heating, kinetic energy and surface energy terms are negligible. (6) At any instant, thermodynamic equilibrium conditions exist at the bubble interface. (7) The bubble may be lumped thermally. This assumption is reasonable because of low heat capacity and high thermal diffusivity in the gas phase. (8) The initial pressure distribution is uniform and the same throughout both phases. The static-pressure change due to bubble motion is considered negligible. (9) The inertia force induced by the radial motion of the bubble surface is sufficiently small to be neglected. Thus, the growth or collapse of the bubble is controlled entirely by heat and mass transfer in the liquid phase and mass transport in the gas phase. This is a good assumption particularly because the study is focused on the growth or collapse of a bubble starting with a finite initial sizeo (10) The bubble moves at a velocity sufficiently large so that a full internal circulation known as the Hill~s vortex prevails within the gas phase. This implies that the stress and velocity are continuous across the interface without forming any surface active material film (30). It is convenient to analyze the problem in a frame of reference fixed at the center of the bubble as shown in Fig. 3-1, i.e., a spherical co-ordinate system (r, 0, f)

-15r U p(t) I e Figure 3-1. Uniformly growing spherical bubble moving with velocity Uoo(t)

-16moving with the bubble. As pointed out by Lighthill (45 ) for an incompressible fluid, the frames of moving and stationary references are equivalent if they are in a uniform, though not necessary constant, translational motion to one another. This would be true in this case. By the adaption of a coordinate system attached to the moving bubble, the governing equations are obtained in two regions of physical space. A.1 Liquid Region (R<r < c ) Continuity equation: p?: —: - ( 3-1) Momentum equation in Q direction. at + 4~i33'aU J~3 21W r =y if 1AP>;3 rt S/R ) (3-2) Momentum equation in r directions'Dt 2rag ar < tA( 2g — r3 33) Energy equation 3T /T+ T 7( T (3-4) -—. v?~ T- V - Mass equation ^ + $3/+,sF" DVT (3-5)

-17The appropriate boundary conditions are u(r, e, o) = v(r, o0)=) o (3-6a) (R, eat) —'(s, t), ](R,,t)- it) (3-6b) u9 (OY, t ) 6- ( (tj), Trt)=r6(t 27z,)C (3-6c)'^A(R. a, t)r-^?:r(R. fag t )< (3-6d) u (oG,, t) dJ fJ (o 6, t) -- Uoo(t) (3-6e) Tf,o)= Tco (3-7a) (Rn, I} = 7(t) (3-7b) rT(an At^^~~)- ~Tc (3-7c) T(r, A t)- Tr et.; t ) (3-7d) X(r,, o)= X (3-8a) (R,, t)- -Xs(t) (-8b)

-18X(oo,bt)=X00 (3-8c) X(r, i, t ) X() et.t ) (3-8d) where / = r2O9(4 / _ f 2 A.2 Gas Phase (O0 r- R) IX~ ~X~i1Li _-' D' X (3-9) fX'.: ~)X x' - (-) Xi e,( o) - X O) R' (% a, t) A= )Z5 (Gt ) (3-10b) X (o,,t) _ o (3-lOc) ar X(,e t)- Xo(#S Qt;e tv) o-o(3-10d) A 3 Interface Boundary Conditions (r=R) The continuity of momentum, mass and energy fluxes at the bubble surface yields, Continuity of mass. (V- s )'s( D'')- x), (3- 1) (^'? — =(- z")? (3-12)

-19Continuity of energy~ (, -R)P,~' T) (3-13) (^)fB- -^ )3 Continuity of momentum in radial direction: p R P 4p-+ (8- V - 4 ) P ~av +aZstv \ _ (3-14) -7(kr r /'i;vr r) In addition, the assumption of the thermodynamic-equilibrium conditions at bubble surface yields: From the dew curve in Fig. 3-2 Xs-, - Ts) (3-15) From the boiling curve in Fig, 3-2 wrl; =(f)L (/;-oo )T-7 (3-16) From the inviscid-flow theory, it is found that the potential flow distribution H-U~(t)S~In (l+S 7 ) (3-17) r-. U(t() cos e (l- + 3-) may satisfy the governing equations (3-1), (3-2) and (3-3) and the boundary conditions (3-6b), (3-6c) and (3-6e)o If sufficient evidence can be found that the fluid flow over the bubble may be approximated as an irrotational flow, then

-20Dew Curve I TOO I D i / ^\ y —Ts gIrJ ^-I Ei H O 2 a, S F —- IBoiling Curve X0Xo X XsX Xs 1.0 MOLEFRACTION X Figure 3-2. Phase - equilibrium diagram of binarycomponent systems

-21equations (3-17) and (3-18) would represent the velocity field around the bubble. This approximation is reasonable as was revealed by the following observation: The tangential shear stress at the bubble surface may be expressed as /In /,ia3 o* 4L (3-19) on the liquid side and /r( l?,o-t) + r,( U (3-20) on the gas side. Since the viscosity of the gas is very small compared to that of the liquid and the velocity gradientin — if...''J 1- for the Hill vortex flow is of the order of unity, the tangential shear stress on the gas side of the bubble surface is very small, i.e., e(KO, 6) —0 On the other hand, with the substitution of equations (3-17) and (3-18), the tangential shear stress in the liquid at the bubble surface may be expressed as r^e-, At) - V&-4.^ s/o' or " _ 9. /2B (3-21) The equation indicates that er would be negligibly small if the Reynolds number for the flow around the bubble is relatively large. In other words, so long as the Reynolds number for the flow around the bubble is relatively large

-22the viscosity of both the liquid and gas play an unimportant role even in the vicinity of the bubble. Thus, the flcw around the bubble may be approximated as an irrotational flow. In addition, it is found that the velocity field as expressed by equations (3-17) and (3-18) satisfy the boundary condition (3-6d). A further argumrent on the approximation of the irrotational flow around a moving bubble is presented in detail in Appendix I. An examination reveals that equations (3-17) and (3-18) also satisfy the initial conditions (3-6a) since the bubble is initially at rest and with zero growth rate. The interfacial condition (3-12) for the case of >~ t simply yields the expression for the radial velocity component of the liquid flow at the bubble surface as ^?s24~(%R,,z)- Rit) <(3-22) which is satisfied by equation (3-18). Now equations (3-17) and (3-18) representing the velocity field around the bubble, are substituted into equations (3-4) and (3-5). The solutions of the resulting equations subjected to the appropriate boundary and interfacial conditions represent the liquid temperature and concentration fields around the bubble. It is rather unfortunate that the analytical solution to the resulting equations is still beyond the reach of the presently available techniques excluding the use of a digital computer. However, from the physical reasoningone obtains an important conclusion that the viscous flow around

-23(This page is blank)

-24the bubble induced by its ascending motion and growth or collapse may be approximated as a potential flow. In other words, around the surface of the moving bubble there exists a hydrodynamic boundary layer:in which the velocity components may be approximated by those of a potential flow, as expressed by equations (3-17) and (3-18). Now, effort will be directed toward Solving the energy equation (3-4) and the mass diffusion equations (3-5) and (3-9), It is convenient to analyze the problem of energy and mass transport in a frame of reference fixed on the bubble surface, i.eo, a coordinate system moving with the bubble surface. Equations (3-4), (3-5) and (3-9) in the spherical coordinates fixed at the center of the moving bubble is now transformed to the orthogonal coordinates (x,y) fixed at the forward stagnation point on the bubble surface, X measures the distance along the bubble surface and y measures the distance normal to the bubble in the outward direction for the liquid phase or in the inward direction for the gas phase. With the boundary-layer simplification, equations (3-4), (3-5) and (3-9) in the new coordinates reduce to at T r i T Ad T (3a23) for heat transfer in the liquid phase, at e a3 v^ " u 2 (3-24)

-25for mass transfer in the liquid phase, and DX' ~.,'' D Js X v+3X (/3 -D)2.#4[(k+) (3-25) for mass transfer in the gaseous phase, respectively. U and V correspond respectively to u and v expressed by equations (3-17) and (3-18) in the new coordinates. U' and V' are still unknown quantities. Due to the shift of reference systems, a uniform but time-varying velocity R(t) which is equal to the growth rate of the bubble is superimposed on the normal component of the velocity field around the bubble surface. Thus the normal velocity V becomes zero on the bubble surface and -R(t) at a distance from the bubble in the liquid phase. This leads us to approximate V in the form of id ^)^:- (3-26) which satisfies the conditions at the bubble surface and at a distance from the bubble. Or, in a simpler form, V may be expressed as V(t)-=-2 2 + O(f-3) (3-27) with a deviation from equation (3-26) by the order of --- For a thin hydrodynamic boundary layer over the bubble,the order of may be considered to be very small. Now, the velocity component U along the bubble surface is assumed to take the form of

-26U (t)-K UU Ct) (3-28) in which K is the proportionality constant to be determined by-equating the total heat fluxes over the entire bubble surface obtained for the potential flow distribution and the flow field expressed by equations (3-27) and (3-28). The method of determining K is presented in detail in Appendix II. The appropriate boundary conditions t6 be satisfied are obtained from equatiohs (0-7), (3-8) and (3-10) as, o t,. O )-t, oo, t ) - T(o,, )'- T -r(X to, t3- i() t) (3-29) for equation (3-23); x((, a o) =- X) (r, o, t)-=X o, ( t )- )'. X (C, ot) = X," () f((eot) ati n (3-) -30) for equation (3-24)9 and X't8,o)=XO, Xito )=X, X Y -, X(f) X'(x, i, t)= X(-x, it) (3-31)

-27for equation (3-25). In addition, two interfacial (or matching) conditions to be simultaneously satisfied are obtained from equations (3-11) and (3-13) as: The energy balance at the bubble interface requires 74 C^ ) ^ ^S 7-,). d3 (3-32) The mass balance at the bubble interface requires 4z d ([~ ~ - Y ) ~~-S:- z)PnsI - E D f' d - 0 for boiling in a binary liquid mixture, and 77, E(1-Xs)pR=n - 2nRDf Ja|XT d (3-34) for injection cooling. The thermodynamic equilibrium condition between the temperature and concentration at the bubble interface yields i - c21 __ XQOO /(3-35) T - T XS- Xo for injection cooling; and x5- x~-"o'x' (' - Tx ) ~(3-36a) and Xs^Xco l ($j (7a - Tooos ) (3-36b)

-28for boiling in a binary liquid mixture. B. MECiHANISM OF BUBBLE DYNAMICS B.1 Single-Component System In order to understand the basic mechanism in the dynamics of the moving bubble, the differential equation (3-23) and its appropriate boundary and interfacial conditions (3-29) and (3-32) are partially non-dimensionalized as follows. With the definitions of tJt, ^ ^,' /_, u v,, -T:-s~ TD' -~ 0 JX ) (3-37) -.P' (To-TS?..' d —- d ~ 2 in which U() K Uooo and is the thickness of the time-dependent thermal boundary layer, equations (3-23), (3-29) and (3-32) may be rewritten, following simplification, as -i ). z +^ (4t+ i )- (3-38) G(o, )=.. o'(,;t).0 O,= (YjO)= (3-39) aff ss/'R )~ (3-40)

-29The combination of equations (3-26) and (3-40) yields Zi — i& )o l( i (3-41) It must be noted that the magnitude of () within the thermal boundary layer is 0(C) and * j ab> y,) 2 a0e o(/) In order to satisfy two of the boundary conditions it is ao necessary to retain the u term in equation (3-38). An examination of equation (3-38) with the substitution of equation (3-41) reveals that the following regions are of particular interest. Region a: For =0( ), the transient term is 0(1). Region b: For =)(/-3} the radial convective term is 0(1) Region ca For Jd =O/?O2) the axial and radial convective terms are of the same order. Region d. For S=O(Pe-) the axial convective term is 0(1), Now these regions will be examined separately as follows, Regions a and b correspond to the initial stage of the bubble growth in which diffusion and/or radial convection are the dominant mechanisms. In the initial stage of the bubble growth one defines Th^e equ ati. (3-42) Then equation (3-38) may be rearranged as t^/^ -J^ )it ^ ^ ( )- (3.43)

-30It may be shown from an examination of equation (3-43) that only the second term on the left-hand side and the term on the right-hand side of the equation remain when s << I and P <~ /.This indicates that the diffusion process is the dominant mechanism in region a. However, if = 0 (O ) but 1 t / so that 2(< 1, then both diffusion and radial convection become important but the contribution of the axial convection is still negligible. The contribution of the axial convection depends on the magnitude of the parameter. Regions c and d correspond to the intermediate and asymptotic stages of the bubble growth. By defining r- Uo v'^ - - ef SI 0 (3-44) equation (3-38) may be rewritten as'dO __ / _ 2_ - (3-45) Equation (3-45)) indicates that when >, the radial convection becomes negligibleo This corresponds to region d. However, when O-.0/) both the radial and axial convection are of the same order, which corresponds to case c, This can be summarized in Table 3-1.

-31TABLE III-1 MECHANISMS OF BUBBLE DYNAMICS CONTROLLED BY THERMAL DIFFUSION Stage of Dominant Region Situation Bubble Growth Mechanism a 2<0~ I initial diffusion b 2wt, </ initial diffusion and radial convection c =, / intermediate diffusion, axial convection and radial convection d fr' > / asymptotic diffusion and axial convection (i) For regions a and b, the analysis of the bubble dynamics of translating bubbles reduces to that of stationary bubbles (ii) In regions c and d, the axial convection effect is important. B.2 Binarv-Component Systems Since the energy equation for a binary-component system is the same as that of a single-component system, the parameter &r is still a criterion in determining the mechanism of the bubble dynamics. With the definition of: t, c-, 4 J,, a-,,: T- %L, X=ZaZsa), O ).'LO =o — c, %J~oft), L>M-Js) -T A, T' —,-X-= x; D

-32the interfacial conditions (3-32), (3-34) and (3-35) for injection cooling may be rewritten as M2r =if 1 (3-47) x Xu4 a/) < (3-48) +x'0X0 JC9 P) A ff for a linear relationship -- = T between the temperature and (3-49) concentration at the bubble surface As was discussed in the previous section III-Bo,, when the bubble dynamics is controlled by heat transport, equation (3-47) indicates /, T,) are 0(1). Hence, according to equation (3-49) X is 0(1). Therefore, one finds from equation (3-48) that -or, %) 0(/) Hence, if 9 I L < _ Z would be vanishingly smiall and mass diffusion has very little effect. The bubble dynamics is then heat transport controlledo Similarly if (3 J / t:he bubble dynamics will be mass transport controlled. In reference to III-Bol one may now realize that the bubble dynamics in binary systems can be classified into the same four domains as in the case of a single-component system. Each domain may be further divided into three subdomains depending on the magnitude of the parameter (3 For example, Table III-2 presents three subdomains of the domain a.

-33TABLE III-2 MECHANISMS OF BUBBLE DYNAMICS IN BINARY-COMPONENT SYSTEMS Region Key parameter Stage Dominant Mechanism 2 al < < small time energy diffusion a2 2Z</ (3> j small time mass diffusion a3 <| | small time both mass and 4[s )' \ energy diffusion The other three domains b, c and d can be divided in the same way as the domain a into three subdomains, B.2 Boiling in Binary-Component system The parameters which determine the mechanism of the bubble dynamics'are obtained from the governing differential equations and interfacial conditions as r___ e ) g p, I _- (3-50) It should be noted that L appears in the parameters rather than LU as in the case of injection cooling analysis. Due to the large density JA P' ratio, is generally much smaller as 4JL~ y compared to 9 L for injection cooling. However, and I play an identical role as ~ in injection cooling, that is, they indicate the relative importance between the mass:&energy transport to the bubble growth. In addition, it is disclosed in equations (3-50) that do and 0( which are the gradients of the dew or boiling curve in the thermodynamic equilibrium diagram and the interface

-34relative volatility ratio ( Xo,- Xoo ) are,also among the physical properties involved in the criteria.

-35TA I- - w iTB LJ/ To LL!T 1 j 0 Xo X 1.0 MOLEFRACTION X Figure 3-3. Ihase - equilibrium diagram of noncondensing gas-cryogenic liquid system.

CHAPTER IV SOLUTIONS FOR DYNAMICS OF MOVING BUBBLES Two asymptotic solutions, valid for small and large times, are obtained in the following to predict the bubble growth rate in single- and binary-component systems. A. SINGLE-COMPONENT SYSTEM Aol Small Time Solution The small time solution, which signifies the solution for small times after the start of the bubble motion is equivalent to the case << where QJ is the velocity component along the bubble surface as expressed by equation (3-28), and U is its time-derivative. With the definition of; t Too.(.'h /-.. ^ ^ --- ) (4-1)!~, l o where ~'= - ).f (4-2a) o''- "' —2 -'X, t=e,(4-2b)'. — —'-i ti~'r-, ~,!9 ( e'"-' (4-2c) the energy equation (3-23) and its appropriate boundary conditions reduce (3-29) to -36

-37a? ^ ~- 2 2.t ) + j-/ -32 ) -- - (, - ) t, (4-3) and O(oA'.-A,) )= O t(O, n,)- / (4-4) respectively.. The temperature function 9 may be expanded as 3(,~',:; ) - so^ (~ ) l:/) +; Hi /(~) rz e;. t ^,/?)<y4 - - - (4-5) In applying the method of power series, the series concerned must be convergent. For the temperature function as expressed by equation (4-5), it requires that the tn and,E be small and are in a decreasing sequence. In an isothermal system Levich (44), Moore (48,49) and Chao (16.) obtained a relationship between the drag coefficient an-d the Reynolds number of a steadily ascending air bubble as CD oC Re This indicates that if the change in the flow velocity is sufficiently small compared to the velocity i tself, then the translatory velocity of the gas bubble may be directly proportional to the squate of the bubble radius or uoa d 7R

-38On the other hand, if the change in the flow velocity is very large compared to the velocity itself, Miyagi (47) proposed U00Q/- 2fet The relationship was obtained experimentally and theoretically by the concept of added mass, For small times the bubble velocity as approximated by equation obtained from Miyagi (47) gives u. (+ e' )(e -/) i 2 eec,, ct, -{(i e O - /+ e'Ct. ew -4e( + e* As t --— ^ O, thpe first equation gives J = - 4S 0 (t3) Consequently,;,= oct (.. O't,),,o-O( )" -—. and the parameters F-2 are in a decreasing sequenceo As to the asymptotic nature of (?B), there is a certain uncertainty because the bubble behavior is still unknown at this time. However, in case the bubble grows from a finite initial size by a diffusion mechanism, this growth may be expressed as - w eR, t+ QQot(/+ Q,t a, t- -— ) where.0/ a j^ - are constants.

-39Then it can be found that as t — O ^- ~(t4.?7v- ~ ) S -—, t- O2t ) Although the parameters Cn do no.t form an asymptotic sequence of decreasing order, it will be seen later that the coefficients of the tn terms in the temperature funct ion 0 decrease very rapidly. Now, equation (4-5) is substituted into equations (4-3) and (4-4) followed by collecting terms multiplied by the various powers and products of the'r, and r It yields. Zeroth order solution (4-6) First order soluticn - C0 ta t i 1/t 2 4 1/i_:-o —,. +e.,,,,f -.*;- o n -Lgt-~tWF- + 61 a~rs E

-40Second order solution o) -~ Ot4- n 7o +2tt-at -. O (4-8) ~t, -^-^^!^ ^ Zo/) 2 <2 @DQ,>o..-2 /i i>.. f!.2 l0 2 00)- /, eo()= o; (4) and all other functions at 0=0 and LO equal to zero. -o+, I?/3 +~f, The universal functions C^ are numerically integrated 6;564~ 4 f-i2 --- -wh appropriate bou)y _,lo)~ a, BCao)- O; (4-9)

-41One is now reminded thiat the absolute value of the coefficients for the t2 terms such as 0(o):jQ()D 68(o) —- decrease like 0.5731, 0.1857, 0.0458. Depending upon the accuracy sought for the solution, the following three approximations are presented. Zercoth order approximation which is accurate to 0(1) Ir2, s 0 & (4-lla) t.irst order approximation which is accurate to 0(g2) rf- = ~^ o^o70 s (4-llb) Seccond order approximation which is accurate to 0( Z) ( 0-0 ~7? r) (4-11c) in which Z is defined as - or 4d( A.2. Larg.e time soluti'on The large time solution which signifies the solution for large times from the beginning of the transient behavior is equivalent to the case U<K /. The flow and temperature field at large times may be considered to be in a quasi-steady state. Thus, at large times the rate of molecular diffusion is balanced by the rate of convection. By defining T,^~-r (7' (s,, --- ) (413) IS" *00

-42__- (4-13a) r U: cZ (-J2/J-\t. (4-13b) It is interesting to note here that the parameters Zt defined for the small-time solution are in an identical form with the parameters j? defined for the large-time solution since the relation of /, OX 2 is.imposed for large timeso By expanding J in the form of -" -It"(R. at +,At -+- - )* orU -Z E CQa( + Q~, t t - --- )'::It can be r:eadily seen that unless U-> 0 as - > 0o t ~3 form a sequence of decreasing order as t — from equation (4-13b), since only a finite range of 2C ioe., from 0 to 7C R is to be considered in the analysis. The energy equation (3-23) and its appropriate boundary conditions (3-29) may be rewritten as 5-t2 o I -225 + 49 (,-3+S, A) (4-14) 37

-43f(o)00), 0o, /0 (4-15) With the substitution of the temperature function 0 expanded-in the form of ) (7y a; )- QC) +II'/, - - -(4-16) into equations (4-14) and (4-15) one obtains Zeroth order approximation -y-.%'- lo' ~ (4-17) First order approximation 0, o-;j 2 o+o -E' ~ r"/ e' -i -9 + + -2''o (4-18) Second order approximation (j,) - 2'....' -:^ o 2.o- } - o 2.9....-2 9 t 2 / / 1%-,6 +%^^= (4-19) ~~~~~z~~~~~~~~- ~~~~~~~~~~~F,~ ~ ~ ~ wf

-44Third order approximation 3 (J) -AQe I-go 2 -WDt3, O = (4-20) 2 ^^ ^ 2 FF >^F^ f subjected to the appropriate boundary conditions Gj(o)= 1, ~oC D) -~ (4-21) and all other functions at =0 and to are zero. At large times the energy flux varies significantly over the bubble surface. Therefore, a surface integration is required to evaluate the total heat flux over the entire bubble surface. With the substitution of O CO ) obtained by the numerical integration of equations (4-17) through (4-21) into equation (3-32) one obtains, ~(0..64.. + o, /443 r~ _oo.''4a9r ~ o1o//d~ r 0 G)(fi4)(3 4-, 0J- 1 -'o / o o ~::. 14__ ~) _ o ~,0 / 1? /' o ~ o W} ~. r ~( 4 - 2 2 ) In case a quasi-steady state is assumed to exist at large times, the relation U - holds. Then equation (4-22) through the substitution of c>jpisd ~,n =- -rnr \ l^p -- 0?,J?7 and.' 7 SP \' I i \ I6 ^ -I\ 1 fP4

-45becomes'2. S,, 2 /:o/o6+ 0.3/,L5 -/97+ o. O.O'2 t 0.3688. - 0.6/ - 4uU6~~ 676Uj4 # (4-23) o,7+3.) 17 2 Just for the sa'ke Of compari'son Wiith the quasi-steady behavior, if r.aid U were considered to be independent, then equation (4-22) may be rearranged, as -Q.14-32At o.o64 j U -- )d (4d-24) (-o34a-^ ~.,32S^ * —- (4-25) i/ T2/ 0,. 2

-46Both equations (4-23) and (4-25) converge very rapidly as ( —t2 2 or when the axial convection becomes a dominant mechanism. Three approximate solutions are obtained from equation (4-23) for the bubble behavior at large times: First approximation which is accurate to 0(1) (4-26a).2 Z/ / Second approximation which is accurate to ( &) i^ - o o. S (4-26b) Third approximation which is accurate to 0( Z6j) ^ 1 - 1.010J + 0.3 ^ ^lsL 2 U/^ i/f^'.(4-26c) 2 0I 3'8 u u~~~~~3 B, BINARY-COMPONENT SYSTEM INJECTION COOLING CASE The problem of bubble dynamics for injection-cooling case may be simplified by the following physical argument: Noting that the diffusion coefficient of a gas is much larger than that of a liquid, the concentration distribution in the gas phase may become uniform more rapidly than in the liquid phase. Therefore, one can assume that the bubble growth at large times is controlled by heat transport in the liquid phase while the concentration distribution in the gaseous phase is uniform at a saturation value. However, during small

-47times, the bubble behavior is controlled by both the mass diffusion in the gas phase and the energy diffusion,, in the liquid phase. For such case the boundary layer concept cannot be applied to the gas phase since the diffusion of, gas toward the. bubble center is very, rapid. The problem i.s analyzed for two cases; A Bubble behavior at small times and - > 0. By an integral technique, Larsen (42): obtained the solution for the case at finite e, that is, both mass and energy. diffusion but no convective ef:fiects are presented.o Here, consideration is given to the case only energy difffusion controls the bubble dynamics B. Bubble behavior at- large times during which the gas concentration inside the bubble distributes uniformly. B31. Small-time solution with 6 -> 0 As —' 0, the gas concentration is considered uniformly distributed and equal to that at a saturation state X~, Hence the bubble dynamics is governed only by the energy equation (3-23) subject to the conditions of thermodynamic equilibrium, and energy and mass balance at the bubble surface, Due to the timewise variation in the gas concentration during the bubble growth the temperature at the bubble surface is time-dependent. Consequently the fluid temperature T must also depend on Ts and its time derivatives Ts T'S o By defining -hT0T - Q( o 0/ oJ8"- ) (4-27) [& 1

-48where r1 —~~~ ~~~~ ^ ~~(4-28a). TS''' (4-28d) ^;=~R.-)'P':~/.- L)2)' [, --- (4-28b) The asymptot-ic nature of the parameters, and ad,has. been already discussed' in IV-Aolo If T varies with time in the form of Tis Q,+Va +'{j t — (4-29) where ( is a constant and of 0(1) and f l> 2, it can be shown that the parameters are in a decreasing ~Sc.Now WI U (4-27= sequence. Now with the substitution of equations (427) and (4a28) the energy equation (3-23) and its appropriate boundary conditions (3-29) can be rewritten as -I"~~~S Q.0S-) 2where is a constant and of 0(1) 3and

-49and Q(~0o^) finy bn)=l, Deoo^ /J S )= 0(4-31) respectively. With the introduction of the temperature function 0 expanded in the power series as e(i * tn wiS~ in) = (3 ~/?) t+Z ^5) (a) ^ ^i %(v)-+Z^/; ^ ".<(1l~t — (4-32) into equations (4-30) and (4-31) followed by collecting terms, one obtains- a set of ordinary differential equations. The differential equations.for 90 of the zeroth order approximation, 9Qr and by of the first order approximation and Q 0>B and Q i, of the second approximation are identical with equations (4-6), (4-7) and (4-8) respectively. The equations and appropriate boundary conditions for the functions related to the parameters n are First order approximation (4-33)' -K29 + + +20-6,0 Second order approximation involving. 3 (2A 2 - ~' 2dr,:3 0tS 2 6Sg S,300 0k t 8 2 -*?6J)2(3 +090 t26o (s2

-50%(fo) _ g(O ) - 0, —- (4-35) The boundary conditions are that all functions at 7 - O and 0o are zero. The solution for the bubble dynamics of the injection cooling case may be obtained in a similar way as the case of a single component system except the expression corresponding to equation (4-1.0) for the single-component system becomes much more complicatedo Therefore, only the solution up to the first-order approximation will be presented here. Since f = O(*t ) 0(t4)) the terims with the parameter Fo are neglected. The substitution of equation (4-32) and the numerical value of the first derivatives of the functions tnO at = 0 into equation (3-32) produces 0.3eo jo pi +-? a s (4-3-6) Since the gas concentration is uniformly distributed inside the bubble, equation (3-34) reduces to y o. -Xc^ _ S__ X___ )(4-37) X*0- /'_ With the substitution of /)' T obtained from a linear relaticn ship between the interfacial temperature and concentration, the parameter )0, defined as _O — ~ e, (4-38) -/-"XcJ

-51and the parameter ) defined by equation (4-37) to eliminate T) Tj? and X in equation (4-36), one obtains *^ os~l O -^ rs" ) (4.39) "J5, + l.SOfI-Xj~ -ow9 ~(-^~j-l 2ly4 ( The bubble growth predicted by equation (4-39) is accurate to o4-~ ) or 0 ) For the case of fanift'e 3, Larsen (42) obtained the bubble growth equation by an integral technique as 3~~~~~~~~-1 * 2(I-X,)V7'' )__ where ( 3 Lu Bo2 Large-time solution At large times the bubble growth is governed by the energy transport in the liquid phase subject to the conditions of thermodynamical equilibrium and energy balance at the interface, With the introduction of i7 -e(? -) ) (4-41) TS o V)= Li <JC (4-42a)

-52^~ 10^' fJ U t2 - (V / S )) ~ (4-42b) TS -' T s — - Ito (4=42c) the energy transport equation (3-23) and its boundary conditions (3-29) can be rewritten as (-lSt^%r c(r, f — [ -l and2a2 ana {' &ao (o 2,, =, 0 O ( ^ )- <. (4-44) respectively, By substituting the temperature function (4-45) into equations (4-43) and (4-44), one finds that the....,. functions d0,; and 0; have to satisfy equations (4-17), (4-18) through (4-21), and (4-33) through (4-35), respectively In a similar fashion as IV B-l, the substitution of equation (4-45) and On(~) into equation (3-32), one obtains the solution for the bubble growth up to the firstorder approximation.

-533.O3-2 37a p', ~ -- (4-46) With the aid of equations (3-49), (4-37) and (4-38) to eliminate:. T*, T* and X*, in equation (4-46), it yields the bubble growth, equation as /l/o (1-X,){,- 1) ^ -/ -- v —--—,,...................... (4-47),rt q W(/-Xx ) xa,- X. e- ) _3- /I) which is accurate up to 0)^ ) or ( C. BINARY-COMPONENT SYSTEM- BOILING Scriven (61) was among the first who analyzed the bubble dynamics problem for nucleate boiling in binarycomponent systems. The growth of bubble in such a system is governed by the simultaneous transport of energy and mass in the liquid except during the initial stage of the bubble growth in which the sutface tension and inertia effects are important. Consider a binary liquid mixture which is in thermodynamic equilibrium at points 1 and 1 in Fig. 3.2. The liquid temperature is Too, and the concentrations are X > and X0 in the liquid and vapor phases, respectively. Suppose at zero time the system is superheated by A1t- ipO O. Then as time goes on the equilibrium condition assumed existing at the bubble interface will shift along the equilibrium lines from points 1 and 1v respectively to points 2.and 2V which are at the temperature Ts and

-54concentrations X and X. respectively. Thus the interfacial temperatures and concentrations on the vapor as well as liquid sides are time-wisely varying. The concentration and temperature gradients in the vapor phase may be considered negligible compared with those in the liquid phase, since the diffusivitiesare much larger in the vapor phase than those in the liquid phase. Therefore, the bubble dynamics is considered governed by the energy and mass transport in the liquid side, The solutions of the energy equation (3-23) for the iiquid phase subject to the boundary conditions (3-29) are obtained in the previous section as equation (4-32) for small times and equation (4-45) for large times. In the following only the mass equation (3-24) for the liquid phase subject to the boundary conditions (3-30) will be first solved for small and large times. Effort will then be directed toward finding the expressions for the bubble growth at both small and large times. Co 1 Small time solution With the definition of Xs_ X0. X-c- X(XlT 23)?2' ) X 0" g, — g< 25 — -) (4-48) where'- J~ V(4-49a) 2.,, 2 RA L) U 2 2(tJJ(_)_- - (449b)

-55FO= r0 S ) 4)X- (4-49c) _u., * u a-xO "x~ -x~' S —-49d If the mass transport equation (3-24) and its appropriate boundary conditions (3-30) may be expressed as,x,+ xr(X 2 -2) - 2;+- 2x_., ) ),, —? -2;;) 2'; 3(~ 4 _ 4) a, + ) -_, (4-50)' ) - - - and:X(o,,,, )- ),,,,?,,nJn ): ~ (4-51) respectively~ Now the concentration function X is expanded into the Taylorgs power series as X(?.,, A ) = o(f)) +z z;(' )+- Z;,,( ) (4-52) +S; Xg,), ) The substitution of equation (4-52) into equations (4-50) and (4-51) followed by collecting terms yields a...

-56set of the ordinary differential equations for the functions Xo, X/ 2-( -." Since the energy and mass equations are identical in the form and the temperapture and concentration functions are expanded in an identical form, the differential equations which define the functions XCn are also identical with those for the functions 0 in equation (4-32) provided that ri is replaced by Sn. Therefore, the solutions for a' a'nd En are the sameo With the substitution of the numerical values of (0o) and Xn(o the interfacial conditions (3-32) and (3-33) may be expressed as 0.t 4 0.10' 0.' 0. -= Tr (4-53) 06 + o.4x/o x" % ~'"r:~-F,,, aO ~S - /o. 0'4ott2D t =S(SI 4 (4-54) T/-T' The interfacial equilibrium conditions(3-36) are rewritten as =- a, (l-"T') (4-55) X = Jo (X- Tjl ) (4-56) The important parameters in equations (4-53) and (4-54) are the parameters C and p. 0 and p. are analogous to 3 in the injection-cooling case and indicate the effect of mass diffusion, As Qj ->O or equivalently

-57Lu- >a, the process is predominantly heat-transfer controlled. On the other hand, as O -, > Do or equivalently L- >O j the mechanism is predominantly mass-transfer controlled. T*, the interfacial temperature in dimensionless form, may be obtained by combining equations (4-53), (4-54) and (4-56) as 7- __ _p___ / (o~s o0&o8 #?) (4-57) The combination of equations (4-53) and (4-57) produces the bubble-growth equation as r2 ~ olf- -1~ o -- os t)7 T </ (4-58) o " - ~ ( -- S+ The approptiate initial conditions are r(o ) =/, O (4-59a) T ~/ ('~ ) ~=-. ~/, 0 ~(4-59b) T(o) may be determined by substituting equation (4-59a) into equation (4-57)o The bubble behavior predicted by equation (4-58) is accurate to the largest value of 0(-s); 0(T, ) ~o/or ) C 2 LargEe-time solution For large times, the concentration function is defined as K.. 4 4 4?,, —,' S'o $,, —-) (4-60) xs-xop

-58where r = ij- (4-61a),R u /- ffiJY -- 4(4-61b) ~"=Xe''w^u' <XO r (4-61c) or in the Taylor s power series as,- = i I) X~tG; t^. ) +2 ~,~,/?2)t —- (4-62) where the functions X n are identical with the functions fn in equation (4-45) At large times, convective effects become important, and the bubble behavior deviates remarkedly from that of a stationary bubbleo In addition, the mass and energy fluxes vary considerably over the bubble surface. Therefore, the surface integration of the fluxes over the entire bubble surface is required in taking the interfacial balance of energy and mass. Taking the surface integral of the flux terms and substituting the numerical values of the transport functions, the interfacial energy and mass balance conditions (3-32) and (3-33) yield 1/01t/ /3 32 JT.' 0 = t r TP (4-63) /o/Sr9 +.'/3,,2 p *-. - -. -^^ r- ___ (4-64)

-59The combination of equations (4-63), (4-64), (4-55) and (4-56) produces T7 / ~+ o, _ /f.ol/9+o7or (4-65) 2(3 r "~'J -/$3 (o, —1,1-~3 (4-66) subject to the initial conditions which correspond to the final valusof-the small-time solution (4-58)o An examination of equations (4-65) and (4-66) discloses that the bubble growth at large times depends on three parameters T'/ (3 and p,. ~ determines the effect of the axial convection, while go and / indicate the relative importance of the mass to energy transport.

CHAPTER V RESULTS AND DISCUSSION Equations (4-6) through (4-9), (4-17) through (4-21) and (4-33) through (4-35) for the universal functions fi or equivalently Xq were integrated by means of an IBM 7090 digital computer. The Q / and X* functions associated with the temperature:and concentration respectively are plotted in Figso 5-1 through 5-6. The fundamental requirement in the analysis is that the temperature and concentration functions, 0 and X, expressed in the form of the asymptotic series as equations (4-5), (4-16), (4-32), (4-45), (4-52) and (4-62) or equivalently, the bubble dynamic equations (4-10), (4-11), (4-23), (4-25), (4-26), (4-39), (4-47), (4-58) and (4-66) converge in the specific domains of time. The speed and radius of convergence of these asymptotic series may be determined by comparing the results thus obtained with the existing exact solution and/or experimental results. The radius of convergence is especially important since it determines the possibility of joining two asymptotic- solutions obtained for small and large times, A DYNAMICS OF MOVING BUBBLE IN SINGLE-COMPONENT SYSTEM As was discussed in Chapter III, the dynamics of a moving bubble at small times after the transient is initiated is dominated mainly by diffusion and/or radial convection. -60

1.0 0.8 0.6 0.4 \!'~ o\ -e XTTO X.O. c. 02 -> 717 c f 1.0 / -- 2..0 — 3.0.4.0 -0.2 -0.4 er X ro MOW 0 r/,~, -0.6 Figure 5-1. Universal functions (9o0 tbo 0 O 0, and XX ) associated with temperature and concentiation fields at small and large times

0.4 o &l0o.Xt, xtoo / \ 0.3 / 0.2 - \ > 0.1 // \ Figure 5-2. Univ l f n (. associated with temperature and concentration ^ \\ / \ -.-0.3 / / / o 0o~i r r. Figure 5-2. Universal function (_C -0-r/ -.0 associated with temperature and concentration fields at small and large times.

0.4 0.3 0.2 \ C0.1 e ^ \ C 0' —~-~- - I -~Lc~c. ~ I 9 q?0 ^^^ ^^ 0 ^ Q ---- - ~^^ ^- 4.0 ^ ^.= O. s'06 xsl - 0.2 I V ^f' Xs -0.3 e- )o o Figure 5-3. Universal functions (0 )( g0 ) associated with temperature and foncent at iot fields at small and large times.

0.4 0.3 0.2 e x B' X 0~~~e 0.3 / n X8 ~~~igr 5-4 1JieIljnt~l - - * ^ ^ "r ~ ro^.0^ 4.0 c-0 s.1 -0.2 Figure 5-4. -Universal functions ated with temperature an and? &S"I ) associa th temperature and concentration fields at small and large timeS.

1.0 0.8 0.6,o,X 0.4 C 0.2 xO /0 0 C -0.4 e40^^r10^" )XC —^-^^ —3.0- 4.0 )XC -0.2:\.< --— ^ /\ ^~80' X80 - 0.4 et o'^o0 -0.6 Figure 5-5. Functions( 6o, 9Os l\oAo )XA / XX~ for the case of u = v2 at large times. 2.O st larg times

0.4 0.3 0.2 ]=.i-j e El "E X, j —^yr ^"^^o^^-^^ ~ r 0~~~~~~~~~~~~~~~~~~~~~~~~~00'E O~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~a Cb O 7 / y' -,xt. -0.2- / 0 i 0 e 04 e~0~0, x4,00, -0.3 Figure 5-6. ( 0^ XX,,eQX)Q Xf,^ for the case of u = r2 at large ti s.6

-67That is)the axial convection effect is less significant, or corresponding to the case of T<< I The exact solution of Florschuetz and Chao (21) was employed to compare with the results obtained by equations (4-11). Equation (4-lla) which yields r= / A S <(5-1) upon the integration indicates only the diffusion effect; Equations (4-llb) and (4-11c) include the effect of radial convectiono Fig. 5-7 shows that the series is rapidly converging and have a rather large radius of convergence. It is instructive to examine the behavior of the bubble dynamic equation (4-11) as — W0. When -- 0 corresponding to the case of stationary bubble, equation (5-1) gives ( r-/ ) r = const. While as s — o and T becomes sufficiently large, the bubble tends to behave like T T constant. For the latter, one finds, a / ='2' and io 0,o Consequently the governing equation reduces to 3 di, ~6,d 4- O 2 d d which gives r; J- j ( 5-2) This behavior is also seen from Fig. (5-7), in which the bubble growth becomes a function of square-root time as 7 or r-> oo Another way, we may consider that equation (5-1) corresponds to the solution for the dynamics of a plane interface while

3.0 0 -- 0 2.8- l o a. l2.nd Order /' It 2.6 r h.. l 2.4( O 0 24 Order 1. D 2.0. Oth Order r=l- r (s Cn 1.86- w. 1.4 z LUJ 1.2 1.0 0 2 4 6 8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2x102 DIMENSIONLESS TIME rs = 4JA2a /T R2 t Figure 5-7. Convergence of bubble-dynamics solution for single-component system at small times.

-69equation (5-2) corresponds to the solution for the dynamics of a spherical interface. For the former case there is no relative convective effect to accompany the motion of the phase front. At large times in which Zs 7 ~,2 the bubble dynamics is dominated by the diffusion and axial convection. Equations (4-26:a) and (4-26b) can be integrated to yield 2 (r2- ) ) go//(v- z-, ) (5-3a) and r (- ). o - ) =/,o/ (/ - 0d) (5-3b) respectively, in which'?7 and' P7 are the time and bubble radius at the matching point of the small-and large-time solutions. For the case of? > b Z let ) and Y be the solutions of equations (5-3a) and (5-3b) respectively. Then T and 2are related as ir, \MHA _ X0fi k? p t I ]+ (5-4) Equation (5-4) shows that ars P/ as Since -— iI = O -5). one finds that as' 0! the series becomes rapidly convergent, Fig. 5-8 shows the convergence of the solution up to the -second-order approximations for several values of F at the matching point of the small- and large-time solutionso It is revealed from the figure that if f = 0,1 A 0,2, the solution up to the second-order approximation can give sufficiently accurate predictiono From equation (5-3a) and (5-3b), one can see that ro (I|+tt) or o o( c0+~i~)Y (5-5)

2.2 0o I "~~~/ o: / / j., 2.0- 2 2. / / B, -_ - -Zeroth-Order Approximation u(,8- Q / / / ^,./ First-Order Approximation 1.8 / /^^ <' ___.Second-Order Approximation w 1.6 m /'~' Also I //sn -cmcoponent at,t s ^'/^^ y^' ^^^^^ 1~~~~A 1.18 2.0 0.38 0 o ^ ^' B 1.184.00.19 z ^ ^ C I1. 4.00.24 1 0DI 1.14 10 0.08 5I I I 2 0 1.0 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 38 DIMENSIONLESS TIME, rt =(U0/R) t Figure 5-8. Convergence of Dubble-dynamics solution for single-component at large times.

-71This confirms that the coordinates ~ are indeed the asymptotic coordinates as ZtO —-~ if U can be expressed as t 1a t So long as ( as shown in equation (5-5), ^ can be arranged in a very rapidly decreasing sequence as Fig. 5-9 shows the solutions in the entire time domain for various values of the parameter &, As ---. D the matching point shifts to the right, indicating that the growth solution of a stationary bubble expresses that of a moving bubble for longer time, But as — > 00, the matching point shifts towards the origin, indicating that the axial convection begins earlier to dominate the bubble growth. It is obvious from Fig. 5-9 that the bubble growth behavior is completely different in these two time domains. Florschuetz and Chao (21) studied experimentally the bubble collapse in a single-component system (vapor bubble in subcooled water) under a reduced gravitational condition. Some theoretical predictions for the collapse of a stationary bubble obtained by means of the thin thermal boundary layer concept deviate from the experimental results, This deviation is attributed to the slight translatory motion of the bubbles which was observed during the free-fall. It is therefore, attempted to compare the present theory for the moving bubble with their experiments (21), The results are summarized in Table 5-1.

y:10.0 / 8.0 4.5 3.0 (p=0.09) (0.09) (0.10) (0.11) 0 Exact Solution 0 a 3.0 2.8 2.6 2.4 U ~2.2 ~m 2.0 0 1.8 cO |J 1.6 1.4 co LU 1.2 1.0 o.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 28 3.0 3.2 34 3.6 DIMENSIONLESSTIME, rs=4J a/7r Ro2 ) t Figure 5-9. Growth of moving bubbles in single-component system.

-73TABLE 5-1 RESULTS FOR EXPERIMENTAL BUBBLE COLLAPSE IN SINGLE-COMPONENT SYSTEM (21) Run Ja Approx. Pe -, iR 0o.s A 26,4 0,3 1.7 2.9 0,6 0,2 B 25.1 0.8 4.5 4.7 0.7 0.1 C 15,8 1.0 5,6 5.4 0.7 01o D 14.8 Normal 44.0 14.8 0.9 0,1 Gravity It should be noted that the bubble velocity U0o was measured when 0 = 0.8, If the plane-interface solution (5-1) is still valid at r = 0.8 j r may be obtained as The matching time and radius,,' and r1 were determined under the conditions that the parameter,o for large-time is equal or less t-hah 0,2 and the growth curve near the imatching time is smootho Professor Chao indicated in a personal communicatioh that the dimensionless bubble velocity UL 1 o In other wb0ds, t he bubble velocity was almost constant during the collapseo Therefore, the use of equation (4-25) is appropriate. The equation up to the second-order approximation for U = 1 is then integrated to give Run A~ r O -o.0 C/(r- 0o.)- 0&6(- o. ) Run Bo 2 o -o ifr-o.)-1 / (a- a /)

-74Run Co ri o. - o.// (a)r-o (-o f-Oo/ ) Run D: r. o39 0 (r-0 )-/(/(O —O/) The results are presented graphically in Fig. 5-10, which shows good agreement with the experiments (21), The effect of the bubble shape on the dynamic behavior is an important factor to be considered, since most bubbles deviate from a spherical shape. Consider a potential flow distribution around an ellipsoid. In order to use a uniform flow distribution in the present analysis, it is first examined what correction should be made to K in equation (3-28) for the spherical bubble. This will be examined under steady state as was proven for the spherical case, By defining W = a/b as the ratio of the major and minor axis, the potential function f is obtained as f=>?4rcooSo+ (WJ7 Y / ) (5-6) Where 3 (W;$ - 1-'w/' /) By employing a method similar to the analysis of the spherical case, Boussinesq(12) shows that, for thin boundary layer, total heat flux Qe over the entire surface of an ellipsoid is Q = S8 -r f 0 su-7) ivalent radius Re for the ellipsoid is where the equivalent radius R for the ellipsoid is eq

_ Experimental Results * Large - Time Solution up to the First- Order Approximation for Constant Velocity With U=l 1.0 [..-x Large- Time Solution up to the First - Order Approximation for the Case of u= r o.8 \ I oc \ 10. CollExact Solution For Exact Stationary Bubble (7y 0) A m (x(= 2.9) UI) LLJ -j.2 B co ^ D (D 4.7) ~ ( 7^~< =14.8) 0 -- 0.2.4.6.8 1.0.2.4.6 DIMENSIONLESS TIME, rs (4JAa /7r RO) t Figure 5-10. Collapse of moving bubble in singlecomponent system

-763 a/b For a sphere of radius R, equation (5-7) reduces to 0Qua - rT-T s6 <)(5-8) The analytical results obtained for a spherical bubble may be applied to an ellipsoidal one by the use of the equivalent radius R or equivalently using eq' L- 8 4i ___ (5-8a) TC 31B W in equation (3-28). The numerical values of K which is related to Kepg as k= TT ke (5-8b) are presented in Table 5-2 for several ellipsoidal bubbles of different W, TABLE 5-2 K AND KeLL o FOR ELLIPSOIDAL BUBBLES W 1o0 102 1.4 1l6 1,8 290 2.5 3.0 KeLL 1,0 1.07 1.15 1.240 1.33 1,41 1,62 1,83 K 2.52 2 72 2.93 3.16 3.39 3.60 4.10 4.66 Since the forward stagnation region of an ellipsoidal bubble widens as W increases, the uniform flow around the bubble would also increase correspondingly.

-77 - The experimental results obtained by Langston and Eustis (40) are now compared with the present analysis for the collapse of ellipsoidal bubbles initiated by a step increase in the system pressure. Since the saturation temperature corresponding to the final pressure was much lower than that of the vapor inside the bubble, the bubble (40) collapsed due to the loss of the internal: energy in the vapor rather than due to condensation. However, the present analysis may still be applied-to this situation provided t t T that the Ja number is replaced by the parameter - - Because of If 1, it is recessary to consider the solution up to the first-order approximation. However, the interface temperature varies since the gas temperature decreases with time, equation (4 -46) should be employed. It was reported in reference (40) that the relationship Tr= eiA:7t exists between T- TooTJ where is a constant. Therefore, the bubble-dyna-mics equation (4-46) up to the first-order approximation may be written as ~ d/1 -e; -,/ o r A -,' The prediction agrees very well with the experiments as may be seen from Fig. 5-11. B DYNAMICS OF A MOVING BUBBLE IN A BINARY COMPONENT SYSTEM. —--- INJECTION COOLING Equation (4-39) which describes the bubble growth at small times corresponding to (3->0 is integrated numerically.

1.00 I Solution up to First-Order - -o- - Experiments [40].98-, 0 Ao A.94 W <:: L. C 1.58 87.ii D.90 w w(n~~~~~~~~~ ~A 2.56 263 J z.88 B 1.54 145 Zn C 1.58 87 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.22.4 2.6 2.8 3.0 32 34 36 38 4.0 DIMENSIONLESS TIME, rt (UaDo /Ro) t Figure 5-11. Collapse of ellipsoidal bubbles.

-79The results are presented in Fig. 5-12 for a comparison with Larsen's solution (42) obtained by an integral technique. Equation (4-40) by Larsen is presented graphically in Fig, 5-13 for X = 0.01 and Xoo = 0.90 to demonstrate the effects of the parameter Q at small times. It is seen from the figure that as Q increases, the bubble growth decreasesrapidly. This is because the gas inside the bubble is nonsoluble to the surrounding liquid and consequently the bubble growth becomes dominated by mass diffusion in the bubble as increases. At large times where the effect of the axial convection in the liquid becomes important and the concentration of the gas phase inside the bubble reaches a practically uniform condition, equation (4-47) is integrated. For Xo = 0,01 and Xoo = 0~90 the effects of q, the bubble radius at the matching point of the two asymptotic solutions, and are illustrated in Pigs. 5-14 and 5-15o They indicate that the bubble growth becomes slower as' increases. The effects of the initial and final compositions, Xoo and Xo on the bubble growth are illustrated in Figso 5-16 and 5-17. For 0 = 0.2, Fig. 5-16 shows that the effects of Xo is entirely negligible while Figo 5-17 demonstrates a large effect due to X oo Two asymptotic solutions for small and large times are matched at two different times qx ='0O56 and q4 = 0o36 as illustrated in Pigo 5-18o It is disclosed from the figure that the matching point of the two solutions seems to affect the bubble dynamics very little as time elapses,

1.0 Lorsen E421.9 _8 m u _ Present Analysis' 7 C).6<[ cr r'5 Q3 Om.4 00 n3.3 IJ z.2 o ~2 0 Lz -1 - o0 0 1.0 2.0 3.0 4.0 x 102 DIMENSIONLESS TIME, r (4 Ja /T R )t Figure 5-12. Growth of moving bubble in injection cooling at small times for (3 —->

1.0 CO 0- 0.9i 0lJ Cn 50 2. DIMENSIONLESS TIME, r 3.0 4.OxI-2 or rlME- i -i n ea i 2R 4.0hlO Figure 5-13. V oo gb f fb e o fn' o n t h e g r o w t h o f m o v i n g bubble in injection cooling at small times For X= 0.01 and X0 0.90.

7 1.0,8:=0.20.9 ".8 C' D.7 cr.6 - ta.4w,3 3 -J.5 z o.2 C) z ~ * 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 DIMENSIONLESS TIME, Tt U /Ro t Figure 5-14. Effects of r on the growth of moving bubble in injection cooling at large times. For Xo = 0.01, Xoo 0.90 and qC = 0.436.

0 26 4, ^^ ^S^" --- =.56,.2 CO) LU LiJQD08 ~ ~ ~ ~ ~ ~ ~ ~ ~ - W I4I I (n 0 1.0 2.0 3.0 4.0 5.0 6.0 DIMENSIONLESS TIME, Tt = U0/Ro t Figure 5-15. Effects of 4 and qV on the growth of moving bubble in injection Pooling at large times. )~Cn~For Xo =. 0.90. 3 o.2 0 1.0 2.0 3.0 4.0 5.0 6.0 DIMENSIONLESS TIME., Tt R t Figure 5-15. Effects of r and on the growth of moving bubble in injection dooling at large times. For X0 = 0.01 and Xo, = 0.90.

.91. ".8 Cx X0.7 o~x Xo=0.03 o x x 0.02 <.6 x o 0.01 J 5.5 0 co x 3.4 in 00 (CD w.3 -J z o.2 z LmJ 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 20 2.2 2.4 2.6 2.8 3.0 32 34 DIMENSIONLESS TIME, rt = Uo/RO t Figure 5-16. Effects of X0 on the growth of moving bubble in injection cooling at large times. For Xoo = 0.90 and - = 0.2.

1.0 09- 9 X( —-' —0.90 0). TO.9;-0 *w -.II 0.2 I I 5 ~J 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 DIMENSIONLESS TIME, rt = UO/Ro t Figure 5-17. Effects of X oo on the growth of moving bubble in injection cooling at large times. For X0 = 0.003 and t = 0.2.

-._ Large-Time Solution ___.Experiment [421 Small-Time Solution ~l~'l/3= 7 7.8 [42] / -y. — = 1.0 mi 1-0,, ~~In-8- ~ ~~~ q9-a-0.36 _.2/y.8 - -.. 0L ——'-"1 cf / = 0-.6 I I _J1.4.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.236 cnfor X= 0.003 and X o = 0.56 r 0.2 oO 0.2.4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 DIMENSIONLESS TIME, =(lo x a/R 2)t Figure 5~18. Growth of moving bubble in injection cooling for Xo = 0.003 and X~o = 0.856.

-87The theoretical prediction agrees well with the experiments (42) over a considerable range of timeo C, DYNAMICS OF A MOVING BUBBLE IN A BOILING BINARY-COMPONENT SYSTEM At small times during which both diffusion and radial convection but not axial convection control the bubble growth, the present analysis would yield an identical result as a stationary bubble. However, for stationary bubbles the result obtained in the present analysis is different from that of Scriven (61) for stationary bubble in that the bubble may grow from a finite size, The bubble-dynamics equation (4-58) is integrated numerically for several values of o and /, The results are plotted in Fig, 5-19 for the bubble growth history and Figo 5-20 for the bubble-surface temperature, Figure 5-19 shows that as Jo and (3/ decrease, the bubble grows more rapidly. When (3 and G, approaches zero, the bubble dynamics becomes heat-transport controlled and approaches that of single-component system, The bubblesurface temperature as shown in Fig. 5-20 remains almost constant with time for large values of (0 and (9/ But as (3 and P3 decrease, the temperature T* tends to increase with time and asymptotically approaches unity at large times, Figures 5-21 and 5-22 illustrate the effects of o, and! respectively on the surface temperature and growth of moving bubble, P0 and (3/ play an identical role on bubble growth as. for the injection cooling case. An increase in either 3 or (3: signifies more mass-diffusion effect on the

o-0 =0o ^ / o0.1 // 0.3 2.4 / / 2.2 0.5 2. o~~~, ~~~~0.8 WD~~~~~~~ 1.48~~. Lm 1.4 -j1.4 0: I0 1.0 0 1.0 2.0 3.0 4.0 5.0 DIMENSIONLESS TIME,, =(4J, a /r R0) t Figure 5-19. Growth of moving bubble in boiling binarycomponent system at small times.

I.o - - - -- - -- o- - -:,o 8.8~~~~~~ ~0.3 8 ^"~~0.5'3~~~~~~~0 < 0.8,,.4 1.0 % - I1.5 1.8.2 2.4 _3.5 0 1.0 2.0 3.0 4.0 5.0 DIMENSIONLESS TIME, s =(4J 2a/ 0 R ) t Figure 5-20. Surface temperature ot moving bubble in boiling binary-component system at small time.

8.50 8 s.80= 0.5.= 40 1: —--------- 1.0 ".30 ___.20 -— 2.0 gr'.6.8 0.5P —o~ LU-i ~~~~~~~~~~~~~~~1.0 IJ 0' 5 ~ — D 1.4 C'.....~:t Ls'^^ or U) 2~~~~~~~~~~~~~~~~~.0 lmes to ft~ - 22.0 z 1.2 0 LU LaI'.0 DIMENSIONLESS TIM 2. 0 ip rs (4da, n I 2 5.0 Figure 0-21. e o growth of mov0 Ofl the Surface temperature and iof movng bubble in boiling binary. component system at mes o tsma. I times for (3.u

98 I/3= o0.0o 1^.50 - H.40 0.5 ~~~~~~.30 >1.0 1.5 - c/i.~~2.0.20 - 5 1.8, 1.6 w::) 1./4 w 2.0 DIMENSIONLESS TIME, T< (4J^/ /t R2 ) t 1o t 1f m vn.2b l n b i in i y.0 Figure 5-22. Effects of 9 n the surface temperature and growth of moving bubble in boiling binary-. component system at small times fox = 2.0.

-92growth and consequently lower growth rate. It is interesting to investigate how the gradients 0,o (/ in the phaseequilibrium diagram and the relative volatility (Xo, - X5) affect the bubble growth. As o(, or (cd -Oj)/d o increases for constant-4,the relative volatility between the liquid and gas sides of the bubble surface would cause the effect of mass transport to increase. Figure 5-23 illustrates the bubble-surface temperature, phase-equilibrium diagram, (3 and (3/ against the concentration for the ethyl glycol-water solution,, ( and (/ are calculated based on the phase-equilibrium diagram of the solution. The bubble-surface temipaeture T* isobtained from equation (4-53) and (4-54) when the interfacial conditions are at steady state, that is, T* = 0 and X* = 0 The figure shows that T* reaches the minimum value at approximately X - 0,95 where the bubble growth rate is minimum as analytically predicted by Scriven (61) and confirmed experimentally by Benjamin and Westwater (6)0 This is due to the characteristics of the phase-equilibrium diagram which determines the values of (o and Q/ As shown in Figo 5-23, I is maximum at X = 0,95. The effects of d on the growth of moving bubble in boiling binary-component system are illustrated in Fig. 5-24 for the case of (3o = (3/ 005 It is disclosed that the bubble growth rate increases with an increase in,.

-93U) 1.0.6 I-.4 200 0 w 180- / e 160 - 40, 140 -, 2 120 ao 100.8 -4.0 0.2.4.6.8 1.0 ETHYLGLYCOL IN MASS FRACTION Figure 5-23. Bubble surface temp rature, phase equilibrium diagram, and d I for moving bubble in boiling eth ylglycol-water solution.

2.4 =10.0o 8.0 4.0 2.0 2.2 0 0 2.0 1.8 i -1.6 CD f/ in:3 D () 1.4 Cf) LLd w -j z 0 Cn 1.2 z 1.0 0 1.0 2.0 3.0 4.0 5.0 DIMENSIONLESS TIME, r =(4Js a/7a R ) t Figure 5-24. Effects of c on the growth of moving bubble in boiling binar y-component system for O, = 0.5, and (= 0.5.

CHAPTER VI CONCLUSION The dynamics of a single moving spherical bubble in both single- and binary-component systems was analyzed for the case in which the mechanism is controlled only by the transport of energy and/or mass. With the application of the boundary layer simplification and the approximation of the velocity field around the bubble, the solutin of the transport equations of heat and mass was obtained asymptotically both for small and large times by employing the coordinate perturbation method. In the case of the singlecomponent system, the small-time solution is characterized by the dominance of diffusion and/or radial convection effects, while the large-time solution is subjected to the dominance of the diffusion and axial convection effects. The total heat transfer across the bubble surface is obtained in terms of an asymptotic series of the coordinates, f2 and F. It is then substituted into the energy balance condition at the bubble surface. The growth or collapse rate of the moving bubble is determined by the combination of the resulting equation and the appropriate bubble velocity. The prediction is accurate to the perturbation order retained for the heat flux solution. The solutions for both the heat flux and bubble dynamics at small and large times are found to converge very rapidly. -95

-96It is disclosed that the parameter r =,' ef / Jd which shows the relative importance of the axial convection to the radial convection controls the dynamics of the moving bubble in a single-component system, Whether the elapsed time is "small time" or "large time" is determined by the relative magnitude between ZS and I /2 e If s>. /2 the time corresponds to "large time"; if Zs< /&2, the time is designated as "small time". The bubble behavior is entirely different in these two time domains, However, the bubble behavior in these two time domains may be successfully joined at a certain time interval. The analytical results agree very well with the experiments (21, 40) for both spherical and ellipsoidal bubbles o For binary-component systems, the dynamics problem is complicated by the time dependency of the interfacial temperature and'concentrationO This difficulty was overcome through the use of the coordinate perturbation method. The transport equations of mass and energy are solved in terms of coordinates which include Ts(t) and X,(t) in addition to U s (t) and R(t). The bubble dynamics is determined by substituting the fluxes into the interfacial mass and energy balance conditions and then combining with thermodynamic equilibrium condition together with the appropriate relation of the translating velocity. Due to the complication, the solutions up to the first -order approximation are obtained for the mass and heat fluxes as well as the bubble dynamicso

-97In the case of injection cooling the effect of the axial ccnvection on the dynamics of moving bubble is negligible at small times. Therefore, the integral solution of Larsen (42) may be usedO At large times, the bubble dynamics is controlled only by heat transport in the liquid, since the gas concentration inside the bubble becomes nearly uniform due to its large diffusivityo However, the interfacial temperature and concentration vary with time following the thermodynamic equilibrium relationship. In addition to r, the parameter J< /U, which signifies the relative importance between the mass and. energy transport, affects the dynamics of moving bubble in injection cooling. It is disclosed that the bubble growth is controlled by the mass transport for large values of (3 and by the energy transport for small values of (3 The parameter' plays an identical role in the injection-cooling case as that in a single-component system, The analytical results give good qualitative agreement with the experiment (42)o For a boiling binary-component system, the analysis is performed for the growth or collapse of a moving bubble starting from a finite initial size. The bubble dynamics in this case is controlled by the simultaneous mass and energy transport in the liquid sideo The relative importance of the mass to energy transport is determined by.two parameters ~ ___,and f J P _ ___ Q, ddd\. Ja P X-,X oP / Oo —do, 0 ~"a d, In the 3 and [3 there is the common factor JW/fU. jr/P which corresponds to Q in the injection-cooling ease, In addition, the gradient of the phase equilibrium curve -oa dff

-98and the relative volatility (X - XI ) are included in the two parameters. In general,..the interfaci-al temperature and concentration vary with time. However, for large values of (3 and (3, j they remain almost constant as in the case of a stationary bubble (61).

APPENDIX I FLOW DISTRIBUTION AROUND A BUBBLE The physical model to be analyzed consists of a spherical bubble growing at the rate of R(t) in a free stream of an incompressible liquid with velocity U00 (t)o Following the analysis of Moore (49), we will show that the irrotational flow distribution is exact everywhere to O(Re.t), it is assumed that the velocity field within the boundary layer can be expressed as / U - A= - + / lp - + -y P = P + P where u, v and p are mean velocity components and pressure u' vV and pI are the time varying perturbation order velocity components and pressure in the irrotational flow field expressed by equations (3-17) and (3-18). From continuity equation (3-1) it is found a/= 00 However, from equation (3-21)'_ _ +. ti I _/ -99

-100Since -= 0 (S ) = 0(), = (/ t 0( ) ) therefore a -DH ) 2 In order that the viscosity effect may be of the same order as convective terms in equations (3-2), and (3-3), s3 O (R' ),) Therefore, U =- 0( ), ( ), =O0(/ ), o (/e ) or, the irrotational flow distribution of equations (3-17) and (3-18) is valid to 0 (1C )

APPENDIX II DETERMINATION OF THE CONSTANT K The interfacial boundary conditions by allowing an arbitrary time dependence on Uoo(t) and R(t) determine the unknowns Uoo(t) and R(t)o This scheme of analysis suggests that with regards to transport equations we are interested only in the total heat flux values over the bubble surface. Consider the energy equation uAT +:DT T.. (A-1) subject to the boundary conditions -i(c, oo~ )- T(co, ) =7T, ( 0 ) IS 7(A-2) where u and v are the potential flow distribution and may be expressed as vi- 3 ~ Z cos- - By defining 6 -2 ( l^ 3 _ _ _ _ _ _ _ -101

-102the equation and boundary conditions may be reduced to Q dO d2Q._ = J 2 d? d 7' e(o) I, 9 (00)=0 Total heat flux over the entire bubble surface can be obtained as 77 3 JL I -'(;Tm- /- e (A-3) Now one assumes that the flow over the sphere is uniform and may be expressed as K'U k Uli (and - 0) where K is a constant to be determined. For this case th-e total heat flux over the entire bubble surface can be obtained as = (T -T) ^ Pk2 (A-4) K is determined by the comparison of equations (A-4), (A-3) as, ___ _- (A-5) Consider the unstey e y Consider the unsteady energy energy equation T 2 T T 2 d- 2T (^A-6) subject to the initial and boundary conditions

-103= T(, oo,t )= T(o,, t )= T)(x,,) = Too Tcx, o,t )= Ts where u and v are the potential flow distribution written as u=- 3U, s(t) sir, =- 3C -(t ) c At small times, one defines d cR 00'' RL I (A-7) L,/~lj,., It \2 f), A ^, p,, S' -U00 A O where 3 /q which varies between ( O'-7 ) over the sphere. It is noted that UJf+)=-cD)as -; — 0 and therefore'tO(ct' o,0=o, o ) -) =o0(t$ —- ast- 0 With -the temperature fu-nction Q expanded as T-'To 0) the total heat flux over the entire bubble-surface can be obtained.as 2 I / l, T -T, j ()+ d6 ) ( -] 5-'si Since o0(o) cn (O) —- are constants and sin(3 and cos(3 are orthogonal in the interval ( 0, )T ), the contribution of I( term to the total heat flux is entirely negligible. Therefore, the effect of the sinusoidal velocity variation on the total heat flux must come from the dC, term which is O(t )

-104or 0() as may, be seen in Chapter IV. Consequently the value of K determined for a steady state may also be applied to the unsteady state within the error of O(t3) at small time s At large times the temperature field may be considered to be in quasi-steady state, that is the zeroth-order approximation represents the steady-state solution. The deviation of the unsteady behavior from the steady one is mainly due to the delay in the boundary-layer response to the variation of the free stream. For potential flow the velocity field is established everywhere by pressure wave propagation. Therefore, the potential flow field is established instantaneously. So the delay may come only from the response of the thermal boundary layer. The order of magnitude of such delay can be established by the following physical reasoning. /( is known as the diffusion time and indicates the order of magnitude of the time required for the energy (generally the vorticity) to diffuse through the thickness S of the boundary layer. The characteristic time due to flow unsteadiness is u Therefore, it is expected that the ratio of the two characteristic times indicates the order of magnitude of the delay in the boundary layer response. One now defines the ratio as = The boundary layer thickness q is wellknown for a flat plate and stagnation region~ S _ d for a flat plate UCQ

-105z —- ~ {for a stagnation region u0o From the definition of the similarity variable 7 for steady state solution', the energy boundary layer thickness may be found to be / ~ " l3d' 23/ Uoo dk1(3 Consequently If may take the following forms depending on the geometry., _ (aO ~P for a flat plate ( C iOL for a stagnation region " C-aI, 3_, e _o 3R jo 2 if for a sphere. It is important to note that f is oi_~bR for the three different geometries, provided that 3 is 0(1), Therefore, the large-time behavior deviates from its steady one by O( 0 o R The conclusion may also be derived from the following analysis, One defines U-3 U00) "2 R a0 3 sir y =_ U S [ 0 Q 3/.U1's~n(@ ) B U 3)'u,sLr( in which the parameters,& are in an asymptotic sequence at large times for a large class of function U o (t).

-106Equation (A-6) may then be rewritten as,,o (,^ o[ _s' I~ ^ 8,3I 5df2O _.~,;~ 3/S?9 ~ 7 The solution may be expressed in the power series as OCi( 7, F) - ) C9,p) t, + - - - (A.-8) where the functions 9 er mast satisfy the equations?, dqo d(2' - o lIe.,,,,:d% d9' 2 4 (I4 co-aJ sr ) 4 C) Ie4j 9 DQ F~,u'v3t 2 i7 Sk3?( 22 subject to the conditions Q (0) = 1, o( O ) = 0, 6(~8"p)- (S)= 60156 )D= C, The function Q is the zeroth-order approximation of the quasi-steady solution and is the steady state solution of equation (A-l)o Equation for Qo may be solved by the iterative method. The results are presented in Fig. A-lo The total heat flux over the entire bubble surface up to the first-order approximation c-an bq.expressed as __, 3YJs'V Y2 9Wo.. ark(T,~r,)4 tJO R;3Fy UGooF3 4 -'nR(Tz>-5)T R Pe )2 oS72 The integral of the equation was performed graphically in Fig, A-2.

%E (77R) =6 e810(0o,o) pEo(T,O) O 00-1% C OI. 0=0 ^ d10%* ** ^^ 0 0 75 X ft ^ ---- -— x —- — x —-------- -—,75.6 ce (?Z —-<r^' ^ ^ -950 --------- - -o -— ~ — ^- 90 ^.4 0100 110~ 0A 11n~ 40 *2 1200 750 X~~~1 1200 110 0 0 1.0 2.0 3.0 4.0 F7 Figure A-i. Funct ions e0Q(,9), 0ea,')1/os0., o) and ((, 0O'O,O)for equation (A6) II EC)~~~~~~~~~~

0.3 0- ('/) 0 0.2 1/2 -O. I - Ir o -0. 00 200 40~ 600 800 1000 1200 1400 1600 180 1, DEGREES Figure A-2. Functions O (01 0) and. ) 5;/~. -

-109If the velocity field is approximated by b — kb L and orfO, one obtains the total heat flux as = -7 TOO — T-,LA o/-J 0 From the comparison of the last two equations one finds that the result obtained from the latter expression deviates about 2% from that given by the former expression. It was shown in Chapter V that at a matching point where f is a maximum, eO.,2. Since the zeroth- and first-order approximations of the total heat flux are of the order of unity and E, respectively, two per cent error in the first-order approximation would contribute about 0,4% error to the total heat flux or equivalently to the bubble dynamics, Therefore, it may be concluded that the uniformflow approximation M=kU o and f= 0 may be employed for the determining of the bubble dynamics, in the steady as well as unsteady thermal and concentration fields,

-110APPENDI1X III COMPUTER PROGRAM S THIS PROGRAM COMPUTES THE UNIVERSAL FUNCTIONS FOR TRANSPORT EuLATIONS FUR ENERGY AND MASS. THE FUNCTIONS ARE VALID FORi buTH SMALL AND LARGL TIMES IF THE RELATIlON BETWEEN U_ LE RADIUS AND T.RANSLATJRY SPELD IS TO BE INDEPENDENT. THE COMPUTATION IS BASED ON ThE RUNGE-KUTTA METHCD. T.E.. 1__ I..$$ZEROTH _O.RD.E.R_.,s$ SOLUT I ON T.E.2 IS..$$ to ORDER'$ SOLUTION T E 3 I $ $ ORDER _ $$__ SOL U'_.0 N T.E.4 IS..$'C ORDER $$ SOLUTION T_.E.,5 _ IS,._$ O....... ORDER $$ SOLUTI G O T E.6 IS..$$ o ORDER $$ SOLUT I O T._E.5.7 IS_ $$. ORDER *$ SOLUTIOi T.E.8 IS.$ $ a ORDER $S SOLUTION T.E 9 IS,$$......... ORDER...SUOLUTION T.E.10 IS..$$' ORDER 5$ S3OLUTION T.E,....9IS..$$, R, r ORDER 2. SOLUtT I 0'l T.E.12 IS..S t' ORDER S$ LUTITCI0, T.E.13 IS.. $S $ ORDER US SOLuTIO' T.E.14 IS..B$ gO ORD!)R.5 SL JTI,.Pt T.E.16 IS..)$$ RSE "CORDR S )LTIiT T O E. 1 6 I S O.. rO S $ o JG R- D Ei S L;T I ),i T.E.17 IS.. S,&) lRDE R e SLuL TI'Ji T.E.18 IS eO$ (6')CRDER mS SuLUTlIO, T.E.19 IS.. S s, ORD-)ER bi S',JL TI',T.E.20 IS..$$ ORDER D. SOvLUT I'i FOR SMALL TIMES FOR LARGE TIiMES ~N ~} U fei)~ri()/ 9 ~} Ul...............(_ * (; I ~ -7 A ^TnU/) ) AT TV' AT \\J),7 Rj } 0, L) ( ~ jf —- la U & j j i ij) a (

-111$COMPILE MAD,PRINT OBJECTEXECUTE,DUMP MAD (03 JAN 1966 VERSION) PROGRAM LISTINJ...... DIMENSION Yt2),F(2),Q(2) INTEGER J I N,M,T,ZET,TEND,TMAX1,TMAX2 INTEGER GRS I Z',G RS I Z 1 G RS IZ 2 VECTOR VALUES DEL=$1H,F6.4,3F15.8*$ VECTOR VALUES EPS=$IH,1F 10.4*$ DIMENSION P(121),YI(121),Y2(121),F2(121),R(55) DIMENSION A( 3509,AD),B(3509,BD),C(3509,CD) VECTOR VALUES AD=2,1,120 VECTOR VALUES D=2-, 1,120 VECTOR VALUES CD=2,1,120 READ AND PRINT DATA EXECUTE SETRKD.(2,Y( 1),F(1),QX,STEP) T=1 CO=CO PRINT COMMENT $ FOR SMALL TIME BINARY SYSTEM $ TMAX =TMAX2 POTEN=O. GRS I Z1=GRSIZ GRSIZ2=GRSIZ1 SWST GRSIZ=GRSIZI............. SW6 ZET=O SW10 ZET=ZET+1 WHENEVER ((-1).P.ZET).L.0,YX=YDP+ZET-1 WHENEVER ((-1).P.ZET).G.0,YX=YDP-ZET+1 R(O)=1.5 R( 1)=1. N=1 SW2 Q-=CF*(R(N)/R (N-1 )-1) WHENEVER Q.G.87.20,Q=87.20 WHENEVER Q.L.(-87.20),Q=-87.20 DX=-DY*((2.3504/(EXP.(Q)-EXP.(-Q)))+1) DY=DX N=N+1 YX=YX+DY Y(2)=YX WHENEVER T.E.1 Y(1)=1. OTHERWISE Y(1)=0.'ZND OF CONDI IONAL J=O 1=0 I=O X=O START J=J+1 P(J)=X Yi(J)=Y( 1) Y2 ( J ) =Y(2 ) Y2(J)=Y(2) F2(J )=F(2) SW1 F(1)=Y(2) WHENEVER T.E.1, F( 2)=-0. 5*X*Y(2) WHENEVER T.E.2, F(2)=-.5*X*Y(2)+Y(1)-2.*X*B(1,J) WHENEVER T.E.3, F(2)= —.5*X*Y(2)+2.*Y(1) + A(2,J) WHENEVER T..E.4, F(2)=-.5*X*Y(2)+3.*Y(1) + A(3,J)

-112WHENEVER T.E.6, F(2)=-.5*X*Y(2)+Y(1)-..5*X*B(1,J) WHENEVER T.E.7, F(2)=-.5*X*Y(2)+'2.*Y(1)+A(6,J) WHENEVER T.E.10, F(2)=-.5*X*Y(2)+2.*Y(1)-2.*X*B(2,J)-A(2,J) WHENEVER T.E.11, F(2)=-.5*X*Y(2)+3.*Y(1) -X*b(3,J)+A(10,J) 0 -.5*A(3,J) WHENEVER T.E.12, F(2)=-.5*X*Y(2)+4.*Y(1) + A(11,J) *2. WHENEVER T.E.13, F(2)=-.5*X*Y(2)+2.*Y(1)+.5*B(6,J)*X-2.*A(6,J) WHENEVER T.E.14, F(2)=-.5*X*Y(2)+3.*Y(1)+.25*B(7,J)-1.5*A(7,J) 0 +A(13,J) WHENEVER T —.E.15, F(2)=-.5*X* Y( 2) +2.*Y(1) +'..5'*X*B(2,J )-A( 2,J)-2* 0 X*B(6,J) WHENEVER T.E.16, F(2)=-3.5*X*Y(2)+Y(1)+A(lJ) WHENEVER T.E.17, F(2)=-0.5*X*Y(2) +2.*Y(1)+A(16,J) WHENEVER T.E.18, F(2)=-0.5*X*Y(2)+2.*Y(1) WHENEVER T.E.19, F(2)=-O.5*X*Y(2)+2.*Y(1)-2.*X*B(16,J)+A(2,J) o - 1.,,,J) WHENEVER T.E.20,F(2)=-0.5*X*Y*Y2)+2.Y(1)+A(6,J)+0.5*X*B(16,J) S=RKDEQ.(0) WHENEVER S.E.1.0,TRANSFER TO SWI WHENEVER J.E.GRSIZ,TRANSFER TO SW3 TRANSFER TO START SW3 WHENEVER.A3S.Y(1).l,.0.0o05.AND. N.L.NMAX R(N)=Y(1) -RANSFER TO SW2 OTHERWISE TRANSFER TO SW15 END OF CONDITIONAL SW15 WHENEVER.ABS.Y(2).G..'01.AND.ZET.L.ZMAX,TRANSFER TO SW10 TRANSFER TO PTPT PTPT WHENEVER ZET.E.ZMAX,PRINT COMMENT $ EQ. CAN NOT BE MADE 0 SATISFY dOUNDARY CONDITIONS *$ PRINT RESULTS T PRINT RESULTS ZET WHENEVER N.E.NMAX,TRANSFER TO SW14 PRINT COMMENT $ R(M) $ THROUGH SW5, FOR M=2,1,M.G.N SW5 PRIIT FORMAT FPS,R(M) THRl)UGH SW4, FOR I=1,1,I.G.GRSIZ A(T, I)=Y1(I) 6(T,I)=Y2(I) SW4 C(T,I)=F2(I) THROUGH SW7,FOR I=1,4,I.G.GRSIZ SW7 PRINT FORMAT DEL,P(I),Yl(I),Y2(I),F2(I) SSS r=T+ WHENEVER T.E.5.OR.T.E.4.OR.T.E.8.OR.T.E.9,TRANSFER TO SSS WHENEVER T.G.GTEID,TRANSFER Or SW14 TRANSFER TO SW6 SW14 END OF PROGRAM THE FOLLOWING NAMES HAVE OCCURRED ONLY ONCE IN THIS PROGRAM. COMPILATION WILL CON[INUE. SWST *020 CF *028 GRSIZ2 *019 POTEN *017 STEP *012 TEND *094 TMAX1 ___ _*016 TMAX2 *016

-113THIS PROGRAM COMPUTES THE BUBBLE GROWTH OR BUBBLE COLLAPSE FOR BINARY SYSTEMS. THE COMPUTATION IS BASED ON THE $RUNGE-KUTTA$ METHOD. T.E.1...7.$$$ SMALL TIME SOLUTIONS FOR BINARY BOILING SYSTEM $$$ T*E.2.....$$$ LARGE TIME SOLUTIONS FOR BINARY BOILING SYSTEM $$$ T.E.3.....$$$ SMALL TIME SOLUTIONS FOR INJECTION COOLING SYSTEMS $$$ T.E.4.....$$$ LARGE TIME SOLUTIONS FOR INJECTION COOLING SYSTEMS $$$

-114$COMPILE MAD,PRINT OBJECT,EXECUTE,DUMP,PUNCH OBJECT MAD (03 JAN 1966 VERSION) PROGRAM LISTING........ DIMENSION Y(3),F(3),Q(3) INTEGER J,T,MAMB, TMAXTMIN,JM,J1,MO,MOO VECTOR VALUES DELl=$lH,F8.4,4F15.8*$ VECTOR VALUES DEL2=$1H,F8.4,5F15.8*$ READ READ AND PRINT DATA WHENEVER T.E.3, EXECUTE SETRKD. ( 1,Y(1),F(1),Q,X,STEP) WHENEVER T.E.4, EXECUTE SETRKD.(1,Y(1) F(1),QX,STEP) WHENEVER T.E.1, EXECUTE SETRKD.(2,Y(1),F(1),Q,X,STEP) WHENEVER T.E.2, EXECUTE SETRKD.(2,Y(1), F(1)Q,X,STEP) JM=JM TMAX=TMIN XSO=SQRT.(XINI) MO=10 X001=1.-XOO R003=( 1.0-XO)/XO01 BATO=( 1.-XO )/(XOO-XO) COFO=.5759*.5631/BETO COF1=0.5759*0. 2858/BETO COF2=. 5759*2815/BETO DOOFO=RAMO*BETO DOF1=DOFO/1.0159 DOF2=1. 3932/DOFO A00=(1.0 + BETO+ BET1)/BETC ROO=R003.P..33333 AGAIN J=O BAT=0.5631( 1.-BATO*0.5) BET=(1.0-BET1 )0.5631 RINI=1.0 + 2.*BET*XSO RANI=1.0 + 2.*BAT*XSO PRINT RESULTS ROO,BET,RINI,BAT,RANI WHENEVER T.E.1 Y( 1 )=RINI Y(2)=BET/XSO OR WHENEVER T.E.3.OR.T.E.4 Y( 1 )=RSTA OR WHENEVER T.E.2 Y(1)=RSTA Y(2)=VSTA END OF CONDITIONAL X=XINI MB=O START J=J+1 SW1 XSQ=SQRT.( X) WHENEVER X.L..005 STEP=O.0005 OTHERWISE STEP=STAP END OF CONDITIONAL WHENEVER T.E.2 YSQ=SQRT.(Y(1)) Y2=Y(1)*Y(1) Y32=Y(1)*YSQ AO1=1.059*Y2 + 0.702*Y(2)

-115TST=AOO-AO1/(DOF1*Y32 *Y(2)) F(1)=Y(2) F(2)=Y(2)*Y(2)*DOF1/YSQ*(-TST*AO1+RAMO*Y(2)*Y32+ 0 DOF2*(0.5079-1.143*Y (2)/Y2 )/YSQ) OR WHENEVER T.E.1 BCON=Y.(.2) *XSQ BCON2=BCON*BCON RIR=Y(2)/Y( 1) RIRT=RIR*X A01=0.5631 + 0.5710*RIRT TST=AOO -AO1/(BETO*BCON) F( 1)=Y(2) F ( 2) =( COF 1*RIRT-COF2+BCON2-AO.L.*TST*BCON)/COFO-*Y (2)/X OR WHENEVER T.E.3.OR.T.E.4 Y2=Y(1)*Y( 1) Y3=Y2*Y(1) YY3=R003/Y3 TST= XOOI-(YY3-1.0) WHENEVER T.E.3 F(1)=O.5631*TST/(XSQ + 1.728*XO01*X*YY3/Y(1)- 0.5716*X/Y(1) 0 *TST) OTHERWISE F(1)=0.5631*TST/XSQ END OF CONDITIONAL END OF CONDITIONAL RREL=(Y(1)-1.O)/(ROO-1.0) SW2 S=RKDEQ.(O) WHENEVER S.E.1.OTRANSFER TO SW1 WHENEVER T.E.2.OR.T.E.1 WHENEVER J.G.GRSIZ PRINT RESULTS TX,Y(1),Y(2),F(2),TST,RREL TRANSFER TO SW5 END OF CONDITIONAL OTHERWI SE WHENEVER Y(1).G.RMAX.OR.Y(1).L.RMIN.OR.J.G.GRSIZ SW6 PRINT RESULTS TXY(1),Y(2),TSTRREL TRANSFER TO SW5 END OF CONDITIONAL END OF CONDITIONAL J1=J WHENEVER X.L..01 MA=MOO OTHERWISE MA=MO END OF CONDITIONAL WHENEVER J1/MA.E.MBTRANSFER TO START MB=J1/MA WHENEVER T.E.4,PRINT FORMAT DEL1,X,Y(1),F(1),TSTRREL WHENEVER T.E.3,PRINT FORMAT DEL1,X,Y(L),F(1),TSTRREL WHENEVER T.E.1,PRINT FORMAT DEL2, X,Y(1),Y(2),F(2),TST,RREL WHENEVER T.E.2,PRINT FORMAT DEL2, X,Y(1),Y(2),F(2),TSTRREL TRANSFER TO START SW5 TRANSFER TO READ END OF PROGRAM THE FOLLOWING NAMES HAVE OCCURRED ONLY ONCE IN THIS PROGRAM. COMPILATION WILL CONTINUE. AGAIN *025 SW2 *078

REFERENCES 1. Acrivos, A., "Mass Transfer in Laminar Boundary Layer Flow with Finite Interfacial Velocity." J. A.I.Ch E., 1960, p. 410 2. Alen, H.S., "The Motion of a Sphere in a Viscous Fluid." Phil..ag., Vol. 50, 1900, p. 323 3. Baid, M.H.I. and Davidson, J.F., "Gas Absorption by Large Rising Bubbles." Chem. Eng. Science, Vol. 17, 1962, pp. 87-93 4. Barlow, E.J. and Langlois, W.E., "Diffusion of Gas From a Liquid into an Expanding Bubble." J.I.B.M., Vol. 6, 1962, p. 393 5. Beek, W.J. and Kramers, H., "Mass Transfer with a Change in Interfacial Area." Chem. Eng. Science, Vol. 17, 1962. pp. 909-921 6. Benjamin, J.,E. and Westwater, J.W., "Bubble Growth in Nucleate Boiling of a Binary Mixture." International Development in Heat Transfer, Part 2,.1962, pp. 212-223 7. Birkoff, G., et. al., "Spherical Bubble Growth." Physics of Fluid, Vol. 1, No. 3, 1958, p. 201 8. Blassius, H., Dissertation, Gottingen, 1907 9. Blassius, H., "Grenzschichten in Flussigkeiten mit Kleiner Reibung." Zeit. Math. u. Phys., Vol. 56, 1908, pp. 20-37 10. Bond, W.N., "Bubbles and Droplets and Stokes Law." Phil. Mag., Vol. 4, 1927, pp. 889-898 11. Bond, W.N. and Newton, D.A., "Bubbles, Drops and Stokes Law." Phil. Mag., Vol. 5, 1928, p. 794 12. Boussinesq, M.J., "Calcul du Pouvoir Refroidissant des Courants Fluids." J. Math., Vol. 1, 1905, p. 285 13. Boussinesq, M.J., "Vittesse de la Chute Lente Devenue Uniform, Goutte Liquide Spherique, dans un Fluide Visquex de Poids Specifique Moindre." Compt. rend., Vol. 156, 1913, p. 1124 14. Chambre, P.L., "On the Dynamics of Phase Growth," Quart. J. Mech. & App. Math., Vol. 9, 1956, p. 224 -116.

-11715. Chang, I.D., "Navier-Stokes Solutions at Large Distances from a Finite Body." J. Math. Mech., Vol. 10, 1961, pp. 811-876 16. Chao, B.T., "Motion of Spherical Gas Bubbles in a Viscous Liquid at Large Reynolds Number." Physics of Fluid, Vol. 3, 1962, pp. 69-79 17. Datta, R.I., Napier, D.H. and Newitt, D.M., "The Properties and Behavior of Gas Bubbles Formed at a Circular Orifice." Trans. Brit. Inst. Chem. Eng., Vol. 2b, 1950, p. 14 18. Davidson, J.F. and Schuller, O.G., "Bubble Formation at an Orifice in a Viscous Liquid." Trans. Brit. Inst. Chem. Eng. 1960, pp. 144-154 19. Davis, R. and Taylor, G.I., "The Mechanics of Rising Bubbles." Proc. Royal Soc. A., Vol. 200, 1950, pp. 375-390 20. Dergadebidian, I., "Observations on Bubble Growth in Various Superheated Liquid," J. Fluid Mech., Vol. 9, 1960, p. 39 21. Florschuetz, L.W. and Chao, B.T., "On the Mechanics of Vapor Bubble Collapse — A Theoretical and Experimental Analysis." Univ. of Illinois, Mech. Engr. Dept., 1963 22. Forster, H.K., "On Conduction of Heat into a Growing Vapor Bubble." J. App. Physics, Vol. 25, No. 8, 1954. p. 1067 23. Forster, H.K., "Diffusion in a Moving Medium with Time Dependent Boundaries." J. A.I.Ch.E., Vol. 2, 1957, p. 535 24. Fisch, S., "Steady-State Diffusion into a Streaming Sphere at Low Re. No." J. Chem. Phys., 1954, pp. 123-135 25. Frossling, N., "Evaporation, Heat Transfer and Velocity Distribution in Two Dimensional and Rotationally Symmetrical Laminar boundary Layer Flow." NACA T. M. 1432, 1958 26. Garner, F.H. and Hammerton, D., "Circulation Inside Gas Bubbles." Chem. Lng. Science, Vol. 3, 1954, p. 1 27. Gaywood, J. and Birkhoff, G., "Fluid Flow Pattern." J. Appl. Phys., Vol. 20, 1949, pp. 646-659 28. Goldstein, S. and Rosenhead, L., "Boundary Layer Growth." Proc. Cambridge Phil. Soc., Vol. 32, 1936, pp. 392-401

-11829. Groothius, H. and Kramers, H., "Gas Absorption by Single Drops during Formation." Chem. Eng. Science, Vol. 4, 1955, pp. 17-25 30. Haberman, W.L. and Morton, H.K., "An Experimental Investigation of the Drag and Shape of Air Bubbles Rising in Various Liquid." David Taylor Model Basin, Rept. No. 802 31. Hadamand, J., "Mouvement Permanent Lent d'une Sphere Liquide et Visquese dans un Liquide Visqueux. " Comptes Rendes, Vol. 152, 1911, p. 1735 32. Han, C.Y. and Griffin, P., "The Mechanism of Heat Transfer in Nuclear Pool Boiling." M.I.T. Report, 1962 33. Hartunian, R.A. and Sears, W.R., "On the Instability of Small Gas Bubbles Moving Uniformly in Various Liquids." J. Fluid Mech., Vol. 3, 1957, p. 27 34. Higbie, R., "The Rate of Absorption of a Pure Gas into a Still Liquid During Short Periods of Exposure." Trans. A.I.Ch.E., Vol. 31, 1935, pp. 365-389 35. Howarth, L., "On Calculation of the Steady Flow in the Boundary Layer Near the Surface of a Cylinder in a Stream." Aeron. Res. Coun., Rep. & Mem., No. 1632, 1935 36. Hsu, C-J, "Heat Transfer to Liquid Metals Flowing Past Spheres and Elliptical-Rod Bundles." Int. J. Heat and Mass Transfer, Vol. 8, 1965, pp. 303-315 37. King, L.V., "On the Convection of Free Heat from Small Cylinders in a Stream of Fluid: Determination of the Convection Constants of Small Platinum Wires with Applications to Hot-Wire Anemometry." Proc. Roy. Soc., Vol. 90, 1914, p. 563 38. Kronig, R. and Brink, J.C., "On the Theory of Extraction from Falling Droplets." J. APP. Science Res., Ser. A, Vol. 2, 1950, p. 142 39. Kronig, R, Van Der Veen, B. and Ijzerman, P., "On the Theory of Extraction from Falling Droplets II." J. App. Science. Res., Ser. A, Vol. 3, 1952, pp. 103-110 40. Langston, L.S. and Bustis, R.H., "Heat Transfer from a Non-Spherical Bubble Rising in an Isothermal Liquid." Department of Mechanical Engineering, Stanford University, 1964

-11941. Larsen, P.S. et al., "Cooling of Cryogenic Liquids by Gas Injection." Advances in Cryogenic Engineering, Plenum Press, Vol. 8, 1963 42. Larsen, P.S., "The Dynamics of Gas-Vapor Bubbles in Binary Systems." Pb.D. Thesis, University of Michigan, 1963 43. Leonard, J.H. and Houghton, G., "Effect of Mass Transfer on the Velocity of Rise of Bubbles in Water." Nature, London, 1961, p. 687 44. Levich, V.G., Pysicochemical Hydrodynamics. English Translation, Prentice Hall Inc., 1962 45. Lighthill, M.J., "The Response of Laminar Skin Friction and Heat Transfer to Fluctuations in the Stream Velocity." Proc. Roy. Soc., Vol. A 224, 1954, pp. 1-23 46. Lunnon, R.G., "Fluid Resistance to Moving Spheres." Proc. Roy. Soc., A. 1926, p. 303 47. Miyagi, 0., "The Motion of a Bubble Rising in Water," Tech. Repts, Tohoku Imperial Univ., Vol. 5 1925, p. 135 48. Moore, D.W., "The Rise of a Gas Bubble in a Viscous Liquid." J. Fluid Mech., Vol. 6, 1959, pp. 113-130 49. Moore, D.W., "The Boundary Layer on a Spherical Gas Bubble." J. Fluid Mech., Vol. 16, 1963, pp. 161-175 50. Moore, F.K., "Unsteady Laminar Boundary-Layer Flow." NACT TN 2471, 1952 51. Moore, F.K., "The Unsteady Laminar Boundary Layer of a Wedge and a Related Three-Dimensional Problem." Heat Transfer and Fluid Mech. Inst., 1957, p. 397 52. Nielson, A.E., "Diffusion Controlled Growth of a Moving Sphere. The Kinetics of Crystal Growth in a Pottassium Perchlorate." J. Chem. Phys., Vol. 65, 1961, pp. 46-49 53. Plesset, M.S. and Zwick, S.A., "A Non-Steady Heat Diffusion Problem with Spherical Symmetry." J. App. Phyl 1952, p. 95 54. Pouhlhausen K., "Zur naherungsweisen Intergration der Differutialgleichung der Laminaren Grenzschicht " Z.A.M.M., Vol. 1, 1921, p. 252

-12055. Ramshaw, C. and Thornton, J.O., "Droplet Circulation and Interfacial Disturbances in Gas Liquid System." Nature, London, Vol. 184, 1959, p. 719 56. Rayleigh, Lord, "On the Pressure Developed in a Liquid during the Collapse of a Spherical Cavity." Phil. Mag., Vol. 34, 1917, p. 94 57. Rosenberg, B., "The Drag and Shape of Air Bubbles Moving in Liquids." David Taylor Model Basin Rep't 727, 1953 58. Rybczynski, W., "Uber die fortshrietende Bewegung einer flussigeis Kugel in einen zahen Medium." Bull. Acad. Sci. Crac. 1911, pp. 40-46 59. Saffman, P.G., "On the Rise of Small Air Bubbles in Water." J. Fluid Mech., Vol. 1, 1956, p. 249 60. Schlichting, H., Chapter XI, Boundary Layer Theory, McGraw-Hill Book Company, 1960 61. Scriven, L.E., "On the Dynamics of Phase Growth." Chem. Eng. Science, Vol. 10, 1959, pp. 1-13 62. Sibulkin, M., "Heat Transfer near the Forward Stagnation Point of a Body of Revolution." J. Aero. Sco., 1952, p. 570 63. Sparrow, E.M. and Gregg, J.L., "Non-Steady SurfaceTemperature Effects on Forced Convection Heat Transfer.' J. Aero. Sci., Vol. 24, 1957, pp. 776-777 64. Stokes, Sir G., Mathematical and Physical Papers, _, _...._. _ Cambridge Univ. Press 65. Stuart, J.T., Chapter VII, Laminar Boundary Layers, edit. by L. Rosenhead, Oxford Press, 1963 66. Tadaki, S. and Maeda, S., "Shape and Velocity of Single Air Bubbles Rising in Various Liquid." Kagaku Kogaku, Japan, Vol. 25, 1961, pp. 254-264 67. Tokuda, N., "A New Perturbation Scheme by Generalized Coordinate Expansions." (forthcoming), 1966 68. Tokuda, N. and Yang, W-J., "Unsteady Stagnation Point Heat Transfer Due to Arbitrary Timewise Variant Free Stream Velocity." Proc. 3rd Int. Heat Trans. Conf., pp. 223-230

-12169. Van Dyke, M., Perturbation Methods in Fluid Mechanics, Academic Press, 164 70. Yang, K.T., "Unsteady Laminar Boundary Layers in an Incompressible Stagnation Flow." J. App. Mech., Vol. 25, 1958, pp. 421-427 71. Yang, W-J and Clark, J.A., "On the Application of the Source Theory to the Solution of Problems Involving Phase Change. Part I. Growth and Collapse of Bubbles." J. Heat Transfer, Trans. ASME, Series C, Vol. 86, No. 1 1963 72. Young, N.O., Goldstein, J.S. and Block, M.J., "The Motion of Bubbles in a Vertical Temperature Gradient." J. Fluid Mech., Vol. 6, 1959, pp. 321-349 73. Zuber, N. and Forster, H.K., "Dynamics of Vapor Bubbles and Boiling Heat Transfer." J. AIChE, Vol. 1, 1955, p. 531.

UNIVERSITY OF MICHIGAN 3 9015 03525 2306