Direct numerical simuiation of bubble-cluster's dynamic characteristics

689

2008,20(6):689-695

DIRECT NUMERICAL SIMUIATION OF BUBBLE-CLUSTER’S DYNAMIC CHARACTERISTICS* CHEN Hui School of Power and Energy, Northwestern Polytechnical University, Xi’an 710072, China, E-mail: [email protected] LI Sheng-cai, ZUO Zhi-gang, LI Shuai School of Engineering, University of Warwick, Coventry CV4 7AL, UK

(Received August 29, 2008, Revised October 8, 2008)

Abstract: A Direct Numerical Simulation (DNS) for understanding the dynamic response of bubble cluster to pulses of pressure perturbations has been studied by using a front-tracking method. The results show that owing to high nonlinearity, the bubble shape and volume oscillations caused by passing by pressure wave will be transformed into an in-phase volumetric oscillation of whole bubble cluster at a particular low-frequency. The value of the frequency is independent of the pulse excitations but the characteristics of the bubble cluster such as its bubble size, bulk void fraction and its spacial distribution etc. It is believed that this study provides important information for us to understand the coupling mechanism of cavitation cloud involved in cavitation resonance, a phenomenon noticed by one of the authors more than two decades ago. Key words: numerical simulation, cavitation, dynamics

1. Introduction  Cavitation induced pressure oscillations at low frequency is an important phenomenon that affects operation and life of hydraulic systems[1]. More than two decades ago, Li et al.[2-4] firstly noticed a phenomenon named as cavitation resonance on the UM Venturi. That is, a particular component of these pressure fluctuations would be manifested and significantly magnified with attributions similar to resonance, revealing as system instability. Later on, more studies were performed by the authors and other researchers. A macro-mechanism has been proposed that the cavitation resonance is induced through the coupling of the cavitation cluster (cloud) with the liquid phase of the flow system. However, the coupling mechanism is not fully understood[1]. Extensive researches have been done for understanding the dynamics of cavitation cluster 

* Project supported by the EPSRC’s Warwick-IMRC of United Kingdom (Grant Nos. R.ESCM.9001, R.ESCM.9219). Biography: CHEN Hui (1970-), Male, Ph. D. Candidate

(cloud). Experimental work can be traced back to the investigation by Campbell and Pitcher[5], who measured the propagation velocity of shock wave in bubbly liquids. Theoretically, Van Wijngaarden[6] derived an analytical model for pressure wave propagating through liquid/bubble mixture. Mørch[7] explored the collapse process in three different forms of bubble cluster as spherical, cylindrical and one-dimensional cavity cluster. Based on linearised Rayleigh-Plesset equation, d'Agostino and Brennen[8] studied the dynamics of spherical cavitation cloud. Similar approach was also employed by Omta’s investigation[9]. A vertical shock tube was used by Beylich and Gulhan[10] to study the wave structures in a bubbly liquid. Kameda et al. [11] found that the discrepancy between the experimental results and the numerical results based on the Wijngaarden’s model was attributed to the non-uniformity of bubble distribution. Several numerical methods[12,13] were successfully used to capture the interface between gas and liquid in some cases of two-phase flows. From experimental investigations on hydrofoils, Reisman et al. [14] observed two types of shock waves, having

690

global and local structure respectively, in collapsing cavitation clouds. Considering a fully nonlinear effect, Wang and Brennen[15] found that interaction parameter was crucial for the dynamics and noise of cavitation cloud. However, we are still not equipped with enough knowledge to understand the coupling mechanism, in particular it involves high void-fraction cloud. The first step towards the full understanding of the coupling mechanism is to understand how the cloud of high void-fraction respond to pressure perturbation at individual bubble level, i.e., at the microscopic level. The front-tracking method developed by Tryggvason et al. [16] has been adapted to simulate the dynamic response of a 3-D bubbly cluster. 2. Mathematical model The equations governing the unsteady motion of viscous, immiscible two-phase fluid for the whole domain are the Navier-Stokes equations. The conservative form of the motion equation is:

w  U u + ’< U uu = ’p + wt ’< P  ’u + ’T u + U g + Fı

outer liquid respectively. A fast Poisson solver is employed to solve the indicator function.

’ 2 I = ’
where G  x is the gradient of the indicator function at every stationary grid point.

G  x = ¦ D  x  x (l) n l 's  l

where n  l is the unit vector normal to an interface element l , whose surface area is represented by l 's  . D is a distribution function that determines what fraction of the quantity at every interface point should transfer to its adjacent stationary grid point. We employ an equation in our code which was introduced by Peskin[17]:





D x  x  l =  4h l

where u is the velocity vector, U , p and P are the density, pressure and viscosity respectively, g is the gravitational acceleration, Fı is the surface tension force which is zero everywhere except at the interface. The compressibility of liquid is neglected and its continuity equation is:

(6)

l

if xi  xi (1)

(5)



S

– «¬1+ cos 2h  x  x »¼ , ª

3

3

l

i

º

i

i =1

 2h, i = 1,3

(7a)



(7b)

D x  x  l = 0 , otherwise

where h is the size of fixed mesh. Therefore the surface tension force at each grid point can be expressed:





Fı  x = ¦ D x  x  FV  l

l

l

(8)

where FV  is surface tension force of an interface point. l

’< u = 0

(2)

The bubbles are compressible with pressure inside varies isothermally. An indicator function I  x is employed to indicate the location either in the bubbles ( I  x = 1 ) or in the outer liquid ( I  x = 0 ). The values of local density and viscosity at each grid are thus determined directly as:

U  x = U o +  Ub  U o I  x

(3)

FV  = V k  n 's  l

l

where the subscripts b and o indicate bubble

(4) and

l

(9)

where V is the surface tension coefficient, and k  l is the curvature of an interface element. Similarly, velocity at an interface point is interpolated from the velocities of its neighbour grid points.





u = ¦ D x  x  ui l

i

P  x = P o +  Pb  P o I  x

l

l

(10)

The new location of interface is therefore updated by integrating

691

 l

dx = u l dt

(11)

surface suddenly rises by 'p . Total 32 spherical bubbles of different sizes are initially assumed to be distributed with a void fraction from E = 11% to E = 27% (Fig.2). The gas pressure inside the bubble is assumed to obey the isothermal law.

3. Numerical methods In this code, all parameters such as bubble size, time, density, viscosity are dimensionless numbers. Among them are Ub / U o and Pb / Po . The others are reduced based on the Morton number, M = g Po4 / UoV 3 and the Eötvös number,

Eo = Uo gde2 / V . The whole computation domain is a bubbly liquid filling in a 1×1×4 3-D rectangular column. It is discretized by a stationary structured mesh of 34×34×128 grid points. A set of dynamic triangular mesh (Fig.1) is used to represent the front, through which information is transferred between the inside of bubble and the outer liquid. When the front deforms, some grid points will be added, deleted or redistributed to maintain appropriate geometry ratio and adequate resolution. In order to avoid numerical instabilities at the fronts owing to steep gradients of properties, the interface is set to have a small thickness in the order of the fixed mesh size.

Fig.2 Boundary conditions of the computational domain

Fig.1 Interface meshes

The boundary conditions on the transverse sections are taken to be periodic. The bottom of the domain is set as wall condition. The pressure applying on the top

Fig.3 The flow chart of present DNS code

The momentum equation is calculated by second-

692

order difference scheme for the spatial variables and second-order time-integration method for the time advection. The flow chart of present DNS code is shown in Fig.3. Firstly, the fronts are constructed according to the initial conditions. Once the indicator function has been determined, the density, viscosity and surface tension force can be updated accordingly. Then a predictor-corrector scheme is employed to calculate the velocity field. And the pressure equation is solved by a Successive Over Relaxation (SOR) method. Finally, the fronts are therefore reconstructed. This procedure will iterate until an accuracy criterion is reached. Fig.6 Velocity vector at different time

The pressure variations against time are obtained at different locations along the route of wave propagation. Based on these data, the mean speed of wave can be estimated. According to Fig.7, the mean speed of shock wave is about 14.6 m/s.

Fig.4 Unsteady pressure distribution on the diagonal section

4. Results and discussions Figure 4 shows the unsteady pressure distribution on the diagonal section at different dimensionless time. It is found that once the pressure at the top face is suddenly jumped by 'p , a pressure wave is generated and begins propagating through the bubbly liquid. In addition, waves with higher order of frequencies decay much more quickly. When this wave passes by, bubbles deform and move responding to the variations of local pressure and velocity, referring to Fig.5.

Fig.7 Fig.5 Deformation and movement of bubbles at different time

The velocity vectors at different time are also obtained, referring to Fig.6.

Local pressure fluctuations subject to different step impulses

Based on the Van Wijngaarden’s model[18] , the sound speed C0 is:

693

C02 =

J p0 dp = dU Ul E 1  E

(12)

For this case, i.e. the void fraction E = 26.8% , the sound speed C0 is about 16.8 m/s. Considering the fact that the strong nonlinearity of the bubble dynamics involved in our case, i.e., neither the bubble oscillations nor the exerting pressures are small enough for linear assumptions, the resultant wave speed should be slower than that of its counterpart linear case (i.e., the sound speed). The fact that the wave speed estimated from our DNS study is lower than sound speed actually indicates a good reflection of the dynamic nature of this bubble cluster in response to pressure excitations of finite values. The presence of bubbles significantly reduces the wave speed, making it much lower than the sound speed either in pure water (1480 m/s) or in pure air (332 m/s). Referring to Fig.7, the results also indicate that for the same void fraction, the wave speed remains unchanged, which is independent of the excitation characteristics, i.e., their amplitude ( 'p ) and duration ( 'T ). However, by examining Fig.7, it can be seen that after a few cycles of bubble oscillations, pressures at different locations along the wave propagation route begin to oscillate in phase at a particular frequency. The amplitudes of this particular oscillation depend on the input energy of excitation and they decay gradually owing to the viscosity of liquid. It is suggested that this is a collectively volumetric oscillation of individual bubbles, attributing to the global feature of the system, referring to Fig.8. Obviously the bubbly cluster of lower void fraction has a higher frequency and vice versa. The modes with higher frequencies decay more quickly.

ZB2 = 3J

For an isentropic gas bubble, the angular frequency of individual bubble ZB is[18]:

(13)

Neglecting surface tension, its natural frequency is

fB =

1 2SR0

3J p0

Ul

(14)

For our case (the non-dimensional radius of bubble R0 = 0.2 ), f B | 10000Hz . It is much higher than the oscillation frequency we observed from DNS results, referring to Fig.8. Therefore, it is obvious that this mode of low-frequency oscillation observed is not attributed to individual bubble oscillations but volumetric oscillation of the whole cluster.

Fig.9

Fig.8 Volume oscillations for different void fractions

p0 2V  2 Ul R0 Ul R03

Oscillations at different locations by two sinusoidal excitations

In order to examine the influence of impulse on this volumetric oscillation, we change the excitation from step-rise to two sinusoidal styles: 'p = 0.5sin Zt shown in Fig.9(a) and 'p = 0.5sin Zt shown in Fig.9(b). The durations of these two impulses are equal ( 'T = 1 ). Despite of different excitations, the resultant frequencies of the volumetric oscillations are all the same, also referring to Fig.7, whereas the phase angle of oscillation is subjected to the impulse.

694

Finally, the effect of spatial distribution of bubbles is examined. For the first simulation, 16 bubbles with radius R0 = 0.2 and other 16 bubbles with radius R0 = 0.15 are distributed alternately. For the second case, the bubbles of radius R0 = 0.2 are located in the upper domain and the others in the lower domain. And, for the third one, the upper and the lower are swapped from the second case. For all these simulations, the void fraction remains unchanged as E = 19.1% . The results are shown in Fig.10. It is clearly demonstrated that the frequency of induced volumetric oscillations not only depends on bulk void fraction, but also on spatial distribution of void.

become oscillating in phase at a particular low-frequency, leading to a purely volumetric oscillation of whole cluster. Its frequency is independent of excitations, but dominated by the characteristics of the cluster, such as its bulk void fraction and its spatial distribution. (4) Despite of the fact that it is still premature by using this study to explain the micro mechanism of cavitation resonance, this DNS does provide some important clues about the dynamic features of cavitation cluster (cloud) of high void fraction involved in cavitation resonance. (5) Therefore, we believe that this numerical study is a step towards a better understanding of the microscopic nature of cavitation resonance. Acknowledgments The authors are grateful to the EPSRC’S Watwick-IMRC for financial supported. The authors also wish to express appreciation to Professor Tryggvason G. of Worcester Polytechnic Institute and Professor Delale C. F. of Istanbul Technical University for providing their original code of front-tracking method. References

Fig.10 Volumetric oscillations for different bubble distributions but with the same bulk void fraction

From the viewpoint of energy transfer, the energy input from the perturbation is firstly fed into the oscillations in both modes, i.e., shape and volume oscillations, of individual bubbles through the passing by wave. The energy transfers from shape oscillation to volumetric oscillation, leading to an in-phase volumetric oscillation of whole cluster. The transfer from the strong shape oscillation of bubbles to the volumetric oscillation of whole cluster is a reflection of the nonlinearity of this phenomenon. Consequently its frequency spectrum reveals a tendency of shifting towards lower frequencies. 5. Conclusions (1) The present DNS using the front-tracking method is relatively effective and accurate to simulate the dynamic response of bubble clusters with high void fraction, which is closer to real flow situation. (2) A pressure impulse, applied on the top surface, would trigger a pressure wave propagating through the whole domain, owing to the non-linearity whose mean speed is lower than the sound speed in this bubbly fluid. (3) Through bubble interaction, local pressures

[1]

LI S. C. Cavitation of hydraulic machinery[M]. London: Imperial College Press, 2001. [2] Li S. C., ZHANG Y. J. and HAMMITT F. G. Investigation of low-frequency pressure fluctuation associated with cavitating venturi flow[R]. Report No. UMICH 014571-64-I, Ann Arbor, USA: University of Michigan, 1983. [3] LI S. C. Pressure fluctuations in cavitating draft-tube flows[C]. ASME Winter Annual Meeting. Anaheim, California, USA, 1992, FED136: 1-6. [4] LI S. C., ZUO Z. G. and LIU S. H. et al. Cavitation resonance[J]. ASME Journal of Fluid Engineering, 2008, 130(3): 31302(PP1-7). [5] CAMPBELL I. J. PITCHER A. S. Shock waves in a liquid containing gas bubbles[J]. Proc. R. Soc. London Ser. A, 1958, 243(1235): 534-545. [6] Van WIJNGAARDEN L. On the equations of motion for mixtures of liquid and gas bubbles[J]. J. Fluid Mech., 1968, 33: 465-474. [7] MɎRCH K. A. Energy considerations on the collapse of cavity cluster[J]. Appl. Sci. Res., 1982, 38(1): 313-321. [8] D'AGOSTINO L., BRENNEN C. E. On the acoustical dynamics of bubble clouds[C]. ASME Cavitation and Multiphase Flow Forum. Toronto, Canada, 1990, 2: 72-75. [9] OMTA R. Oscillations of a cloud of bubbles of small and not so small amplitude[J]. J. Acoust. Soc. Am. 1987, 82(3): 1018-1033. [10] BEYLICH A. E., GÜLHAN A. On the structure of nonlinear waves in liquids with gas bubbles[J]. Phys. Fluids A, 1990, 2(8): 1412-1428. [11] KAMEDA M., SHIMAURA N. and HIGASHINO F. et

695

al. Shock waves in a uniform bubbly flow[J]. Phys. Fluids, 1998, 10(10): 2661-2668. [12] WANG Zhi-dong, WANG De-guan. Free surface reconstruction with VOF methods[J]. Journal of Hydrodynamics, Ser. A, 2003, 18(1): 52-62 (in Chinese). [13] DONG Yu-hong. The nonlinear behavior of interface between two-phase shear flow with large density ratios[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(5): 587-592. [14] REISMAN G. E., WANG Yi-Chun and BRENNEN C. E. Observations of shock waves in cloud cavitation[J]. J. Fluid Mech., 1998, 355(1): 255-283.

[15]

WANG Yi-Chun, BRENNEN C. E. Numerical computation of shock wave in a spherical cloud of cavitation bubbles, ASME Journal of Fluid Engineering, 1999, 121(4): 872-880. [16] TRYGGVASON G., BUNNER B. and ESMAEELI A. et al. A front tracking method for the computing of multi-phase flow[J]. J. Comput. Phys., 2001, 169(2): 708-759. [17] PESKIN C. S. Numerical analysis of blood flow in the heart[J]. J. Comput. Phys., 1977, 25(3): 220-252. [18] LAUTERBORN W. Cavitation and inhomogeneities in underwater acoustics[M]. New York, USA: Springer-Verlag, 1980, 127-140.