NUMERICAL SIMULATION OF SPHERICAL, CYLINDRICAL AND AXIAL BUBBLE CLOUDS COLLAPSE*

MAHDI Miralam, EBRAHIMI Reza, SHAMS Mehrzad

Faculty of Mechanical Engineering, K.N.Toosi University of Technology, Tehran, Iran, E-mail: miralam_mahdi@yahoo.com

(Received August 28, 2011, Revised May 1, 2012)

NUMERICAL SIMULATION OF SPHERICAL, CYLINDRICAL AND AXIAL BUBBLE CLOUDS COLLAPSE*

MAHDI Miralam, EBRAHIMI Reza, SHAMS Mehrzad

Faculty of Mechanical Engineering, K.N.Toosi University of Technology, Tehran, Iran, E-mail: miralam_mahdi@yahoo.com

(Received August 28, 2011, Revised May 1, 2012)

The nonlinear dynamics of a spherical, cylindrical and axial cloud of cavitation bubbles were numerically simulated in order to learn more about the physical phenomena occurring in the cloud cavitation. The simulations employed the fully nonlinear continuum mixture equations coupled with the Gilmore equation for the dynamics of bubbles by considering the compressibility of liquid. A set of the Navier-Stokes equations was solved for the gas inside a spherical bubble, considering heat transfer through the gas inside the bubble and the liquid layer. The flow field around the cylindrical and axial cloud was obtained by solving the Navier-Stokes equations using a finite volume method and a dynamic layering mesh scheme. The calculated strength of shock wave in the liquid around the cloud was of the order of 1×106Pa and the propagation of this shock wave lasted for 10 μs. The conducted investigations illustrate that the shock wave propagates before the cloud has completely collapsed. A good agreement with experimental data was observed.

cloud cavitation, shock wave, spherical cloud, cylindrical cloud, axial cloud

Introduction

Cloud cavitation is a periodic phenomenon occurring in a body which involves formation and collapse of sheet cavities. Experimental studies of cloud cavitation have shown that noise, erosion and vibration are mainly caused by the shock wave emission due to cloud cavitation collapses[1,2].

The earliest analytical study of the dynamics of cloud bubbles were conducted by Van Wijngaarden[3]. Van Wijngaarden used the continuum mixture models in the study of the behavior of a collapsing layer of bubbly fluid next to a flat wall. Agostino and Brennen[4]analyzed the linearized dynamics of a spherical cloud of bubbles using a continuum mixture model coupled with the Rayleigh-Plesset equation for the dynamics of the bubbles. Morch[5,6]and Hanson et al.[7]studied the collapsing cloud based on the energy transfer theory. They speculated that the collapse of acloud proceeded layer by layer from the outside layer of the cavity cloud towards its inner part. Zhang et al.[8]used the idea proposed by Hanson to simulate cloud cavitation collapses in a Venturi tube. Wang and Brennen[9]investigated the fully nonlinear dynamics of spherical bubbles and found that a bubbly shock wave developed as a part of the nonlinear collapse of the bubble cloud. Mahdi et al.[10,11]demonstrated that both thermal conduction and liquid compressibility has great effects on bubble dynamics and shock wave strength due to bubble collapses, but radiation heat transfer effect is negligible. Wu and Lu[12,13]used the Homogenous Equilibrium Model (HEM) for simulating laminar cavitating flows around the 2-D hydrofoil and 3-D axis-symmetric body. The so-called homogeneous equilibrium model is based on a single fluid approach, in which the relative motion between liquid and vapor is neglected. Fu et al.[14]have calculated the cavitating flow around a revolution body with the HEM approach of void fraction transport equation. Reisman et al.[15]measured very large impulsive pressures on the suction surface of an oscillating hydrofoil experiencing cloud cavitations and demonstrated that these pressure pulses were associated with the propagation of the bubbly shock waves. Konno etal.[16]showed that cloud cavitation collapses in spherical, cylindrical and axial features. Duration of collapse was the order of 10 to 100 microseconds and the maximum pressure generated by this collapse related to the collapse type. In their experiment, the peaks of impulsive force did not meet the instant of final collapsing, but were often some 10 to 100 microseconds earlier than the final collapsing.

In the present paper, with considering liquid compressibility, heat transfer through the gas inside the bubble and the liquid layer, the nonlinear dynamics of spherical, cylindrical and axial cloud of cavitation bubbles were studied. For a spherical cloud, the time-dependent boundary condition at the surface of cloud can be obtained using the Bernoulli equation in the spherically symmetric incompressible liquid outside the cloud. A CFD program is developed to compute pressure field on the surface of cylindrical and axial clouds.

Fig.1 Schematic of physical model

1. Physical model

Bubble cloud consists of a limited number of cavitation bubbles. The bubble cloud may be spherical (Fig.1(a)), cylindrical or axial (Fig.1(b)). This cloud is surrounded by the liquid, and compared with the compressibility of the cloud, the compressibility of the liquid is neglected. First, η bubbles per unit liquid volume inside the cloud are uniformly distributed. The mass transfer between two phases of liquid and gas is neglected and no coalescence or breakup occurs in the bubbles. Therefore, the value of η is always constant and uniform within the bubble cloud.

Applying an order of magnitude analysis[17]indicates that, for the present flows, relative motion of two phases can be neglected. However, a homogeneous continuum bubbly mixture model can be used for analyzing the problem. The spherical bubbles are uniformly distributed with equal initial radius within the bubble cloud and an amount of air and vapor exists within the bubbles. Over time, the bubble radius R(r,t) changes as a function of time and radial coordinate r relative to the bubble cloud center. The radius of the spherical cloud A(t) is time dependent. B(t) and 2H(t) are the radius and height of the cylindrical (axial) cloud, respectively. When the cloud behaves cylindrically, H is constant and B(t) changes with time and, in the axial cloud, B is constant and H(t) changes with time. At the moment t=0, the cloud and liquid are in equilibrium. Over time, the far-field driving pressure of the liquid, P∞(t), changes. In this paper, the cloud response to P∞(t) is investigated.

2. Mathematical model

2.1 Single bubble

The Gilmore equation describes the behavior of a spherical bubble within a static, compressible and inviscid liquid subjected to a sinusoidal wave. As no gravity or other asymmetrical disturbing effects are considered, the equation describing the evolution of bubble radius as a function of time is[18]

The dots in Eq.(1) refer to the first and second order time derivations. C andbH are the speed of sound and the liquid enthalpy at the interface between the gas-filled bubble and liquid, respectively.

where P0is the ambient pressure and ρl= 1000kg/m3

is the liquid density. The constant coefficientsandare considered 3.5×108Pa and 6.25, respectively, based on the experimental NIST data for water pressure and density. The pressure at the interface in the liquid is expressed as

where μl=0.001Ns/m2and S=0.0725N/m are the surface tension and dynamic viscosity of the liquid, respectively.

The gas pressure inside the bubble,bP, can be obtained by solving the mass and momentum conservation equation for the gas inside the bubble[15]

where0bP is the gas pressure at the bubble center and r is radial coordinate from the bubble center. By decomposing the density into center- and radial-dependent parts

a set of solutions for the mass conservation equation may be obtained which includes

and

The constant a is related to the gas mass inside a bubble, m, by

where P0＇=P0+2S/R0, and Pb0and Tb0are pressure and temperature at the bubble center, respectively. A temperature profile of the gas inside the bubble can be obtained from the energy equation with the pressure profile given in Eq.(5)[19]

where T∞=300K and Tblare the ambient temperature and bubble liquid interface temperature, respectively, A*and B*are the coefficients in the temperature-dependent gas conductivity with a form like as kg=A*T+B*andη*=(R/δ)(kl/B*), where kl=0.4W/mK is liquid conductivity, A*=5.528× 10-5J/msK2and B*=1.165×10-2J/msK and δ is the thermal boundary layer thickness.

The temperature at the bubble wall can be obtained easily using Eq.(9)

At the bubble center, the state equation for an ideal gas is Pb0R3/Tb0=const . Therefore,

where γ is the specific heat ratio of gas.

The temperature distribution in the liquid layer adjacent to the bubble wall, which is important for determining the heat transfer through the bubble wall, is assumed to be quadratic, that is,

Such a second-order curve satisfies the following boundary conditions

The energy conservation equation for the liquid under the influence of bubble wall motion is expressed by

where αlis the thermal diffusivity, and the radial velocity of liquid due to bubble motion can be obtained from mass conservation for an incompressible liquid ?（urer）=0, that is,

with the temperature profile and boundary conditions given in Eqs.(12) and (13), and the velocity profile in the liquid (Eq.(15)), the above equation becomes[18]

2.2 Cloud bubble

The basic approach is the fully nonlinear solutions for cloud used by Wang and Brennen[9]for the spherical cloud dynamics. The continuity and momentum equations for the bubbly mixture (neglecting the contributions of gravity and viscosity) are

wheremP andmiu are the mixture pressure and velocity, respectively,mix is a spatial coordinate, R is bubble radius, t is time and η is the population of bubbles per unit volume of liquid.

It is assumed that the ratio of liquid to vapor density is sufficiently large so that the volume of evaporated or condensed liquid is negligible. It is also assumed that bubbles are neither created nor destroyed. Thus, η is constant and uniform in the cloud and it is related to void fraction, α, by where α is the total volume of bubbles per unit volume of the mixture.

3. Flow field around the bubble cloud

By solving the one-dimensional continuity and momentum equations in the incompressible liquid flow outside the spherical cloud, the boundary condition on the surface of the cloud is obtained analytically

A two-dimensional code is developed to solve the Navier-Stokes equations around the cylindrical and axial bubble clouds. The surrounding cloud flow is assumed to be incompressible, laminar and unsteady. The governing equations are

where r and z are the radial and axial coordinates, u and w are the liquid velocity in the radial and axial directions, andlP is the liquid pressure.

4. Boundary conditions and initial conditions

All bubbles inside the cloud are assumed to have the same initial size. Therefore, the following initial conditions should be applied:

At the cloud center because of geometric symmetry

In Eq.(20) and Fig.1(b), P∞(t ) is the changes of the far-field liquid pressure experienced by the bubble cloud in the cavitating flow passing around a body. A sinusoidal form is chosen for P∞(t )

For a bubbly cloud which passes with the velocity U∞around a body with the size L, tGis the time period for which the cloud experiences pressure changes. CPmin=(P-P0)/0.5ρlis the minimum pressure coefficient of the flow around the body. P0is the ambient pressure which can be calculated knowing the cavitation number, σ, gas saturated pressure, Pν, and using the following relation

Figure 2 shows the changes of pressure coefficient of the far-field liquid from the bubble cloud for CPmin=-0.7, σ=0.4 and tG=5ms in terms of time. Figure 3 demonstrates the changes of spherical bubble radius by R0=100μm subjected to this pressure change using the Gilmore equation. Moreover, the results related to using the Rayleigh-Plesset equation employed in Ref.[20] are also given in this figure. Changes of bubble radius are exactly the same for both modes from the first step of growth until just before the first bubble collapse, however, after the bubble collapse, the changes of bubble radius become completely different. For the situation simulated using the Gilmore equation, the bubble radius damps after the fast collapse around the initial radius, which is due to the damping effects caused by liquid compressibility.

5. Numerical method

5.1 Inside the bubble cloud

The solving method of bubble cloud governing equations is based on the Lagrangian coordinate system. By the Lagrangian integrating of the continuity and momentum equations (Eqs.(17) and (18)) based on the mixture velocity, the changes in the radius of bubbles and pressure distribution inside the bubble cloud are obtained.

Fig.2 Variation of far field driving pressure with time

Since the governing equations for tbe bubble cloud have been written in the Eulerian form, first, they should be rewritten in Lagrangian coordinates. In the Lagrangian coordinate system, the density of a mixture material element,0(,)rtρ, is related to its initial density through the following relation

In the above equation, J is the Jacobian of the coordinate transformation from Lagrangian to Eulerian and is defined in the following way

where the value of β for spherical, cylindrical and axial modes are 2, 1 and 0, respectively.

By substituting the value of J from Eq.(29) into Eq.(28) and integrating it, the value of0(,)rrt is given as

By differentiating the above relation, the mixture velocity0(,)urt is obtained as

The pressure distribution within the bubble cloud, Pm(r0,t), is obtained by integrating the momentum Eq.(18)

whereSr is A, B and H for spherical, cylindrical and axial cloud respectively.

These equations form the basis of the present method for solving the flow field inside the cloud. The values of R(r0,t+Δt), ?R(r0,t+Δt)/?t , T0(r0,t+ Δt), Tl(r0,t+Δt) and δ(r0,t+Δt) were calculated by the numerical solving of Eqs.(1), (9), (11) and (16) using the Runge-Kutta fourth-order method. Equations (30) and (31) were used for calculating r(r0,t+ Δt ) and um(r0,t+Δt), respectively. The pressure Pm(r0,t+Δt) was calculated using Eq.(32). These set of nonlinear equations is solved by using the Newton-Raphson method. The pressure on the cloud surface is calculated by using Eq.(20) for spherical cloud. Eq.(32) was coupled with the Navier-Stokes equation for axial or cylindrical cloud.

5.2 Outside the cylindrical and axial clouds

The finite volume method was applied to transform the partial differential equation to algebraic relations. Staggered grid and upwind scheme were used for the discretization of differential equations. A segregated solver with SIMPLER algorithm was employed as the velocity-pressure coupling. The discretized equations were solved using the Conjugated Gradient Solver (CGS)[21]. The dynamic layering mesh scheme was used for modeling the moving boundary.

6. Results and discussion

Bubble cloud including the spherical air bubbles with the initial radius R0=100μm flows with velocity U∞=10m/s in the 20oC water (ρl= 1000kg/m3, μl=0.001Ns/m2and S=0.0725N/m ) which experiences pressure changes P∞(t) is displayed in Fig.2. The initial void fraction within the cloud is α0=5%. The initial radius of the spherical cloud was considered A0=100R0. The radius and height of cylindrical and axial cloud is B0=H0=87R0to ensure equal bubble volume cloud.

Fig.4 Time history of dimensionless of bubble radius at six different positions in the spherical cloud. Parameters used are, Cpmin=-0.7, tG=5ms , σ=0.4, α0= 5%, R0=100μm and A0=100R0

Fig.5 Bubble size and pressure distributions in the spherical cloud. Parameters as in Fig.4

Fig.6 Bubble center temperature distributions in the spherical cloud. Parameters as in Fig.4

Fig.7 Time history of dimensionless of bubble size at six different positions in the cylindrical cloud for H0=B0= 87R0. Other parameters as in Fig.4.

Fig.8 Time history of dimensionless of bubble size at six different positions in the axial cloud for H0=B0=87R0. Other parameters as in Fig.4.

Fig.9 Time history of dimensionless of cloud size. Parameters as in Fig.4

Fig.10 Calculated cloud wall velocity for the clouds shown inFig.9

Fig.11 Time-dependent pressure coefficient at clouds wall shown in Fig.9

Fig.12 Pressure distribution in pure liquid around the spherical cloud as shown in Fig.4 at collapse point

Fig.13 Pressure field around the cylindrical cloud as shown in Fig.5 at different instants

Figure 4 shows dimensionless bubble radius inside the spherical cloud in terms of time for six nodes in the Lagrangian coordinates from the cloud surface r0=A0to near the cloud center r0=0.1A0. Similar to the single bubble, the radius of bubbles started to increase with the decrease in the liquid pressure around the cloud. The bubbles close to the cloud surface, experienced liquid pressure earlier, therefore, they grew with faster gradient. The farther bubble location from thecloudsurface,theslowerthegradientofthe bubble growth. The bubbles on the cloud surface started to collapse with pressure recovery. Outer bubbles, as a shield, protected the inner bubbles from the liquid’s pressure increase, thus, inner bubbles collapsed and reached a bigger size later. The collapse of bubbles on the cloud surface led to more severe collapse of the bubbles close to them. Therefore, the pressure of bubble cloud was higher at the location of the collapse of these bubbles. As shown in Fig.5, at the location of bubble collapse, a shock wave was formed and propagated toward the cloud center.

Figure 5(a) shows the pressure distribution and bubble radius within the bubble cloud at the moment t=8.121 ms. At this moment, the shock wave was located at the distance r =6.2×10-3m from the cloud center with strength of 4.5×105Pa. Over time, the wave got closer to the cloud center and its strength increased in a way that at t=8.2261ms, the wave reached r =3.9×10-3m. At this time, its strength increased to 5.7×106Pa (Fig.5(b)). Figure 6 demonstrates the temperature distribution within the bubble cloud at t =8.1218ms and t =8.2261ms . By shock wave passing through the cloud, heat wave was also formed.

Figure 7 shows the radius changes of bubbles within the cylindrical bubble cloud. The radius changing process of the bubbles within the cloud was similar to that of the spherical bubble cloud. For the cylindrical bubble cloud, the maximum radius of bubbles was approximately equal. Figure 8 illustrates a similar trend in the axial bubble cloud. The maximum radius of the bubbles got smaller by getting closer to the cloud center and, for the spherical cloud, and vice versa.

Fig.14 Pressure field around the axial cloud as shown in Fig.6 at different instants

The size changes of spherical, cylindrical and axial bubble clouds are given in Fig.9. The size changing process of the cloud was approximately the same for all the modes. Also, the size changing of the bubble cloud was similar to the size changing of the single bubble. By decreasing the liquid pressure, the cloud grew and reached its maximum value. By increasing the pressure for the second time, the cloud collapsed, oscillated around its initial size and damped. In similar condition the maximum radius of bubble cloud occurs approximately at one moment and the size of axial bubble cloud was more than that of two others. The collapse of cylindrical and spherical cloud happened approximately at the same time, however, it happened later for the axial one.

Velocity changing process of the bubble cloud wall is given in terms of time in Fig.10. Cloud wall velocity increased at first and reached its maximum value. The maximum velocity of the axial cloud was higher than that of two others. When the cloud collapsed for the first time, the size of its wall velocity increased and reached its maximum value. Then, the velocity decreased with a fast gradient. In this situation, the amount of pressure on the cloud surface suddenly increased and led to the shock wave propagation in the liquid around the cloud.

Figure 11 shows pressure coefficient on the cloud surface at different times. At the moment of cloud collapse, sudden pressure increase happened on the cloud surface. The pressure coefficient in the axial cloud was about 920 and, for spherical and cylindrical clouds, it was 770 and 550, respectively. Pressure distribution in the flow field around the spherical cloud at the moment of collapse is given in Fig.12. The pressure on the cloud surface was 2.75×107Pa (equal to the pressure coefficient of 550) and decreased to 2×106Pa at the 0.15 m distance from the cloud center, which is in agreement with the experimental data given in Ref.[22].

Figure 13 demonstrates pressure distribution and stream lines around cylindrical cloud at different times. Cloud wall velocity and its position relative to the cloud center were determined in each figure. In Figs.13(a)-13(d), the cylindrical cloud grows and itsradius increases to 1.6×10-3m. Cloud collapse started approximately from 5.53 ms. After this time, the cloud wall velocity became negative and the cloud size started to decrease (Figs.13(e)-13(f)). At 8.31 ms, pressure coefficient in cloud wall reached 920 and a shock wave propagated form the cloud. At this moment, the cloud size did not reach its minimum value yet and the cloud wall velocity was -2.8 m/s and its size was decreasing. Konno et al.[16]confirmed this issue by conducting experimental investigations on the collapse of cloud cavitation on the hydrofoil. Bubble cloud reached its minimum size at 8.322 ms (Fig.13(i)) and, after that, its size became larger for the second time (Figs.13(j) and 13(k)). The processes of pressure and velocity changes around the axial cloud were similar to those of the cylindrical cloud. A similar physics could be observed for axial cloud which is presented in Fig.14.

[1] DULAR M. Relationship between cavitation structures and cavitation damage[J]. Wear, 2004, 257: 1176-1184.

[2] WANG Y. C., BRENNEN C. E. Shock wave development in the collapse of a cloud of bubbles[C]. Cavitation and Multiphase Flow. New York: American Society of Mechanical Engineers, 1994, 194: 15-19.

[3] Van WJINGARDEN L. On the collective collapse of a large number of gas bubbles in water[C]. Proceeding 11th International Congress on Applied Mechanics. Berlin: Spring-verlag, 1964.

[4] AGOSTINO L., BRENNEN C. E. Linearized dynamic of spherical bubble clouds[J]. Journal of Fluid Mechanics, 1989, 199: 2126-2134.

[5] MORCH K. A. On the collapse of cavity cluster in flow cavitation[C]. Proceeding First International Conference on Cavitation and Inhomogenieties in Underwater Acoustic. Springer, 1980, 4: 95-100.

[6] MORCH K. A. Energy considerations on the collapse of cavity cluster[J]. Applied Scientific Research, 1982, 38: 313.

[7] HANSON I., KEDRINSKII V. K. and MORCH K. A. On the dynamics of cavity cluster[J]. Journal of Physical D: Applied Physics, 1981, 15(9): 1725-1734.

[8] ZHANG Xiaodong, FU Yong and LI Zhiyi et al. The numerical simulation of collapse pressure and boundary of the cavity cloud in venturi[J]. Chinese Journal of Chemical Engineering, 2009, 17(6): 896-903.

[9] WANG Y. C., BRENNEN C. E. Numerical computation of shock waves in spherical cloud of cavitation bubble[J]. Journal of Fluids Engineering, 1999, 121(4): 872-880.

[10] MAHDI M., SHAMS M. and EBRAHIMI R. Effects of heat transfer on the strength of shock waves emitted upon spherical bubble collapse[J]. International Journal of Numerical Methods for Heat and Fluid Flow, 2010, 20(4): 372-391.

[11] MAHDI M., EBRAHIMI R. and SHAMS M. Numerical analysis of the effects of radiation heat transfer and ionization energy loss on the cavitation bubble’s dynamics[J]. Physics Letters A, 2011, 375(24): 2348-2361.

[12] WU Lei, LU Chuan-jing and XUE Lei-ping. An approach in modeling two-dimensional partially cavitating flow[J]. Journal of Hydrodynamics, Ser. B, 2002, 14(1): 45-51.

[13] WU Lei, LU Chuan-jing and LI Jie et al. Numerical simulations of 2D periodic unsteady cavitating flows[J]. Journal of Hydrodynamics, Ser. B, 2006, 18(3): 341-344.

[14] FU Hui-ping, LU Chuan-jing and WU Lei. Research on characteristics of flow around cavitating body of revolution[J]. Journal of Hydrodynamics, Ser. A, 2005, 20(1): 84-89(in Chinese).

[15] REISMAN G. E., WANG Y. and BRENNEN C. E. Observation of shock wave in cloud cavitation[J]. Journal of Fluid Mechanics, 1998, 355: 255-283.

[16] KONNO A., YAMAGUCHI H. and KATO H. et al. On the collapsing behavior of cavitation bubble clusters[C]. Proceedings Fourth of the International Symposium On Cavitation. Pasadena, USA, 2001.

[17] BRENNEN C. E. Fundamentals of multiphase flow[M]. New York: Cambridge University Press, 2005.

[18] GILMOR F. R. The growth of collapse of a spherical bubble in a compressible liquid[R]. Hydrodynamics Laboratory Report 26-4, Pasadena, USA: California Institute of Technology, 1952.

[19] KIM K. Y., BYUN K.-T. and KWAK H.-Y. Temperature and pressure fields due to collapsing bubble under ultrasound[J]. Chemical Engineering Journal, 2007, 132(1-3): 125-135.

[20] WANG Y. C. Shock waves in bubbly cavitating flows[D]. Ph. D. Thesis, Pasadena, USA: California Institute of Technology, 1996.

[21] KIM C. J., RO S. T. Efficient and robust matrix solver for the pressure correction equations in two-and threedimensional fluid flow problems[J]. Journal of Numerical Heat Transfer, 1995, 27(4): 355-369.

[22] REISMAN G. E., BRENNEN C. E. Shock waves measurements in cloud cavitation[C]. 21st International Symposium on Shock Waves. Great Keppel Island, Australia, 1997.

10.1016/S1001-6058(11)60279-5

* Biography: MAHDI Miralam (1980-), Male, Ph. D., Assistant Professor

SHAMS Mehrzad,

E-mail: shams@kntu.ac.ir