Numerical Simulation of Bubble Cluster Induced Flow by Three-Dimensional Vortex-in-Cell Method.

Bin Chen1 State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China e-mail: [email protected]

Zhiwei Wang State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China; Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Xueyuan Avenue 1068, Shenzhen 518035, China

Tomomi Uchiyama EcoTopia Science Institute, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan

Numerical Simulation of Bubble Cluster Induced Flow by Three-Dimensional Vortex-in-Cell Method The behavior of air bubble clusters rising in water and the induced flow field are numerically studied using a three-dimensional two-way coupling algorithm based on a vortexin-cell (VIC) method. In this method, vortex elements are convected in the Lagrangian frame and the liquid velocity field is solved from the Poisson equation of potential on the Eulerian grid. Two-way coupling is implemented by introducing a vorticity source term induced by the gradient of void fraction. Present simulation results are favorably compared with the measured results of bubble plume, which verifies the validity of the proposed VIC method. The rising of a single bubble cluster as well as two tandem bubble clusters are simulated. The mechanism of the aggregation effect in the rising process of bubble cluster is revealed and the transient processes of the generation, rising, strengthening, and separation of a vortex ring structure with bubble clusters are illustrated and analyzed in detail. Due to the aggregation, the average rising velocity increases with void fraction and is larger than the terminal rising velocity of single bubble. For the two tandem bubble cluster cases, the aggregation effect is stronger for smaller initial cluster distance, and both the strength of the induced vortex structure and the average bubble rising velocity are larger. For the 20 mm cluster distance case, the peak velocity of the lower cluster is about 2.7 times that of the terminal velocity of the single bubble and the peak average velocity of two clusters is about 2 times larger. While for the 30 mm cluster distance case, both the peak velocity of the lower cluster and two clusters are about 1.7 times that of the terminal velocity of the single bubble. [DOI: 10.1115/1.4026968] Keywords: vortex-in-cell method, two-way coupling, bubble cluster, aggregation effect, vortical flow, gas-liquid two-phase flow

1

Introduction

Gas-liquid two-phase flow exists widely in various engineering processes. The behaviors of bubbles in gas-liquid two-phase flow are often surprising and unexpected because they suffer various forces [1]. Bubbles rising in water will induce ambient liquid flow and the dynamic changed bubble distribution induces a different flow field. At the same time, the induced flow field increases the inhomogeneity of the bubble behavior and distribution. The interaction between bubbles and ambient flow field often lead to a turbulent flow and make the flow highly complicated. These interactions have attracted many researchers’ attention not only in scientific research but also in engineering application, while a comprehensive knowledge of the interaction mechanisms between the two phases is still lacking [2]. The case of bubble clusters rising in quiescent water due to buoyancy is used in the present work as a reasonable option to investigate the interaction between phases in gas-liquid two-phase flow since in this way bubble behavior and bubble-fluid interactions can be studied without considering the effect of local flow conditions. Bubbles in a cluster released into quiescent water will rise due to buoyancy, inducing an ambient liquid flow, which in turn will also affect the behavior of the bubbles in the cluster. With the development of computational fluid dynamics techniques, numerical simulation of this process is expected to provide valuable information on the interaction mechanisms between the 1 Corresponding author. Contributed by the Fluids Engineering Division of ASME for publication in the JOURNAL OF FLUIDS ENGINEERING. Manuscript received July 8, 2013; final manuscript received February 24, 2014; published online May 12, 2014. Assoc. Editor: Zhongquan Charlie Zheng.

Journal of Fluids Engineering

two phases [3]. The vortex method is one of the effective numerical methods to reveal the interaction mechanisms [4]. In the vortex method, the flow domains where vorticity is nonzero are initially discretized into a number of vortex elements, whose motion can be obtained by summing up the contribution of all other elements using the Biot–Savart law. Subsequently, the flow field can be obtained through the induced velocity of all the vortex elements. The evolution of the flow field is determined by the motion of vortex elements, so the instantaneous large-scale coherent structures can be readily identified. The major advantages of the traditional discrete vortex method are the automatic adaptability of the computational elements and low numerical dissipation. In the past several decades, vortex methods have been applied to simulate various typical fields such as density stratified fluids [5], jet [6], combustion [7], flow around cylinders [8–10], vortex rings [11,12] and plume [13,14]. However, when using the Biot–Savart equation, the computational time will increase drastically as the number of vortex elements N is increasing since the required operations are of order O(N2). Therefore, the method is rendered impractical for highresolution three-dimensional simulations because of the tremendous computational time requirements. Different from the pure Lagrangian traditional vortex methods, the vortex-in-cell (VIC) method is a mixed Eulerian–Lagrangian algorithm, which combines the best features of the Eulerian and Lagrangian methods [15]. In VIC methods the vortex elements are tracked in the Lagrangian frame while the velocity field is calculated on Eulerian grids after distribute the vorticity carried by the vortex elements to the grid nodes. The computational cost of VIC methods increases linearly with the number of vortex elements. When the number of the vortex elements is N and the cell number is M, the

C 2014 by ASME Copyright V

AUGUST 2014, Vol. 136 / 081301-1

Downloaded From: http://fluidsengineering.asmedigitalcollection.asme.org/ on 05/20/2014 Terms of Use: http://asme.org/terms

computational complexity can be reduced to O(N þ Mlog2M) [16–18]. Furthermore, there are many efficient solvers available for the Poisson equation on a regular Eulerian mesh. Therefore, VIC methods are preferred because they are conceptually more efficient and have been proved to be hundreds or thousands times faster than direct integration from the Biot–Savart equation. Christiansen [19] was the first to apply the VIC method to a fluid flow problem. Cheng et al. [20] simulate the twodimensional viscous incompressible flows over a bluff body using a hybrid vortex method, in which the diffusion-vortex method is used in the region near the body surface to solve the Navier–Stokes equations, while the vortex-in-cell method is used in the exterior domain. Abdolhosseini and Milane [21] proposed a two-dimensional VIC method for the simulation of a spatially growing high Reynolds number mixing layer. In the method the vorticity from the point vortex is distributed to the four neighboring nodes surrounding each vortex using an area weighting scheme and the velocity on the nodes is calculated using a central difference scheme. Two fractional steps were used to solve the convection and viscous diffusion, respectively. Liu [22] proposed a vortex particle-in-cell method for the simulation of threedimensional unsteady incompressible viscous flow. In the method the projection of the vortex strengths onto the mesh is based on volume interpolation and the change in vorticity due to diffusion is also computed on the Eulerian mesh and projected back to the particles. Milane and Abdolhosseini [15] developed a vortex filament method based on the VIC method for the investigation of a spatially growing uniformly sheared flow. In the method the Lagrangian vortex filaments are initially subjected to a threedimensional random perturbation to simulate an initial fluctuating velocity field with nearly isotropic conditions. Uchiyama and Naruse [23] proposed a two-dimensional VIC method for the simulation of gas-particle two-phase free turbulent flow. In the method the velocity vector field was represented as the summation of the gradient of a scalar potential and the rotation of a vector potential, both of which are solved on the Eulerian grid. The Lagrangian vortex methods were originally conceived for the simulation of inviscid unsteady flows. Chorin [24] was the first to study slightly viscous flow using a random walk vortex method. Then Leonard [16] indicates that the random walk vortex method requires a large number of particles compared to the Reynolds number to achieve accurate results of viscous diffusion. Later several approaches for viscous flows have been proposed to simulate the viscous diffusion [25]. The approaches can be classified into two types; in one type, vortex elements are transported while their strength remains fixed, such as the random walk [24–26] and diffusion velocity methods [27,28]. While in another type, the strength of vortex elements is changed, such as the core-spreading method [29] and particle strength exchange (PSE) method [30,31]. The PSE method is a widely used algorithm in which the diffusion equation is converted into integral-differential equations and then discretized in space by approximating the integral using a quadrature rule. This method has been successfully applied to several typical flows [32–36]. By introducing bubble/particle motion equations, vortex methods are also applicable to two-phase flow analysis. Crowe et al. [37], who were the first to work on two-phase flows by vortex methods, studied the dispersion of particles in turbulent flows. Their methodology assumed that there were no collisions between particles and only one-way coupling was considered. Brecht and Ferrante [38] studied the evolution of one and two inviscid buoyant bubbles by VIC method. Chen and Marshall [39] proposed a two-dimensional vortex method to simulate the two-phase interaction in a two-phase flow with heavy particles. Two-way coupling was used in their method by introducing a particle-induced vorticity source term. Walther and Koumoutsakos [40] presented a three-dimensional vortex method for the simulation of particleladen flows. In their simulation the flow field was computed using Lagrangian vortex elements advected with the local velocity and viscous diffusion was considered by the particle strength 081301-2 / Vol. 136, AUGUST 2014

exchange algorithm. Yang et al. [41] studied the two-way interaction between small bubbles and the large-scale coherent structures in plane turbulent mixing layers using a discrete vortex model. They investigated the effect of large-scale vortex structures on the dispersion and condensation of bubbles and found that larger-size bubbles modify the vortex strength significantly [42]. Nevertheless, the two-way coupling between the two phases in their method was incomplete because the effect of void fraction was not considered [43]. Therefore, this method is only valid in flow with very low void fraction (typically about 1%) and the development of a two-way coupling VIC method for gas-liquid flow with larger void fraction is of particular scientific interest. The authors have also developed three-dimensional VIC methods for the simulation of bubble plume [13] and the particle-laden gas flow [44], as well as the interaction between vortex ring and bubble plume [14]. In their simulation the mass conservation equation was independently derived and a smoothing function was used to calculate the bubble volume fraction. Hence, the effect of void fraction was fully considered and they claimed the void fraction can be up to 18% in the simulation. Overall, significant developments have been achieved in the vortex methods simulation of two-phase flow in recent years while the studies are still far from maturity. There are still many issues unsolved, for example, how to introduce the vorticity source term to fully consider the two-way coupling between both phases, how to accelerate the computation for the two-phase flow with high void fraction, etc. The motivation of this paper is to develop a two-way coupling algorithm based on the VIC method, in which the interactions between the two phases can be fully considered with high computational efficiency, for the simulation of the rising of bubble clusters in quiescent water. The remainder of this paper is organized as follows: In Sec. 2, a two-way coupling numerical algorithm based on the VIC method for the problem is introduced. The numerical model is validated and the simulation conditions are described in Sec. 3. The simulation results are presented in Sec. 4. The behavior of a single bubble cluster as well as two tandem bubble clusters and the flow field induced by the rising bubbles are presented and the effect of void fraction and the distance between two bubble clusters on the induced flow are analyzed. The conclusion is given in Sec. 5.

2

Numerical Algorithm

2.1 Conservation Equations for Two-Phase Mixtures. The following assumptions are employed in the simulation: (1) Both phases are incompressible and no phase change occurs. (2) The mass and momentum of the gas phase are very small and negligible compared with those of the liquid phase. (3) The bubbles maintain their spherical shape. (4) Interactions among individual bubbles are neglected and neither fragmentation nor coalescence occurs. Based on these assumptions, the conservation equations for the mass and momentum of the bubbly flow are expressed as @al þ r ðal ul Þ ¼ 0 dt Dul 1 ¼ rp þ tl r2 ul þ al g al ql Dt al þ ag ¼ 1

(1) (2) (3)

where al and ag are liquid and gas fraction, respectively, ul is the velocity of the liquid phase, ql is the density of liquid phase, p is the pressure, tl is the kinematic viscosity of liquid phase, and g is the acceleration of gravity. According to Helmholtz’s theorem, any vector field can be represented as the summation of an Transactions of the ASME


the direction of the vorticity strength vector. With the discretization of vorticity field, the vorticity equation (8) can be solved in Lagrangian frame. In the VIC method, the vorticity of the liquid phase x is recorded on a regular Eulerian mesh, which is given by the summation of the contribution induced by all vortex elements xðxÞ ¼

N X cn x xn f r3 r n¼1

f ðeÞ ¼

Fig. 1 Interpolation of velocity from the neighboring mesh nodes

irrotational (curl-free) vector field and a solenoidal (divergencefree) vector field. Thus, the liquid velocity is written as ul ¼ uP þ uR

(4)

where uP ¼ r/ is the irrotational velocity vector, and uR ¼ r w is the solenoidal velocity vector; / and w are scalar and vector potential, respectively. So we can obtain ul ¼ r/ þ r w

(5)

Substituting Eq. (5) into Eq. (1) the following equation is obtained: @al þ r ðal r/Þ þ r w ral ¼ 0 @t

(7)

where x is the vorticity of the liquid phase Taking the curl of Eq. (2), the vorticity equation for the bubbly flow is derived. For the three dimensional calculation, it is expressed as @x tl 1 Dul þ r ðul xÞ ¼ x ðrul Þ þ r2 x þ ral g @t al al Dt (8) where the second term of the LHS is vortex convection and the three terms of the RHS are vortex stretching, viscous diffusion, and the vorticity source term, respectively. 2.2 Implementation of Vortex-in-Cell Method. In order to solve the vorticity equation (8), the vorticity field will be discretized into spherical blobs and the discretization of the vortex element proposed for single-phase flow simulation is employed here. Initially, there are no vortex elements in the flow field and new vortex elements will be generated due to the rising bubbles. Each vortex element has a viscous core and its own strength. As shown in Fig. 1, each spherical blob is shown as a short small cylinder to illustrate its position and vorticity strength simultaneously. The position of the spherical blob is the geometric center of the corresponding small cylinder, and the obliquity of the cylinder denotes Journal of Fluids Engineering

þ 1Þ7=2

(10)

where N is the number of vortex elements, c and r are the strength and core radius of each vortex element, respectively, and f(e) is the high order smoothing core function proposed by Winckelmans and Leonard [29]. On the Eulerian mesh, the vector potential w of the liquid phase will be determined from the Poisson equation (7) by the successive overrelaxation method. After solving the scalar potential / from Eq. (6), taking into account bubble void fraction, the velocity field of the liquid phase can be obtained by Eq. (5). 2.3 The Motion and Redistribution of Vortex Elements. After discretization, the locations of vortex elements are modified to account for convection, based on their velocity obtained by interpolation of the liquid velocity on the eight neighboring nodes dxn ¼ ul dt

(11)

while the strength of vortex elements are changed to account for vortex stretching [40] dcn ¼ rul cn dt

(12)

(6)

Taking the curl of Eq. (5) and considering r w ¼ 0, the vector Poisson equation is derived r2 w ¼ x

15 8pðe2

(9)

The motion of the vortex elements is governed by the local liquid velocity. The computational cell, where a vortex element is located, is divided into eight child cubes by the position of the vortex element. The local velocity of the vortex element can be interpolated from the velocity at the eight nodes of the cell. As can be seen in Fig. 1, the weight of each node depends on the volume defined by the corresponding diagonal. The second term of the RHS of Eq. (8) is calculated in the Lagrangian frame. For the calculation of the viscous diffusion term, a commonly used model in vortex methods is the so-called particle strength exchange scheme introduced by Degond and Mas-Gallic [30], in which the Laplacian operator is replaced by an equivalent integral one. In the present work the strength of the vortex elements is updated using the algorithm proposed by Winckelmans and Leonard [31] to account for viscous diffusion effect, expressed as Nv dcn 2tl X ðVn cm Vm cn Þgr ðxn xm Þ ¼ 2 dt r m¼1

(13)

where V is the volume of vortex element and gr is a regularization function, which is expressed as 1 j xj (14) gr ðxÞ ¼ 3 g r r gðeÞ ¼

1d f ðeÞ e de

(15)

As to the last term of the RHS of Eq. (8), the vorticity source term generated by the gradient of void fraction is considered as AUGUST 2014, Vol. 136 / 081301-3


the results using the two interpolation functions, which may indicate the use of the function in the redistribution process is accurate enough to obtain exact results. An example in one dimension is shown in Fig. 2. Initially there is a vortex element (the solid circle) and its strength is redistributed to neighboring nodes of grid cells according to the M04 function. Then the vortex element is replaced by a new set of vortex elements released from the nodes of the grid (the hollow circles). In three dimensions a tensorial product of the M04 interpolation function is used as smooth function. So the circulation in node (xi, yj, zk) is the summation of the contribution of every vortex element and it can be expressed as

Fig. 2 Redistribution of vortex elements using the M04 function

the change of circulation C in a computational cell when employing the Reynolds transport theorem ð DC @x ¼ þ r ðul xÞ dV (16) Dt V @t Substituting the right hand side of Eq. (8) into Eq. (16), the change of circulation in a cell can be expressed as ð DC 1 Dul ¼ ral ðg Þ dV (17) Dt Dt V al We can assume the change of circulation in a computational cell is DC in a time step Dt. In the case that there is no vortex element in the cell, a new vortex element will be generated at the center of the cell and the circulation of the new generated vortex element is the same as the change of the circulation of the cell. While in the case that the number of vortex elements in a cell is n, the change of the circulation for each vortex element during Dt is supposed to be added on every vortex element in the cell equally and expressed as DC/n. The strength of the vortex element c is equal to the product of its circulation and a length vector [45] cn ¼ xn Vn ¼ Cn ln

(18)

In this work both one-way coupling and two-way coupling algorithms are used to simulate the interaction between liquid flow and gas bubbles. In the one-way coupling algorithm, the bubble behavior is governed by the bubble motion equation and the bubble-induced vorticity source term is not used. So there is no effect of the gas bubbles to the liquid-phase flow. While in the two-way coupling algorithm, besides the influence of the liquid flow field on bubble motion, the bubble-induced vorticity source term is also used and new vortex elements are generated to account for the effect of the gas phase on the liquid phase. In the VIC method, the flow field is represented by individual vortex elements. With time elapse, there will be a distortion on the distribution of vortex elements induced by the flow map. The movement of vortex elements is not uniform, resulting in a distorted local vorticity field. The redistribution of vortex elements can be used as an optimization technique to improve the accuracy of the simulation. In the redistribution process the interpolation function is very important to the accurate and efficient assignment of vortex element and the M04 interpolation function proposed by Monaghan [46] is used in this work. This moment-conserving kernel function has been proved to introduce minimal numerical dissipation in simulations. A higher order interpolation function, the M6000 interpolation function introduced by Bergdorf and Koumoutsakos [47], was also used as a validation. There are little differences between 081301-4 / Vol. 136, AUGUST 2014

Ci;j;k

x x y y z z j n i n k n ¼ cn Wc Wc Wc Dx Dy Dz n¼1 Nv X

Wc ðeÞ ¼

8 > 1:5jej3 2:5jej2 þ1 > >
> > :

(19)

0:5ð2 jejÞ ð1 jejÞ 1 jej 2

(20)

jej > 2

0

where (xn, yn, zn) is the position of the nth vortex element, Ci,j,k is the circulation in node (xi, yj, zk), and Wc ðeÞ is the M04 smoothing interpolation function. However, there is still the dilemma how often to apply the redistribution [35]. On the one hand, it is better to redistribute frequently to best preserve overlap between neighboring vortex elements. On the other hand, one does not want to apply it too frequently to keep the redistribution-induced dissipation [48] much smaller than the viscosity-induced one since higher order redistribution introduces a hyperviscous type of dissipation. In this work, vortex elements were redistributed every ten time steps. 2.4 Simulation of Bubble Motion. The motion of bubbles is expressed by the following equation [43]: dug 1 þ CV Dul 1 3CD ¼ ur jur j dt b þ CV Dt b þ CV 4d b1 CL g ur ðr ul Þ þ b þ CV b þ CV

(21)

where ug is the bubble velocity and ur is the relative velocity (ur ¼ ug – ul); d is the bubble diameter; ql and qg are the liquid and bubble density, respectively; b is the density ratio; and b ¼ qg/ql. CV, CL, and CD are the added mass coefficient, lift coefficient, and drag coefficient, respectively. Both CV and CL are 0.5, and CD is given by the following equation [49]: 8 24 3 > > ð1 þ ReÞ Re 0:01 > > Re 16 > > < 24 CD ¼ ð1 þ 0:1315Reð0:820:05 log10 ReÞ Þ 0:01 < Re 20 > Re > > > > > : 24 ð1 þ 0:1935Re0:6305 Þ 20 < Re 260 Re (22) where Re ¼ db jur j=tl is the local Reynolds number based on velocity difference between the two phases. The calculation of void fraction is very important in two-phase flow, while a notable disadvantage of the VIC method in comparison to the direct integration method is the numerical diffusion caused by the grid. The continuity of void fraction is helpful to improve the accuracy of simulation. In this paper, the following expression was used in the calculation of void fraction to counteract the grid effects: Transactions of the ASME


Table 1 Comparison of computational complexity for VIC and Biot–Savart law Number of operations Number of vortex elements 1 104 1 104 1 104 5 104 5 104 5 104 1 105 1 105 3 105 3 105

ai;j;k

Cell number

Biot–Savart law

Vortex-in-cell

Acceleration rate

40 40 80 60 60 120 80 80 160 40 40 80 60 60 120 80 80 160 60 60 120 80 80 160 60 60 120 80 80 160

1.0 108 1.0 108 1.0 108 2.5 109 2.5 109 2.5 109 1.0 1010 1.0 1010 9.0 1010 9.0 1010

6.6 105 2.4 106 6.2 106 7.0 105 2.5 106 6.2 106 2.5 106 6.3 106 2.7 106 6.5 106

151 41 16 3559 1006 403 3953 1600 32,895 13,947

x x y y z z pdb3 j n i n k n ¼ Wa Wa Wa 6DxDyDz Dx Dy Dz n¼1 Nb X

(23) 8 2 2 > > < 0:5ð1:5 þ jejÞ 1:5ð0:5 þ jejÞ jej < 0:5 Wa ðeÞ ¼ 0:5ð1:5 jejÞ2 0:5 jej 1:5 > > : 0 jej > 1:5 (24) where ai,j,k is the local void fraction in grid cell (i,j,k) and Wa ðeÞ is the M3 smoothing interpolation function. 2.5 Summary of the Two-Way Coupling VIC Method. Now we can summarize the numerical procedure in our method. When the flow field ul of the liquid phase at time t is known, the velocity of the liquid phase at next time step t þ Dt is computed by the following procedure: (1) Calculate bubble motion from Eq. (21) to get the new bubble position. (2) Calculate the void fraction from Eq. (23) and obtain the liquid volume fraction from Eq. (3). (3) Calculate the convection of the vortex elements from Eq. (11) to update their positions. (4) Calculate the stretching and the effect of viscous diffusion on the strength of the vortex elements from Eq. (12) and Eq. (13), respectively. (5) Calculate the change of circulation in a cell from Eq. (17) to consider the effect of vorticity induced by the change of void fraction and change the strength of vortex elements in the cell. (6) Redistribute vortex elements every tent time steps. (7) Calculate the vorticity field of the liquid phase on the nodes of the Eulerian grid by summing the contribution of all vortex elements from Eq. (9). (8) Calculate vector potential w from Eq. (7). (9) Calculate scalar potential / from Eq. (6). (10) Obtain the new liquid velocity field from Eq. (5). In this method the effect of void fraction is considered in Eq. (1) and Eq. (2) while the method proposed by Yang et al. [41] solves the equations with al ¼ 1 and ag ¼ 0. Consequently, our method is able to take into account the volumetric fraction of each phase and, thus, can handle flows with larger void fractions. The limitation of the present method lies in the assumptions mentioned in Sec. 2.1. Since we assume that bubbles maintain their spherical shape, it is only valid in the cases of bubbles with small diameter. It is also assumed that neither break nor coalescence occurs and the interaction among bubbles is neglected, yielding limited density of the bubble population and, thus, the local void fraction. Journal of Fluids Engineering

In the traditional vortex method, one needs to sum up the contribution of all other vortex elements by the Biot–Savart law to obtain the induced velocity in each time step. The computational complexity is O(N2), where N is the number of the vortex elements in the flow field. However, in the present VIC method, the liquid velocity can be obtained by solving the vorticity Poisson equation (6) and equation (7) by the successive overrelaxation method on an Eulerian grid. When the number of the vortex elements is N and the cell number is M, the computational complexity can be reduced to O(N þ Mlog2M) [16–18]. As shown in Table 1, the computation complexity for the VIC method is several orders of magnitude smaller than that for direct integration by Biot–Savart law. The acceleration rate increases rapidly with vortex element number while decreasing with the increase of cell number. In typical cases in the present simulation, there are 3 105 vortex elements in the flow field and the cell number is 80 80 160; hence, the calculation using the present VIC method is about 10,000 times faster than that using direct integration by the Biot–Savart law.

3

Simulation Conditions

In this paper, the behavior of clusters of air bubbles rising in a water tank is studied by numerical simulation. Clusters of thousands of small gas bubbles are released free in a water tank. The bubbles rise due to buoyancy, setting in motion the surrounding liquid. Then the motion and distribution of bubbles are also affected by the flow field originally induced by themselves. In this simulation the size of the water tank is 200 mm 200 mm 400 mm. Four thousand bubbles randomly distributed in a spherical bubble cluster are released free from the bottom of the tank. As the bubble clusters are rising, the average bubble position zb will increase with time. At a first stage, the numerical simulation was performed with different grids to establish the independence of the solution to the grid resolution. Figure 3 shows the evolution of average bubble position for different grid densities. For the 20 20 40 grid, the average bubble position zb increases almost linearly with time, indicating that the resolution of the grid is too small. For larger grid densities of 60 60 120 and 80 80 160, the two curves on the figure almost coincide with each other, which indicates grid-independent solutions. Therefore, the water tank is resolved using 80 80 160 cubic cells in the present simulation, with cubic cells of edge size 2.5 mm in all directions. The core radius of the vortex elements is the same as the cell width, ensuring the overlap among neighboring vortex elements [16]. In the computation, the time t is normalized by the diameter of the cluster D and the terminal rising velocity of a single bubble UT and the normalized time t* can be expressed as t ¼

tUT D

(25)

AUGUST 2014, Vol. 136 / 081301-5


Fig. 5 Comparison between present simulation and experimental observation of bubble plume (i: experimental observation of Alam and Arakeri [51]; ii: present simulation) Fig. 3 Average position of bubbles in z direction at different grid resolutions

Fig. 4

Terminal rising velocity of single bubble

Fig. 6 Variation of plume centerline velocity on gas flow rate at 50 mm above plume source

Since there is an obvious difference between the velocity of bubble and liquid-phase flow, a subcycle is introduced to adjust the time step [8]. The liquid-phase field is computed only once at each time step, while bubble motion is determined at several substeps. A satisfactory precision can be obtained in the solution as long as the subcycle time step is small enough, and the computational efficiency is improved at the same time. The time step Dt* is 0.01 for water flow and 0.0002 for bubble motion in all cases. The terminal velocity of a single bubble rising in water is very important and there is vast literature on this subject. Zheng and Yapa [50] studied the terminal velocity of rigid particles with spherical shapes. They found that for very low Re number the rising velocity of a single bubble is described by the well-known Stokes’ law. For larger Re numbers (750 < Re < 3.5 105), the drag coefficient can be assumed to have a constant value of 0.5 and this approximation results in another expression of terminal velocity. Distinguishing these ranges with the use of a critical diameter the terminal velocity may be expressed as 8 2 gd ðql qg Þ > > > b > < 18tl ql UT ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 8gdb ðql qg Þ > > > : 3ql 081301-6 / Vol. 136, AUGUST 2014

db < dc (26) db dc

where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2l ql 3 dc ¼ 9:52 gðql qg Þ Figure 4 shows the comparison of the terminal velocity of a single bubble with different diameters between the present simulation and Zheng and Yapa’s correlation [50]. Both one-way coupling and two-way coupling options have been used in the computation. We can see that the simulation results of one-way coupling and two-way coupling are similar for small bubbles, while the result of two-way coupling is much more close to the analytical result of Zheng and Yapa for larger bubbles. The difference between the present simulation result and the correlation is acceptable in all cases, indicating that the present two-way coupling model can be used to simulate the gas-liquid two-phase flow satisfactorily. The motivation of the present paper is to investigate the rising of bubble clusters in water by the two-way coupling algorithm described in Sec. 2. However, since it is quite difficult to investigate the behavior of one or two bubble clusters, a continuous bubble plume was simulated by the two-way coupling model proposed in the manuscript and compared with experimental data. Transactions of the ASME


Fig. 7

Initial distribution of the bubble clusters

Alam and Arakeri [51] experimentally investigated the plane laminar bubble plume using flow visualization. Figure 5 shows the comparison of the plume structure at two different gas flow rates between the present simulation and experimental observation of Alam and Arakeri [51]. In the investigation they found that initially the bubble core remains compact and appears straight, while there is a sudden loss of the compactness of the bubbly core at a certain distance from the source. At higher gas flow rates the com-

Fig. 8

pactness is lost at a shorter height, accompanied by formation of discrete alternate vortices and discrete vortices do not form at smaller gas flow rates. As can be seen from the comparison, the sinuous instability and the discrete alternate vortices of plane bubble plumes were captured preferably using the present two-way coupling model, which indicates the complicated interaction between the two phases can be simulated successfully with our two-way coupling model.

Bubble distribution of the rising of single bubble cluster (ag 5 0:0135, QD=UT 5 0:1)


AUGUST 2014, Vol. 136 / 081301-7


Fig. 9 Liquid Q value field of the rising of single bubble cluster (ag 5 0:0135,QD=UT 5 0:1)

Alam and Arakeri [51] also measured the mean plume centerline velocity at a height of 50 mm above the plume source. They found that the plume remains laminar for the range of gas flow rates studied and the centerline velocity increases linearly with the gas flow rate. Figure 6 shows the comparison of the variation of plume centerline velocity between the present simulation and the experimental measurement of Alam and Arakeri [51]. Though the present simulation result is smaller than the experimental result, the linearly increased trend is preferably simulated. As can be seen from the comparisons, the good agreement with the experimental observation can undoubtedly validate our two-way coupling model. The initial position of the bubbles for two test configurations studied in this work, namely, a single bubble cluster and two tandem bubble clusters, are shown in Fig. 7. The diameter of the bubbles is 0.15 mm and the bubble cluster is spherical with a diameter of 10 mm. In both cases the bubble clusters are initially released on the central axis of the computational region at t* ¼ 0. The distance between the center of bubble cluster and bottom of the computational region is 25 mm. For the configuration of a single bubble cluster, two cases are considered with 2000 or 4000 bubbles in the cluster and the local void fraction ag, which is defined as the ratio of the total volume of bubbles to the initial volume of the bubble cluster, is 0.00675 or 0.0135, respectively. While for the cases of two tandem bubble clusters, there are 2000 bubbles in each bubble cluster and the local void fraction is 0.00675. The initial position of the lower bubble cluster is the same as in the single bubble cluster cases and the distance between two bubble clusters varies from 20 mm to 30 mm. 081301-8 / Vol. 136, AUGUST 2014

4

Results and Discussion

4.1 Single Bubble Cluster. First, a single bubble cluster that contains small gas bubbles rising in water is simulated. Initially, 4000 bubbles are randomly distributed in the cluster and released free. Figure 8 shows the evolution of bubbles in the cluster at successive time instances and in Fig. 9 the induced flow field is illustrated by the Q field. In the figure, the initial void fraction in the bubble cluster is 0.0135. The Q criterion, introduced by Hunt et al. [52], has become one of the most widely used local vortexidentification criterion in recent years. The Q value is given by 1 (27) Q ¼ ðkXk kSkÞ 2 where S is the symmetric component of the velocity gradient tensor with Sij ¼ 1=2ð@ui =@xj þ @uj =@xi Þ while X is the antisymmetric component with Xij ¼ 1=2ð@ui =@xj @uj =@xi Þ. The Q value is normalized by the initial bubble cluster diameter and the terminal rising velocity of a single bubble. As can be seen from the figure, initially the bubble cluster is spherical as it is released free. The bubble cluster starts to rise due to the buoyant force. The motion of the bubbles induces ambient liquid flow, and thus, a small vortex ring structure is generated immediately after the release of bubbles (t* ¼ 0.1). Because of the entrainment of the vortex ring structure, bubbles in the center of the cluster rise more quickly than those in the periphery and the bubble cluster deforms into convex gradually while rising. At the same time, the small vortex ring led by the rising bubble cluster is strengthened Transactions of the ASME


Fig. 10 Liquid velocity field of the rising of single bubble cluster on the x 5 0.1 m plane (ag 5 0:0135)

by the rising bubbles very quickly. Simultaneously, both the shape of the bubble cluster and distribution of the bubbles in cluster are changed by the vortex ring structure. Both the vortex ring and the bubble cluster expand significantly at about t* ¼ 3.0. Then the vortex ring structure is unable to catch up with the bubble cluster, and the distance between the vortex ring and the bubble cluster gets larger and larger. A series of secondary vortex rings is generated around the primary vortex ring until about t* ¼ 5.6. After that the secondary vortex ring structures are gradually weakened by the ambient water, and they disappear at about t* ¼ 13.0, while the distance between the vortex ring structure and the bubble cluster increases further. The rate of deformation of the bubble cluster decreases and the vortex ring structure rises more and more slowly after that. The strength of the primary vortex ring decreases gradually until its final fading. At t* ¼ 28.0, the strength of the vortex ring is so low that the vortex ring structure can only partially be visualized. Figures 10 and 11 show the snapshots of the bubble distribution and flow field on the central vertical plane (x ¼ 0.1 m plane). In the figure the initial void fraction in the bubble cluster is 0.0135. Figure 10 depicts the liquid velocity field and Fig. 11 the liquid vorticity field, in which the vorticity is normalized by the terminal Journal of Fluids Engineering

velocity of a single bubble and the initial diameter of the bubble cluster. Bubbles in a thin slice are superimposed on the figure to demonstrate the relative position between the bubbles and the vortex structure. The center of the slice is positioned on the x ¼ 0.1 m plane and the thickness of the slice is 0.1D. Two almost symmetrical vortex structures with opposite rotation can be identified in the flow field. Considering that this plane is the central vertical cross section in the flow filed, the two symmetrical vortices indicate that a vortex ring structure has been formed. The vortex ring structure is generated at the very beginning of the rising of the bubble cluster, and it is, subsequently, strengthened by the rising bubbles. The upward velocity at the center of the vortex ring becomes largest at about t* ¼ 0.6, which is due to the bubbles initially being released free and so the acceleration is very large. With the rising of bubbles, the bubble cluster expands quickly horizontally and on the central vertical plane (x ¼ 0.1 m) it deforms into an “arch,” which indicates the bubble cluster is whirled up by the vortex ring structure induced by the rising of themselves. The vortex ring is led by the rising bubble cluster and keeps on rising and expanding in diameter. There is also an obvious difference in the rising velocity of the vortex ring and the bubble cluster. The bubble cluster starts to deviate from the vortex ring gradually, while the vortex AUGUST 2014, Vol. 136 / 081301-9


Fig. 11 Liquid vorticity field of the rising of single bubble cluster on the x 5 0.1 m plane (ag 5 0:0135, xx 5 xx D=UT )

Fig. 12 Effect of void fraction on flow field of the rising of single bubble cluster on the x 5 0.1 m plane (xx 5 xx D=UT , t* 5 16.0)

081301-10 / Vol. 136, AUGUST 2014

Transactions of the ASME


Fig. 13 Comparison of the evolution of average bubble rising velocity of single bubble cluster

ring is strengthened by the cluster until about t* ¼ 3.0. The distance between the vortex ring and the bubble cluster becomes larger and larger with time. Bubbles continue to rise quickly due to buoyancy while the vortex ring structure is left behind, travel-

ing slowly upward by its self-induced velocity. Hence, the strength of the vortex ring structure is weakened by the ambient flow gradually. The rising velocity of vortex ring decreases significantly, while the liquid velocity in the peripheral area increases gradually. The acceleration effect of mutual interaction becomes weaker and weaker as the distance between the bubble cluster and vortex ring gets larger and larger. Therefore, the bubble cluster expands and rises slower, which induces a distinct upward velocity of the fluid in the wake. The upward velocity of the fluid is larger at the edge of the bubble cluster and the bubble distribution in the slice is like an “M.” The effect of void fraction on the bubble distribution and induced flow field is then considered. Figure 12 shows the liquid velocity field and vorticity field for different void fractions on the x ¼ 0.1 m plane at t* ¼ 16.0. In Figs. 12(a) and 12(b), the void fractions are 0.0135 and 0.00675, respectively. It can be seen from the figure that both the positions of the vortex structure and bubble cluster are a little higher in the case of the larger void fraction than that of the smaller one, which indicates that both the induced liquid velocity and the rising velocity of the bubble cluster are larger for the larger void fraction case. The width of the vortex structure and the size of vortex core are also larger for the

Fig. 14 Liquid velocity field of the rising of two tandem bubble clusters on the x 5 0.1 m plane (cluster distance: 30 mm)


AUGUST 2014, Vol. 136 / 081301-11


Fig. 15 Liquid vorticity field of the rising of two tandem bubble clusters on the x 5 0.1 m plane (cluster distance: 30 mm, xx 5 xx D=UT )

larger void fraction case. The bubble slice is flatter for the case of ag ¼ 0.00675, forming an arch, while the distribution for the case of ag ¼ 0.0135 is more like an “M.” In the case of ag ¼ 0.0135 the bubbles in the slice are mainly distributed on the sides while in the ag ¼ 0.00675 case the bubbles are distributed uniformly. Figure 13 shows the comparison of the evolution of average bubble rising velocity for different void fractions. The velocities are normalized by the terminal rising velocities of a single bubble. As can be seen from the figure, in both cases, three stages can be identified in the evolution of the average bubble velocity after the release of bubbles. In the first stage, the bubble cluster rises due to buoyancy and induces a vortex ring structure, which is strengthened by the rising bubble cluster. At the same time, the bubbles are also intensively accelerated by the vortex ring structure, so there is a straight climb of the average bubble velocity from zero to above the rising velocity of a single bubble. In both cases the average bubble velocity reaches a maximum value at about t* ¼ 0.6, which is larger for the larger void fraction case. Then the bubble cluster is gradually separated from the vortex ring structure since there is an obvious difference between the rising velocity of the bubble cluster and the vortex ring. The acceleration effect of the vortex ring becomes weaker because of the deviation and so there is a rapid decrease on the average rising velocity until about t* ¼ 3.0. As the distance between the bubble cluster and the vortex ring gets larger and larger, the acceleration effect of the vortex 081301-12 / Vol. 136, AUGUST 2014

ring gets weaker and weaker and becomes negligible at the last stages. Therefore, the average rising velocities decrease and tend to become equal to the terminal velocity of a single bubble. While in all stages the average bubble velocities are larger than the terminal velocity of a single bubble, which indicates that there is an aggregation effect in the rising process of bubble cluster and bubbles are accelerated by the induced flow all the time. 4.2 Two Tandem Bubble Clusters. In this section the rising of two tandem bubble clusters is studied, in which the rising of the lower bubble cluster will be affected by the induced field of the upper bubble cluster. Figures 14 and 15 show the snapshots of bubble distribution and flow field on the x ¼ 0.1 m plane of the rising of two tandem bubble clusters. Figure 14 depicts the liquid velocity field and Fig. 15 the liquid vorticity field. In each bubble cluster there are 2000 bubbles with 0.15 mm diameter, and the void fraction in both bubble clusters is 0.00675. The initial distance between the two bubble clusters is 30 mm. Bubbles in a thin slice are also superimposed on the figure to demonstrate the relative position between bubbles and vortex structure. As the bubble clusters are rising, two vortex rings are generated and the distance between the two bubble clusters gets smaller and smaller with time. Both bubble clusters deform significantly and the rising of the lower bubble cluster is strongly affected by the induced flow Transactions of the ASME


Fig. 16 Merging of the two vortex ring structures in the rising process of two tandem bubble clusters (cluster distance: 30 mm, QD=UT 5 0:3)

Fig. 17 Effect of initial cluster distance on the flow field on the x 5 0.1 m plane on the rising process of two tandem bubble clusters (cluster distance: 20 mm, xx 5 xx D=UT )

of the upper bubble cluster. At t* ¼ 4.7, the upper bubble slice deforms into an arch, which is similar to that of a single bubble cluster, while the lower bubble slice deforms into an “n” because the induced velocity is larger at the center of the vortex ring. Subsequently, the two vortex rings start to merge and a new vortex ring is generated at about t* ¼ 6.5. The lower bubble slice also deforms gradually into a smaller arch. At t* ¼ 9.0, the two bubble clusters are almost merged and the bubbles in the upper cluster are distributed in a wider range. The acceleration effect becomes weaker and weaker as the distance between the bubble cluster and vortex ring gets larger and larger, and therefore, the bubble cluster expands and rises slowly after that. Since initially there are two tandem bubble clusters in the flow field, each bubble cluster will induce one vortex ring structure. With the rising of the clusters, the two vortex ring structures will Journal of Fluids Engineering

gradually merge. As shown in Fig. 16(a), two tandem vortex rings are generated and strengthened by each bubble cluster. The interactions between the two vortex rings also lead to the expansion of the upper vortex ring and the shrinkage of the lower one, as shown in Fig. 16(b) where the upper vortex ring is apparently larger than the lower. As a result, the upper vortex ring decelerates and the lower vortex ring speeds up, catching up and merging with the upper. Finally, one new vortex ring is formed as shown in Fig. 16(c). The effect of initial distance between the two bubble clusters on bubble distribution and induced flow field is investigated in Figs. 17 and 18. The figures show the liquid vorticity field on x ¼ 0.1 m plane and the thin bubble slice is also superimposed on the figure. In Figs. 17 and 18 the distances between the two clusters are 20 mm and 30 mm, respectively. As can be seen from the AUGUST 2014, Vol. 136 / 081301-13


Fig. 18 Effect of initial cluster distance on the flow field on the x 5 0.1 m plane on the rising process of two tandem bubble clusters (cluster distance: 30 mm, xx 5 xx D=UT )

Fig. 19 Comparison of the evolution of average bubble rising velocity of two tandem bubble clusters

figure, in both cases the two vortex rings merge, although the bubble behavior and the flow fields induced by the rising bubble clusters are quite different. In the cluster distance 20 mm case, the distance between the two bubble slices is small, so the aggregation effect is stronger and the two vortex rings merge quickly after the release of the bubbles. While in the 30 mm case, the aggregation effect is stronger and the two vortex rings are almost fully developed before they approach each other. In the first stage, both bubble slices deform into an arch in both cases, while the range of bubble distribution is larger for the larger distance case. The upward velocity is the largest in the central area and the lower bubble slice deforms into an inverted “V” in Fig. 17 and an “n” in Fig. 18, as the two bubble clusters come close. The distributions of the two bubble slices are also different after the generation of the new vortex rings. In Fig. 17 the two bubble slices mix quickly and deform into an arch, while in Fig. 18 the lower bubble cluster cannot catch up with the upper one and the two bubble slices can 081301-14 / Vol. 136, AUGUST 2014

still be distinguished. The upper bubble slice is like an arch and the lower bubble slice is like another arch with smaller curvature radius distributed in a narrower range. The extent of the vortex structure is much larger for the larger distance case in all stages. Figure 19 shows the comparison of the evolution of average bubble rising velocity for the two different initial cluster distances. In Figs. 19(a) and 19(b) the distances between the two clusters are 20 mm and 30 mm, respectively. In Fig. 19(a), the evolution of average rising velocity of both clusters is similar to that for the case of a single bubble cluster. The same three stages of evolution of the average velocity can also be identified, that is, the straight climb stage, the rapid decrease stage, and the slow asymptotic decrease stage. Moreover, the lower cluster and the upper cluster reach different maximum rising velocities at different times. The maximum rising velocity for the lower cluster is much larger since its motion is strongly affected and accelerated by the flow field induced by the upper bubble cluster. There is Transactions of the ASME


only one peak value of the average velocity in Fig. 19(a), while there are two peak values in Fig. 19(b). In Fig. 19(b), the initial distance between two bubble clusters is larger, so the lower bubble cluster will catch up with the upper bubble cluster later. The two peaks are observed in the evolution of the average rising velocity of the lower bubble cluster. The first peak is similar to the case of a single bubble cluster. The rising velocities almost coincide in this stage, indicating that the interaction between the two clusters is negligible at this time. Later on, the lower bubble cluster is strongly affected by the flow field induced by the upper bubble cluster. The acceleration of the lower bubble cluster due to its interaction with the upper vortex ring, results in the second peak, which is even higher than the first one. In both cases the rising velocity is larger than the terminal velocity of a single bubble, which means that bubbles are accelerated by the flow induced by themselves in all cases. The aggregation effect can be clearly identified from the peak velocities in the figure. For the 20 mm cluster distance case, the peak velocity of the lower cluster is about 2.7 times of the terminal velocity of the single bubble and the peak average velocity of two clusters is about two times larger, while for the 30 mm cluster distance case, both the peak velocity of the lower cluster and two clusters are about 1.7 times of the terminal velocity of the single bubble.

5

Conclusions

The behavior of air bubble clusters rising in water and the induced flow fields are numerically studied using a threedimensional two-way coupling algorithm based on the VIC method. The small gas bubbles in the clusters rise due to the buoyant force, setting in motion ambient liquid. Thus, a small vortex ring structure is generated immediately after the release of each cluster. The vortex rings are led by the bubble cluster and expand with them. At the same time, the bubbles are accelerated by the vortex ring structure because of the aggregation effect, and the shape and distribution of the bubble cluster are also changed significantly under its influence. The mechanism of the aggregation effect in the rising process of bubble cluster is revealed and the transient processes of the generation, rising, strengthening, and separation of vortex ring structure with bubble clusters are illustrated and analyzed in detail. The effect of void fraction and distance between two bubble clusters are studied and the corresponding bubble behavior and flow field are presented. The simulation clarified the interaction between the bubble motion and the induced water flow satisfactorily. The main conclusions of the study are summarized as follows: (1) A cluster of bubbles will be accelerated by the vortex ring structures induced by the rising of themselves and so the average rising velocity is larger than the terminal rising velocity of a single bubble. (2) Bubbles in the center of the cluster rise more quickly than those in the periphery because of the acceleration of vortex ring structure. (3) The initial local void fraction affects the peak average rising velocity significantly, while there is only slightly effect on the stable average rising velocity. (4) The merging of two tandem vortex rings will occur with the rising of tandem bubble clusters. Due to the interaction between the two merging vortex rings, bubbles in the lower cluster rise more quickly while expand less in horizontal direction. (5) The average bubble rising velocity is larger for the smaller cluster distance case.

Acknowledgment This study was supported by the Natural Science Foundation of China (51176152, 51121092), the International Science & Technology Journal of Fluids Engineering

Cooperation Plan of Shaanxi Province (2013KW30-05) and the Fundamental Research Funds for the Central Universities.

References [1] Mudde, R. F., 2005, “Gravity-Driven Bubbly Flows,” Annu. Rev. Fluid Mech., 37, pp. 393–423. [2] Cheng, M., Hua, J. S., and Lou, J., 2010, “Simulation of Bubble-Bubble Interaction Using a Lattice Boltzmann Method,” Comput. Fluids, 39(2), pp. 260–270. [3] Chiang, T. P., Sheu, W. H., and Hwang, R. R., 1997, “Three-Dimensional Vortex Dynamics in a Shear-Driven Rectangular Cavity,” Int. J. Comput. Fluid D., 8(3), pp. 201–214. [4] Cottet, G. H., and Poncet, P., 2003, “Advances in Direct Numerical Simulations of Three-Dimensional Wall-Bounded Flows by Vortex-in-Cell Methods,” J. Comput. Phys., 193(1), pp. 136–158. [5] Sohn, S. I., and Hwang, W., 2005, “Numerical Simulations of Vortex Sheet Evolution in Stratified Shear Flow,” J. Phys. Soc. Jpn., 74(5), pp. 1472–1478. [6] Borthwick, A. G. L., and Barber, R. W., 1992, “Numerical-Simulation of JetForced Flow in a Circular Reservoir Using Discrete and Random Vortex Methods,” Int. J. Numer. Meth. Fl., 14(12), pp. 1453–1472. [7] Ogami, Y., and Fukumoto, K., 2010, “Simulation of Combustion by Vortex Method,” Comput. Fluids, 39(4), pp. 592–603. [8] Chen, B., Wang, C., Wang, Z. W., and Guo, L. J., 2009, “Investigation of GasSolid Two-Phase Flow Across Circular Cylinders With Discrete Vortex Method,” Appl. Therm. Eng., 29(8–9), pp. 1457–1466. [9] Akbari, M. H., and Price, S. J., 2010, “Wake Synchronization for a Pair of Staggered Cylinders in Cross-Flow With the Upstream Cylinder in Transverse Oscillation,” Eng. Appl. Comp. Fluid, 4(1), pp. 58–70. [10] Mustto, A. A., and Bodstein, G. C. R., 2011, “Subgrid-Scale Modeling of Turbulent Flow Around Circular Cylinder by Mesh-Free Vortex Method,” Eng. Appl. Comp. Fluid, 5(2), pp. 259–275. [11] Chen, B., and Wang, Z. W., 2009, “Simulation of Interaction Between Two Vortex Rings,” Prog. Comput. Fluid Dyn., 9(3–5), pp. 292–299. [12] Wang, Z. W., and Chen, B., 2012, “Numerical Investigation of the Evolution of Elliptic Vortex Ring,” Prog. Comput. Fluid Dyn., 12(1), pp. 19–26. [13] Uchiyama, T., and Matsumura, S., 2010, “Three-Dimensional Vortex Method for the Simulation of Bubbly Flow,” ASME J. Fluids Eng., 132(10), p. 101402. [14] Wang, Z. W., Uchiyama, T., and Chen, B., 2013, “Numerical Simulation of the Interaction Between Vortex Ring and Bubble Plume,” Appl. Math. Model., 37(24), pp. 10007–10026. [15] Milane, R. E., and Abdolhosseini, R., 2004, “Development of a ThreeDimensional Vortex-in-Cell Method for a Spatially Growing Uniformly Sheared Flow,” Int. J. Comput. Fluid D., 18(1), pp. 47–69. [16] Leonard, A., 1980, “Vortex Methods for Flow Simulation,” J. Comput. Phys., 37(3), pp. 289–335. [17] Greengard, L., and Rokhlin, V., 1987, “A Fast Algorithm for Particle Simulations,” J. Comput. Phys., 135(2), pp. 280–292. [18] Qian, L., and Vezza, M., 2001, “A Vorticity-Based Method for Incompressible Unsteady Viscous Flows,” J. Comput. Phys., 172(2), pp. 515–542. [19] Christiansen, J. P., 1973, “Numerical Simulation of Hydrodynamics by the Method of Point Vortices,” J. Comput. Phys., 13(3), pp. 363–379. [20] Cheng, M., Chew, Y. T., and Luo, S. C., 1997, “A Hybrid Vortex Method for Flows Over a Bluff Body,” Int. J. Numer. Meth. Fl., 24(3), pp. 501–523. [21] Abdolhosseini, R., and Milane, R. E., 2000, “On the Effect of Vortex Grid Density in the Vortex-in-Cell Simulation of Mixing Layers,” Int. J. Comput. Fluid D., 13(2), pp. 161–183. [22] Liu, C. H., 2001, “A Three-Dimensional Vortex Particle-in-Cell Method for Vortex Motions in the Vicinity of a Wall,” Int. J. Numer. Meth. Fl., 37(5), pp. 501–523. [23] Uchiyama, T., and Naruse, M., 2004, “Numerical Simulation for Gas-Particle Two-Phase Free Turbulent Flow Based on Vortex in Cell Method,” Powder Technol., 142(2–3), pp. 193–208. [24] Chorin, A. J., 1973, “Numerical Study of Slightly Viscous Flow,” J. Fluid Mech., 57(4), pp. 785–796. [25] Cottet, G. H., and Koumoutsakos, P., 2000, Vortex Methods: Theory and Practice, Cambridge University Press, Cambridge, UK. [26] Creuse, E., Giovannini, A., and Mortazavi, I., 2009, “Vortex Simulation of Active Control Strategies for Transitional Backward-Facing Step Flows,” Comput. Fluids, 38(7), pp. 1348–1360. [27] Beaudoin, A., Huberson, S., and Rivoalen, E., 2003, “Simulation of Anisotropic Diffusion by Means of a Diffusion Velocity Method,” J. Comput. Phys., 186(1), pp. 122–135. [28] Milane, R. E., 2004, “Large Eddy Simulation (2D) Using Diffusion–Velocity Method and Vortex-in-Cell,” Int. J. Numer. Meth. Fluids, 44(8), pp. 837–860. [29] Cottet, G. H., Koumoutsakos, P., and Salihi, M. L. O., 2003, “Vortex Methods With Spatially Varying Cores,” J. Comput. Phys., 162(1), pp. 164–185. [30] Degond, P., and Mas-Gallic, S., 1989, “The Weighted Particle Method for Convection-Diffusion Equations,” Math. Comput., 53(188), pp. 485–507. [31] Winckelmans, G. S., and Leonard, A., 1993, “Contributions to Vortex Particle Methods for the Computation of Three-Dimensional Incompressible Unsteady Flows,” J. Comput. Phys., 109(2), pp. 247–273. [32] Koumoutsakos, P., Leonard, A., and Pepin, F., 1994, “Boundary Conditions for Viscous Vortex Methods,” J. Comput. Phys., 113(1), pp. 52–61. [33] Ploumhans, P., and Winckelmans, G. S., 2000, “Vortex Methods for HighResolution Simulations of Viscous Flow Past Bluff Bodies of General Geometry,” J. Comput. Phys., 165(2), pp. 354–406.

AUGUST 2014, Vol. 136 / 081301-15


[34] Eldredge, J. D., Leonard, A., and Colonius, T., 2002, “A General Deterministic Treatment of Derivatives in Particle Methods,” J. Comput. Phys., 180(2), pp. 686–709. [35] Ploumhans, P., Winckelmans, G. S., Salmon, J. K., Leonard, A., and Warren, M. S., 2002, “Vortex Methods for Direct Numerical Simulation of ThreeDimensional Bluff Body Flows: Application to the Sphere at Re ¼ 300, 500, and 1000,” J. Comput. Phys., 178(2), pp. 427–463. [36] Schrader, B., Reboux, S., and Sbalzarini, I. F., 2010, “Discretization Correction of General Integral PSE Operators for Particle Methods,” J. Comput. Phys., 229(11), pp. 4159–4182. [37] Crowe, C. T., Gore, R. A., and Troutt, T. R., 1985, “Particle Dispersion by Coherent Structures in Free Shear Flows,” Part. Sci. Tech., 3(3–4), pp. 149–158. [38] Brecht, S. H., and Ferrante, J. R., 1989, “Vortex-in-Cell Simulations of Buoyant Bubbles in Three Dimensions,” Phys. Fluids A, 1(7), pp. 1166–1191. [39] Chen, H., and Marshall, J. S., 1999, “A Lagrangian Vorticity Method for TwoPhase Particulate Flows With Two-Way Phase Coupling,” J. Comput. Phys., 148(1), pp. 169–198. [40] Walther, J. H., and Koumoutsakos, P., 2001, “Three-Dimensional Vortex Methods for Particle-Laden Flows With Two-Way Coupling,” J. Comput. Phys., 167(1), pp. 39–71. [41] Yang, X. G., Thomas, N. H., Guo, L. J., and Hou, Y., 2002, Two-Way Coupled Bubble Laden Mixing Layer,” Chem. Eng. Sci., 57(4), pp. 555–564. [42] Yang, X. G., Huang, X. B., Rielly, C., Zouaoui, Z., and Guo, L. J., 2005, “The Modifications of a Plane Mixing Layer by Condensed Bubbles,” Chem. Eng. Sci., 60(24), pp. 6876–6886.

081301-16 / Vol. 136, AUGUST 2014

[43] Uchiyama, T., and Degawa, T., 2006, “Numerical Simulation for Gas-Liquid Two-Phase Free Turbulent Flow Based on Vortex in Cell Method,” JSME Int. J. B-Fluid. T., 49(4), pp. 1008–1015. [44] Uchiyama, T., 2013, “Numerical Simulation of Particle-Laden Gas Flow by Vortex in Cell Method,” Powder Tech., 235, pp. 376–385. [45] Uchiyama, T., and Yagami, H., 2008, “Numerical Simulation for the Collision Between a Vortex Ring and Solid Particles,” Powder Technol., 188(1–2), pp. 73–80. [46] Monaghan, J. J., 1985, “Extrapolating B Splines for Interpolation,” J. Comput. Phys., 60(2), pp. 253–262. [47] Bergdorf, M., and Koumoutsakos, P., 2006, “A Lagrangian Particle-Wavelet Method,” Multiscale Model. Sim., 5(3), pp. 980–995. [48] Koumoutsakos, P., 1997, “Inviscid Axisymmetrization of an Elliptical Vortex,” J. Comput. Phys., 138(2), pp. 821–857. [49] Clift, R., Grace, J. R., and Weber, M. E., 1978, Bubbles, Drops and Particles, Academic Press, New York. [50] Zheng, L., and Yapa, P. D., 2000, “Buoyant Velocity of Spherical and Nonspherical Bubbles/Droplets,” ASCE J. Hydraul. Eng., 126(11), pp. 852–854. [51] Alam, M., and Arakeri, V. H., 1993, “Observations on Transition in Plane Bubble Plumes,” J. Fluid Mech., 254, pp. 363–375. [52] Hunt, J. C. R., Wray, A. A., and Moin, P., 1988, “Eddies, Stream, and Convergence Zones in Turbulent Flows,” Center for Turbulence Research, Stanford University, Report No. CTR-S88, pp. 193–208.

Transactions of the ASME


Numerical simulation of glottal flow.

A method for computer simulation of ultrasound Doppler color flow images--I. Theory and numerical method.

Numerical simulation of supersonic gap flow.

Numerical Simulation of Dry Granular Flow Impacting a Rigid Wall Using the Discrete Element Method.

Flow studies in canine artery bifurcations using a numerical simulation method.

Numerical simulation of bubble induced shear in membrane bioreactors: effects of mixed liquor rheology and membrane configuration.

Modified Numerical Simulation Model of Blood Flow in Bend.

Numerical simulation of endocytosis: Viscous flow driven by membranes with non-uniformly distributed curvature-inducing molecules.

Numerical simulation of aerosol particle transport by oscillating flow in respiratory airways.

A Numerical Simulation of Cell Separation by Simplified Asymmetric Pinched Flow Fractionation.

Evaluation of Egorov's improved separate solution method for determination of low selectivity coefficients by numerical simulation.

Numerical Method of Characteristics for One-Dimensional Blood Flow.

Numerical simulation of pressure therapy glove by using Finite Element Method.

Numerical simulation of pressure waves in the cochlea induced by a microwave pulse.

Cell membrane deformation and bioeffects produced by tandem bubble-induced jetting flow.

Gas flow induced by ultrasonic cavitation bubble clouds and surface capillary wave.

Numerical simulation of pipeline deformation caused by rockfall impact.

Numerical modeling of bubble dynamics in viscoelastic media with relaxation.

Numerical Simulation of the blood flow behavior in the circle of Willis.

Hemodynamic study of TCPC using in vivo and in vitro 4D Flow MRI and numerical simulation.

Numerical simulation and clinical implications of stenosis in coronary blood flow.

Lattice Boltzmann numerical simulation and experimental research of dynamic flow in an expansion-contraction microchannel.

Hemodynamic numerical simulation and analysis of oscillatory blood flow in growing aneurysms.

Numerical simulation of steady flow in a two-dimensional total artificial heart model.