A numerical study of the early-stage dynamics of a bubble cluster ＊

Ya-zhen Shi, Kai Luo, Xiao-peng Chen, , Dai-jin Li

1. School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, China

2. School of Mechanics, Civil Engineering and Architecture, Northwestern Polytechnical University, Xi’an 710129,China

Abstract: The dynamics of multiple cavitating bubbles is numerically simulated, with the ambient pressure lower than the saturated vapor pressure, using a pseudopotential lattice Boltzmann method (LBM) coupled with the Carnahan-Starling equation of state.Dual-bubble and multi-bubble systems are tested, and the method for the bubble cluster is validated. It is found that the bubble can either grow or collapse in the early stage, depending on the configuration of the bubble cluster, characterized by the bubble number, the inter-bubble distance and the initial radii. In the induced flow, the bubbles are mutually interacted. Scaling relations of the interaction are proposed according to the numerical results. With consideration of the interactions, the simplified Rayleigh-Plesset equations(RPEs) for multiple bubbles can describe the evolution of the bubbles approximately. The results may serve as the basis for improved cavitation models.

Key words: Cavitation, bubble dynamics, interaction, lattice Boltzmann method (LBM), Rayleigh-Plesset equation (RPE)

Introduction

The cavitation involves a liquid-to-vapor phase transition due to the low ambient pressure to form small vapor cavities. It occurs in various hydraulic machineries, such as pump, propeller, artificial heart and ultrasonic cleaner[1]. It is found that the snapping shrimp can induce cavitation for predation, the maximum swimming speed of the dolphin is limited by the occurrence of cavitation[2], and the coastal erosion has a special cavitating mechanism. The pioneering work dates back to Besant and Lord Rayleigh,who explored the dynamics of a bubble in the fluidic environment. The analyses were further extended with considering the effects of the surface tension by Plesset. Hitherto, the Rayleigh-Plesset equation (RPE) is one of the most widely used equations in hydrodynamics. The theories provide support for the development of modern technologies,including HIFU[3], sonochemical reactions[4]and so forth.

In real constitutions, large or huge number of bubbles may exist in the liquid phase. The bubbles are in sophisticated motions, such as expanding,oscillation and translation. The motions are coupled with the flow field, to induce extra influences on the behaviors of the bubbles[5]. Mettin et al.[6]employed Keller-Miksis equations, coupled with the bubble pressure emission terms, to study the interaction between bubbles in a strong acoustic field. Ida[7]found the suppression effect of the cavitation inception due to the interactions. It is suggested that a larger bubble can prominently decrease the cavitation threshold pressure of a smaller bubble nearby. On the side of experimental study, Bremond et al.[8]explored the expansion and the collapse of the bubbles confined on a patterned plate. The evolution of the inner bubbles is shielded by the outer ones, and an extended RPE is validated. In the recent years, the study is extended to more specific fields for the purpose of the bubble manipulation[9-11]. Meanwhile, the numerical simulation was developed as a powerful tool in exploring the bubble dynamics[12]. Zhang et al.[13]used a volume of fluid (VOF) method to simulate the deformation, the growth and the interaction of bubbles in the rising process. The effects of the fluid properties and the bubble volume were analyzed. Ma et al.[14]employed a Eulerian- Lagrangian approach to study the dynamics of a bubble cloud exposed in a sinusoidal pressure field. The cloud was shown to have a nonlinear resonance to the excitation.

Lattice Boltzmann methods (LBMs) were developed as an efficient method for the flow simulation[15]. Compared with the conventional numerical methods, the LBM does not need anad hocmodel to track the free surface in multiphase flows. A typical multiphase LBM model is the Shan-Chen model. High density ratio can be achieved by incorporating different equations of state (EOS) of the real gas. Sukop and Or[16]employed the Shan-Chen LBM to simulate the bubble dynamics. The growth of a natural cavitating bubble in the quiescent liquid and the shear flow was simulated by Chen et al.[17-18]. It is shown that the shear flow seldom influences the bubble growth. The stress-induced cavitation was captured by using a multirange Shan-Chen LBM[19],and it is found that the cavitation number is not a good predictor of the onset of the cavitation.

1. Numerical method and validations

1.1 Multiphase lattice Boltzmann method

The LBM is based on the classical particle kinetic theories. By discretizing time and space according to the lattice model (see Fig. 1), the governing equation can be expressed as

wherefi(i= 1,2,…) denotes the density population of thei-th microscopic velocity,ei. The left-hand-side (LHS) of the equation indicates a streaming process: the particles jump fromxtox+eiΔtin a time interval Δt, while the right-hand-side (RHS) shows the gain/loss of the populations due to the collision (the first term) and the external force (the second term). In the equation, the Bhatnagar-Gross-Krook (BGK) collision model is applied, which implies that the density distribution approaches an equilibrium one BGK after the collision,andτis a relaxation time. The equilibrium distribution is assumed to be Maxwellian, and for the D3Q19 cubic model (Fig. 1), it is expressed as

In the expression, the weightsωiare 1/18 fori= 1,2,… ,6, 1/36 fori=7,8,…,18, and 1/3 for the rest particles (i=0), respectively. Furthermore, the local macroscopic velocity and the density of the fluid can be estimated by using the moments offi:. The Navier-Stokes equations can be recovered by applying the Chapman-Enskog analysis

where the forcing term is neglected. Meanwhile, it can be deduced that the kinematic viscosity[17].

Fig. 1 A schematic diagram of three-dimensional and nineteen velocity (D3Q19) lattice model

According to the mean field theories[20], the phase separation can be explained by the concept of the inter-molecular force, which can be introduced into the LB equation by using a properly designed forcing term. Following Shan and Chen[20], we set the force as

whereGdenotes the interacting strength (among the particles), andψ(x,t) is the effective density. To achieve a high density ratio in the simulations,ψis determined according a certain EOS of the fluid[21]

wherep≡p(ρ) for an isothermal process. Then,various EOS’s can be used. Among various EOS’s,the Carnahan-Starling (C-S) EOS enjoys a better performance with a large liquid-vapor density ratio, a wide stable temperature range and a lower spurious current[21]. The C-S EOS is expressed as

Since the bubble dynamics is our main concern,the fundamentals and the details of the method are not described here and can be found in our previous papers[17-18]and other papers[16,20].

1.2 Modified RPE’s

With the spherical bubble assumption, the RPE is applied to describe the dynamics of a cavitation bubble in a space filled with liquid[22]

wherer∞?Ris the location of the far field boundary with the pressureP∞,Ris the instantaneous bubble radius.lρ,lμandSare the liquid density, viscosity and interfacial tension,respectively. The subscriptBdenotes “bubble”, and the over-dot represents time derivatives. The first term on the RHS is the pressure difference, the second term represents the effect of the viscosity forces, and the third term shows the surface tension effects,respectively.

As can be shown, the RPE describes the liquid phase flow. The bubbles in a bubble cluster would affect each other, if they are in a certain motion, and an interacting force will be induced. If only the volumetric motions of the spherical bubbles with fixed locations are considered and modelled as the monopoles in the flow field, the pressure disturbed by theith bubbles can be expressed as

whereDiis the distance between the bubble and a probe point[6]. The local increment of the pressure in a multi-bubble system is ∑pi. Adding this extra pressure to Eq. (8), the modified RPE becomes

In the later analyses, the above equations are solved by the Runge-Kutta method, and the initial condition matches the numerical setup. At last, it should be known that the equations can be improved if we consider more details, such as the time dependentDi(t). However, according to our preliminary tests,these details could be ignored in the early-stage evolution.

1.3 The numerical setups and validations

The dimension of our computational domain is 2 000×2 000×1 000 (Fig. 2). Owing to the symmetry of the computation results, the symmetric boundary condition is imposed on the bottom face of the domain,while a fixed pressure is imposed on other faces by using the non-equilibrium extrapolation method[23].The applied pressure is set to be lower than the saturated vapor pressurePL, and higher than a critical one,Pcr=S/2R. In the computations, the dimensionless temperatureand kinematic viscosityν=0.1667(τ=1) are adopted,respectively. The surface tensionσis 1.38×10-2.

Fig. 2 A schematic diagram of computational domain

Two cases are simulated to validate our numerical method: a single bubble’s expansion(R0=20) and two identical bubbles withR10=R20=20,D0=50. As is known, given a proper interfacial tension, a critical bubble radius is obtained, above which the bubble will grow, and it will collapse otherwise[24]. As shown in Fig. 3, the results are consistent with the theoretical predictions. It indicates that the evolution of the bubble dynamics and the interaction between the bubbles can be well captured with the LBM.Meanwhile, Fig. 3 also shows that the interacting pressure suppresses the growth of the bubble. On the other hand, notingin the present work,the convergence of the calculation is also implied by the results, where the ratio of the interface thickness(ε～4) to the bubble radius varies. The parameter corresponds to a bubble (R=27μm) with ΔP=34 200 Pa in the water.

Fig. 3 Comparison between the LBM simulations and the theoretical predictions(Eqs.(8)and (10)). t is time in lattice unit

2. Numerical results and discussions

There are usually many bubbles (a bubble cluster)in the liquid phase in real cavitating flows. The interactions among the bubbles should surely be considered[7]. In this section, various bubble systems are simulated and analyzed.

2.1 Dynamics of two unequal bubbles in quiescent liquid

We first consider a dual bubble system, with the initial radii of. Both of them will grow, if they are set in the domain individually, as the case in Section 1.3.However, when the bubbles are set in the domain with the distanceD0=50, the small bubble collapses (see Fig. 4). Meanwhile the large bubble grows nonuniformly, and a tip pointing to the small one appears gradually. The phenomenon was also reported by Quinto-Su and Ohl[25]. Equation (9) (the right-hand side) shows that the neighboring bubble exerts additional pressure on the target one. The force depends on the inter-bubble distance and their growing speeds.

The evolution of the small bubble is shown in Fig.5 with variousD0. The results show that the simulations are consistent with the theory again. It is revealed that the small bubble grows similarly as the individual one as the inter-bubbles distance is large. With the decrease ofD0, the growth of the small bubble is suppressed drastically and it may collapse due to the inter-bubble actions. It is interesting to note that the small bubble undergoes a non-monotone growth asD0～100. It expands slightly at first. AfterDdecreases andD0increases a little, the additional pressure increases to compress the small bubble quickly.

Fig. 4 The evolution of two bubbles with R L 0=25, R S 0=20.Both of them would grow, if they are put in the domain individually.The subscripts S and L denote small and large bubbles, respectively

Fig. 5(a) The evolution of the small bubbles with different D0.The symbols (single bubble result) indicate the simulating results, and the curves are calculated from Eq. (10).

A semi-qualitative analysis of the suppression is conducted as follows. The evolution of the smaller bubble in the presented system is influenced not only by the ambient pressure, but also by the growth of the neighbouring one. The driven pressure on the smaller bubble is estimated as (see Eq. (10))

Fig. 5(b) Comparison of the small bubble

By combining it with the critical radiusR*, a dynamic critical radius,R*(t), is obtained. The bubble will expand when the instantaneous radiusRS(t) is greater thanR*(t), and otherwise it will collapse. Figure 5 shows that, forD0=50,R*(t) is always greater thanRS(t) and the smaller bubble collapses. ForD0=100,in a short period of time, which leads to a slight expansion.Nevertheless,R*(t) surpassesRS(t) aftert=300,and the bubble collapses eventually. ForD0=350,D(t) is so large that the expansion of the large bubble hardly affects the small one. Then, the small bubble keeps expanding. So, the behaviours of the bubbles can be qualitatively evaluated by the critical pressure[7].

2.2 Dynamics of multi-bubble systems

The presented dual-bubble analysis shows that the lattice Boltzmann method predicts the bubble evolutions correctly, and the mechanism of the interbubble action is briefly explored. More sophisticated cases are presented by using “multi-bubble” systems.Figure 6 shows the spatial arrangements of the multi-bubble systems. There is always a smaller bubble surrounded by some equal-sized larger ones in the systems. We aim to understand the influences of the following parameters: the bubble number, the topology of the system, the bubble size and the inter-bubbles distance, from the simulations. As shown in Fig. 6,D0is identical in the systems and

2.2.1Effects of bubble number and topology of the bubbles

The influences of the bubble number (NB) and the topology of the large bubbles are studied in this subsection. Figure 7 shows the collapsing time (t*,defined in the inset) of the small bubble and two features can be observed.t*decreases more quickly by increasingNBto the extent thatt*∝NB-0.3. On the other hand, for lowerNB, the topology of the bubbles impacts the shrink of the central small bubble.One can get a hint from the growth of equal-sized bubbles in Section 2.1. The growth of the large bubbles is influenced by their interaction, which is influenced by the distance between themselves. Such effects are shown in Fig. 8, where the growths of the large bubbles in two 3-bubble systems are compared.Smaller inter-large-bubble distance leads to a stronger interaction, and the large bubble grows slower.Consequently, the small bubble is compressed less.On the contrary, highNBsystems with different topologies have closer inter-large-bubble distances,and the influence of the topology is weak (see models 5a, 5b in Fig. 6).

Fig. 6 (Color online) Topologies of multi -bubble systems. The initial radius RL 0=25 for the larger bubble(light blue), and R S 0=20 (yellow) for the smaller one. The figure in the model name indicates the bubble number

Fig. 7 NB-dependence of collapsing time ( t*), which is defined in the inset. The circles show the numerical results,and the dashed line shows a slope of -0.3. D 0 =100 is adopted

Fig. 8 The large bubble growth in models 3a and 3b with D0=100. The inter-large-bubble distances are 200 and 100 for models 3a and 3b, respectively

The interacting pressures, exerted by the large bubbles on the small one, are presented in Fig. 9 for a further investigation. As shown in Fig. 8, the large bubble grows almost with a constant speed (t＞400).In fact, the bubble growth speed is estimated aswith some cavitating models[26-28](although not very precisely). Meanwhile, the estimation also leads to a collapsing time of. Hereinafter, we consider the impact of the interacting pressure on the collapsing time. It is convenient and reasonable to focus on the equal-sized large bubbles and neglect the influence of the small one. Eq. (10) is modified with the unknown functionRL

Fig. 9 The disturbed part of pressure exerted by a large bubble on the small bubble . DIda denotes the effective interbubble distance for large bubbles (Eq. (13)). The solid curve indicates the simplified relation pi ∝DIda/(DIda +RL ), and the dashed one is a fitting curve of pi ∝

2.2.2Effects of D0and RL0

We adjustD0andRL0in models 7 and 9 to explore the influences of the inter-bubble distance and the bubble size. The corresponding collapsing timesare shown in Fig. 10. Relationsandare obtained from the simulations. We can normally anticipate that for a largerD0or a smaller, a longer collapsing time is obtained according to Eqs. (9) and (11), due to a lowerpi.

Fig. 10 Collapsing time of the small bubble as functions. The results for models 7 and 9 are presented, with RS0=20

The presented analyses demonstrate that the bubble cloud can be regarded as a whole. The bubbles in the cloud grow with interactions. The cloud evolution should be considered when we focus on the evolution of a single bubble.

3. Concluding remarks

The evolution of multi-bubble systems is simulated by using a pseudopotential multiphase lattice Boltzmann method. The collapse of the central bubble is investigated particularly, although it will grow if it is put in an environment with the applied ambient pressure. The interaction among the bubbles is analyzed in the framework of the Rayleigh-Plesset theory. Empirical scaling relations are proposed associated with the numerical results. The geometrical parameters and the topology of the bubble system are studied. The concluding remarks are summarized as follows:

(1) With the LBM, the bubble dynamics can well be captured when the bubble is surrounded by many other bubbles. The additional pressure exerted by the other bubbles may compress the target bubble,although it is located in an environment with a low pressure.

(2) The evolution of the bubbles is well predicted by using the modified Rayleigh-Plesset equations, if they are equal-sized. The influences of the bubble size and the inter-bubble distance can be well predicted in the early stage. For the large bubble, we may haveand, respectively. That further leads toscalings, correspondingly.

(3) Given the number density, it is easy to determine the strength of the driven pressure on the central bubble, which includes the environmentaland interactingcomponents. The collapsing time for the bubble is estimated as

The present study implies a strong nonlinearity of the bubble cloud evolution. Based on the numerical and theoretical results, scaling laws are proposed,which implies a possibility to build a simple new cavitation model with the consideration of the bubble interaction. The studies of more general situations are desirable, such as the bubbles with distributed sizes.The study also lays solid fundamentals to explore an improved cavitation model with considerations of the bubble interactions.