A boundary element method for calculating the shape and velocity of two-dimensional long bubble in stagnant and flowing liquid

Engineering Analysis with Boundary Elements 30 (2006) 539–552 www.elsevier.com/locate/enganabound

A boundary element method for calculating the shape and velocity of two-dimensional long bubble in stagnant and flowing liquid H. Ha-Ngoc a,*, J. Fabre b b

a Institute of Mechanics, 264 Doi Can, Hanoi, Vietnam Institut de Me´canique des Fluides, Alle´e du Pr. Camille Soula, 31400 Toulouse, France

Received 13 July 2005; accepted 15 February 2006

Abstract A numerical method based on a boundary element method (BEM) has been developed for computing the velocity and the shape of long bubbles moving steadily in stagnant and flowing liquid in 2D case: plane and axisymmetrical. The flow is assumed to be inviscid and incompressible. The method consists in solving simultaneously a Poisson equation characterizing the flow and an equation for bubble shape in the form of a functionaldifferential equation resulting from both Bernoulli equation and the jump conditions at the interface. The Poisson equation is solved by a BEM with an iterative loop for nonlinear source term while the system of nonlinear algebraic equations obtained by discretizing the equation on the interface is solved by the Powell’s hybrid algorithm. The bubble shape and velocity are obtained as a part of the solution. The problem of multiple solutions is investigated numerically and the maximum velocity criterion is used for selecting the physical solution. The results obtained by the simulation are in good agreement with the experimental and numerical results of previous studies. q 2006 Elsevier Ltd. All rights reserved. Keywords: Long bubble; Boundary elements; Two-phase flow

1. Introduction The motion of long bubbles in two-phase flow of gas and liquid occurs in variety of industrial systems including geothermal wells, oil and gas wells and pipelines, steam generating boilers, etc. For some operating conditions, the gas is mainly contained in long bubbles. In vertical pipe, they have a cylindrical bullet-like shape and they are frequently referred to as Taylor bubbles since the early work of Davies and Taylor [1]. In horizontal or inclined flow, the bubbles are asymmetrical and they are sometimes referred to as Benjamin bubbles [2]. The flow characteristics depend then critically on their shape and velocity. Theoretical and experimental investigations of these long bubbles in circular tube or rectangular channel has been a subject of great interest for about a halfcentury. Experiments have put in evidence the effects of fluid properties and duct inclination on the shape and velocity of these bubbles [3–5]. In general, Taylor bubbles are * Corresponding author. Tel.: C84 4 7628807; fax: C84 4 8333039. E-mail address: [email protected] (H. Ha-Ngoc).

0955-7997/$ - see front matter q 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.02.007

axisymmetrical and Benjamin bubbles in narrow rectangular channel can be considered as plane bubbles so in the purpose of mathematical modeling they have 2D character. In stagnant liquid and for the inertial regime corresponding to low viscosity, the translational velocity V can be expressed as V pffiffiffiffiffiffi Z Fr gD

(1)

where the Froude number Fr, in which D is the tube diameter or the channel width and g the gravitational acceleration, is constant. In the case of low surface tension, the experiments showed that FrZ0.35 for circular tube [3,6–8], and FrZ0.23 for 2D channel [9,10] for vertical case. The influence of surface tension and inclination on bubble motion was investigated experimentally by many authors such as Zukoski [3], Bendiksen [11] and Hasan and Kabir [12]. The results showed that the velocity is a decreasing function of surface tension. Moreover, it has a ‘peculiar’ behavior with the inclination. The velocity increases as the inclination declines from 908 (vertical position) to about 30–458, then it decreases when the inclination goes to 08 (horizontal position). In flowing liquid, the velocity of a bubble must depend on the velocity of the liquid flowing upstream as well as on the

540

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

translation due to buoyancy. Nicklin et al. [7] suggested from their experiments that the bubble velocity V can be expressed as V Z C0 U C VN

(2)

where U is the average liquid velocity, VN is the bubble velocity in stagnant liquid, and C0 is an empirical coefficient that accounts for the velocity distribution in the liquid ahead the bubble tip. This coefficient is thus sensitive to the flow regime in the flowing liquid. Nicklin et al. [7] have shown that for upward vertical flow, C0 depends on the liquid Reynolds number ReL (ReLZUD/nL) in decreasing from about 1.9 for ReLZ100 to an asymptotic value of about 1.2 for ReLO8000. The bubble velocity was observed to be independent of the viscosity for a fairly wide range of low viscosities. According to White and Beardmore [8] and Wallis [13], its effects can be negligible when NfO300, where Nf is the inverse viscosity number defined as NfZ(gD3)1/2/nL, or according to Zukoski when ReBO100, where ReB is the bubble Reynolds number defined as ReBZVD/nL. Neglecting the viscosity is justified by observing that the growth of the boundary layer created at the wall and at the bubble surface, will be slow, and also prevented by the acceleration of the liquid around the bubble nose. These results suggested that in the case of low viscosity potential flow theory could be applied. Simplified theoretical approaches assuming inviscid fluid were undertaken long before systematic experimental results were available. Dumitrescu [14] and David and Taylor considered the potential flow around an axisymmetrical bubble rising in a vertical tube. They found that the approximate theoretical solution of the stream function is given as a truncated series of Bessel function. With a constant-pressure condition at the interface around the bubble vertex, they found FrZ0.351 and 0.328, respectively, when surface tension was ignored. Collins et al. [15], Bendiksen [16] and Nickens and Yannitel [17] have extended the method to the case of flowing liquid and nonzero surface tension. Although these approaches could not provide the exact bubble shape, they agreed reasonably well with the experimental results [3,7,11] for the bubble velocity. For the case of a plane bubble, the theoretical analyses using potential flow have been taken by Birkhoff and Carter [18], Garabedian [19], Collins [10] and Vanden-Broe¨ck [20,21] for vertical flow with zero surface tension. Coue¨t and Strumolo [22] have generalized the technique employed by Birkhoff and Carter and Vanden-Broe¨ck to include the effects of surface tension and tube inclination. The method was based on conformal mapping and could apply only for the case of stagnant liquid or irrotational flow. Garabedian [19] was the first to point out the existence of multiple solutions. According to his conclusion, the solution would exist for each value of the Froude number less than a critical value Frc. He made a stability analysis based on energy approach and concluded that the maximum bubble velocity corresponding to the critical Froude number Frc should be accepted as a physical solution. Vanden-Broe¨ck and Coue¨t and Strumolo have shown that for

nonzero surface tension, a countable number of solutions exist and they used the same criterion to isolate the physical solution. For the case of horizontal channel, Benjamin has developed a method that uses a global momentum balance to determine the bubble velocity and film thickness without solving the flow field. Ha-Ngoc [23] has extended the method for the case of flowing liquid. However, the method can be applied only for zero surface tension. Numerical solutions of the full Navier–Stokes equations based on the finite difference method have been obtained for the motion of long bubbles in tube [24–27]. They take into account the effects of viscosity and surface tension. However, they are usually costly in computer time and they are limited to 2D flows: 3D calculations for asymmetrical bubbles are not yet available. It must be stressed in passing that Mao and Dukler [24,25] observed in their calculations that multiple solutions can be found. The authors proposed a simple geometric criterion to single out the physical solution: the bubble tip has to be spherical. BEM techniques have been applied to a wide variety of free surface problems including rising bubbles [28–30], liquid jets [31–33] and droplets [34,35]. One important advantage of the BEM compared to other methods is that nonlinearity of the gas/ liquid interface condition can be satisfied more accurately. Indeed, once the BEM solution has been found, the interfacial velocity is obtained directly from solution in terms of normal and tangential derivatives of the velocity potential or stream function. Moreover, with the BEM the optimal placement of the nodes on the interface provides a mean for maximizing accuracy of surface curvature calculations. Consequently, it gives more accurate capture of the interface shape. In this paper, a method based on the boundary element method has been developed to solve the equations for irrotational and rotational inviscid flow of a liquid past a bubble at large Reynolds number, in the 2D case, i.e. for plane and axisymmetrical bubbles. Differently from other works on rising bubbles using BEM [28–30], this work deal with both irrotational and rotational liquid flow, so instead of velocity potential, the stream function is to be used. In general, the boundary elements are then used in an iterative way for the Poisson equation of the stream function with nonlinear source term. This equation is solved together with an equation at the interface so that the bubble shape and velocity can be determined as a part of the solution. In Sections 2 and 3, the mathematical formulations of the problem and the solution procedure are described. Afterwards, the problem of multiple solutions and the convergence of the numerical procedure with respect to the bubble length and to the discretization are investigated. The results for the bubble velocity and shape in stagnant and flowing liquid are presented in Section 5 in comparison with experimental, theoretical and other numerical results. 2. Mathematical formulation of the problem For the purpose of clarity, the mathematical formulation for the plane case is presented in detail in the following paragraph.

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

The equations for the axisymmetrical case that can be obtained in a similar way will be presented afterwards more briefly. 2.1. Plane case Let us consider a long bubble moving with a constant velocity V relative to the channel. For a bubble whose length is greater than, say twice the channel width, experiments have shown that its motion is independent of its length. Thus, from a mathematical point of view the bubble length can be considered as infinite. The geometry as well as the main variables of the problem are presented in Fig. 1. The moving coordinate system attached to the bubble is used to formulate the problem. In this system, the bubble is at rest and the flow is steady provides that the bubble has reached a constant velocity. The coordinate system is chosen so that the x-axis coincides with the channel symmetry axis and the y-axis passes through the stagnation point S on the bubble. The channel is inclined from the horizontal by an angle a. In general, the bubble is asymmetric so that the stagnation point does not coincide with the bubble vertex and does not lie on the channel symmetry axis but is at a distance e from it (Fig. 1a). Far upstream the bubble, the liquid is assumed to flow as if it was alone in the channel. Its velocity profile is symmetrical with a maximal value uc on the centreline. As pointed out by Dumitrescu [14] and Davies and Taylor [1] under the inviscid assumption, the flow is irrotational for a bubble moving in still liquid. However, for flowing liquid, the flow is rotational and the distribution of vorticity far upstream the bubble corresponds to that of the liquid flowing alone in the channel. For large enough Reynolds number, the inviscid assumption is expected to apply in the bubble nose region. Indeed, the diffusion of vorticity by the viscosity is a rather slow process at large ReL. For a bubble moving in still liquid, it creates thin boundary layers that develop slowly at the bubble surface and at the wall. For the case of a bubble in a flowing liquid, the main role of viscosity is to develop the flow and to fix the vorticity ahead of the bubble nose [15]. The theoretical basic for such a framework has been well established by Lamb [36]

541

and Batchelor [37] and it has been applied successfully in various problems, e.g. the flow in axial compressors and turbines [38]. Collins et al. [15] and Bendiksen [16] used the theory successfully for the prediction of velocity of bubble rising in flowing liquid in vertical tubes. 2.2. Equation for the stream function The steady 2D flow of an inviscid, incompressible fluid can be characterized by stream function j(x, y) in terms of which the velocity components express as ux Z

vj ; vy

vj uy ZK vx

(3)

where ux, uy are velocity components in x- and y-directions, respectively. In 2D flow, there is no vorticity stretching and the only component of vorticity is determined by: uZ

vuy vux K : vx vy

(4)

Replacing Eq. (3) into Eq. (4) we obtain: v 2 j v2 j C 2 ZKu: vx2 vy

(5)

In a steady inviscid flow, the vorticity is simply transported by the flow. It remains constant along each streamline and can be expressed as a function of j alone: u Z f ðjÞ

(6)

We obtain then the equation for the stream function j in form of a Poisson equation: v 2 j v2 j C 2 ZKf ðjÞ: vx2 vy

(7)

Eq. (7) will be solved with the boundary conditions at the channel walls, at the interface and the flow conditions far upstream and downstream. The boundary conditions at the channel walls and the function f(j) will be determined by a given velocity profile far upstream. 2.2.1. Velocity profiles and function uZf(j) The function f(j) is obtained from the conditions imposed upstream at xZKN. The velocity distribution ua in the fixed coordinates is indeed known there and assumed to depend only on y. Hence, in the moving coordinates attached to the bubble, the velocity at xZKN is: uðKN;yÞ Z V Kua ðyÞ:

(8)

By choosing the arbitrary constant, j(KN, 0)Z0, the stream function and the vorticity at xZKN can be expressed as: ðy Fig. 1. Coordinate system for a gas bubble in a channel: (a) fixed coordinates attached to the channel wall; (b) moving coordinates attached to the bubble.

jðKN;yÞ Z uðKN; tÞdt 0

(9)

542

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

duðKN;yÞ : uðKN;yÞ ZK dy

(10)

It is generally not possible to eliminate the coordinate y between Eqs. (9) and (10) to get an explicit form of Eq. (6). However, by taking into account Eq. (8) this relation can be put in the following parametric form that is valid in the whole flow domain 8 ðY > > > > > < jðYÞ Z VY K ua ðtÞdt; (11) o > > dua ðYÞ > > > : uðYÞ Z dY ; where Y is a parameter that has the same dimension as y and that belongs to the segment [Ka,Ca]. In the case of a bubble moving in stagnant liquid, the velocity profile far upstream is uniform and has a value equal to the bubble velocity. The flow is irrotational, therefore, we have: u Z f ðjÞ Z 0:

(12)

In the case of flowing liquid, relation (6) can be obtained by solving (11) numerically (Section 3.1). 2.2.2. Boundary conditions The boundary conditions for the stream function j can be established as follows: (a) The channel walls are streamlines along which j is constant. At the channel walls, we have thus the Dirichlet conditions: jðx;aÞ Z jðKN;aÞ

jðx;KaÞ Z jðKN;KaÞ

(13)

In general, the values j(KN, a), j(KN, Ka) can be determined by the integral (9). However, for a symmetrical velocity profile of mean velocity U these values can be expressed simply by: jðKN;aÞ Z aðV KUÞ

jðKN;KaÞ Z aðU KVÞ

(14)

(b) The bubble surface is also a streamline on which the stream function is constant. So at the interface, we have also the Dirichlet type condition jjint Z Jcte

(15)

where the constant Jcte is the value of the stream function at the interface. In general, this value is unknown and it will be looked for as a part of the solution.

2.2.3. Bernoulli equation at the interface We assume that the gas density is small with respect to the liquid one so that the movement of the gas can be neglected and the pressure in the bubble considered as constant. Suppose that the interface is defined by xZX(y) the Bernoulli equation between the stagnation point and any point [X(y), y] of the interface can be written 1 2 p p u KgXðyÞsin aKgðy C eÞcos a C Z 0 (17) 2 r r qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where uZ u2x C u2y is the absolute velocity of the liquid and subscript 0 refers to the stagnation point. The pressure of the liquid at the interface is related to the gas pressure pB by the Kelvin–Laplace equation: pB Z p C sK:

(18)

By using Eq. (18) the pressure can be eliminated from Eq. (17) and we get 1 2 s u KgXðyÞsin aKgðy C eÞcos a C ðK0 KKÞ Z 0 2 r

(19)

where s is the surface tension and K, K0 are the curvature of the interface at the point of interest and at the stagnation point, respectively. The curvature at a point [X(y), y] can be determined by: KðyÞ Z

X 00 ðyÞ : f1 C ½X 0 ðyÞ
2 g3=2

(20)

2.2.4. Equations in dimensionless form pffiffiffiffiffiffiffiffi To obtain dimensionless equations we use 2ga and a, the channel haft-width, as velocity and length scales. The dimensionless quantities are defined as follows: x xZ ; a

y hZ ; a

j j Z pffiffiffiffiffiffiffiffiffiffi ; 2ga3

ux ffi; ux Z pffiffiffiffiffiffiffi 2ga ua u Z pffiffiffiffiffiffiffiffi : 2ga

uy uy Z pffiffiffiffiffiffiffiffi ; 2ga (21)

For convenience, the stars will be dropped whenever there is no risk of confusion. The dimensionless equations introduce two important dimensionless numbers, the bubble and Froude number Fr and the dimensionless surface tension S defined by: V Fr Z pffiffiffiffiffiffiffiffi ; 2ga

SZ

s : rga2

(22)

The Poisson equation (7) becomes: (c) At upstream or downstream infinity, it may be assumed that the velocity component uy is zero. It yields the following Neumann condition: lim

x/GN

vjðx;yÞ Z 0: vx

(16)

v2 j v2 j C 2 ZKf ðjÞ: vh vx2

(23)

The boundary conditions (13) at the channel walls become:

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

(24) jðx;K1Þ Z jðKN;K1Þ ZKFr C U

ðs

ðs

jðx;1Þ Z jðKN;1Þ Z FrKU

sin a cosðqðtÞÞdt C cos a sinðqðtÞÞdtKS½KðsÞKKð0Þ
0

At the interface and at infinity Eqs. (15) and (16) take the form: jðzðhÞ;hÞ Z Jcte ;

543

0

Z u2 ðsÞ (30)

(25) By taking the derivative with respect to s, we obtain:

vjðx;hÞ Z 0: lim x/GN vx

(26)

sin½a C qðsÞ
KSK 0 ðsÞ Z 2uðsÞu 0 ðsÞ:

(27)

Note that the velocity at the interface can be expressed in terms of the stream function j by    vj  (32) u Z  ; vn

The Bernoulli Eq. (19) becomes u2 K½zðhÞsin a C ðh C eÞcos a
C SðK0 KKÞ Z 0 with KðhÞ Z

z 00 ðhÞ : f1 C ½z 0 ðhÞ
2 g3=2

(28)

2.2.5. Bernoulli equation in curvilinear coordinate The interface coordinates (xI, hI) can be expressed as functions of the arc length s of the interface starting with sZ0 at the stagnation point S. Let q be the angle of the tangent at the interface with respect to the channel wall (Fig. 2) then we have ðs

xI ðsÞ Z cos½qðtÞ
dt; 0

ðs

hI ðsÞ ZKe C sin½qðtÞ
dt:

(29)

0

where s belongs to the interval ]KN,CN[ and q to [0, p]. In the curvilinear coordinates, Eq. (27) can be written as follows:

(31)

where n is the vector normal to the interface. Moreover, the curvature of the interface defined by Eq. (28) can be expressed by: KðsÞ ZKq 0 ðsÞ:

(33)

Thus, by introducing the stream function j, Eq. (31) can be written:  0 vj vj ðsÞ sin½a C qðsÞ
KSq 00 ðsÞ Z 2 ðsÞ: (34) vn vn This equation is a functional-differential type equation for q(s) because the right-hand side determined by the Poisson equation (23) is a function of the function q(s). For short, this equation will be referred to as the interface equation. Eq. (34) is a functional-differential equation of second order for q. We expect that on the haft-line [0, N[ it will have a unique solution with two conditions at infinity sZN. From a mathematical point of view, the problem seems to be difficult. Nevertheless, a similar problem has been solved by Vied and Bajalieva [39]. Further theoretical investigations are behind the scope of this article. We assume the existence and uniqueness of the solution of Eq. (34) to investigate the problem numerically. 2.4.6. Boundary conditions for interface equation and condition for bubble velocity We assume that the interface is smooth, parallel to the channel wall and that it does not oscillate at infinity. Thus, the angle of its tangent equals to 0 at sZCN and to p at sZKN (Fig. 2). Furthermore, it can be expected that the curvature and its derivative are equal to zero. Thus we have: lim qðsÞ Z 0;

s/CN

lim qðsÞ Z p:

s/KN

(35)

Taking into account Eq. (33), the conditions that the curvature and its derivative vanish are equivalent to: lim q 0 ðsÞ Z 0

(36)

lim q 0 ðsÞ Z 0

(37)

s/KN

Fig. 2. Definition of curvilinear coordinates.

s/CN

544

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

and 00

lim q ðsÞ Z 0:

s/GN

(38)

By combining Eqs. (35) and (38) with Eq. (34) we obtain:  0 vj vj ðsÞ ðsÞ (39) sin a ZK lim 2 s/KN vn vn  0 vj vj ðsÞ ðsÞ: sin a Z lim 2 s/CN vn vn

(40)

The conditions (39) and (40) are the results of combining the purely geometrical conditions (35) and (38) and the dynamic equation (34). We will use these conditions instead of Eqs. (35) and (38) in our numerical procedure and we will see further that they can make a good convergence. We distinguish two different geometrical cases of the problem: Case (i). The case when the bubble leans against the upper channel wall. This case will cover the situations when the channel inclination lies between 08 (horizontal) and a large angle (a!708 according to Maneri and Zuber [5]). In this case, the interface is part of the streamline that belongs to the upper wall. The value Jcte can be, therefore, determined from the velocity profile given upstream. The point of intersection between the interface and the wall is a stagnation point. For the interface equation, two boundary conditions are being given at infinity sZCN. The bubble velocity can be determined from a given angle q of the interface at sZ0, i.e. qð0Þ Z q0 ;

(41)

where q0 is the given contact angle of the triple-line at the stagnation point. Case (ii). The case when the bubble does not lean against the channel wall. This case will cover the situations when the channel is vertical or nearly vertical (aO808 according to Maneri and Zuber [5]). In contrast to the first case, the interface is an unknown streamline and the value Jcte must be calculated simultaneously with the interface shape and the bubble velocity. The interface equation can be solved for two branches ]KN, 0] and [0,CN[ with two boundary conditions given at each infinity sZKN and CN. The position of the stagnation point (the distance from channel axis e) should be determined by the condition of vanishing velocity. By taking account Eq. (32) this condition is equivalent to: vj ð0Þ Z 0: vn

(42)

The conditions of continuity of both the slope angle q and the curvature K are used to determine the value of the stream function at the interface Jcte and the bubble velocity Fr. They are: lim qðsÞ Z lim qðsÞ

(43)

lim KðsÞ Z lim KðsÞ

(44)

s/0C

s/0C

s/0K

s/0K

For the vertical case, the bubble does not lean against the channel walls but it can be assumed symmetrical. Hence, the interface is the continuation of the streamline that belongs to the axis of the channel, so the problem can be also considered as a special case (i) with q(0)Zp/2. 2.3. Axisymmetrical case The cylindrical coordinates with longitudinal coordinate x and radial coordinate y are used to formulate the problem. This case is always considered as case (i) with slope angle q(0)Z p/2 at the stagnation pffiffiffiffiffiffi point. The dimensionless equations are obtained with gD and RZD/2 as velocity and length scales where D is tube diameter. 2.3.1. Equation for the stream function v2 j v2 j 1 vj ZKh2 f ðjÞ: C 2K h vh vh vx2

(45)

where f(j)Zu/h is determined by the velocity distribution given upstream and u is the only nonzero vorticity component in axisymmetrical flow. This equation can be written as a Poisson equation in the plane (x, h) with new source term depending also on the derivative of the stream function: v2 j v2 j 1 vj Kh2 f ðjÞ: C 2Z 2 h vh vh vx

(46)

2.3.2. Interface equation cosðqÞ C SK 0 Z 2uI uI0

(47)

with uI, the absolute velocity at the interface is given by    1 vj  ðsÞ (48) uI ðsÞ Z  hI ðsÞ vn and the curvature derivative is given by K 0 ðsÞ ZKq 00 ðsÞK

q 0 ðsÞhI ðsÞ C cos½qðsÞ
sin½qðsÞ
: h2I ðsÞ

(49)

2.3.3. Boundary conditions The boundary conditions for the stream function and the interface equation can be obtained similarly as in the plane case for the case (i). 2.4. Summary We seek two functions, the stream function j(x, h) and the interface slope angle function q(s). For the case (i), the only constant to find is the Froude number Fr, but for the case (ii) three constants must be found: the Froude number Fr, the stream function value at the interface Jcte and the position of stagnation point e. The equations for the stream function and the interface slope angle for both plane and axisymmetrical cases can be written in the following general form:

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

– Poisson equation in plane (x, h)   v2 j v2 j vj C 2 Z P h;j; vh v2 x v h

where (50)

(51)

where L is a second order differential operator and = is a functional determined by derivation of the velocity square at the interface. The boundary conditions and the conditions for the unknown constants can be obtained by the considerations given above. 3. Solution method We seek a numerical solution for the two functions j and q and three constants Fr, Jcte and e satisfying the set of nonlinear equations (50) and (51) and the boundary conditions specific to each case (i) or (ii) as described in Section 2. For the numerical solution, the conditions at infinity are replaced by the same conditions applied at constant x far enough from the origin so that they do not influence the solution. The Poisson equation (50) is solved by a BEM [40] with an iterative loop for the nonlinear source term. Linear elements are used to approximate the interface and other boundaries. For Eq. (51), a finite difference scheme is used and the resulting set of nonlinear algebraic equations is solved by the Powell’s hybrid algorithm [41]. We describe shortly the numerical procedure in the following paragraphs. More details can be found in Ha-Ngoc [23]. 3.1. Determination of uZf(j) and limitations The numerical solution of Eq. (11) that gives the vorticity uZf(j) can be done provided that the function j(Y) varies monotonically within [Ka,Ca]. For each j, we can determine a Y uniquely by solving the first equation of the set (11). Then the value of u can be obtained by replacing Y in the second equation. The condition for the function j(Y) to be monotonic requires that the velocity in the moving frame given by Eq. (8) does not change its sign. In the case of co-current flow (uaO0), this happens when the bubble velocity is greater than the maximum velocity, i.e. when VOuc. In the case of countercurrent flow (ua!0), j(Y) is monotonic when the bubble velocity is positive, i.e. VO0. By increasing the mean velocity ua, the bubble velocity may become less than the maximum velocity uc (V!uc) in co-current flow or negative (V!0) in counter-current flow. Hence, for a given velocity profile, the problem can be solved if the mean velocity lies between two values for which we have VZuc and VZ0. In dimensionless formulation, the problem can be solved if the liquid Froude number FrL is such that C FrK C ! FrL ! FrC

U FrC C Z pffiffiffiffiffiffiffiffi 2ag

(53)

for co-current flow and when VZuc,

– Functional-differential equation Lðs;q;q 0 ;q 00 ;a;SÞ Z =ðs;q;Fr;Jcte ;eÞ

545

(52)

U FrK C Z pffiffiffiffiffiffiffiffi 2ag

(54)

for counter-current flow and when VZ0. K The numbers FrC C , FrC depend on the shape of the velocity profile, the geometry and the liquid properties. They may be determined from simulation results. 3.2. Numerical procedure for the interface equation and for the determination of unknown constants We assume that a branch of interface (that begins from the stagnation point) has a length of l chosen long enough so that it does not affect the solution. The interface length is kept unchanged during the simulation. Let the interface branch be discretized into N points S1(sZ0), S2,., SN(sZl) and q1, q2,., qN be values of the function q at these points (qIZq(SI)). The interface is approximated simply by a polygon formed by these N points. The interface is constructed from the angles q1, q2,., qN so that the distance between two consecutive points does not change during the simulation. Thus, the following algorithm is used for the interface reconstruction 0 1 8 > > q C q I IK 1 > AdsIK1 > xI Z xIK1 C cos@ > > 2 < 0 1 (55) for I O 1 > > q C q > I IK 1 > AdsIK1 > h Z hIK1 C sin@ > : I 2 where dsIK1 is the distance between SI and SIK1 and x1Z0, h1ZKe. In general, the solution of the Poisson equation (50) with the interface defined by S1, S2,., SN depend on q1, q2,., qN and on three constants Fr, Jcte and e. Solving Eq. (50) by a BEM give the values of the normal derivatives of j at the interface. Thus, the value of the functional = can be calculated at point SI and we can write: =ðs;q;Fr;Jcte ;eÞ Z GðSI ;q1 ;.;qN ;Fr;Jcte ;eÞ Z GI ðq1 ;.;qN ;Fr;Jcte ;eÞ

(56)

Thus, Eq. (51) at point SI becomes: Lðs;q;q 0 ;q 00 ;a;SÞ Z LðSI ;q1 ;.;qN ;a;SÞ Z LI ðq1 ;.;qN ;a;SÞ Z GI ðq1 ;.;qN ;Fr;Jcte ;eÞ

(57)

The two conditions at infinity (sZl) (37) and (40) can be written: C1ðq1 ;.;qN Þ Z 0

(58)

C2ðq1 ;.;qN ;Fr;Jcte ;eÞ Z 0

(59)

546

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

The forms of LI, C1 and C2 depend on the approximation formula used for the derivative terms in Eqs. (57)–(59). Any finite-difference formula of second degree or higher can be used for the approximation. Eq. (57) yields NK2 nonlinear equations for IZ2, 3,., NK1 that form together the conditions (58) and (59) a set of N equations for N unknown quantities q1, q2,., qN. However, as mentioned before, there are still unknown constants to be found. For the case (i), the condition (41) is used for finding the only unknown constant: the Froude number Fr or the bubble velocity. In this case, we can obtain NC1 equations for NC1 unknown quantities q1, q2,., qN and Fr. The resulting set of nonlinear equations is solved by Powell’s hybrid method. For the case (ii), a similar procedure can be applied for which we can obtain NC3 equations for NC3 unknown quantities q1, q2,., qN, Fr, Jcte and e. Here, three conditions at the stagnation point are used: they are given by Eqs. (42)–(44). Finally, the numerical procedure can be described as follows: Step 1. Suppose q1, q2,., qN, Fr, Jcte and e known in the kth iteration loop. Step 2. Solve the Poisson equation (50) in the domain determined by the interface constructed by q1, q2,., qN, e and the channel walls to obtain the values of GI(IZ2, 3,., NK1). Step 3. Solve the set of nonlinear equations to obtain q1, q2,., qN, Fr, Jcte and e for the (kC1)th iteration loop. The convergence is reached when k kC1 maxðjqkC1 KqkI j;jFr kC1 KFr k j;jJkC1 I cte KJcte j;je

Kek j; IZ 1;2;.;NÞ! 31

(60)

where k is the index of the iteration loop and 31/1 is a convergence tolerance. 3.3. Solution of a Poisson equation by BEM The Poisson equation (50) solved by a BEM with linear boundary elements. When the source term is a constant or a function of the domain coordinates only, the method can be found in classical references [40]. But when the source term depends also on the function and on its derivatives, a special iterative procedure should be applied. This procedure is as follows. The iteration loop begins with PZ0. Suppose that the boundary is divided into N elements with NeZ(Ne1CNe2) nodal points where Ne1, Ne2 are the number of nodal points on the boundary with the Neumann condition G1 and the Dirichlet condition G2, respectively (note that for linear elements NZNe). The solution with PZ0 gives Ne1 values q01 ;q02 ;.; q0Ne1 of vj/vn on G1 and Ne2 values j01 ;j02 ;.;j0Ne2 of j on G2. This allows calculating all the values of j in the domain, referred to as j0.

The first approximation for P is determined by       vj vj0 Z P0 C 6 P y;j0 ; KP0 P1 y;j; vy vy

(61)

where 0!6 %1 is a relaxation factor. Next, the Poisson equation is solved with PZP1 which allows to determine the values q11 ;q12 ;.;q1Ne1 of vj/vn on G1 and the values j11 ;j12 ;.;j1Ne2 of j on G2. Then, the second approximation for P, P2 can be determined. In general, the kth approximation for P is determined in an iterative way by:       vj vjkK1 Pk y;j; Z PkK1 C 6 P y;jkK1 ; (62) KPkK1 vy vy The procedure is repeated until the following condition is satisfied Kjki j;jqkC1 Kqkj j; maxfjjkC1 i j i Z 1;2;.; Ne2 ; j Z 1;2;.; Ne1 g! 32

(63)

where 32/1 is a convergence tolerance. 4. Sensibility analysis In this section, we first address the question of the multiple solutions. It is shown numerically that there exists a discrete set of different solutions for which the bubble velocity has an upper limit: this maximum velocity solution will be chosen as the observable solution. Next, the influence of the bubble length and the discretization size on the solution is investigated. These are numerically established so that they do not affect the solution. 4.1. Problem of multiple solutions Using the numerical procedure described in Section 3, the calculation starts with an initial bubble shape and guess values of Fr0 for case (i) or Fr0 and J0cte for case (ii). The initial bubble shape used in our simulation is determined from the formula:    h Ce 2 zðhÞ ZKA ln 1K K1% h% e; (64) d 0!d!1Ce for the case (i) and from 8 0 0 12 1 > > > h > > KA1 ln@1K@ A A for 0%h%1; 0!d1 !1; > > d1 < zðhÞ Z 0 0 12 1 > > > > >KA2 ln@1K@ h A A for K1%h%0; 0!d2 !1: > > : d2 (65) for case (ii). Some initial bubble shapes are presented in Fig. 3. It can be pointed out that the parameters d, d1, d2 decide the size of the

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

Fig. 3. Bubble shape with FrZ0.225 and some initial bubble shapes.

liquid film while A, A1, A2 determine the magnitude of the interface curvature at the stagnation point. The initial position of the stagnation point for the case (ii) is assumed to be at the coordinate origin (i.e. e0Z0), but it is adjusted during the calculation from the condition of vanishing velocity (42). For some values of A1, A2 or d1, d2 we can start from an asymmetrical initial bubble shape. The numerical program has been tested intensively for various initial bubble shapes and various guess values of Fr0 and J0cte with convergence tolerances 31Z10K6 and 32Z10K3. The results showed that for both cases (i) and (ii) we may get different solutions depending on the initial conditions. To search all the possible solutions, we applied the procedure similar to that of Vanden-Broe¨ck [21] and Coue¨t and Strumolo [22]. For case (i), the problem is solved with a given Froude number and the two conditions (37) and (40) at infinity, but without the condition (41) that gives the contact angle at the stagnation point. The numerical results show that in such conditions the solution does not depend on the initial bubble shape. Thus, for each given Fr we obtain one specific value of the contact angle at the stagnation point q1Zq(0). Therefore, the condition of the given contact angle (41) can be used for determining Fr. The numerical results show that the angle q1 is an oscillating function for Fr smaller than certain critical value Frc, then it decreases monotonically with FrOFrc. As a consequence, it may exist a countable (possibly infinite) number of solutions of Fr!Frc for which we have q1Zq0. And according to Garabedian’s maximal velocity criterion [19], the solution corresponding to the maximal Fr should be the physically observable solution. To search the possible solutions for case (ii), we solve the interface equation for two separated branches starting from the stagnation point with conditions (37) and (40) at infinity for sZCN or (36) and (39) for sZKN. The resulting solution possesses two different interface slope angles at the stagnation point. The conditions for the interface being smooth at the stagnation point requires the equality of the slope angle and of its first derivative (i.e. the curvature): conditions (43) and (44). These two conditions are used for determining Fr and Jcte. Similarly to case (i), for each given value of Fr the interface equation is solved for two branches. The two constants Jcte and e are determined simultaneously by means of the two conditions (42) and (44). The difference between the two angles of the interface at the stagnation point is a function of Fr: Dq(Fr)Zq(0C)Kq(0K). It can be shown that the function Dq(Fr) is an oscillating function around zero for Fr smaller than a certain critical value Frc, and then it decreases

547

monotonically with FrOFrc. Therefore, it may exist a countable (possibly infinite) number of solutions of Fr!Frc for which we have Dq(Fr)Z0. The good solution for Fr can be selected according to the Garabedian’s maximal velocity criterion [19]. To illustrate the considerations presented above, let us consider the bubble motion in vertical channel and stagnant liquid. For vertical flow, the problem can be considered as case (i) when the bubble is assumed symmetrical, and case (ii) when the symmetry is not required. In the first case, the stagnation point lies on the channel axis, then JcteZ0 and q(0)Zq0Zp/2. Upstream the conditions at infinity are placed at a distance of 4 (i.e. 4 times the channel haft-width) from the origin of the coordinate system. Downstream they are at a distance such that the interface length from the stagnation point is equal to 4. Because the interface curvature may be large near the bubble tip fine meshes are required there. The interface is thus discretized into non-uniform meshes, the mesh step ds changing from dsminZ0.025 near the bubble tip to dsmaxZ 0.1 far downstream. We have intensively tested various initial bubble shapes determined by (64) and (65) and various Fr0 for case (i) and Fr0, J0cte for case (ii). As a result, we obtained three different solutions approximately identical for both cases (i) and (ii) with FrZ0.225, 0.135 and 0.067 when SZ0.01. For very large ranges of d and A, for example, for all AO0, d!0.9 and Fr0O 0.10, the solution converges to FrZ0.225. The second solution with FrZ0.135 is found with AZ1.2, dZ0.99 and Fr0Z0.14 and the third solution with FrZ0.067, with AZ1.2, dZ0.99 and Fr0Z0.07. The bubble shapes corresponding to the three solutions are shown in Fig. 4. The comparison between the solution with FrZ0.225 and some of the initial bubble shapes are presented in Fig. 3. When the bubble is assumed symmetrical, the angle at the stagnation point q(0) is illustrated in Fig. 5 for Fr varying from 0.02 to 0.34 for SZ0.001, 0.01, 0.1. We can see that q(0) oscillates around its expected value of 908 for Fr!0.225 and then decreases monotonically for FrO0.225. It can be noted that, as S approaches zero, both the wavelength and the oscillation amplitude decrease and the density of solutions could cover the interval 0!Fr!0.225. For SZ0.01, we can find again the three solutions FrZ0.225, 0.135 and 0.067 for which q1Z908. According to the criterion of maximal velocity, the solution with FrZ0.225 is the physical solution for SZ0.01.

Fig. 4. Bubble shapes for three solutions in the case of vertical channel with SZ0.01.

548

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

Fig. 5. Variation of angle q1 as a function of Fr for vertical channel with SZ 0.001, 0.01 and 0.1.

Fig. 7. Influence of bubble length on the solution for vertical bubble.

The code was tested for different bubble lengths with lZ3, 4, 5 and different surface tensions with SZ0.001, 0.01 and 0.1 for a symmetrical bubble in a vertical channel. The discretization was the same for all cases with dsminZ0.025

and dsmaxZ0.1. In all of our simulations, we placed the upstream boundary at a distance luZl from the origin of coordinate system. The results are presented in Fig. 7. For this particular case, the bubble velocity Fr is not very sensitive to the truncation of the flow domain provided the downstream boundary is located farther than three (lZluO3). Therefore, it is sufficient to place the boundary conditions at infinity at a distance from the bubble nose greater than, say, three times the channel haft-width. The bubble length lZ4 is what was chosen for the vertical case. Similarly the results for the horizontal bubble are shown in Fig. 8. The simulations have been carried out for different bubble lengths with lZ3, 4, 5, 8 and different surface tensions, SZ0.001, 0.01 and 0.1. The discretization is the same for all cases with dsminZ0.025 and dsmaxZ0.2. There is a small difference between the solutions for lZ3 and 4, but for lR5 the solutions are not sensitive to the bubble length. We chose lZ5

Fig. 6. Variation of the function Dq(Fr) for SZ0.01 when the flow symmetry is not assumed.

Fig. 8. Influence of bubble length on the solution for horizontal bubble.

For the case when the flow symmetry is not assumed, for each Fr fixed we have found a symmetrical solution with JcteZ0 and eZ0 independent on the initial bubble shape. That suggests that only a symmetrical solution should exist in this case. The variation of the function Dq(Fr)C908 in comparison with q1(Fr) is shown in Fig. 6. We can see that the three solutions with FrZ0.225, 0.135 and 0.067 are also found. Again the symmetrical solution with FrZ0.225, JcteZ0 and eZ0 is the physical solution for this case. 4.2. Influence of bubble length

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

Fig. 9. Influence of discretization on the solution for horizontal bubble.

in our simulations for both the horizontal and the inclined geometries. 4.3. Influence of discretization The code is tested for different discretization sizes for dsminZ0.1, 0.05, 0.025, 0.0125 with lZluZ4 for the vertical case and lZluZ5 for the horizontal case. The results have shown that the solutions are more sensitive to the discretization for the small value of S (S!0.01). In the vertical case, the simulation did not converge for SZ0.001 with dsminZ0.1, while it converges well for a smaller discretization. In the horizontal case, it can be noted that there is a considerable variation of the bubble velocity for SZ0.001, it change from 0.4748 for dsminZ0.1 to 0.4898 for dsminZ0.0125. However, for the discretization smaller than 0.05, the difference between the solutions can be negligible. Therefore, the discretization with dsminZ0.025 can be used for the calculations. The variation of Froude number as function of discretization size is presented in Fig. 9 for horizontal case.

549

Fig. 10. Variation of Fr as a function of surface tension in the case of axisymmetrical bubble.

about 0.225 for the plane case and 0.345 for the axisymmetrical case. These values are in close agreement with the experimental and theoretical results of Zukoski [3], Collins [10,15], Griffith [9], Garabedian [19], Coue¨t and Strumolo [22] and Bendiksen [16], For SO0.1, the velocity decreases monotonically when the surface tension increases and the comparison with other results is quite good. The bubble shape (Fig. 12) obtained for SZ0.001 is in excellent agreement with the experimental results for both the plane and the axisymmetrical cases [42,43]. In addition, the values of the curvature radius at the bubble nose for SZ0.001 are in very good agreement with the numerical results of Vanden-Broe¨ck [20,21] and the experimental results of Collins [10] (0.64 for the plane case and 0.72 for the axisymmetrical case). It should be note that although we applied different

5. Results In this section, we present the simulation results for bubbles moving in both stagnant and flowing liquid. The influence of surface tension, inclination and velocity profile is investigated. The bubble velocity and shape are compared against experimental, theoretical and previous numerical results. 5.1. Bubble shape and velocity in stagnant liquid The results for an axisymmetrical and plane bubble in vertical channel with surface tension parameter changing from SZ0.001 to 0.7 are presented in Figs. 10 and 11 (aZ908). We can see that for S!0.1 the velocity remains almost constant

Fig. 11. The variation of Fr with surface tension for three inclination angles.

550

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

Fig. 12. The bubble shape with SZ0.001 (a) plane case; (b) axisymmetrical case.

asymmetrical initial bubble shapes from (65) the solution has always converged to the symmetrical one. When the channel is slightly deviated from the vertical, the bubble may be asymmetrical although the liquid still drains in two liquid films around the bubble. We look for the existence of such a pattern when the channel is deviated from the vertical. According to Maneri and Zuber [5] this bubble pattern may occur for inclination lying between 90 (i.e. the vertical) and 608. But for inclination between 80 and 608, the flow is characterized by the instability of the upper liquid film: they referred to it as a transition regime. With the algorithm described above for the case (ii), our code worked successfully for an inclination between 90 and 608. We observed that the velocity increases very slowly when the inclination changes from 90 to 808 (Fig. 14). This feature was also put in evidence by Zukoski [3] in his experimental works. At these inclinations the bubble velocity varies from 0.225 to 0.244 for SZ0.01. The calculated non-symmetrical bubble shapes with different stagnation point positions are presented in Fig. 13. In horizontal flow, the bubble leans against the upper channel wall. We look for the existence of such a pattern for an inclined channel whose inclination lies between aZ0 (i.e. horizontal) and 808: this includes the transition range defined by Maneri and Zuber [5]. The algorithm for case (i)

Fig. 13. Non-symmetrical bubble shapes in inclined channel.

is used to solve the problem. Like in the vertical case, the bubble velocity decreases as S increases for all inclinations. A plot of Fr versus S for inclination aZ0, 45 and 908 is given in Fig. 11. Incidentally, we found a mistake in the presentation of the results of Coue¨t and Strumolo [22]: if our results are in close agreement with theirs for vertical flow, a good agreement is obtained only when the value of S is doubled for inclined flow. After this correction our results agree well with that of Coue¨t and Strumolo [22] for small S. However, it remains for larger S a small discrepancy that might be due to the way the curvature of the interface is approximated. In Fig. 14, we summarized the dependence of the bubble velocity Fr on the inclination angle and we compared our numerical predictions to the experimental data of Maneri [4] and to the numerical results of Coue¨t and Strumolo [22] over the full range of inclination. The agreement is quite acceptable for a!608. This is true also near the vertical for aO808 even if some discrepancy is observed between experimental and numerical results: however, the 2D shape is not guaranteed in experiments. However, for the transition regime (608!a! 808) our results clearly point out the existence of two different solutions obtained for cases (i) and (ii) and both differs from their solution. 5.2. Bubble shape and velocity in flowing liquid Unlike in a stagnant liquid when we have to deal with irrotational flow and to solve Laplace equation, in flowing liquid we have to deal with rotational flow and to solve a nonlinear Poisson equation (50). For solving Poisson equation by a BEM, we have to take area integrals over the whole flow domain, so that it is necessary to discretize the flow domain. Because this domain changes during the calculation we have to use a special adaptive triangulation procedure. At each iteration step, the flow domain is divided into four regions: a region near the stagnation point, a region around the bubble

Fig. 14. Variation of Fr with inclination angle.

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

Fig. 15. Typical triangulation of the flow domain: (a) overview; (b) in detail near bubble nose.

tip, an upstream region and a liquid film region. In each region, we can use a simple triangulation process. Typical triangulation of the flow domain is illustrated in Fig. 15. The velocity of the bubble is expected to follow the Nicklin correlation (2). In general, the coefficient C0 depends on the velocity profile far ahead the bubble tip. It depends also on surface tension and inclination. We consider two cases of velocity profile: a profile of a laminar fully developed flow and a profile that mimics the mean velocity distribution of a fully developed turbulent flow. Although the liquid is assumed inviscid the viscosity has still a particular role in shaping the velocity profile. The liquid Reynolds number ReL characterizes the velocity profile imposed far ahead of the bubble. For large though moderated ReL, the flow far upstream is laminar and the velocity has a parabolic distribution. For very large ReL, the flow may be turbulent and the velocity distribution is a function of ReL. For axisymmetrical flow in tubes, in the cylindrical coordinate system the following correlation was used by Collins et al. [15] and Bendiksen [16]) uðhÞ Z uc ð1Khr 2 Kð1KgÞh2n Þ

(66)

where uc is the maximum velocity in the tube axis and g and n are two parameter that depend on the Reynolds number and that can be determined from the wall friction law. According to Bendiksen [16], by using the Prandtl friction law pffiffiffiffiffiffiffiffi 1=fw Z 3:5 log10 ðReL ÞK2:6

(67)

551

Fig. 16. Coefficient C0 for an axisymmetrical bubble moving in flowing liquid in vertical tube.



log10 ðReL ÞK0:743  g K 1K n Z ðgK1Þ log10 ðReL Þ C 0:31 2

K1 K1

(69)

Note that Eq. (66) includes the parabolic distribution by setting gZ1. The simulation results for C0 is presented in Fig. 16 for an axisymmetrical bubble moving in vertical tube for both laminar and turbulent velocity profile. The agreement between the results is quite good. 6. Conclusions The numerical method presented in this paper has been shown to reliably predict the velocity and shape of a long bubble rising in stagnant or flowing liquid of low viscosity. The investigation of multiple solutions of the problem shows that the supposition on the bubble shape at infinity (no oscillation of the interface) in fact is to limit the number of the solution to a discrete set. Moreover, the maximal velocity criterion proposed by Garabedian for selecting a unique physically observable bubble velocity was shown to be indispensable for the problem; otherwise, the stability analysis of solutions should be performed. One important advantage of the method is that the bubble shape can be obtained with precision and all boundary conditions on the interface should be satisfied exactly on the interface. The simulation results are in good agreement with both previous theoretical predictions and experimental observations.

where fw is the friction factor, the following values of g and n yield:

Acknowledgements

7:5 gZ 4:12 C 4:95ðlog10 ðReL ÞK0:743Þ

The work reported in this paper was supported by a grant from ToTal and partly by Vietnam Programme for Basic Research.

(68)

552

H. Ha-Ngoc, J. Fabre / Engineering Analysis with Boundary Elements 30 (2006) 539–552

References [1] Davies RM, Taylor G. The mechanics of large bubbles rising through extended liquids and through liquids in tubes. Proc R Soc Ser A 1950;200: 375–90. [2] Benjamin TB. Gravity currents and related phenomena. J Fluid Mech 1968;31:209–48. [3] Zukoski EE. Influence of viscosity, surface tension, and inclination angle on motion of long bubbles in closed tubes. J Fluid Mech 1968;25:821–37. [4] Maneri C. The motion of plane bubbles in inclined ducts PhD Thesis, Polytechnic Institute of Brooklyn; 1970. [5] Maneri C, Zuber N. An experimental study of plane bubbles rising at inclination. Int J Multiphase Flow 1974;1:623–45. [6] Harmathy TZ. Velocity of large drops and bubbles in media of infinite or restricted extent. AIChE J 1960;6:281–8. [7] Nicklin J, Wilkes JO, Davidson JF. Two phase flow in vertical tubes. Trans Inst Chem Eng 1962;40:61–8. [8] White ET, Beardmore RH. The velocity of rise of single cylindrical air bubbles through liquids contained in vertical tubes. Chem Eng Sci 1962; 17:351–61. [9] Griffith P. The prediction of low-quality boiling voids; 1963 [ASME Paper No. 63-HT-20]. [10] Collins R. A simple model of the plane gas bubble in a finite liquid. J Fluid Mech 1965;22:763–71. [11] Bendiksen KH. An experimental investigation of the motion of long bubbles in inclined tubes. Int J Multiphase Flow 1984;10:467–83. [12] Hasan R, Kabir CS. Predicting multiphase flow behavior in a deviated well SPE 15449, 61st annual technical conference, New Orleans. LA 5–8 October; 1986. [13] Wallis GB. One dimensional two-phase flow. New York: MacGraw-Hill; 1969. [14] Dumitrescu DT. Stro¨mung an einer Luftblase im Senkrechten Rohr. Z Angew Math Mech 1943;23:139. [15] Collins R, De FF, Davidson JF, Harrison D. The motion of large gas bubble rising through liquid flowing in a tube. J Fluid Mech 1978;89: 497–514. [16] Bendiksen KH. On the motion of long bubbles in vertical tubes. Int J Multiphase Flow 1985;11:797–812. [17] Nickens HV, Yannitell DW. The effects of surface tension and viscosity on the rise velocity of a large gas bubble in close, vertical liquid-filled tube. Int J Multiphase Flow 1987;13:57–69. [18] Birkhoff G, Carter D. Rising plane bubbles. J Math Mech 1957;6:769–79. [19] Garabedian PR. On steady-state bubbles generated by Taylor instability. Proc R Soc A 1957;241:423–31. [20] Vanden-Broe¨ck JM. Bubbles rising in a tube and jets falling from a nozzle. Phys Fluids 1984;27:1090–3. [21] Vanden-Broe¨ck JM. Rising bubble in two-dimensional tube with surface tension. Phys Fluids 1984;27:2604–7. [22] Coue¨t B, Strumolo GS. The effects of surface tension and tube inclination on a two-dimensional rising bubble. J Fluid Mech 1987;184:1–14. [23] Ha-Ngoc H. Etude the´orique et nume´rique du mouvement de poches de gaz en canal et en tube. Thesis, Institut National Polytechnique de Toulouse, France; 2003.

[24] Mao Z-S, Dukler AE. The motion of Taylor bubbles in vertical tubes I. A numerical simulation for the shape and rise velocity of Taylor bubbles in stagnant and flowing liquids. J Comput Phys 1990;91: 132–60. [25] Mao Z-S, Dukler AE. The motion of Taylor bubbles in vertical tubes II. Experimental data and simulations for laminar and turbulent flow. Chem Eng Sci 1991;46:2055–64. [26] Bugg JD, Mack K, Rezkallah KS. A numerical model of Taylor bubbles rising though stagnant liquids in vertical tubes. Int J Multiphase Flow 1998;24:271–81. [27] Benkenida A. De´veloppement et validation d’une me´thode de simulation d’e´coulement diphasique sans reconstruction d’interface Application a` la dynamique des bulles de Taylor. Toulouse, France: The`se de doctorat, Institut National Polytechnique; 1999. [28] Blake JR, Taib BB, Doherty G. Transient cavities near boundaries. Part 1. Rigid boundaries. J Fluid Mech 1986;170:479–97. [29] Chakrabarti R, Harris PJ, Verma A. On the interaction of an explosion bubble with a fixed rigid structure. Int J Numer Methods Fluids 1999;29: 389–96. [30] Harris PJ. An investigation into the use of the boundary integral method to model the motion of a single gas or vapour bubble in a liquid. Eng Anal Bound Elem 2004;28:325–32. [31] Heister SD. Boundary element methods for two-fluid free surface flows. Eng Anal Bound Elem 1997;19:309–17. [32] Kamiya N, Nakayama K. Prediction of free surface of die swell using the boundary element method. Comput Struct 1993;46:387–95. [33] Yoon SS, Heister SD. A fully non-linear model for atomization of highspeed jets. Eng Anal Bound Elem 2004;28:345–57. [34] Lundgren TT, Mansour NN. Oscillations of drops in zero gravity with weak viscous effects. J Fluid Mech 1988;194:479–510. [35] Murray IF, Heister SD. On a droplet’s response to acoustic excitation. Int J Multiphase Flow; 25, 531–50. [36] Lamb H. Hydrodynamics. Cambridge: Cambridge University Press; 1932. [37] Batchelor GK. An introduction to fluid dynamics. Cambridge: Cambridge University Press; 1967. [38] Hawthorne WR. In: Sovran G, editor. Fluid mechanics of internal flow. Amsterdam: Elsevier; 1967. [39] Vied IA, Bajalieva SS. A boundary problem for a system of second degree integral-differential equations of Volter type. Investigations on integral-differential equations, vol. 13. Frunze: Ilim Press; 1980 [in Russian]. [40] Brebbia CA. The Boundary element method for engineers. London: Pentech Press; 1978. [41] Powell MJ-D. A hybrid method for nonlinear algebraic equations. In: Rabinowitz P, editor. Numerical methods for nonlinear algebraic Equation. London: Gordon and Breach; 1970. [42] Brown RAS. The mechanics of large gas bubbles in tube. I. Bubble velocities in stagnant liquid. Can J Chem Eng 1965;43:217–23. [43] DeJesus JD, Ahmad W, Kawaji M. Experimental study of flow structure in vertical slug flow Proceedings of second int. conf. multiphase flow, Kyoto, 3–7 April; 1995.