Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00003429/00001
## Material Information- Title:
- A computational approach for two-dimensional, gas-liquid flows with surface tension and viscosity
- Creator:
- Khunatorn, Yottana
- Publication Date:
- 1996
- Language:
- English
- Physical Description:
- viii, 74 leaves : illustrations ; 29 cm
## Subjects- Subjects / Keywords:
- Surface tension ( lcsh )
Two-phase flow ( lcsh ) Surface tension ( fast ) Two-phase flow ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaves 73-74).
- General Note:
- Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Mechanical Engineering.
- Statement of Responsibility:
- by Yottana Khunatorn.
## Record Information- Source Institution:
- |University of Colorado Denver
- Holding Location:
- |Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 37757065 ( OCLC )
ocm37757065 - Classification:
- LD1190.E55 1996m .K48 ( lcc )
## Auraria Membership |

Full Text |

A COMPUTATIONAL APPROACH FOR
TWO-DIMENSIONAL, GAS-LIQUID FLOWS WITH SURFACE TENSION AND VISCOSITY by Yottana Khunatom B.Eng., Chiang Mai University, 1992 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical Engineering 1996 This thesis for the Master of Science degree by Yottana Khunatom has been approved William H. Clohessy V3 U Date Khunatom, Yottana (M.S., Mechanical Engineering) A Computational Approach for Two-dimensional, Gas-Liquid Flows with Surface Tension and Viscosity Thesis directed by Assistant Professor Sam W. Welch ABSTRACT Two-phase fluid flows including surface tension effects are studied using a numerical method. An existing finite difference method with staggered-grid and interface tracking are used as our basic numerical method. The surface tension force is expressed in terms of the curvature of the interface by considering the surface tension as a force that acts on the interface. The curvature is found by geometrical method. The discrete form of the surface tension force is found by interpolation to the computational cell. The viscous terms are calculated with consideration of the discontinuity in viscosity across the interface. Simulations of two-dimensional buoyancy- driven bubbles are presented as a numerical result with surface tension and viscosity. In addition, an oscillating bubble and the Rayleigh-Taylor instability, which are inviscid problems, are attempted to test the applicability of our numerical approach. The numerical results from the buoyancy-driven bubbles concur very well with the literature. The Rayleigh-Taylor instability simulation gave reasonable results in both unstable and stable wave number regions. The results from oscillating bubble are accurate only for the first half cycle and tend to deteriorate afterward. This, we believe, is due to the numerical dissipation attendant with the donoring convection scheme. This work not only shows the accuracy of our approach, but also illustrates interesting fluid mechanical behavior in the simulations presented. 111 This abstract accurately represents the content of candidates thesis. I recommend its publication. Signed Sam W. Weld IV CONTENTS List of Figures......................................................vi Acknowledgments....................................................viii Chapter 1. Introduction.......................................................1 1.1 Motivation.....................................................1 1.2 Scope of Study.................................................2 2. Numerical Method...................................................4 2.1 The Governing Equations........................................4 2.2 The Numerical Procedure........................................6 2.3 The Finite Difference Method...................................6 2.4 The Gauss-Seidel Iteration.....................................9 2.5 Tracking the Interface........................................12 2.6 Density Calculation.......................................... 16 2.7 Viscous Stress Calculation....................................17 2.8 Surface Tension and Curvature.................................22 2.8.1 Surface Phenomena..........................................22 2.8.2 The Approximation of Surface Tension Effects...............26 2.8.3 The Calculation of Curvature...............................30 2.8.4 The Numerical Method...................................... 32 3. Numerical Results.................................................34 3.1 Buoyancy Driven Bubbles...................................... 34 3.1.1 Theoretical Background.....................................34 3.1.2 Numerical Result...........................................38 3.2 Oscillating Bubble............................................45 3.2.1 Theoretical Background.....................................45 3.2.2 Numerical Result...........................................49 3.3 Rayleigh-T aylor Instability..................................51 3.3.1 Theoretical Background.....................................51 3.3.2 Numerical Result...........................................59 3.4 Conclusions and Recommendations...............................72 Bibliography ........................................................73 v LIST OF FIGURES Figure Page 2.1. Location of all parameters on a staggered cell..................7 2.2. Geometry in interface cell.....................................15 2.3. Mass fluxed to the right cell by velocity u....................17 2.4. Void and viscosity index.......................................19 2.5. Interface displaced through distance dQ, with pressure pi and p2 in fluids 1 and 2 respectively.............23 2.6. The surface tension force at the interface on each side of the control volume..........................................26 2.7. Geometry of curvature calculation..............................31 3.1 Graces graphical correlation..................................37 3.2 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 1..................................40 3.3 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 2..................................41 3.4 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 3..................................42 3.5 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 4..................................43 3.6 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 5..................................44 vi 3.7 Oscillating bubble at start, half period, and full period.....50 3.8 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 1 at time 0.15 sec...........61 3.9 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 2 at time 0.15 sec...........62 3.10 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 3 at time 0.15 sec...........63 3.11 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 4 at time 0.15 sec...........64 3.12 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 1 at time 0.25 sec...........65 3.13 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 2 at time 0.25 sec...........66 3.14 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 3 at time 0.25 sec...........67 3.15 Interface shape, velocity and pressure field Rayleigh-Taylor instability case 4 at time 0.25 sec...........68 3.16 Interface shape, velocity field Rayleigh-Taylor instability case 1 at time 0.45 sec...........69 3.17 Interface shape, velocity field Rayleigh-Taylor instability case 2 at time 0.45 sec...........70 3.18 Interface shape, velocity field Rayleigh-Taylor instability case 3 at time 0.45 sec...........71 vu ACKNOWLEDGMENTS I would like to express my sincerest thanks Sam W. Welch for the thesis topic and his guidance. I would also like to thank to my family for their love and encouragement. Finally, I would also like to thank Piradee Sograi and all my friends for their support and encouragement. vui 1. Introduction 1.1 Motivation Fluid interfacial motion induced by surface tension plays a fundamental role in many natural and industrial phenomena. Studies of capillarity, low-gravity fluid flow, hydrodynamic stability, surfactant behavior, cavitation, and droplet dynamics in clouds and in fuel sprays used in internal combustion engines are all examples of phenomenon associated with surface tension. Surface tension is caused by the following mechanism. The molecules of a substance usually attract one another. These attraction forces are different for different substances. Hence at the interface surface that separates the two media, there is a surface force because of the unsymmetrical attraction force. The motion of fluids with interfaces is a moving boundary problem. Generally, moving boundary problems are complicated and difficult to solve analytically. These problems not only deal with the transfer of mass, momentum and heat, but they must also deal with the development of the interfacial shape and dynamics. Due to the complexity and nonlinearity of these problems, there are a number of approaches using analytical, 1 experimental, and numerical methods to understand the physics of .these flows. This thesis describes a numerical method with which to study a class of these flows. 1.2 Scope of Study In this work, we focus on using numerical models in studying the effect of surface tension on two-phase fluid flows in two dimensions. The calculations are done under the assumption of an isothermal system. Previous to this work, a finite difference method with a staggered-grid (MAC type grid) was developed as the basic numerical scheme. An interface tracking method based on the Volume-of-Fluid (VOF) method was implemented [1]. In this work we have added numerical models for surface tension and viscosity for gas-liquid flows. Surface tension is treated as a surface force that acts on the interface. In our method it is necessary to calculate the curvature of the interface. This curvature calculation is accomplished by fitting the osculating circle to the interface at the point where the curvature is needed. Then the discrete form of the surface tension force is obtained by simple interpolation onto the computational grid. The sample simulations presented in this work are 2 The buoyancy driven bubble Oscillating bubble Rayleigh-Taylor instability The first of these simulations is for viscous fluids. Therefore, a simple approach to treating the spatial dependent of viscosity is implemented. The second and third simulations consider inviscid fluids. All three simulations have the surface tension as an important parameter that affects the physics. 3 2. Numerical Method In this chapter we present a numerical method which may be used to study two-phase flows in which viscosity and surface tension are important. The method is based on the volume of fluid (VOF) method that is widely used for treating free surface flows. The existing interface tracking method of Mortensen, Trapp, and Welch [1] is extended to include surface tension and discontinuous viscosity. First, we present the governing equations in differential form including a surface tension force. Then we review the pre-existing flow solver and interface tracking algorithm. Lastly, we present the addition of the viscous stress terms and the surface tension force. 2.1 The Governing Equations In this work, we consider two dimensional, transient, incompressible, viscous, and immiscible fluids. The governing equations of mass and momentum conservation with the surface tension force added are Conservation of Mass du dv dx + dy (2.1) 4 Conservation of Momentum x-direction da du da -----1- u------1- v dt dx dy 1 dP 1 d da\ Id p dx p dx\ dx) p dy da chi' --- + -- dy dx) (2.2) y-direction dv dv dv at + + v 5^ 1 8P 1 8 P dy p dy ( dv) 1 8 f du Sv'j 2P V dy) + P dy ^ dy dx) + 9y (2.3) where w and v are the components of the velocity in x and y direction respectively. The density is denoted by p, and P refers to the pressure. The body force components in each direction are gx and gy. The second and the third terms on the right hand side of the momentum equations are the viscous stress terms. Note that the viscosity is a function of position. Note also that the surface tension force components (Fa and Fsy ) are to be modeled as volume forces [2], 5 2.2 The Numerical Procedure The basic finite difference method is a variant of the traditional MAC method of Harlow and Welch [3], The velocity components and pressure are defined on an Eulerian, staggered grid. At each time step, finite difference approximations of the governing equations are constructed. The resulting linear system of equations is solved by employing a Gauss-Siedel iteration. An interface tracking procedure has been added and this thesis describes the addition of surface tension and discontinuous viscosity. A description of the numerical method is as follows: Approximate new time velocities are determined by taking an explicit step in the momentum equation. The surface tension force and the viscous stresses are treated explicitly during this step. Using Gauss-Seidel iterations on the mass conservation and momentum equations solve for the new time pressure and the exact new time velocity. Determine the new time density from the mass flux which flows into or out of each cell with the new time velocity. Track the new interface by the subgrid counting procedure. 2.3 The Finite Difference Method As described above, the finite difference discretizations are based on an Eulerian mesh. A staggered (MAC) grid is employed. In the staggered grid, the velocity components are calculated for the points that lie on the faces 6 of the cell, while all other variables are calculated at the cell center. This is shown in Fig. 2.1. Location of all parameters on a stagger cell Figure 2.1 The finite difference approximations of the governing equations are 2n+1 l>3 - Uj- Ax, hi vn+1 hj AVj (2.4) At = -(fuxx) (FUXY) + gx - avg Pi+l,j Pi,j Ax + (fsx).'. + (viscxx + viscxy) Paver Pa avg (2.5) 7 1 v"+1 V. . 1,3 At where = -(FUYX) (FUYY) + gy - avg Pj,j+i Pi,j AyiH 1 1 + ----(fsy) + ----(VISCYX + VISCYY) Pavg Pavg (2.6) FUXX = 2Axjti + U"w) + - ui.j) + u) + )_ Fray = + "*0+ h*i.i|Kj ui.j) + i.l) + W|"lH i.i) (2.7) FUYX = 2Ax, L *3+s + + kj.!^ V*Ul) - vul) FUYY = 2A>V * + rM*>) + kj^lKj - + ^.j) + V) (2.8) 8 The cells are numbered by indices i and j. These indices indicate cell center position along the horizontal and vertical directions, respectively. Cell boundary positions are labeled with half-integer values of the indices. The dimensions of the rectangular cells are Ax and Ay. In contrast to the other variables, which use the same index as the cell numbering, the velocity components, u and v, are indexed by the same index number but at the boundary of the cell. The indexing scheme is shown in Fig. 2.1. In addition to the space index subscripts, we use a superscript to number the time cycle. For example, u?*1 designates the horizontal velocity at the time t = (n+J)At, in which At is the time increment per cycle. When the cycle number superscript is omitted, it is assumed that its value is n, corresponding to the time t = (n)At. The average density, pavg, is used in the difference equations due to the presence of two-phase cells. The calculation of this average density is discussed in the following section. Because of discontinuity of the viscosity across the interface, our differential approximation of the viscous terms requires special treatment. This will be discussed in detail in a later section. 2.4 The Gauss-Seidel Iteration The Gauss-Seidel iteration is a simple but efficient point-iterative procedure for a large system of equations. Convergence is related to the condition of diagonal dominance of the matrix of coefficients [4], The procedure for a general system of algebraic equations is 9 Make an initial guess for all unknowns. Solve each equation for the unknown whose coefficient is on the diagonal, using guessed values initially and the most recently computed values thereafter for the other unknowns in each equation. Repeat iteratively the solution of the equations in this manner until changes in the unknowns become small. The exact new time velocities must satisfy the momentum equations with the new time pressures. The discrete momentum equations in the x- direction, Eq.(2.5), may be written u;1 = ultj At(FUXX) At(FUXY) + gx - At avgr pn+l ____ pn+l ci+1.3 *1.1 Ax i+i + (Fsx). + (VISCXX + VISCXY) avgr avg- (2.9) Introducing the pressure increment, AP, we have At Tpn+1 pn+ll _ [*1+1.1 *i,1 J ~ i + 1 At Pavg^1 + i 2 At Pav^X. , [Pi+i,3 Pi,j] [APiVi Api"j+1] (2.10), combining Eq. (2.9) and (2.10) we obtain 10 Un+1 = uexp i/J hi At [a- AP"1] (2.11), where u*xp = Uifj At(FUXX) At(FUXY) + gx - At avg- P P Ax. _ + (psx). + (VTSCXX + VISCXY) avg avg (2.12). Using the same procedure as above in the y-direction, we obtain vn+1 = vexp - 2.3 At [ap,:;.1, ap^;1] (2.13), where ~ At(FUYX) At(FUYY) + gy - At P avg 1 Ay. j J j+I + (psy). (VJSCifX + viscyy) At ^ ''"10 p avg (2.14). And the continuity equation is ^ + ytf Ax, Ay, = 0 (2.15) 11 Substituting Eq.(2.11) and Eq.(2.13) into Eq.(2.15) we obtain <7 - ~ Ax, Ay,- (2.16) In the procedure above, the initial guess of the new time velocities are found by solving the momentum equations with old time pressure values. These velocities are denoted by uexp in Eq.(2.12) and ve!tp in Eq.(2.14). The Gauss-Seidel iteration is used to solve Eq. (2.16) for the pressure increments. This pressure increment then used in Eq.(2.12) and (2.14) to update uexp and vexp, respectively. The iteration is repeated until the pressure increments converge to a sufficiently small value. After convergence we have the new time velocities and pressures. 2.5 Tracking the Interface This section describes the simple subgrid counting procedure presented by Mortensen, Trapp, and Welch [1], The method uses the void 12 fraction to indicate the existence of interface cells. The cell is pure liquid when the void fraction is 1 and pure gas when the void fraction is 0. Two- phases cells exist when the void fraction is between 0 and 1. The void fraction, a, is defined in term of the densities of the two-phases as a P~ P* Pi ~ Pg (2.17) where p is the cell density and pg and pi are the density of the gas and the liquid phases, respectively. Assuming the void fraction is continuous, the gradient of the void fraction is a vector perpendicular to the surfaces of constant void fraction pointing in the direction of maximum void fraction increase. The normal vector at interface in any two-phase cell pointing into the liquid region can be calculated by using an averaged second order gradient of void fraction a* = I- lV A a + i A AJKAxJ^ V2AAx>i(j + |iÂ¥^ av = -3(C. ^ A + ,i The vector is normalized as follows (2.18) (2.19) 13 n. X (2.20a) n. X (2.20b) Since the unit normal vector at the interface in a cell and the volume place the interface at the correct location within the cell. These geometric calculations are done using a simple subgrid counting scheme based on a subgrid mesh in the cell. At the interface cell, the subgrid is constructed. If the dimension on each side of the cell is Ac, the dimension of the subgrid cells is Axsub=Ax-H7sub, where nsiib is the number of subgrid cells in each direction. As shown in Fig. 2.2 the possible location of the interface is along the direction of the normal vector and can be found by comparing the true void fraction with the void fraction calculated from this counting. The true void fraction calculation is discussed in section 2.6. of the cell occupied by gas are known, geometric calculations must be done to 14 Figure 2.2 The distance along the unit normal vector measured from the cell center is defined as S. Dividing the number of subgrid cells on the liquid side of the plane by the total number of subgrid cells in each cell is approximately the void fraction as a function of the interface location S. S can be solved for consistent with the known void fraction by mapping S into a function of cell void fraction and applying the Secant iteration method with this function. As shown in Fig. 2.2 the vector f is a vector from the intersection of the normal vector and the interface to any subgrid cell. The liquid cell is counted if the inner product f h is positive and the cell is counted as vapor cell if the inner product is negative. If a subgrid cell center is within the distance AxSUb+2 of the interface, the unit counter for this subgrid cell is weighted 15 linearly between zero and one. It is zero if the subgrid cell center is exactly a distance Axmb^2 from the interface on the vapor side. It is one if the subgrid cell center is exactly a distance 4x^+2 from the interface on the liquid side. 2.6 Density Calculation The new time density of each cell is calculated by the mass flux across the cell consistent with the new time velocity. The mass flux is determined by calculating the volume of fluid that is fluxed out of the cell during each time step. Shown in Fig. 2.3 is an example mass flux calculation in the x- direction. The positive x-direction velocity will flux the mass on the right face of the cell within the distance uAt of the right face of the cell. The mass flux is calculated as the density in this cell multiplied by the flux volume. In a two- phase cell the mass flux is calculated using the same method but the subgrid counting procedure is employed to find the volume of gas and the volume of liquid. The gas mass flux is then the gas density multiplied by the gas volume and the liquid mass flux is the liquid density multiplied by the liquid volume. 16 Mass in liquid and vapor region fluxed to the right cell by velocity u Figure 2.3 2.7 Viscous Stress Calculation The Navier-Stokes equations with constant viscosity and density are applicable to each single phase. Since these variables are both discontinuous across the interface, we may expect either excessive numerical diffusion or problems with oscillations around the interface if no special treatment is used at the interface. As described in an earlier section, the density at the interface is calculated by the subgrid counting method. The momentum equations, Eq.(2.2) and Eq.(2.3), can be expressed as 17 x-direction xy da du du 1 (do dx + u + v = - + dt dx dy p \ dx dy ) + ~ Fsx + &x (2-21), y-direction dv dv dv 1 + u + v = dt dx dy p dr. yx V dx + (2.22). We are considering isotropic, incompressible, Newtonian fluids, hence the stresses are du ox = -P + 2p _ dv a = -P + 2p dy * xy ~ Tyx ~ M ' dv_ du_ \dx dy> (2.23a), (2.23b) In constructing finite difference approximations for the viscous stress terms in Eq.(2.21) and (2.22), we must take into account the spatial distribution of viscosity. Consider the x-direction momentum control volume shown in Fig.2.4, 18 Hi-ij+i Hij+i Vj j Ui+lj+1 l ^ a Otij+1/2 Jij -1/2 j ;J a V:; 1 +1/2J +lj Ui-ij- IJ 1 Ctij-1/2 Uij-1 Ui+Ij- Void and viscosity index Figure 2.4 the shear stress, txy, acts on the top and bottom faces, therefore we need the viscosity at the top and bottom faces. Similarly the viscous stress, Txx, acts on the left and right faces hence we need the viscosity at those faces. The finite difference approximations for the viscous terms in the momentum equations are 19 VISCXX VISCXY (2.24), VISCYY VISCYX f f \ ( Y\ \ l Ayi+1 J Pi.l l AVj JJ (2.25). Where the indices on the viscosity indicate its special location, as shown in Fig. 2.4. The viscosity at the cell center is calculated by exact void weighting, the viscosity at the right and left junctions in the x-direction momentum cells is calculated by exact void weighting. At the top and bottom faces of the x- direction momentum cells, the viscosity is calculated by average void weighting. A similar procedure is used for the y-direction momentum cells. This may stated mathematically as, 20 (2.26), where htj = + (l ~ aitj)fi9 *j+i, Mi. 1 i+i.:K (2.27), 1 ai,j+l)Mg (2.28), i)/' + - Km),)"* (2.29), >)/J + i1 - (2.30), i).M + i1 - K=.hLK (2.31), i)f+ t1 - (2.32), [a i i) = fa. a + a. i) \ 1+i'J+\)y 2 \ 1,3+* (2.33a), (2.33b), 21 (2.33c), . + a 2 (2.33d). 2.8 Surface Tension and Curvature Flows involving surface tension are the main focus of this work. The including this phenomena into the governing equations. 2.8.1 Surface Phenomena Following [5], we seek to define the surface tension in Laplaces formula, which gives the surface pressure, as a cause of the pressure difference of the jump across the interface when the interface surface is curved. Consider a segment of an interface with area dA cL^dl,, displaced in the normal direction by an infinitesimal virtual displacement, 5Â£. dJ\ and d_L, are the lengths on each side of the element. following section will describe briefly surface phenomena and our method of 22 Interface displaced through distance dwith pressure pi and P2 in fluid 1 and 2 respectively. Figure 2.5 The virtual work done in displacing the area element shown is J (- Pi + P2)5CdA (2.34), where pi and p2 are the pressure of fluid on each side of the interface. If we define 5W is the total virtual work done in displacing the surface we must add the virtual work done in changing the area of the surface. This work is proportional to the change of the area of the surface, 8A. Thus the total virtual work is 5w -J p2)gÂ£dA + odA (2.35), 23 where ct is the fluid surface tension coefficient. If the system is in an equilibrium condition, the virtual work, 8W, is zero. If Ri and R2 are the principal radii of the curvature of the surface. According to the infinitesimal size of the normal virtual displacement, the corresponding virtual increments of length d\ and dL, are (sQR^dl^ and {hC,/R^)dl2, respectively. Then the virtual surface area increment can be expressed as d\ ( 5Cl r 1 + dJL l + Â£ . V dlydlz dl.dl, 1 + 5C + K. R2 5C,dA dlydJ^ (2.36). The total virtual change in area of the interface is 5A = | sc| JL_ Rl (2.37) Substituting Eq.(2.37) into Eq.(2.35) and setting the total virtual work to zero we obtain the equilibrium condition J Sqta p2) a V-Rl >dA = 0 (2.38) 24 Since SÂ£ is arbitrary, the expression in the bracket must be zero: (Pi Pi) CT ^ + (2.39). Eq.(2.39) is Laplaces formula, which relates the surface tension to the fluid pressures and the surface curvature. For the interface between two viscous fluids in motion, Laplaces formula may be generalized to r (Pi ~ P2h = (hj* ~ T2i*K + CT ' 1 1^ + n, (2.40), Since the surface tension coefficient is not always constant. There is another force tangential to the surface, grada. The momentum balance at the interface may be shown to be f Pi \ (2.41). Note that above equation is the projection of the normal and tangent coordinate into the Cartesian coordinate. The last term on both the left hand side and the right hand side are the surface tension terms in the normal and tangential directions of the interface, respectively. This work considers isothermal problems, hence the surface tension coefficient is constant. The 25 last term on the right hand side is zero. Hence, the effect of the surface tension on the momentum at the interface is only the normal direction. 2.8.2 The Approximation of Surface Tension Effects This section describes how the surface tension forces is modeled and applied to the momentum equation. The effect of surface tension is a surface force that acts on the interface, as shown in Fig. 2.6. The control volume form of the momentum equation can be written as: The surface tension force at the interface on each side of the control volume. Figure 2.6 where T is the stress tensor and h is the outward directed unit normal vector to the surface of the control volume. The body force is represented by g and Fsv is the surface tension force. Considering the interface surface as infinitesimally thin, the surface tension force is added into the momentum equation as a line integral. By employing transport theorem [6] f 'FdV = f cJt v O'? Dt + 'Fdivv dV (2.43), where Y can be any scalar, vector, or tensor, and the divergence theorem: j^divvdV = ndS (2.44), the momentum equation, Eq.(2.42), can be rewritten as Dt = jvdiv^dv + \vP$dV + jc^dL (245a)> or l ^ - divT pg dV = Â£FsvdD (2.45b). Consider the right hand side term, which is the surface tension force term. Viewing Fig 2.6, we calculate 27 (2.46), Â£rsvdi = Â£otdi = J crtjdz + J atzdz + J atfdx + J atadx where t is a unit vector in the tangent plane of the surface directed outward on the boundary. Since we consider two dimensional problems, the A components of above equation in the third dimension vanish because tÂ£ and ta are equal in magnitude but opposite in direction. We have tz -Â£j + At, Eq.(2.46)becomes i^L = J afc + tr)dz = j* a(At)dz (2.47). Substituting this into Eq.(2.42) we have Â£ (pv) divf pg dV = J cr(At^jdz (2.48) Our discrete approximation to the control volume equation is (pv) divf pg AxAyAz = cf(aÂ£)Az (2-49). Dividing by AxAyAz we get 28 (2.50). D x (yov) = divT + pg + DtK ' AxAy From the definition of curvature, ic = , where t is the tangent vector dS of the curve, and dS is an increment along the length of the curve. Substituting a discrete approximation to the curvature into Eq.(2.50) we have D , \ (pv) = divT + pg + Dt a(xAs)n AxAy (2.51) The last term on the right hand side is the surface tension term which is added to the momentum equation in this work. Note that we explicitly show the direction of the force is in the direction normal to the interface. cr(xAs)n Surface Tension Force = -L- (2.52) AxAy The x and y component of the surface tension force was included in Eq.(2.2) and Eq. (2.3), where Fsx and FSy are the components of the surface tension force in x and y directions, respectively. We consider flows with constant surface tension. The arc length, AS, and the curvature, k, in Eq.(2.52) are found by a simple geometric method. 29 These parameters influence the accuracy of the surface tension force. Of the two parameters, the curvature is the most difficult to calculate. The next section will describe the curvature calculation. 2.8.3 The Calculation of Curvature This section discusses the calculation of the curvature that is employed in this work. Since by definition, the curvature vector is the change in the unit tangent vector along a curve with respect to changes in the arc length S, along the curve, K dt dS (2.53), and k = \ (2.54), where R is the radius of curvature. The concept our method uses is to solve for the radius of curvature and then the reciprocal of this radius is the curvature, as expressed in Eq. (2.54). Consider Fig 2.7, (x0, yo), (x/, yj), (x2, y2), are points on the interface which define a curve that we want to find the curvature of. 30 Geometry of curvature calculation. Figure 2.7 Let Si be a straight line connecting points (x0, yo) and (xj, yf), and S2 be a straight line connecting points (xi, yj) and (x2, y2). If Li and L2 are the perpendicular bisectors to the line Sj and S2 respectively, the point of intersection of Lj and L2, (xc, yc), is the center of the curvature of the circle that passes through these three points. In effect, we approximate small sections of the interface as circular arcs. These circular arcs are known in the literature as osculating circles [7], 31 Each cell that contains an interface may have two or more neighboring cells that contain interfaces. We have found that selecting neighboring cells which have void fractions closest to 0.5 give us a more accurate curvature calculation. These cells have a longer length of interface. Generally, the longer the interface length, the more accurate the unit normal vector is. The unit normal vector is important in applying the above geometrical method for calculate the curvature. The curvature vector can be calculated by substituting this curvature into Eq. (2.54). 2.8.4 The Numerical Method As described in section 2.3, a staggered grid is used in our finite difference approximation. The surface tension calculations discussed previously provide a cell centered quantity. We interpolate this cell centered quantity to momentum cells by simple averaging. Hence the surface tension terms in Eq. (2.2) and Eq. (2.3) are expressed as = \ ((p4 + (Oj (2.55a), for the x-direction momentum, and ^ = | (W.+ teA) 32 for the y-direction momentum. Here (f^x) and are the components of the surface tension forces in the x direction of the left hand side and the right hand side of the x-momentum control volume, respectively. Similarly, (fsy) , and (^,y j are the components of the surface tension force in the y direction of the top side and the bottom side of the y-momentum control volume, respectively. According to Eq. (2.5) and Eq. (2.6), the numerical expression of Eq. (2.55) are Xf-Lu, + ( normXj'j (2.56a), and (2.56b), where (normx) and (normy) are unit normal vector components of the interface in x and y directions of two-phase cells. 3. Numerical results This chapter presents some background theory and computational results from the three simulations attempted. These simulations were chosen to test key features of the method. The simulation of the buoyancy driven bubbles is a case where the fluids can have arbitrary density difference, and surface tension and viscosity are important parameters that effect the characteristic of the bubbles. A number of bubbles with various Eo (Eotvos number) and M (Morton number) were simulated and the results compare favorably with results in the literature. The second simulation is an oscillating bubble. This inviscid bubble oscillation is driven by surface tension without gravity. The last simulations are of Rayleigh-Taylor instabilities of two inviscid fluids of different density. Details of all simulations will be discussed in the following sections. 3.1 Buoyancy Driven Bubbles 3.1.1 Theoretical Background Bubbles can be considered as fluid particles which occur in many natural phenomena ( e.g. boiling) and industrial processes (e.g. fermentation). The behavior of bubbles differ from that of solid particles in two main aspects: 34 the bubble is deformed by hydrodynamics forces, and an internal flow will result from momentum transfer across the gas-liquid interface. Due to buoyancy forces, unconstrained gas bubbles rise in liquid. The motion of the liquid induces a change in pressure. The inertia forces, pressure, and surface tension lead to deformation of the bubbles. In spite of a number of previous works, there is no satisfactory theory except for the cases of very small deformation at either high or low Reynolds number [10]. However, for engineering purposes, a number of experimental works have been published in terms of dimensionless groups and empirical correlations. In this case of bubbles rising freely in an infinite media, as shown in Fig. 3.1, Grace[8] has summarized available experimental data and presented it as a graphical correlation in term of three dimensionless parameters; Eotvos number Eo gbpdl a Morton number M _ gpAp _2 _3 Pl<* (31), (3.2), and Reynolds number 35 Re Pde U Pi (3.3). Where Ap is the difference between the density of surrounding fluid and bubble fluid. The dimension of the bubble is given as an equivalent diameter, de and the terminal velocity of the bubble is U. The subscript T means the quantity is the property of the surrounding fluid. As expressed in Equation (3.1)-(3.3), the Eotvos number represent the ratio of the buoyancy force to the surface tension. As far as we know, there is no simple interpretation for the Morton number. The Morton number does include surface tension, inertia, and viscosity. The Reynolds number is the bubble Reynolds number. The effects of these dimensionless parameters on the deformation of the bubble is represented by Graces graphic correlation At a constant Morton number, increasing bubble diameter or decreasing surface increased the Eotvos number. Therefore, the shape of bubble changed from spherical to ellipsoidal or wobbling and finally to spherical cap or dimpled. At a constant Eotvos number, increasing the Morton number causes the bubble shape to change from wobbling or spherical cap to ellipsoidal or dimpled and finally to spherical. The Morton number can be increased by increasing surrounding fluid viscosity, by decreasing surrounding fluid density, or by decreasing surface tension. Note that the Reynolds number increases with decreasing Morton number at constant Eotvos number. The Reynolds number increases with an increase of Eotvos number at constant Morton number. The effect of the Reynolds number on the bubble motion will be discussed later. 36 Figure 3.1 Graces graphical correlation 37 According to the previous studies about the effect of Eotvos number and Morton number on the bubble motion, Tomiyama, Zun, Sou, and Sakaguchi [9], they have been presented that the fluctuation motion of the bubble can be terminated by increasing the Morton number or by decreasing the Eotvos number. As mention above, the Reynolds number decrease when the Morton number increase. 3.1.2 Numerical Results A number of bubbles for various Eo and M are simulated to compare to results in the literature. The approximate steady state shape of the bubbles with velocity fields and pressure profiles are shown in Fig. 3.2-3.6. All simulation are done at 65 by 33 grid resolution except Fig. 3.2 and 3.6 which are done at a finer grid resolution, 130 by 66. All calculations are done with the density ratio pj pb =40 and free slip boundaries. The Eotvos number, Morton number, and the viscosity ratio are noted in the caption of each figure. Considering Fig.3.2-3.5, the bubbles in steady state at small Eo ( small bubble or large surface tension) are close to circular shape. This is because at small Eo, the surface tension and viscous forces are much more important than the inertia force. Hence the deformation of the bubble is very small. At higher Eo number the bubbles obtain a steady state at elliptic shape or nearly flat bottom and semicircular top. When Eo number becomes large ( large bubble or small surface tension), the flow separation causes the wake at the rear of the bubbles to expand. As a result, the bubbles become dimple shapes or ellipsoidal/spherical cap shapes, as shown in Fig. 3.6. 38 Considering small Eo and small M bubble (Fig. 3.2), the bubble is in the ellipsoidal region of Graces chart. Note that the bubble has the top part slightly flatter than the bottom part which is because the front stagnation pressure becomes dominant as described in Ryskin and Leals work [10]. Comparing Fig. 3.2 with Fig. 3.3, the pressure change in the liquid around the bubble is higher in the lower Morton number bubble. The bubble flow in the less viscous fluid causes a larger change in the pressure field of the outer fluid. We compare our simulations with Graces graphical correlation by ignoring that Graces chart is based on three-dimensional and infinite boundary flows. A more precise comparison is made with the two-dimension results of Tryggvason and Unverdi [11]. For the same Eotvos number, Morton number, viscosity ratio, and density ratio, our bubble shapes are in agreement. 39 Mo = I0~7 Eo=1.0 Bl = 40 ^-=88 P* Pi Figure 3.2 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 1 40 Mo =10~ 4 Eo=1.0 = 40 = 493 PA Pi Figure 3.3 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 2. 41 Pb n* Figure 3.4 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 3. 42 Mo =10'2 Eo = 10.0 40 = 277 Pi P* Figure 3.5 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 4. 43 0_ Mo = 100 Eo= 104.0 P;/Pi=40 Ji,/fi6=479 Figure 3.6 Bubble shape, velocity and pressure field. Buoyancy-driven bubble case 5 44 3.2 Oscillating Bubble 3.2.1 Theoretical Background This is a good test simulation as there exists an analytic solution with which to compare result. We consider a gas bubble surrounded by liquid in an equilibrium configuration. The flow is invisid and we neglect gravity. According to these conditions, the surface tension force becomes dominant in stress balancing. When the bubble is perturbed from its equilibrium shape the surface tension will tend to restore the bubble to its original shape. An analysis of the vibration of a circular column of gas in surrounding liquid by Sir Horace Lamb [12] applies to this flow. The analytic procedure is as follows. Consider the potential equation Â£l + L^L.JLÂ£l = o dr2 r dr r2 dd2 a solution for the above equation is
(3-4),
k |