- Email: [email protected]
Three-dimensional oscillations of rising bubbles
Engineering Analysis with Boundary Elements 28 (2004) 315–323 www.elsevier.com/locate/enganabound
Three-dimensional oscillations of rising bubbles C. Pozrikidis*,1 Department of Mechanical and Aerospace Engineering; University of California, San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0411, USA Received 8 November 2002; revised 16 February 2003; accepted 30 April 2003
Abstract A boundary-element method is implemented for simulating the three-dimensional shape and volume oscillations of a stationary or moving bubble occupied by a compressible gas and suspended in an effectively inviscid ambient liquid. The disturbance potential due to the bubble is expressed in terms of an interfacial distribution of potential dipoles, amounting to a double-layer potential, and a point source situated inside the bubble, included to account for changes in the bubble volume. The interface is described by a grid of quadratic boundary elements whose nodes are convected either with the interfacial liquid velocity or with the velocity component normal to the interface. The strength density of the double-layer potential is interpolated isoparametrically and then calculated by solving an integral equation of the second kind using an iterative method. Once the solution has been found, the tangential and normal components of the interfacial velocity are computed, respectively, in terms of tangential derivatives of the harmonic and vector potential; the latter is evaluated as a single-layer potential in terms of the strength of the interfacial vortex sheet underlying the double-layer potential. Results of dynamical simulations reveal that rising bubbles oscillate at higher frequencies than stationary bubbles. q 2003 Elsevier Ltd. All rights reserved. Keywords: Potential flow; Bubble oscillation; Vortex methods; Indirect boundary element methods
1. Introduction When a gas or vapour bubble rises due to gravity or translates relative to the surrounding liquid with velocity U at high Reynolds numbers, the free surface deforms in response to the dynamic pressure due to the relative motion; the Reynolds number is defined as Re ; rUa=m; where r is the liquid density, m is the liquid viscosity, and a is the equivalent bubble radius defined by the equation VB ¼ 4pa3 =3; where VB is the bubble volume. The deformed bubble shape is determined by the Weber number, We ; raU 2 =g; expressing the magnitude of the dynamic pressure relative to capillary pressure due to the surface tension g: When We is sufficiently small and gravitational effects are negligible across the bubble diameter, the bubble surface obtains a nearly spheroidal shape of revolution with respect to the direction of the motion. As the Weber number is raised, the bubble obtains a spherical-cap and possibly skirted shape [1 – 3]. Such shapes have been computed by several authors * Tel.: þ1-858-534-6530; fax: þ 1-858-534-7078. E-mail address: [email protected] (C. Pozrikidis). 1 Internet: http://stokes.ucsd.edu/c_pozrikidis 0955-7997/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00080-8
using asymptotic and numerical methods based on boundary-integral formulations [4 – 8]. Because the velocity of a rising bubble is a function of the a priori unknown bubble shape that is determined by the Weber number, it must be either assumed or computed as part of the solution. Laboratory and theoretical studies have shown that when Re and We fall in a certain regime in parametric space, a path instability occurs in which nearly spherical and ellipsoidal bubbles—but not spherical-cap bubbles— stray away from the rectilinear path and exhibit a zigzag or spiraling helical motion [9,10]. Depending on its volume, a bubble may initially rise in the zigzag path and then follow a helical motion, but not vice-versa. Recent experiments have shown that spherical bubbles whose volume falls in a certain range rising in clean water follow a zigzag path, whereas ellipsoidal bubbles follow a spiral path [11]. Although a physical reason for the path instability has not been established with absolute certainty, evidence suggests that shedding of rotational fluid from an attached or detached wake developing behind the bubble is the most probable cause, as first proposed by Saffman [12]. Corroborative evidence is provided by the observation that
316
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
the purity of the liquid determining the mobility of the interface plays an important role in defining the critical Reynolds number for unstable motion. Flow visualization studies by Lunde and Perkins [13] have suggested that the helical path is associated with the instability of a steady wake, whereas the zigzag path is associated with the ejection of hairpin vortices. In an effort to probe the physical origin of the path instability, Hartunian and Sears [14] performed an experimental study and carried out a simplified stability analysis of the motion of spherical bubbles about the steady state. More recently, Meiron [15] and Feng [16] addressed the stability of non-spherical shape and demonstrated that, in contrast to the predictions of the earlier work, a path instability associated with shape oscillations does not arise. Moreover, their analysis revealed that the degeneracy of the frequency of shape oscillations with respect to the azimuthal modes is removed when the unperturbed bubble shape is non-spherical. Deviation from sphericity at the steady state was found to reduce the frequency of oscillation. Lunde and Perkins [13] observed shape oscillations of rising bubbles in the axisymmetric 2-0 mode, where the interface oscillates between oblate and prolate spheroidal shapes, and in the three-dimensional 2-2 mode, where the interface oscillates between two triaxial ellipsoidal shapes whose major axes lie in equatorial planes, and attributed the oscillations to mode excitation due to vortex shedding. The frequency of the oscillations was found to be in good agreement with theoretical analysis and heuristic predictions. These results seem to suggest that vortex shedding is significant so far as to excite normal modes that subsequently evolve as though the ambient liquid were inviscid. Computational studies of finite amplitude axisymmetric or three-dimensional oscillations translating bubbles coupling volume and surface modes have not been carried out. In this paper, we apply an indirect boundary-element developed recently by the author [17] to simulate the three-dimensional oscillations of a bubble translating in an effectively inviscid ambient liquid. In the earlier study, the method was applied to study the oscillations of a stationary bubble suspended in a quiescent fluid. The algorithm is based on the ingenious premise of generalized vortex methods developed by Baker et al. [18 – 20]. In this formulation, the potential flow in the ambient liquid is represented by an interfacial distribution of potential dipoles amounting to a harmonic double layer, and a point source situated inside the bubble included to account for changes in the bubble volume. From a mathematical viewpoint, the presence of the point source is necessary for completing the deficient range of the double-layer potential in the case of exterior flow [21], p. 172. The double-layer representation is completed in a way that preserves linearity by setting the strength of the point source proportional to the surface integral of the density of the double-layer potential. Once the strength
density of the double-layer potential is available, the tangential and normal components of the interfacial velocity are computed, respectively, in terms of tangential derivatives of the harmonic and vector potential. The latter is evaluated as a single-layer potential in terms of the strength of the interfacial vortex sheet underlying the double-layer potential. One important advantage of the generalized vortex method compared to methods based on Green’s third identity is that the density of the doublelayer potential is found by solving an integral equation of the second kind instead of an integral equation of the first kind, and this can be done accurately and economically using iterative methods based on successive substitutions, as will be discussed in Section 2. The numerical algorithm for simulating bubble oscillations in a stationary fluid was discussed by Pozrikidis [17]. For completeness, and to explicitly demonstrate the necessary modifications in the case of a moving bubble or a stationary bubble subjected to an ambient flow, a summary of the algorithm is given in Section 2. Theoretical results on the deformation of rising bubbles are reviewed in Section 3, and results of numerical simulations are presented in Section 4.
2. Numerical method Consider potential flow past or due to the expansion, shrinkage, or oscillations of a gas or vapour bubble occupied by a fluid with negligible density and viscosity, suspended in an effectively inviscid and incompressible ambient fluid. The irrotational velocity field around the bubble, u, can be represented by the gradient of the single-valued potential f; u ¼ 7f: The continuity equation requires that the potential satisfy Laplace’s equation 72 f ¼ 0:
ð2:1Þ
In the mathematical formulation, the harmonic potential f is expressed in terms of (a) the incident field denoted by f1 ; (b) a double-layer integral defined over the bubble surface denoted by D; and (c) a point source situated at a point xs that is located inside the bubble, ð fðx0 Þ ¼ f1 ðx0 Þ þ qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ D
2 sðqÞGðx0 ; xs Þ;
ð2:2Þ
where q is the a priori unknown density of the doublelayer potential, n is the unit vector normal to the bubble surface pointing toward the exterior, sðqÞ is the strength of the point source, and Gðx; x0 Þ is the Green’s or Neumann function of Laplace’s equation in three dimensions, [22], p. 516. If the bubble is suspended in a virtually infinite liquid, Gðx; x0 Þ is the free-space Green’s function; whereas if the bubble is suspended in a semi-infinite fluid bounded by an impenetrable plane
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
wall, Gðx; x0 Þ is the Neumann Green’s function whose gradient vanishes over the wall. To remove the ambiguity in the definition of the strength of the point source while preserving linearity, we set 4p 1=2 ð sðqÞ ; 2bðtÞð Þ q dS; ð2:3Þ SD D where SD is the surface area of D; and b is an arbitrary dimensionless coefficient. Taking the limit of Eq. (2.2) as the evaluation point x0 approaches the free surface D from the exterior, and expressing the limit of the double-layer integral in terms of its principal value denoted by PV; we derive an integral equation of the second kind for q; ðPV 1 qðx0 Þ þ qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ 2 D ð2:4Þ 2 sðqÞGðx0 ; xs Þ;
fðx0 Þ ¼f1 ðx0 Þ þ
where the surface potential on the left-hand side is considered to be known. Rearranging, we derive the standard form ðPV qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ þ sðqÞGðx0 ; xs Þ qðx0 Þ ¼2 2 D þ fðx0 Þ 2 f1 ðx0 Þ : ð2:5Þ The homogeneous integral equation corresponding to Eq. (2.5) is
cðx0 Þ ¼ 22
ðPV D
cðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ þ 2sðcÞGðx0 ; xs Þ; ð2:6Þ
where c is an eigenfunction. The spectrum of eigenvalues of this integral equation is discussed in Ref. [17]. Alternatively, the velocity field associated with the double-layer potential may be expressed as the curl of a vector potential, A, and the velocity may be expressed in the form uðx0 Þ ¼ 70 f1 ðx0 Þ þ 70 £ Aðx0 Þ 2 sðqÞ70 Gðx0 ; xs Þ;
ð2:7Þ
where 70 denotes the gradient with respect to the evaluation point x0 : It can be shown that the vector potential is given by the integral representation ð zðxÞGðx; x0 ÞdSðxÞ; ð2:8Þ Aðx0 Þ ¼ D
where z is the strength of the vortex sheet corresponding to the double-layer given by
z ¼ n £ 7q;
ð2:9Þ
[22], p. 503. The right-hand side of Eq. (2.9) involves only tangential derivatives of the density distribution q alone. According to the general principles of the generalized vortex method, the evolution of the bubble interface is
317
simulated using a boundary-element method according to the following steps: Step 1. Specify the initial shape of the interface, interfacial distribution of the harmonic potential f; bubble pressure, and pressure at infinity. Step 2. Solve the integral Eq. (2.5) for the double-layer density q; as will be discussed later in this section. Step 3. Evaluate the harmonic potential f over the interface using Eqs. (2.2) and (2.3). Step 4. Compute the tangential component of the velocity over the interface, n £ ðu £ nÞ ¼ ðI 2 nnÞ·u; in terms of the surface gradient 7s f; where I is the unit matrix and dij 2 ni nj is the tangential projection matrix. In practice, the surface gradient is computed by solving a linear system of three equations, which arises by requiring that the projection of 7s f in two tangential directions is equal to the corresponding directional derivatives, whereas the projection in the normal direction is zero. Step 5. Compute the strength of the vortex sheet z from Eq. (2.9). Step 6. Compute the interfacial distribution of the vector potential A from Eq. (2.8). Step 7. Compute the normal component of the interfacial velocity nðn·uÞ in terms of (a) the normal component of the incident flow, (b) tangential derivatives of the vector potential A, and (c) the normal component of the point source corresponding to the third term on the right-hand side of Eq. (2.2). Step 8. Advance the position of marker points distributed over the interface either with the liquid velocity or with the normal component. In the first case, the marker points are Lagrangian point particles of the exterior liquid. Step 9. Update the harmonic potential at the interfacial marker points using Bernoulli’s equation, taking into consideration that the pressure of the liquid at the bubble surface is equal to the bubble pressure, pB ; reduced by the capillary pressure, 2gkm ; where g is the surface tension and km is the interface mean curvature [22], p. 489. For a spherical bubble of radius a; km ¼ 1=a: Specifically, the rate of change of the potential following Lagrangian point particles moving with the liquid velocity is given by Df 1 2g p 2c km 2 B ¼ lul2 þ þ g·x; 2 r Dt r
ð2:10Þ
where D=Dt is the material derivative, r is the liquid density, c is Bernoulli’s constant, and g is the acceleration of gravity. On the other hand, if the marker points move with the normal component of the liquid velocity, the rate of change of the potential is given by df 1 2g p 2c km 2 B ¼ lu·nl2 2 lul2 þ þ g·x: 2 r dt r
ð2:11Þ
When gravitational effects are insignificant across the bubble diameter, the last term on the right-hand sides of Eqs. (2.10) and (2.11) is absent. If, in addition, the flow
318
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
vanishes at infinity, Bernoulli’s constant c is equal to the liquid pressure far from the bubble, c ¼ p1 : On the other hand, if the flow far from the bubble executes steady streaming motion with velocity U; c ¼ p1 þ r 12 lUl2 : Step 10. Update the bubble pressure according to the polytropic equation of state pB 2 pV ¼ pB ðt ¼ 0Þ 2 pV
vB ðt ¼ 0Þ vB
l
;
ð2:12Þ
where vB is the bubble volume, pV is the partial vapor pressure of the ambient liquid inside the bubble, assumed to remain constant during the motion, and the exponent l is determined by the prevailing thermodynamic conditions. Under adiabatic conditions, l is equal to the ratio of the specific heats; for an ideal gas, l ¼ 1:4: In the case of a pure gas bubble considered in the present simulations, pV ¼ 0: Step 11. Update the pressure at infinity according to a specified protocol. The solution of the integral Eq. (2.5) required in Step 2 is found by the method of successive substitutions. This involves guessing the distribution of q; computing the righthand side of Eq. (2.5), and then replacing the guessed distribution with the left-hand side of Eq. (2.5). To expedite the computation, except at the initial instant, the guessed distribution is identified with the converged solution at the previous time step. When this is done, only 3– 10 iterations are necessary to obtain the solution accurate to the eighth decimal place. An alternative would be to explicitly construct the influence matrix corresponding to the right-hand side of Eq. (2.5), and then solve the linear system using a direct or iterative method. We found, however, that the above described method of successive substitutions is both more efficient and easier to implement. In the numerical implementation, the bubble surface is discretized into an unstructured grid of six-node curved triangular elements, and all geometrical variables and surface functions are approximated with quadratic functions over each element with respect to local triangle coordinates. The local interpolation functions are given in Ref. [23]. The interfacial grid is generated by successively subdividing the faces of a regular octahedron into descendant elements, and finally projecting the nodes onto the sphere. To improve the accuracy, the normal vector and mean curvature are averaged at each node over the elements sharing the node, and the surface distributions of the normal vector and mean curvature over each element are subsequently described by quadratic interpolation. A numerical instability of unknown origin undermines the longevity of the simulations, as discussed in Ref. [17]. To allow the computations to proceed for an extended period of time, the surface potential is smoothed by expanding it in a Fourier – Legendre series, and then truncating the series to a certain level, as discussed in Ref. [17]. Numerical experimentation has shown that smoothing by spectrum truncation is an effective and
benign method of regularization, as will be discussed in Section 4.
3. Small-deformation theory Consider a gas bubble rising in a quiescent ambient liquid or translating relative to a stationary or moving ambient liquid at high Reynolds numbers, and assume that the bubble is and remains spherical at all times. In a frame of reference moving with the bubble, the fluid far from the bubble executes streaming motion with velocity U. Elementary fluid mechanics shows that, if viscous effects are negligible and the flow is assumed to be irrotational, the velocity field is given by the gradient of the harmonic potential ! 1 a3 ; ð3:1Þ fðxÞ ¼ U·^x 1 þ 2 r3 where a is the bubble radius, x^ ¼ x 2 xs is the distance of the field point x from the bubble center xs ; and r ¼ l^xl [22]. When gravitational effects are insignificant, the pressure far from the bubble is given by p1 ¼ c 2
1 rlUl2 ; 2
ð3:2Þ
where c is Bernoulli’s constant. However, because the pressure over the exterior surface of the bubble varies according to Bernouli’s equation, the interface cannot remain spherical at all times and will deform by an amount that is determined by the Weber number, We ; ralUl2 =g; where g is the surface tension. Deformation may cause the bubble to translate as a whole along the x axis with a velocity that depends on the Weber number. The dimensionless Bernoulli constant d ¼ ac=g plays an important role in determining changes in the bubble volume and pressure due to the motion. Assume, in particular, that the bubble translates along the x axis, and describe the motion in a frame of reference where the spherical bubble is stationary and the incident flow is given by U ¼ 2Ux ex : If the initial bubble shape is spherical and the initial distribution of the potential is given by Eq. (3.1), the bubble will deform to obtain an oblate spheroidal shape, while translating and exhibiting axisymmetric shape oscillations in the 2-0 mode. In real life, the oscillations will eventually be dampened by viscous dissipation, and the bubble will assume a steady oblate spheroidal shape. Moore [24] showed that, when the deviation from sphericity is small, to first order in We, the minor and major bubble axes are given by ax ¼ að1 2 e Þ and as ¼ að1 þ ð1=2Þe Þ; where
e¼
3 We: 16
ð3:3Þ
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
319
The corresponding bubble axes ratio is
x;
as 9 We þ OðWe2 Þ: .1þ 32 ax
ð3:4Þ
Subsequently, Moore [4] derived the improved nonlinear inverse relation qffiffiffiffiffiffiffiffiffi x3 þ x 2 2 2 We . 2 4=3 2 ½ x arcsec x 2 x2 2 12 ; ð3:5Þ x ðx 2 1Þ3 which was shown by Miksis et al. [6] to provide accurate predictions for Weber numbers on the order of unity.
4. Numerical simulations In the numerical simulations, we consider a bubble with a spherical initial shape of radius a; and specify the initial distribution of the surface potential ! 1 a3 ga ðj 2 1Þðj þ 2Þ 1=2 2 eS fðxÞ ¼U·^x 1 þ 2 r3 jþ1 r Tjlml ðcos uÞcosðmwÞ;
ð4:1Þ
where e S is a dimensionless amplitude and Tjlml are the modified Legendre functions [17]. In all simulations discussed in this section, the bubble surface was discretized into 512 elements defined by 1026 surface nodes, the polytropic exponent l in Eq. (2.12) was set to 1.25, the reduced Bernoulli constant d was set to 100, and time integration was carried out by the second-order RungeKutta method with constant time step Dt0 ¼ 0:02: The dimensionless time indicated by a prime is defined as t0 ¼ tðg=ra3 Þ1=2 : Unless specified otherwise, smoothing was effected by truncating the Fourier –Legendre spectrum to the value of 8. Numerical experimentation has shown that, for the interfacial grid considered, this level of truncation strikes an optimal balance between consistency and regularity of the simulated motion. Fig. 1 illustrates the oscillations of a bubble in the axisymmetric 2-0 mode, corresponding to j ¼ 2 and m ¼ 0; for We ¼ 0:25 and e S ¼ 0:25: The initially upward curves represent the bubble axis along the x axis of symmetry, defined as half the maximum bubble size in this direction, reduced by the initial equivalent bubble radius. The initially downward curves represent the bubble axis normal to the x axis, reduced by the initial equivalent bubble radius. In the simulations represented by the solid and hardly visible dotted line, the marker points move with the velocity normal to the interface. In the case of the dotted line, smoothing is not enabled, whereas in the case of the solid line, smoothing is applied after each time step. In the simulations represented by the dashed line and circular symbols, the marker points move with the liquid velocity. In the case of the circular symbols, smoothing is not enabled, whereas in the case of the dashed line, smoothing is
Fig. 1. Axisymmetric oscillations of a bubble in the 2-0 mode, for We ¼ 0:25 and e S ¼ 0:25; computed with different options. (a) Evolution of the reduced x axis, and (b) deformed grid developing when the nodes move with the fluid velocity.
applied after each time step. While the agreement between all simulations is excellent during the initial stage of the motion, only the simulation represented by the solid line survives for an extended period of time. Note, in particular, that because the interfacial marker points are convected toward the rear of the bubble under the influence of the streaming flow, severe grid distortion occurs when the nodes are identified with material point particles moving with the fluid velocity, and the computation breaks down at an early stage of the simulation. Fig. 1(b) depicts a deformed grid at the time when the simulation fails, corresponding to the dashed line in Fig. 1(a). Moore’s [24] small-deformation theory predicts that, in the absence of oscillations, a bubble translating at We ¼ 0:25 has a steady spheroidal shape with minimum axis ax ¼ 0:9531a; maximum axis as ¼ 1:0234a; and axes ratio as =ax ¼ 1:0234; where a is the equivalent bubble radius. Consider this Moore spheroid, and assume that the distribution of the potential over the interface at the initial instant is described by Eq. (4.1) with e S ¼ 0: The heavy solid lines in Fig. 2(a) and (b) show the evolution of the longitudinal and transverse bubble axes, ax and as ; reduced by the initial bubble equivalent radius a0 ; the heavy solid line in Fig. 2(c) shows the evolution of
320
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
Fig. 2. Axisymmetric oscillations of a bubble in the 2-0 mode. The heavy solid lines correspond to the Moore spheroid for We ¼ 0:25: (a) Bubble axis in the x direction, (b) bubble axis in the transverse direction, and (c) bubble volume reduced by the initial value.
the bubble volume. In order to sustain the regularity of the grid for an extended period of time, the marker points are required to move normal to the interface. The numerical results reveal that the bubble exhibits small-amplitude axisymmetric periodic oscillations about the Moore spheroid. To demonstrate the effect of the Weber number on the finite-amplitude motion, we consider the evolution of an
initially spherical bubble, and set the initial distribution of the potential as shown in Eq. (4.1) with e S ¼ 0:15 for the 2-0 mode. The solid and dashed lines in Fig. 2 show results of simulations, respectively, for We ¼ 0 and 0.25. In fact, the dashed line in Fig. 2(a) extends the solid line displayed earlier in Fig. 1(a) to longer times. Consider first the free oscillations of a stationary bubble immersed in a quiescent fluid, We ¼ 0: According to small deformation theory, the x axis of the bubble, reduced by the mean bubble radius, oscillates around the mean value of unity with dimensionless period T 0 ¼ 1:814 and reduced amplitude 0.0946, corresponding to the horizontal and vertical lines in Fig. 2(a), while the transverse axis oscillates with an 1808 phase difference with half the reduced amplitude 0.0473, corresponding to the horizontal and vertical lines in Fig. 2(b). The bubble volume exhibit simultaneous oscillations with dimensionless period T 0 ¼ 0:322; as shown by the vertical lines in Fig. 2(c). These predictions are borne out by the numerical results shown with the solid lines in Fig. 2. Note that, even though the volume mode is not excited at the initial instant, volume oscillations do occur because of nonlinear interactions. Comparing the solid with the solid lines in Fig. 2(a) and (b), we see that raising the Weber number increases the period, and thus reduces the frequency of the shape oscillations, in agreement with theoretical predictions for small-amplitude motion [15,16]. Fig. 3(a) describes the evolution of the x position of the bubble centroid for We ¼ 0:25; corresponding to the dashed lines depicted in Fig. 2. The results show that a bubble that would be stationary in the spherical state, moves opposite to the direction of the oncoming flow. Conversely, a rising bubble slows down as a result of the deformation. The net motion effectively alters the Weber number of the incident flow by an amount that is a priori unknown and must be found as part of the solution. Fig. 3(b) displays the bubble interface as described by the boundary-element grid at time t0 ¼ 3:95; the incident streaming flow is directed from left to right toward the negative direction of the x axis. Similar results are obtained for the three-dimensional 2-2 mode, corresponding to j ¼ 2 and m ¼ 2; where the bubble surface oscillates between two triaxial ellipsoidal shapes whose major axes lie in the equatorial xy or xz plane. Fig. 4(a) –(c) shows the evolution of the bubble axes, defined as half the maximum bubble size in the corresponding direction reduced by the initial bubble radius, for e S ¼ 0:20; and We ¼ 0 (solid lines) or 0.25 (dashed lines); in the second case, the bubble translates along the z axis. Fig. 5(a) displays oscillations in the bubble volume excited by nonlinear interactions. The numerical results are consistent with the predictions of the linear theory for small-amplitude shape and volume oscillations represented by the vertical and horizontal lines in Fig. 4(b) and (c). Note that in the small-amplitude linearized motion, the x axis of the bubble described in Fig. 4(a) remains unchanged.
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
321
Fig. 3. Axisymmetric oscillations of a translating bubble in the 2-0 mode for We ¼ 0:25; corresponding to the simulation represented by the dashed lines in Fig. 2; (a) reduced x position of the bubble centroid, x0 ¼ x=a; and (b) instantaneous bubble surface at time t0 ¼ 3:95; the incident streaming flow is directed from right to left toward the negative direction of the x axis.
The simulation of the 2-2 mode oscillations reveal that, as in the case of the 2-0 mode, raising the Weber number increases the period, and thus reduces the frequency of the shape oscillations. Furthermore, in the case of a translating bubble, oscillations occur about a mean spheroidal shape with two major axes perpendicular to the direction of the motion. Fig. 5(b) shows the evolution of the z position of the bubble centroid for We ¼ 0:25; corresponding to the dashed lines in Fig. 4. The results show that a rising bubble slows down as a result of the deformation. Finally, Fig. 5(c) displays the bubble interface as described by the boundary-element grid at time t0 ¼ 3:95; the incident streaming flow is directed from top to bottom toward the negative direction of the z axis.
5. Discussion The completed double-layer representation combined with Baker’s generalized vortex method has been shown
Fig. 4. Oscillations of the bubble axis in the 2-2 mode for We ¼ 0 (solid lines) or 0.25 (dashed lines).
to provide us with an efficient algorithm for simulating three-dimensional potential flow due to the motion or past compressible bubbles. Numerical instabilities are a serious impediment to the longevity of the simulations over an extended period of time, but smoothing by
322
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
Fig. 5. Oscillations of a translating bubble in the 2-2 mode, corresponding to the simulation represented by the dashed lines in Fig. 4. (a) Oscillations in the bubble volume, (a) reduced z position of the bubble centroid, z0 ¼ z=a; and (c) instantaneous bubble surface at time t0 ¼ 3:98; the incident streaming flow is directed from top to bottom toward the negative direction of the z axis.
spectrum truncation removes numerical irregularities while maintaining the fundamental features of the nonlinear motion. The superiority of the method compared to direct methods based on Green’s third identity lies in its ability to compute the instantaneous flow by solving an integral equation of the second kind, which can be done accurately and economically by the method of successive substitutions. An important aspect of the boundary-element implementation is the option for the boundary-element nodes to be convected normal to the interface with the liquid velocity and tangentially with an arbitrary component, thereby ensuring the regularity of the interfacial grid for an extended period of time. Though regridding is an alternative [25], the introduction of numerical error due to interpolation in parameter space or directly in three-dimensional physical space is a serious concern. The algorithm has been applied with success in this and in the earlier work [17] to simulate bubble oscillations in infinite flow. Although not tested, the method is likely to perform equally well in the case of flow bounded by an impenetrable boundary such as an infinite plane wall. The only necessary modification is the substitution of the free-space Green’s function with a Neumann function appropriate for the geometry of the flow under
consideration. A library of Green and Neumann functions accompany a recent text [26] and is available from the internet site http://bemlib.ucsd.edu
References [1] Wegener PP, Parlange JY. Spherical-cap bubbles. Annu Rev Fluid Mech 1973;4:79–100. [2] Miksis MJ, Vanden-Broeck JM, Keller JB. Rising bubbles. J Fluid Mech 1982;123:31–41. [3] Batchelor GK. The stability of a large gas bubble rising through liquid. J Fluid Mech 1987;184:399 –422. [4] Moore DW. The velocity of rise of distorted gas bubbles in a liquid of small viscosity. J Fluid Mech 1959;23:749 –66. [5] Harper JF. The motion of bubbles and drops through liquids. Adv Appl Mech 1972;12:59–129. [6] Miksis MJ, Vanden-Broeck JM, Keller JB. Axisymmetric bubble or drop in a uniform flow. J Fluid Mech 1981;108:89–100. [7] El Sawi M. Distorted gas bubbles at large Reynolds number. J Fluid Mech 1974;62:163–83. [8] Benjamin TB. Hamiltonian theory for motion of bubbles in an infinite liquid. J Fluid Mech 1987;2:118– 28. [9] Duineveld PC. Rise velocity and shape of bubbles in pure water. J Fluid Mech 1995;292:325 –32. [10] Magnaudet J, Eames I. The motion of high-Reynolds number bubbles in inhomogeneous flows. Annu Rev Fluid Mech 2000;32: 659 –708.
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323 [11] Wu M, Gharib M. Experimental studies on the shape and path instability of small air bubbles rising in clean water. Phys Fluids 2002; 7:L49–L52. [12] Saffman PG. On the rise of small air bubbles in water. J. Fluid Mech 1956;1:249–75. [13] Lunde K, Perkins R. Shape oscillations of rising bubbles. Appl Sci Res 1998;58:387–408. [14] Hartunian RA, Sears WR. On the instability of small gas bubbles moving uniformly in various liquids. J Fluid Mech 1957;3: 27 –47. [15] Meiron DI. On the stability of gas bubbles rising in an inviscid fluid. J Fluid Mech 1989;198:101–14. [16] Feng JQ. The oscillations of a bubble moving in an inviscid fluid. SIAM J Appl Math 1992;52:1–14. [17] Pozrikidis C. Numerical simulation of three-dimensional bubble oscillations by a generalized vortex method. Theor Comput Fluid Dyn 2002;16:133– 50. [18] Baker GR. Generalized vortex method for free-surface flow problems. In Waves in fluid interfaces MRC, Wisconsin: University of Wisconsin; 1983.
323
[19] Baker GR, Meiron DI, Orszag SA. Generalized vortex methods for free-surface flow problems. J. Fluid Mech 1982;123: 477–501. [20] Baker GR, Meiron DI, Orszag SA. Boundary integral methods for axisymmetric and three-dimensional Rayleigh –Taylor instability problems. Physica 1984;12D:19–31. [21] Mikhlin SG. Integral equations and their applications to certain problems in mechanics, mathematical physics and technology. New York: Pergamon Press; 1957. [22] Pozrikidis C. Introduction to theoretical and computational fluid dynamics. Oxford University Press; 1997. [23] Pozrikidis C. Numerical computation in science and engineering. Oxford University Press; 1998. [24] Moore DW. The rise of a gas bubble in a viscous liquid. J Fluid Mech 1959;6:113 –30. [25] Kwak S, Pozrikidis C. Adaptive triangulation of evolving, closed or open surfaces by the advancing-front method. J Comp Phys 1998;145: 61–88. [26] Pozrikidis C. A practical guide to boundary element methods with the software library BEMLIB. Chapman & Hall/CRC Press; 2002.
Three-dimensional oscillations of rising bubbles C. Pozrikidis*,1 Department of Mechanical and Aerospace Engineering; University of California, San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0411, USA Received 8 November 2002; revised 16 February 2003; accepted 30 April 2003
Abstract A boundary-element method is implemented for simulating the three-dimensional shape and volume oscillations of a stationary or moving bubble occupied by a compressible gas and suspended in an effectively inviscid ambient liquid. The disturbance potential due to the bubble is expressed in terms of an interfacial distribution of potential dipoles, amounting to a double-layer potential, and a point source situated inside the bubble, included to account for changes in the bubble volume. The interface is described by a grid of quadratic boundary elements whose nodes are convected either with the interfacial liquid velocity or with the velocity component normal to the interface. The strength density of the double-layer potential is interpolated isoparametrically and then calculated by solving an integral equation of the second kind using an iterative method. Once the solution has been found, the tangential and normal components of the interfacial velocity are computed, respectively, in terms of tangential derivatives of the harmonic and vector potential; the latter is evaluated as a single-layer potential in terms of the strength of the interfacial vortex sheet underlying the double-layer potential. Results of dynamical simulations reveal that rising bubbles oscillate at higher frequencies than stationary bubbles. q 2003 Elsevier Ltd. All rights reserved. Keywords: Potential flow; Bubble oscillation; Vortex methods; Indirect boundary element methods
1. Introduction When a gas or vapour bubble rises due to gravity or translates relative to the surrounding liquid with velocity U at high Reynolds numbers, the free surface deforms in response to the dynamic pressure due to the relative motion; the Reynolds number is defined as Re ; rUa=m; where r is the liquid density, m is the liquid viscosity, and a is the equivalent bubble radius defined by the equation VB ¼ 4pa3 =3; where VB is the bubble volume. The deformed bubble shape is determined by the Weber number, We ; raU 2 =g; expressing the magnitude of the dynamic pressure relative to capillary pressure due to the surface tension g: When We is sufficiently small and gravitational effects are negligible across the bubble diameter, the bubble surface obtains a nearly spheroidal shape of revolution with respect to the direction of the motion. As the Weber number is raised, the bubble obtains a spherical-cap and possibly skirted shape [1 – 3]. Such shapes have been computed by several authors * Tel.: þ1-858-534-6530; fax: þ 1-858-534-7078. E-mail address: [email protected] (C. Pozrikidis). 1 Internet: http://stokes.ucsd.edu/c_pozrikidis 0955-7997/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00080-8
using asymptotic and numerical methods based on boundary-integral formulations [4 – 8]. Because the velocity of a rising bubble is a function of the a priori unknown bubble shape that is determined by the Weber number, it must be either assumed or computed as part of the solution. Laboratory and theoretical studies have shown that when Re and We fall in a certain regime in parametric space, a path instability occurs in which nearly spherical and ellipsoidal bubbles—but not spherical-cap bubbles— stray away from the rectilinear path and exhibit a zigzag or spiraling helical motion [9,10]. Depending on its volume, a bubble may initially rise in the zigzag path and then follow a helical motion, but not vice-versa. Recent experiments have shown that spherical bubbles whose volume falls in a certain range rising in clean water follow a zigzag path, whereas ellipsoidal bubbles follow a spiral path [11]. Although a physical reason for the path instability has not been established with absolute certainty, evidence suggests that shedding of rotational fluid from an attached or detached wake developing behind the bubble is the most probable cause, as first proposed by Saffman [12]. Corroborative evidence is provided by the observation that
316
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
the purity of the liquid determining the mobility of the interface plays an important role in defining the critical Reynolds number for unstable motion. Flow visualization studies by Lunde and Perkins [13] have suggested that the helical path is associated with the instability of a steady wake, whereas the zigzag path is associated with the ejection of hairpin vortices. In an effort to probe the physical origin of the path instability, Hartunian and Sears [14] performed an experimental study and carried out a simplified stability analysis of the motion of spherical bubbles about the steady state. More recently, Meiron [15] and Feng [16] addressed the stability of non-spherical shape and demonstrated that, in contrast to the predictions of the earlier work, a path instability associated with shape oscillations does not arise. Moreover, their analysis revealed that the degeneracy of the frequency of shape oscillations with respect to the azimuthal modes is removed when the unperturbed bubble shape is non-spherical. Deviation from sphericity at the steady state was found to reduce the frequency of oscillation. Lunde and Perkins [13] observed shape oscillations of rising bubbles in the axisymmetric 2-0 mode, where the interface oscillates between oblate and prolate spheroidal shapes, and in the three-dimensional 2-2 mode, where the interface oscillates between two triaxial ellipsoidal shapes whose major axes lie in equatorial planes, and attributed the oscillations to mode excitation due to vortex shedding. The frequency of the oscillations was found to be in good agreement with theoretical analysis and heuristic predictions. These results seem to suggest that vortex shedding is significant so far as to excite normal modes that subsequently evolve as though the ambient liquid were inviscid. Computational studies of finite amplitude axisymmetric or three-dimensional oscillations translating bubbles coupling volume and surface modes have not been carried out. In this paper, we apply an indirect boundary-element developed recently by the author [17] to simulate the three-dimensional oscillations of a bubble translating in an effectively inviscid ambient liquid. In the earlier study, the method was applied to study the oscillations of a stationary bubble suspended in a quiescent fluid. The algorithm is based on the ingenious premise of generalized vortex methods developed by Baker et al. [18 – 20]. In this formulation, the potential flow in the ambient liquid is represented by an interfacial distribution of potential dipoles amounting to a harmonic double layer, and a point source situated inside the bubble included to account for changes in the bubble volume. From a mathematical viewpoint, the presence of the point source is necessary for completing the deficient range of the double-layer potential in the case of exterior flow [21], p. 172. The double-layer representation is completed in a way that preserves linearity by setting the strength of the point source proportional to the surface integral of the density of the double-layer potential. Once the strength
density of the double-layer potential is available, the tangential and normal components of the interfacial velocity are computed, respectively, in terms of tangential derivatives of the harmonic and vector potential. The latter is evaluated as a single-layer potential in terms of the strength of the interfacial vortex sheet underlying the double-layer potential. One important advantage of the generalized vortex method compared to methods based on Green’s third identity is that the density of the doublelayer potential is found by solving an integral equation of the second kind instead of an integral equation of the first kind, and this can be done accurately and economically using iterative methods based on successive substitutions, as will be discussed in Section 2. The numerical algorithm for simulating bubble oscillations in a stationary fluid was discussed by Pozrikidis [17]. For completeness, and to explicitly demonstrate the necessary modifications in the case of a moving bubble or a stationary bubble subjected to an ambient flow, a summary of the algorithm is given in Section 2. Theoretical results on the deformation of rising bubbles are reviewed in Section 3, and results of numerical simulations are presented in Section 4.
2. Numerical method Consider potential flow past or due to the expansion, shrinkage, or oscillations of a gas or vapour bubble occupied by a fluid with negligible density and viscosity, suspended in an effectively inviscid and incompressible ambient fluid. The irrotational velocity field around the bubble, u, can be represented by the gradient of the single-valued potential f; u ¼ 7f: The continuity equation requires that the potential satisfy Laplace’s equation 72 f ¼ 0:
ð2:1Þ
In the mathematical formulation, the harmonic potential f is expressed in terms of (a) the incident field denoted by f1 ; (b) a double-layer integral defined over the bubble surface denoted by D; and (c) a point source situated at a point xs that is located inside the bubble, ð fðx0 Þ ¼ f1 ðx0 Þ þ qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ D
2 sðqÞGðx0 ; xs Þ;
ð2:2Þ
where q is the a priori unknown density of the doublelayer potential, n is the unit vector normal to the bubble surface pointing toward the exterior, sðqÞ is the strength of the point source, and Gðx; x0 Þ is the Green’s or Neumann function of Laplace’s equation in three dimensions, [22], p. 516. If the bubble is suspended in a virtually infinite liquid, Gðx; x0 Þ is the free-space Green’s function; whereas if the bubble is suspended in a semi-infinite fluid bounded by an impenetrable plane
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
wall, Gðx; x0 Þ is the Neumann Green’s function whose gradient vanishes over the wall. To remove the ambiguity in the definition of the strength of the point source while preserving linearity, we set 4p 1=2 ð sðqÞ ; 2bðtÞð Þ q dS; ð2:3Þ SD D where SD is the surface area of D; and b is an arbitrary dimensionless coefficient. Taking the limit of Eq. (2.2) as the evaluation point x0 approaches the free surface D from the exterior, and expressing the limit of the double-layer integral in terms of its principal value denoted by PV; we derive an integral equation of the second kind for q; ðPV 1 qðx0 Þ þ qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ 2 D ð2:4Þ 2 sðqÞGðx0 ; xs Þ;
fðx0 Þ ¼f1 ðx0 Þ þ
where the surface potential on the left-hand side is considered to be known. Rearranging, we derive the standard form ðPV qðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ þ sðqÞGðx0 ; xs Þ qðx0 Þ ¼2 2 D þ fðx0 Þ 2 f1 ðx0 Þ : ð2:5Þ The homogeneous integral equation corresponding to Eq. (2.5) is
cðx0 Þ ¼ 22
ðPV D
cðxÞnðxÞ·7Gðx; x0 ÞdSðxÞ þ 2sðcÞGðx0 ; xs Þ; ð2:6Þ
where c is an eigenfunction. The spectrum of eigenvalues of this integral equation is discussed in Ref. [17]. Alternatively, the velocity field associated with the double-layer potential may be expressed as the curl of a vector potential, A, and the velocity may be expressed in the form uðx0 Þ ¼ 70 f1 ðx0 Þ þ 70 £ Aðx0 Þ 2 sðqÞ70 Gðx0 ; xs Þ;
ð2:7Þ
where 70 denotes the gradient with respect to the evaluation point x0 : It can be shown that the vector potential is given by the integral representation ð zðxÞGðx; x0 ÞdSðxÞ; ð2:8Þ Aðx0 Þ ¼ D
where z is the strength of the vortex sheet corresponding to the double-layer given by
z ¼ n £ 7q;
ð2:9Þ
[22], p. 503. The right-hand side of Eq. (2.9) involves only tangential derivatives of the density distribution q alone. According to the general principles of the generalized vortex method, the evolution of the bubble interface is
317
simulated using a boundary-element method according to the following steps: Step 1. Specify the initial shape of the interface, interfacial distribution of the harmonic potential f; bubble pressure, and pressure at infinity. Step 2. Solve the integral Eq. (2.5) for the double-layer density q; as will be discussed later in this section. Step 3. Evaluate the harmonic potential f over the interface using Eqs. (2.2) and (2.3). Step 4. Compute the tangential component of the velocity over the interface, n £ ðu £ nÞ ¼ ðI 2 nnÞ·u; in terms of the surface gradient 7s f; where I is the unit matrix and dij 2 ni nj is the tangential projection matrix. In practice, the surface gradient is computed by solving a linear system of three equations, which arises by requiring that the projection of 7s f in two tangential directions is equal to the corresponding directional derivatives, whereas the projection in the normal direction is zero. Step 5. Compute the strength of the vortex sheet z from Eq. (2.9). Step 6. Compute the interfacial distribution of the vector potential A from Eq. (2.8). Step 7. Compute the normal component of the interfacial velocity nðn·uÞ in terms of (a) the normal component of the incident flow, (b) tangential derivatives of the vector potential A, and (c) the normal component of the point source corresponding to the third term on the right-hand side of Eq. (2.2). Step 8. Advance the position of marker points distributed over the interface either with the liquid velocity or with the normal component. In the first case, the marker points are Lagrangian point particles of the exterior liquid. Step 9. Update the harmonic potential at the interfacial marker points using Bernoulli’s equation, taking into consideration that the pressure of the liquid at the bubble surface is equal to the bubble pressure, pB ; reduced by the capillary pressure, 2gkm ; where g is the surface tension and km is the interface mean curvature [22], p. 489. For a spherical bubble of radius a; km ¼ 1=a: Specifically, the rate of change of the potential following Lagrangian point particles moving with the liquid velocity is given by Df 1 2g p 2c km 2 B ¼ lul2 þ þ g·x; 2 r Dt r
ð2:10Þ
where D=Dt is the material derivative, r is the liquid density, c is Bernoulli’s constant, and g is the acceleration of gravity. On the other hand, if the marker points move with the normal component of the liquid velocity, the rate of change of the potential is given by df 1 2g p 2c km 2 B ¼ lu·nl2 2 lul2 þ þ g·x: 2 r dt r
ð2:11Þ
When gravitational effects are insignificant across the bubble diameter, the last term on the right-hand sides of Eqs. (2.10) and (2.11) is absent. If, in addition, the flow
318
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
vanishes at infinity, Bernoulli’s constant c is equal to the liquid pressure far from the bubble, c ¼ p1 : On the other hand, if the flow far from the bubble executes steady streaming motion with velocity U; c ¼ p1 þ r 12 lUl2 : Step 10. Update the bubble pressure according to the polytropic equation of state pB 2 pV ¼ pB ðt ¼ 0Þ 2 pV
vB ðt ¼ 0Þ vB
l
;
ð2:12Þ
where vB is the bubble volume, pV is the partial vapor pressure of the ambient liquid inside the bubble, assumed to remain constant during the motion, and the exponent l is determined by the prevailing thermodynamic conditions. Under adiabatic conditions, l is equal to the ratio of the specific heats; for an ideal gas, l ¼ 1:4: In the case of a pure gas bubble considered in the present simulations, pV ¼ 0: Step 11. Update the pressure at infinity according to a specified protocol. The solution of the integral Eq. (2.5) required in Step 2 is found by the method of successive substitutions. This involves guessing the distribution of q; computing the righthand side of Eq. (2.5), and then replacing the guessed distribution with the left-hand side of Eq. (2.5). To expedite the computation, except at the initial instant, the guessed distribution is identified with the converged solution at the previous time step. When this is done, only 3– 10 iterations are necessary to obtain the solution accurate to the eighth decimal place. An alternative would be to explicitly construct the influence matrix corresponding to the right-hand side of Eq. (2.5), and then solve the linear system using a direct or iterative method. We found, however, that the above described method of successive substitutions is both more efficient and easier to implement. In the numerical implementation, the bubble surface is discretized into an unstructured grid of six-node curved triangular elements, and all geometrical variables and surface functions are approximated with quadratic functions over each element with respect to local triangle coordinates. The local interpolation functions are given in Ref. [23]. The interfacial grid is generated by successively subdividing the faces of a regular octahedron into descendant elements, and finally projecting the nodes onto the sphere. To improve the accuracy, the normal vector and mean curvature are averaged at each node over the elements sharing the node, and the surface distributions of the normal vector and mean curvature over each element are subsequently described by quadratic interpolation. A numerical instability of unknown origin undermines the longevity of the simulations, as discussed in Ref. [17]. To allow the computations to proceed for an extended period of time, the surface potential is smoothed by expanding it in a Fourier – Legendre series, and then truncating the series to a certain level, as discussed in Ref. [17]. Numerical experimentation has shown that smoothing by spectrum truncation is an effective and
benign method of regularization, as will be discussed in Section 4.
3. Small-deformation theory Consider a gas bubble rising in a quiescent ambient liquid or translating relative to a stationary or moving ambient liquid at high Reynolds numbers, and assume that the bubble is and remains spherical at all times. In a frame of reference moving with the bubble, the fluid far from the bubble executes streaming motion with velocity U. Elementary fluid mechanics shows that, if viscous effects are negligible and the flow is assumed to be irrotational, the velocity field is given by the gradient of the harmonic potential ! 1 a3 ; ð3:1Þ fðxÞ ¼ U·^x 1 þ 2 r3 where a is the bubble radius, x^ ¼ x 2 xs is the distance of the field point x from the bubble center xs ; and r ¼ l^xl [22]. When gravitational effects are insignificant, the pressure far from the bubble is given by p1 ¼ c 2
1 rlUl2 ; 2
ð3:2Þ
where c is Bernoulli’s constant. However, because the pressure over the exterior surface of the bubble varies according to Bernouli’s equation, the interface cannot remain spherical at all times and will deform by an amount that is determined by the Weber number, We ; ralUl2 =g; where g is the surface tension. Deformation may cause the bubble to translate as a whole along the x axis with a velocity that depends on the Weber number. The dimensionless Bernoulli constant d ¼ ac=g plays an important role in determining changes in the bubble volume and pressure due to the motion. Assume, in particular, that the bubble translates along the x axis, and describe the motion in a frame of reference where the spherical bubble is stationary and the incident flow is given by U ¼ 2Ux ex : If the initial bubble shape is spherical and the initial distribution of the potential is given by Eq. (3.1), the bubble will deform to obtain an oblate spheroidal shape, while translating and exhibiting axisymmetric shape oscillations in the 2-0 mode. In real life, the oscillations will eventually be dampened by viscous dissipation, and the bubble will assume a steady oblate spheroidal shape. Moore [24] showed that, when the deviation from sphericity is small, to first order in We, the minor and major bubble axes are given by ax ¼ að1 2 e Þ and as ¼ að1 þ ð1=2Þe Þ; where
e¼
3 We: 16
ð3:3Þ
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
319
The corresponding bubble axes ratio is
x;
as 9 We þ OðWe2 Þ: .1þ 32 ax
ð3:4Þ
Subsequently, Moore [4] derived the improved nonlinear inverse relation qffiffiffiffiffiffiffiffiffi x3 þ x 2 2 2 We . 2 4=3 2 ½ x arcsec x 2 x2 2 12 ; ð3:5Þ x ðx 2 1Þ3 which was shown by Miksis et al. [6] to provide accurate predictions for Weber numbers on the order of unity.
4. Numerical simulations In the numerical simulations, we consider a bubble with a spherical initial shape of radius a; and specify the initial distribution of the surface potential ! 1 a3 ga ðj 2 1Þðj þ 2Þ 1=2 2 eS fðxÞ ¼U·^x 1 þ 2 r3 jþ1 r Tjlml ðcos uÞcosðmwÞ;
ð4:1Þ
where e S is a dimensionless amplitude and Tjlml are the modified Legendre functions [17]. In all simulations discussed in this section, the bubble surface was discretized into 512 elements defined by 1026 surface nodes, the polytropic exponent l in Eq. (2.12) was set to 1.25, the reduced Bernoulli constant d was set to 100, and time integration was carried out by the second-order RungeKutta method with constant time step Dt0 ¼ 0:02: The dimensionless time indicated by a prime is defined as t0 ¼ tðg=ra3 Þ1=2 : Unless specified otherwise, smoothing was effected by truncating the Fourier –Legendre spectrum to the value of 8. Numerical experimentation has shown that, for the interfacial grid considered, this level of truncation strikes an optimal balance between consistency and regularity of the simulated motion. Fig. 1 illustrates the oscillations of a bubble in the axisymmetric 2-0 mode, corresponding to j ¼ 2 and m ¼ 0; for We ¼ 0:25 and e S ¼ 0:25: The initially upward curves represent the bubble axis along the x axis of symmetry, defined as half the maximum bubble size in this direction, reduced by the initial equivalent bubble radius. The initially downward curves represent the bubble axis normal to the x axis, reduced by the initial equivalent bubble radius. In the simulations represented by the solid and hardly visible dotted line, the marker points move with the velocity normal to the interface. In the case of the dotted line, smoothing is not enabled, whereas in the case of the solid line, smoothing is applied after each time step. In the simulations represented by the dashed line and circular symbols, the marker points move with the liquid velocity. In the case of the circular symbols, smoothing is not enabled, whereas in the case of the dashed line, smoothing is
Fig. 1. Axisymmetric oscillations of a bubble in the 2-0 mode, for We ¼ 0:25 and e S ¼ 0:25; computed with different options. (a) Evolution of the reduced x axis, and (b) deformed grid developing when the nodes move with the fluid velocity.
applied after each time step. While the agreement between all simulations is excellent during the initial stage of the motion, only the simulation represented by the solid line survives for an extended period of time. Note, in particular, that because the interfacial marker points are convected toward the rear of the bubble under the influence of the streaming flow, severe grid distortion occurs when the nodes are identified with material point particles moving with the fluid velocity, and the computation breaks down at an early stage of the simulation. Fig. 1(b) depicts a deformed grid at the time when the simulation fails, corresponding to the dashed line in Fig. 1(a). Moore’s [24] small-deformation theory predicts that, in the absence of oscillations, a bubble translating at We ¼ 0:25 has a steady spheroidal shape with minimum axis ax ¼ 0:9531a; maximum axis as ¼ 1:0234a; and axes ratio as =ax ¼ 1:0234; where a is the equivalent bubble radius. Consider this Moore spheroid, and assume that the distribution of the potential over the interface at the initial instant is described by Eq. (4.1) with e S ¼ 0: The heavy solid lines in Fig. 2(a) and (b) show the evolution of the longitudinal and transverse bubble axes, ax and as ; reduced by the initial bubble equivalent radius a0 ; the heavy solid line in Fig. 2(c) shows the evolution of
320
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
Fig. 2. Axisymmetric oscillations of a bubble in the 2-0 mode. The heavy solid lines correspond to the Moore spheroid for We ¼ 0:25: (a) Bubble axis in the x direction, (b) bubble axis in the transverse direction, and (c) bubble volume reduced by the initial value.
the bubble volume. In order to sustain the regularity of the grid for an extended period of time, the marker points are required to move normal to the interface. The numerical results reveal that the bubble exhibits small-amplitude axisymmetric periodic oscillations about the Moore spheroid. To demonstrate the effect of the Weber number on the finite-amplitude motion, we consider the evolution of an
initially spherical bubble, and set the initial distribution of the potential as shown in Eq. (4.1) with e S ¼ 0:15 for the 2-0 mode. The solid and dashed lines in Fig. 2 show results of simulations, respectively, for We ¼ 0 and 0.25. In fact, the dashed line in Fig. 2(a) extends the solid line displayed earlier in Fig. 1(a) to longer times. Consider first the free oscillations of a stationary bubble immersed in a quiescent fluid, We ¼ 0: According to small deformation theory, the x axis of the bubble, reduced by the mean bubble radius, oscillates around the mean value of unity with dimensionless period T 0 ¼ 1:814 and reduced amplitude 0.0946, corresponding to the horizontal and vertical lines in Fig. 2(a), while the transverse axis oscillates with an 1808 phase difference with half the reduced amplitude 0.0473, corresponding to the horizontal and vertical lines in Fig. 2(b). The bubble volume exhibit simultaneous oscillations with dimensionless period T 0 ¼ 0:322; as shown by the vertical lines in Fig. 2(c). These predictions are borne out by the numerical results shown with the solid lines in Fig. 2. Note that, even though the volume mode is not excited at the initial instant, volume oscillations do occur because of nonlinear interactions. Comparing the solid with the solid lines in Fig. 2(a) and (b), we see that raising the Weber number increases the period, and thus reduces the frequency of the shape oscillations, in agreement with theoretical predictions for small-amplitude motion [15,16]. Fig. 3(a) describes the evolution of the x position of the bubble centroid for We ¼ 0:25; corresponding to the dashed lines depicted in Fig. 2. The results show that a bubble that would be stationary in the spherical state, moves opposite to the direction of the oncoming flow. Conversely, a rising bubble slows down as a result of the deformation. The net motion effectively alters the Weber number of the incident flow by an amount that is a priori unknown and must be found as part of the solution. Fig. 3(b) displays the bubble interface as described by the boundary-element grid at time t0 ¼ 3:95; the incident streaming flow is directed from left to right toward the negative direction of the x axis. Similar results are obtained for the three-dimensional 2-2 mode, corresponding to j ¼ 2 and m ¼ 2; where the bubble surface oscillates between two triaxial ellipsoidal shapes whose major axes lie in the equatorial xy or xz plane. Fig. 4(a) –(c) shows the evolution of the bubble axes, defined as half the maximum bubble size in the corresponding direction reduced by the initial bubble radius, for e S ¼ 0:20; and We ¼ 0 (solid lines) or 0.25 (dashed lines); in the second case, the bubble translates along the z axis. Fig. 5(a) displays oscillations in the bubble volume excited by nonlinear interactions. The numerical results are consistent with the predictions of the linear theory for small-amplitude shape and volume oscillations represented by the vertical and horizontal lines in Fig. 4(b) and (c). Note that in the small-amplitude linearized motion, the x axis of the bubble described in Fig. 4(a) remains unchanged.
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
321
Fig. 3. Axisymmetric oscillations of a translating bubble in the 2-0 mode for We ¼ 0:25; corresponding to the simulation represented by the dashed lines in Fig. 2; (a) reduced x position of the bubble centroid, x0 ¼ x=a; and (b) instantaneous bubble surface at time t0 ¼ 3:95; the incident streaming flow is directed from right to left toward the negative direction of the x axis.
The simulation of the 2-2 mode oscillations reveal that, as in the case of the 2-0 mode, raising the Weber number increases the period, and thus reduces the frequency of the shape oscillations. Furthermore, in the case of a translating bubble, oscillations occur about a mean spheroidal shape with two major axes perpendicular to the direction of the motion. Fig. 5(b) shows the evolution of the z position of the bubble centroid for We ¼ 0:25; corresponding to the dashed lines in Fig. 4. The results show that a rising bubble slows down as a result of the deformation. Finally, Fig. 5(c) displays the bubble interface as described by the boundary-element grid at time t0 ¼ 3:95; the incident streaming flow is directed from top to bottom toward the negative direction of the z axis.
5. Discussion The completed double-layer representation combined with Baker’s generalized vortex method has been shown
Fig. 4. Oscillations of the bubble axis in the 2-2 mode for We ¼ 0 (solid lines) or 0.25 (dashed lines).
to provide us with an efficient algorithm for simulating three-dimensional potential flow due to the motion or past compressible bubbles. Numerical instabilities are a serious impediment to the longevity of the simulations over an extended period of time, but smoothing by
322
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323
Fig. 5. Oscillations of a translating bubble in the 2-2 mode, corresponding to the simulation represented by the dashed lines in Fig. 4. (a) Oscillations in the bubble volume, (a) reduced z position of the bubble centroid, z0 ¼ z=a; and (c) instantaneous bubble surface at time t0 ¼ 3:98; the incident streaming flow is directed from top to bottom toward the negative direction of the z axis.
spectrum truncation removes numerical irregularities while maintaining the fundamental features of the nonlinear motion. The superiority of the method compared to direct methods based on Green’s third identity lies in its ability to compute the instantaneous flow by solving an integral equation of the second kind, which can be done accurately and economically by the method of successive substitutions. An important aspect of the boundary-element implementation is the option for the boundary-element nodes to be convected normal to the interface with the liquid velocity and tangentially with an arbitrary component, thereby ensuring the regularity of the interfacial grid for an extended period of time. Though regridding is an alternative [25], the introduction of numerical error due to interpolation in parameter space or directly in three-dimensional physical space is a serious concern. The algorithm has been applied with success in this and in the earlier work [17] to simulate bubble oscillations in infinite flow. Although not tested, the method is likely to perform equally well in the case of flow bounded by an impenetrable boundary such as an infinite plane wall. The only necessary modification is the substitution of the free-space Green’s function with a Neumann function appropriate for the geometry of the flow under
consideration. A library of Green and Neumann functions accompany a recent text [26] and is available from the internet site http://bemlib.ucsd.edu
References [1] Wegener PP, Parlange JY. Spherical-cap bubbles. Annu Rev Fluid Mech 1973;4:79–100. [2] Miksis MJ, Vanden-Broeck JM, Keller JB. Rising bubbles. J Fluid Mech 1982;123:31–41. [3] Batchelor GK. The stability of a large gas bubble rising through liquid. J Fluid Mech 1987;184:399 –422. [4] Moore DW. The velocity of rise of distorted gas bubbles in a liquid of small viscosity. J Fluid Mech 1959;23:749 –66. [5] Harper JF. The motion of bubbles and drops through liquids. Adv Appl Mech 1972;12:59–129. [6] Miksis MJ, Vanden-Broeck JM, Keller JB. Axisymmetric bubble or drop in a uniform flow. J Fluid Mech 1981;108:89–100. [7] El Sawi M. Distorted gas bubbles at large Reynolds number. J Fluid Mech 1974;62:163–83. [8] Benjamin TB. Hamiltonian theory for motion of bubbles in an infinite liquid. J Fluid Mech 1987;2:118– 28. [9] Duineveld PC. Rise velocity and shape of bubbles in pure water. J Fluid Mech 1995;292:325 –32. [10] Magnaudet J, Eames I. The motion of high-Reynolds number bubbles in inhomogeneous flows. Annu Rev Fluid Mech 2000;32: 659 –708.
C. Pozrikidis / Engineering Analysis with Boundary Elements 28 (2004) 315–323 [11] Wu M, Gharib M. Experimental studies on the shape and path instability of small air bubbles rising in clean water. Phys Fluids 2002; 7:L49–L52. [12] Saffman PG. On the rise of small air bubbles in water. J. Fluid Mech 1956;1:249–75. [13] Lunde K, Perkins R. Shape oscillations of rising bubbles. Appl Sci Res 1998;58:387–408. [14] Hartunian RA, Sears WR. On the instability of small gas bubbles moving uniformly in various liquids. J Fluid Mech 1957;3: 27 –47. [15] Meiron DI. On the stability of gas bubbles rising in an inviscid fluid. J Fluid Mech 1989;198:101–14. [16] Feng JQ. The oscillations of a bubble moving in an inviscid fluid. SIAM J Appl Math 1992;52:1–14. [17] Pozrikidis C. Numerical simulation of three-dimensional bubble oscillations by a generalized vortex method. Theor Comput Fluid Dyn 2002;16:133– 50. [18] Baker GR. Generalized vortex method for free-surface flow problems. In Waves in fluid interfaces MRC, Wisconsin: University of Wisconsin; 1983.
323
[19] Baker GR, Meiron DI, Orszag SA. Generalized vortex methods for free-surface flow problems. J. Fluid Mech 1982;123: 477–501. [20] Baker GR, Meiron DI, Orszag SA. Boundary integral methods for axisymmetric and three-dimensional Rayleigh –Taylor instability problems. Physica 1984;12D:19–31. [21] Mikhlin SG. Integral equations and their applications to certain problems in mechanics, mathematical physics and technology. New York: Pergamon Press; 1957. [22] Pozrikidis C. Introduction to theoretical and computational fluid dynamics. Oxford University Press; 1997. [23] Pozrikidis C. Numerical computation in science and engineering. Oxford University Press; 1998. [24] Moore DW. The rise of a gas bubble in a viscous liquid. J Fluid Mech 1959;6:113 –30. [25] Kwak S, Pozrikidis C. Adaptive triangulation of evolving, closed or open surfaces by the advancing-front method. J Comp Phys 1998;145: 61–88. [26] Pozrikidis C. A practical guide to boundary element methods with the software library BEMLIB. Chapman & Hall/CRC Press; 2002.