Cite this Article
Daemi Mahdi, Taeibi-Rahni Mohammad, Massah Hamidreza. Study of acoustic bubble cluster dynamics using a lattice Boltzmann model. Chinese Physics B, 2014, 24(2): 024302  
Permissions
Study of acoustic bubble cluster dynamics using a lattice Boltzmann model
Daemi Mahdia)†, Taeibi-Rahni Mohammada), Massah Hamidrezab)
Aerospace Engineering Department, Sharif University of Technology, Tehran, Iran
Physics Department, Ferdowsi University, Mashhad, Iran

Corresponding author. E-mail: mdaemi@ae.sharif.edu

Abstract

The search for the development of a reliable mathematical model for understanding bubble dynamics behavior is an ongoing endeavor. A long list of complex phenomena underlies the physics of this problem. In the past decades, the lattice Boltzmann method has emerged as a promising tool to address such complexities. In this regard, we have applied a 121-velocity multiphase lattice Boltzmann model to an asymmetric cluster of bubbles in an acoustic field. A problem as a benchmark is studied to check the consistency and applicability of the model. The problem of interest is to study the deformation and coalescence phenomena in bubble cluster dynamics, as well as the screening effect on an acoustic multi-bubble medium. It has been observed that the LB model is able to simulate the combination of the three aforementioned phenomena for a bubble cluster as a whole and for every individual bubble in the cluster.

Keyword: 43.35.+d; 43.25.+y; 47.55.dp; multiphase lattice Boltzmann model; acoustic field; multi-bubble; bubble cluster dynamics; cavitation
1. Introduction

Bubble cluster dynamics have a long history[1, 2] that is built on many studies. Modeling a detailed bubble behavior is very cumbersome.[3– 6] A bubble usually behaves asymmetrically and its moving boundary has a profound effect on its topology.[7] A comprehensive mathematical model describing the dynamics of a bubble cluster should facilitate its investigation as a single entity, studying each individual bubble, accounting for bubble size distributions and their interaction, and its interaction with external acoustical fields.[8] Clearly, this is not a simple description and, therefore, it does not yield a very attractive theory.[8, 9]

There is no accurate measuring method for parameters such as pressure and density as a function of position for the multi-bubble interior. However, there are reports on the failings of such direct experiments.[8, 10– 12] Hence, resorting to numerical simulation of bubble cluster dynamics is a viable option.[4, 8, 9, 13– 16] A substantial part of current numerical solutions to bubble cluster dynamics in the literature is devoted to macroscopic computational fluid dynamics (MaCFD), in which many aspects of bubble dynamics are examined. The performance of a model for equi-size bubbles was investigated in Ref. [16]. Direct numerical simulation was used to understand the coupling mechanism of cloud.[13] An approximate model was devised to analyze single bubble and bubble cluster interaction.[14] The internal thermal phenomena inside a bubble cluster were considered in Ref. [15]. The growth of a cluster of equi-size bubbles was studied in Ref. [4]. Hu Jing et al. inspected the influence of bubble cluster secondary sound radiation on the dynamics of a bubble cluster in an ultrasound field.[17] Many proposed correction factors were introduced into bubble dynamics equations for proper modeling of bubble– bubble interactions due to their complicated nature.[4, 9] Brennen used equi-size fine arrayed bubble to study cluster dynamics.[3] Wang Yong et al. concluded that in an acoustic field the bubble content plays an important role in multibubble dynamics.[18]

At the continuum level, the approximations that were made for deriving the mathematical models, or the simplifications that were assumed in order to solve them, exhibit inconsistencies between solutions and experimental observations.[4] To overcome these complications in modeling fluid flows, microscopic and mesoscopic views have gained popularity.[19, 20] On a mesoscopic scale, lattice Boltzmann model (LBM) encompasses a greater range of fluid flow phenomena in multiphase cases than macroscopic-based descriptions.[19, 21] LBM can be also thought of as an indirect measure where the experimental shortcomings prevail. Moreover, it renders many advantages compared with other models. Besides its inherent parallelizability, LBM has the ability to define arbitrary geometry, complicated boundaries, and surface interactions for single and multiphase flows. It also features domain scalability and handling multicomponent systems.[22] In recent years, many researchers have enjoyed employing parallel computing in LBM.[23, 24] Alapati et al. indicated that the LBM is well suitable to lending itself to parallel computing.[25] Exploiting the parallelizability capability of LBM has enhanced the method resources.

The present paper is devoted to the utilization of a lattice Boltzmann (LB) model for three-dimensional (3D) single-component multiphase fluid to investigate the bubble cluster dynamics in an acoustic field. We use a 121-velocity Shan– Chen LB model, which seems to be able to address complexities such as deformation, coalescence, and screening effect. To this end, we develop a parallel code using OpenMP on Intel FORTRAN 2013 software in a 24-thread cluster.

The rest of this paper is organized as follows. An LB model for a multiphase medium is presented in Section 2. The results and discussion covering the validation of the multi-phase LBM are given in Section 3, which indicates a study of a symmetric case as a base to evaluate the other two asymmetric cases. Finally, the conclusion is presented in Section 4.

2. Methods and models
2.1. Lattice-Boltzmann model

Even though LBM is a recently developed new model in computational fluid dynamics (CFD) (mesoscopic CFD), it has some promising features to be utilized in numerical simulations. Usually, a finite difference discretized continuum Boltzmann equation accompanied by Bhatnagar, Gross, and Krook (BGK) approximation of collision term is employed to drive a proper LB method.[26, 27]

A Boltzmann– BGK equation that describes the time evolution of the particle distribution function, f(r, e, t), with variable position, r, velocity, e, and time, t, is given as

where a and tr are the force per unit mass acting on the particles and relaxation time, respectively, and f0 is the Maxwell– Boltzmann equilibrium distribution function and is expressed as

where D, m, kB, ρ , u, and T are the space dimension, particle mass, Boltzmann constant, density, macroscopic velocity, and temperature, respectively. The resulting LBM model is

We use a 121-velocity model[28] according to the parameters and their values in Table 1, where ei and ω i represent velocities and weights, respectively. The macroscopic density, ρ , and velocity, u, are calculated from distribution function, fi, at each step:[29]

Hence, to use the LB Eq. (3), an expansion of Eq. (2) is in order. Thus, in order to recover correct hydrodynamics equations at Navier– Stokes level, we, in the isothermal mode, implement a 6th-order Hermitian-based expansion of the Maxwell– Boltzmann equilibrium distribution function[29] as follows:

Table 1. A nine-order accurate Gauss– Hermite quadrature parameters: the directions, magnitudes, ordinates, weights; a fully symmetric set (FS) of points are marked by FS.[28]
2.2. Multiphase LBM model

In our previous work, we showed the formation of a cavitation bubble by our 121-velocity LB model.[29] In that work, we chose a viable and less laborious approach to address our multiphase problem by using the Shan– Chen model.[29] In our Shan– Chen model, the attraction force between a particle and the particles in its first three neighborhoods is described by

and the values of the interaction strength, G, and gi, are adopted from Table 2. In addition, interaction potential, ψ , is given by the following relation:[29]

Here, pressure, p, is calculated from the Peng– Robinson equation of state[30]

where α (T) is obtained from

The variables Rg, ε , T, and Tc represent the gas constant, acentric factor, temperature, and critical temperature, respectively. The parameters, a and b indicate the attraction and repulsion parameters.

Table 2. Interaction strength parameters.
3. Results and discussion

We introduce the aforementioned multiphase LB model for an asymmetric cluster of bubbles in an acoustic field. The acoustic field is imposed by density variations through lattice boundaries.[31] First, we show the validity of our model by solving a single bubble problem via our LBM under the same conditions. The validation process is achieved through comparing the results with those corresponding results obtained by a Cash-karp 45 Runge– Kuttta solution of a Rayleigh– Plesset (single bubble problem) equation. Second, this model is used to investigate the deformation, coalescence of the bubbles and the screening effect of bubble layer(s) in the bubble cluster.

3.1. Validation

In this subsection, we discuss a single bubble dynamics in a standing acoustic field. This bubble, generally described by the time evolution of its radius, exhibits a highly nonlinear behavior. The problem at hand has been studied using various macroscopic theories, bearing different degrees of approximations.[32] However, to the best of our knowledge even now there has been no generally accepted model for the exact description of this phenomenon. The macroscopic models for single bubble dynamics resulted in varieties of the so-called Rayleigh– Plesset (RP) equations. The LBM could be regarded as a mesoscopic alternative to the conventional RP models.

The problem of a gas bubble in a spherical domain, which is exposed to a resonating acoustical field, is discussed. The parameters presented in Table 3 represent a typical sonoluminescence bubble.[33]

Table 3. Properties of medium.

In the present case, our LBM simulation is limited to the isothermal state equation, only. Hence, the comparison will be made solely with the isothermal section of bubble evolution period.[34] The validation compares these results with the corresponding results of the following RP equation:[35]

Here, R, ρ g, ρ l, pb, υ , and σ are the bubble radius, gas density, liquid density, bubble pressure, kinetic viscosity, and surface tension, respectively. External pressure, p∞ , is the sum of environmental pressure, p0, and acoustic pressure, pa. In LBM, the pressure (acoustic) is usually implemented via the density variation at the boundary.[31] In the following relation the applied sound wave is described as a density variation:

where Δ ρ is the density variation at the boundary, and f, t represent frequency and time.

Figure 1 depicts the validity of our LB model in the vicinity of isothermal behavior of the bubble. The existing relative inconsistency between two radius– time curves is due to the difference in the nature between the macroscopic Eq. (9) and that of our LB model span. In Eq. (9), contrary to LBM, a number of oversimplifying approximations are assumed, such as the incompressibility of a liquid, the constant sound speed in a liquid, and negligence of viscosity in liquid and gas.

Fig. 1. Radius of bubble as a function of time; a comparison between the results of LBM and RP model. All parameters are in SI LBM units.

3.2. Asymmetric multi-bubble cluster

We demonstrate the ability for our 121-velocity lattice Boltzmann model[29] to handle intricate phenomena in bubble dynamics. In our study, three different cases were selected to simulate complex geometries and complicated interactions involved in bubble cluster dynamics. The symmetric configuration (8 equi-size bubbles) case serves as a precedent to critique the other two asymmetric clusters. To have a more realistic assessment of the model, and at the same time keep computational cost as low as possible, the asymmetric cases are defined to have a cluster with a minimum number of bubbles.

The geometry of the problem is described by two cubic volumes, one is located inside the other, forming concentric cubes. The facing sides of the concentric cubes are parallel in pairs. The outer cube is a lattice indicating the domain of our computation.

The arrangement of the bubbles in the cluster is done in such a way that each corner, if any, is occupied by a bubble of the cluster. In our lattice, each of the occupying bubbles is laid out by its center on a corner of the inner cube. The geometry of the numerical apparatus remains intact for all of the cases and the same phenomena take part in the dynamics of the problem. We deliberately exclude a bubble from the cluster and a change in the size of a bubble in order to simulate asymmetric cases.

The acoustic field, Eq. (10), is set up by transmitting sound waves to the center of inner cube from all of the lattice planes.

In our study, only the isothermal part of bubble evolution is addressed since this LB model is derived under the assumption of constant temperature.[29]

3.2.1. Symmetrical configuration

In the eight-bubble lattice configuration there is symmetry in each individual bubble and also in the cluster as a whole. Each corner of the inner cube is occupied by one of the eight bubbles of the cluster and each bubble is placed by its center on a corner. All of the eight bubbles in the cluster are of the same size and are assumed to have the same content (Fig. 2(a)).

Fig. 2. Bubble cluster (8 symmetric bubbles) behavior with frequency f = 0.044, at times (a) 0.942-s LBM unit, (b) 14.5-s LBM unit, (c) 16.82-s LBM unit, and (d) 19.14-s LBM unit.

The inception bubble radii are all set at a value less than the critical bubble radius. At the beginning, the bubbles are in contraction mode and they collapse for a short while. Of course, as they sense the negative pressure wave of an acoustic field (sound waves) they start to grow (expansion mode). It is observed that during the growth period, the part of bubble facing the center of the cube grows less than the part facing the lattice planes. This is due to a great difference in the acoustic pressure received by the bubble from the sides. Bubbles are directly irradiated by sound waves from one side, and they receive attenuated acoustic waves from the other side. The attenuation is due to the screening effect because the sound waves are passing through the bubble layer(s). Hence, in the expansion mode, bubbles deform and lose their spherical shapes (Fig. 2(b)) because of the lesser pressure amplitude from one side. When bubbles growth reaches its maximum and the applied pressure is positive, the bubbles start to collapse. They also close in upon each other as they collapse (Fig. 2(c)). The pressure gradient (the force) is due to the bubbles coalescing. The coalescing process of the bubbles occurs after their contact, which follows the collapsing and closing in actions. The spherical symmetry of coalescing process is preserved (Fig. 2(d)). From the observations, the following perceptions are distinguished:

(i) Bubbles experience deformation and size change;

(ii) There exist size and shape distributions;

(iii) No symmetrical applied acoustic pressure is possible for each bubble in a cluster due to the screening effect;

(iv) When the positive pressure wave prevails, the bubbles start to coalesce;

(v) Symmetry of geometry is preserved even in the collapsing, closing in, and coalescing periods; and,

(vi) The contents of coalescing bubbles are stratified. Due to pressure gradient, in every coalescing bubble, the part of less dense contents inclines toward the center of the cluster (Fig. 3).

Fig. 3. Symmetric bubble cluster interior (8 symmetric bubbles) behavior with frequency f = 0.044; slice in x = 100, at time 16.24-s LBM unit.

In this case, using our LB model, the simulated contraction (Fig. 2(a)) and expansion of the bubbles and the cluster (Figs. 2(b) and 2(c)) are simulated. It also shows mutual attraction of bubbles and their collective migration toward the pressure antinode of applied acoustic field. The coalescing process is clearly demonstrated in our simulation. Thus, this symmetric problem is a justified criterion to assess the other two cases.

3.2.2. Asymmetrical configuration 1

The problem geometry remains unchanged and it turns into an asymmetric shape by expelling a bubble from the cluster arrangement (a 7-bubble setup), (see Fig. 4(a)). For this spatial distribution of the bubbles, everything occurs exactly as it happened in case I, except for what is surfaced after the coalescing process has started in our computational domain. In the collapsing cluster, in section where there are no bubbles, the liquid is ruptured, and some smaller bubbles appear (Fig. 4(c)). This is, possibly, because of the pressure influence generated by the bubble interactions with the acoustic field.

Fig. 4. Asymmetric bubble cluster (7 equi-size bubbles) behavior with frequency f = 0.044, at times (a) 1.256-s LBM unit, (b) 17.594-s LBM unit, and (c) 18.222-s LBM unit.

Briefly, the obtained results are similar to those of the previous case. In addition, small bubbles appear in the area where there were no bubbles.

3.2.3. Asymmetrical configuration 2

In this case, the physical configuration of the problem is different from those in the previous cases. Six bubbles of the same size are kept and one is replaced by a smaller size bubble while the eighth corner has no bubble (Fig. 5). The radii of bubbles are set at values less than the critical radius, and the radius of smaller bubble is much less than the other radii. The geometry and rest of the physics of the problem remain the same. Similar results are repeated, again. The tiny bubble does not take part in the coalescing process during the contraction mode (Figs. 6(a) and 6(b)). In addition, it shows a clear nonspherical deformation during the collapse period (Fig. 6(b)). This possibly occurs under the complete influence of pressure distribution caused by the mutual effects of the larger bubbles interaction with the sound field. The smaller bubble growth is in harmony with the expansion of the rest of the cluster (Fig. 6(a)). The rate of deformation and loss of sphericity for the smaller member of the bubble arrays is more intense than the other bubbles (Figs. 6(b) and 7(a)). The tiny bubble becomes vanishingly small as the coalescing is intensified (Fig. 7(b)). Also in the vicinity, a cavitation cloud occurs (Fig. 7(b)). The applied acoustic field affects the density distribution. In the expansion stage, there is no symmetry in this distribution, while the less dense content is distributed on the outer side of the cluster (Fig. 8(a)). On the other hand, in the contraction mode, the bubble contents stratify from outer side to the center of the cluster (Fig. 8(b)). As the cluster contracts, the bubbles closed in upon each other and coalesce, except for the tiny bubble. After a while, the tiny bubble and the surfacing cavitation cloud break away from the pressure antinode (Fig. 9).

Fig. 5. Asymmetric bubble cluster (6 equi-size and one smaller bubbles) at time 0.942-s LBM unit.

Fig. 6. Asymmetric bubble cluster (6 equi-size and one smaller bubbles) behavior with frequency f = 0.044, at times (a) 14.138-s LBM unit, and (b) 18.850-s LBM unit.

Fig. 7. Asymmetric bubble cluster (6 equi-size and one smaller bubbles) behavior with frequency f = 0.044, at times (a) 19.479-s LBM unit, and (b) 19.792-s LBM unit.

Fig. 8. Asymmetric bubble cluster (6 equi-size and one smaller bubbles) behavior with frequency f = 0.044, in times (a) 14.138-s LBM unit, and (b) 18.850-s LBM unit, slice is located at x = 100.

Fig. 9. Asymmetric bubble cluster (6 equi-size and one smaller bubbles) behavior with frequency f = 0.044, at time 20.107-s LBM unit; slice is located at x = 93.

4. Conclusions

Our observations reaffirm that the LBM is a promising mesoscopice numerical tool to address multiphase cases. A great deal of effort has been made to devise and develop schemes to illuminate the capabilities of LBM. We use our devised multiphase LBM to explore the phenomena occurring in a bubble cluster and its bubbles when subjected to an acoustic field. Our simulations are clearly capable of handling deformation, coalescing process, and bubble– bubble plus bubble– sound field interactions. Our 121-velocity multiphase LBM accounted for the macroscopic parameters and the complexities involved in the problem.

Our model well describes the isothermal part of the bubble cluster behavior. Of course, to consider the fast part of the cluster dynamics, our model must be extended for the adiabatic part of the bubble cluster life cycle.

Reference
1 van Wijngaarden L 1968 J. Fluid Mech. 33 465 DOI:10.1017/S002211206800145X [Cited within:1] [JCR: 1.415]
2 Brennen C E 2014 Cavitation and Bubble Dynamics(Proof stage) New York Cambridge University Press [Cited within:1]
3 Brennen C E 2011“An Introduction to Cavitation Fundamentals” Invited Paper, WIMRC Forum 2011, Warwick University UK authors.library.caltech.edu/28373/ [Cited within:2]
4 Sinden D, Stride E and Saffari N 2012 J. Phys. : Conference Series 353 012008 DOI:10.1088/1742-6596/353/1/012008 [Cited within:4]
5 Lahey R T, Taleyarkhan R P and Nigmatulin R I 2007 Nuclear Engineering and Desgin 237 1571 DOI:10.1016/j.nucengdes.2006.12.014 [Cited within:1]
6 Nigmatulin R I 1991 Dynamics of Multiphase Media 1 New York Hemisphere [Cited within:1]
7 Brenner M P 2002 Rev. Mod. Phys. 74 425 DOI:10.1103/RevModPhys.74.425 [Cited within:1] [JCR: 44.982]
8 Nasibullaeva E S and Akhatov I S 2005 Acoust. Phys. 51 705 link.springer.com/article/10.1134/1.2130902 [Cited within:4] [JCR: 0.421]
9 Seo J H, Lele S K and Tryggvason G 2010 Phys. Fluids 22 063302 DOI:10.1063/1.3432503 [Cited within:3] [JCR: 1.942]
10 Kedrinskii V K 2005 Hydrodynamics of Explosion Berlin Springer [Cited within:1]
11 Yoon S W, Crum L A, Prosperetti A and Lu N Q 1991 J. Acoust. Soc. Am. 89 700 DOI:10.1121/1.1894629 [Cited within:1] [JCR: 1.646]
12 Bass A, Putterman S, Merriman B and Ruuth S J 2008 Comput. Phys. 227 2118 DOI:10.1016/j.jcp.2007.10.013 [Cited within:1] [JCR: 3.078]
13 Chen H, Li S C, Zou Z G and Li S 2008 J. Hydrody. 20 689 DOI:10.1016/S1001-6058(09)60003-2 [Cited within:2]
14 Wang C H and Cheng J C 2013 Chin. Phys. B 22 014304 DOI:10.1088/1674-1056/22/1/014304 [Cited within:1] [JCR: 1.148] [CJCR: 1.2429]
15 Matsumoto Y and Yoshizawa S 2005 Int. J. Numer. Methods Fluids 47 591 onlinelibrary.wiley.com/doi/10.1002/fld.833/abstract [Cited within:1] [JCR: 1.093]
16 Shagapov V S, Gimaltdinov I K and Yudin A V 2002 High Temperature 40 256 link.springer.com/article/10.1023%2FA%3A1015211408766 [Cited within:2] [JCR: 0.492]
17 Hu J, Lin S Y, Wang C H and Li J 2013 Acta Phys. Sin. 62 134303 (in Chinese) wulixb.iphy.ac.cn/EN/abstract/abstract54619.shtml [Cited within:1] [JCR: 1.016] [CJCR: 1.691]
18 Wang Y, Lin S Y, Mo R Y and Zhang X L 2013 Acta Phys. Sin. 62 134304 (in Chinese) wulixb.iphy.ac.cn/EN/abstract/abstract54618.shtml [Cited within:1] [JCR: 1.016] [CJCR: 1.691]
19 Chen H, Goldhirsch I and Orszag S A 2008 J. Sci. Comput. 34 87 DOI:10.1007/s10915-007-9159-3 [Cited within:2]
20 Chen H, Teixeira C and Molvig K 1997 Int. J. Mod. Phys. C 8 675 DOI:10.1142/S0129183197000576 [Cited within:1] [JCR: 0.615]
21 Qian Y, d’Humieres D and Lallemand P 1992 Europhys. Lett. 17 479 DOI:10.1209/0295-5075/17/6/001 [Cited within:1] [JCR: 2.26]
22 Mishra S K, Deymier P A, Muralidharan K, Frantziskonis G and Pannala S 2010 Ultrasonics Sonochemistry 17 258 DOI:10.1016/j.ultsonch.2009.05.014 [Cited within:1] [JCR: 3.516]
23 Biferale L, Mantovani F, Pivanti M, Pozzati F, Sbragaglia M and Scagliarini A 2011 Proc. Comput. Sci. 4 994 DOI:10.1016/j.procs.2011.04.105 [Cited within:1]
24 Körner C, Pohl T, Rüde U, Thürey N and Zeiser T 2006 Parallel Lattice Boltzmann Methods for CFD Applications, in Numerical Solution of Partial Differential Equations on Parallel Computers 2006 Berlin Springer [Cited within:1]
25 Alapati S, Kang S and Suh Y K 2009 J. Mech. Sci. Technol. 23 2492 DOI:10.1007/s12206-009-0422-4 [Cited within:1] [JCR: 0.616]
26 He X and Luo L S 1997 Phys. Rev. E 56 6811 DOI:10.1103/PhysRevE.56.6811 [Cited within:1] [JCR: 2.313]
27 He X and Luo L S 1997 Phys. Rev. E 55R6333 DOI:10.1103/PhysRevE.55.R6333 [Cited within:1] [JCR: 2.313]
28 Nie X B, Shan X and Chen H 2008 Europhys. Lett. 81 34005 DOI:10.1209/0295-5075/81/34005 [Cited within:1] [JCR: 2.26]
29 Daemi M, Taeibi-Rahni M and Massah H 2014 Journal of Acoustical Engineering Society of IranAccepted for Publication, 2 in Persian [Cited within:7]
30 Yuan P and Schaefer L 2006 Phys. Fluids 18 042101 DOI:10.1063/1.2187070 [Cited within:1] [JCR: 1.942]
31 Sukop M C and Thorne D T 2006 Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers Berlin, Heidelberg Springer-Verlag [Cited within:2]
32 [Cited within:1]
33 Nigmatulin R I, Akhatov I S, Vakhitova N K and Lahey R T 2000 J. Fluid Mech. 414 47 DOI:10.1017/S0022112000008338 [Cited within:1] [JCR: 1.415]
34 Nigmatulin R I, Akhatov I S, Topolnikov A S, Bolotnova R K, Vakhitova N K, Lahey R T and Taleyarkhan R P 2005 Phys. Fluids 17 107106 scitation.aip.org/content/aip/journal/pof2/17/10/10.1063/1.2104556 [Cited within:1] [JCR: 1.942]
35 Brennen C E 1995 Cavitation and Bubble Dynamics New York Oxford University Press [Cited within:1]