Numerical simulation of bubble flow interactions

316

2009,21(3): 316-332 DOI: 10.1016/S1001-6058(08)60152-3

NUMERICAL SIMULATION OF BUBBLE FLOW INTERACTIONS* CHAHINE Georges L. DYNAFLOW INC., 10621-J Iron Bridge Road, Jessup, Maryland 20794, USA, E- mail: [email protected]

(Received October 18, 2008, Revised December 15, 2008)

Abstract: Bubble flow interaction can be important in many practical engineering applications. For instance, cavitation is a problem of interaction between nuclei and local pressure field variations including turbulent oscillations and large scale pressure variations. Various types of behaviours fundamentally depend on the relative sizes of the nuclei and the length scales of the pressure variations as well as the relative importance of bubble natural periods of oscillation and the characteristic time of the field pressure variations. Similarly, bubbles can significantly affect the performance of lifting devices or propulsors. We present here some fundamental numerical studies of bubble dynamics and deformation, then a practical method using a multi-bubble Surface Averaged Pressure (DF-Multi-SAP©) to simulate cavitation inception and scaling, and connect this with more precise 3-D simulations. This same method is then extended to the study of two-way coupling between a viscous compressible flow and a bubble population in the flow field. Key words: cavitation, bubble flow interaction, numerical simulation

1. Introduction  Study of cavitation inception teaches us that liquids rarely exist under a pure monophasic form and that bubble nuclei are omnipresent. These nuclei are very difficult to eliminate and are always in any industrial liquid at least in very dilute concentrations. However, most engineering studies, and rightfully so, ignore the presence of this very dilute often invisible bubble phase, and consider only the liquid phase to conduct analytical or numerical simulations. This applies to benign situations where pressure variations are not significant and where cavitation, dynamic effects, gas diffusion, and heat transfer do not result in dramatic modification of the flow to warrant inclusion in the modelling of both phases. In this contribution, we are concerned with those conditions where it is important to either evaluate whether cavitation may occur and/or to model the flow in presence of bubbles in significant local concentrations, sizes, or numbers to play a significant role in the flow dynamics. Under these conditions, interaction between the bubbles and 

*Biography: CHAHINE Georges L. (1947-), Male, Ph. D.

the flow can be significant and need to be properly addressed. Three families of situations can be distinguished: (1) Cavitation Inception in uniform or acoustic flow fields. (2) Cavitation Inception in vortical flow fields. (3) Developed cavitation and bubbly flows. In the first two cases, especially if cavitation inception is determined acoustically (i.e., prior to bubbles becoming visible), the interaction between bubbles and liquid flow is “1-way”, i.e., the liquid dynamics affects the bubble dynamics, while bubble effects on the hydrodynamic flow field are negligible. In the third case, bubbles presence and behaviour is influential enough to affect significantly the basic flow field and implementation and use of “2-way” interaction models is warranted. In this contribution, we will describe the work we have conducted at DYNAFLOW over the recent years in order to address the above aspects. This is not intended to be a review of the research field, however, references to other contributions will be done as appropriate. The subdivision of the article follows the

317

three families of situations listed above. 2. Cavitation inception Cavitation and bubble dynamics have been the subject of extensive research since the early works of Besant in 1859[1] and Rayleigh in 1917[2]. A host of papers and articles and several books[3-10] have been devoted to the subject. Various aspects of the bubble dynamics have been considered and included one or several of the physical phenomena at play such as inertia, interface dynamics, gas diffusion, heat transfer, bubble deformation, bubble-bubble interaction, electrical charge effects, magnetic field effects, etc.. Unfortunately, very little of the resulting knowledge has succeeded in crossing from the fundamental “research world” to the “applications world”, and it is uncommon to see bubble dynamics analysis made or bubble dynamics computations conducted for example by propeller or pump designers. This is due to the perceived impracticality of using the methods developed but for “experts”. However, with the tremendous advances in computing resources, these communities now commonly use Navier Stokes solvers and CFD codes to seek better solutions[11-15]. The challenge is thus presently for the cavitation/bubble community to bring its techniques to par with the single phase CFD progress. Consideration of bubble dynamics within a CFD computation should also become common practice. Indeed, the tools already developed are at the reach of all users and should be adopted as more CFD tools to use for advanced design[12,13]. 2.1 Presence of cavitation nuclei A traditional cavitation inception engineering definition is: “a liquid flow experiences cavitation if the local pressure drops locally below the liquid vapour pressure, Pv .” This definition of cavitation inception is only true in static conditions when the liquid is in contact with its vapour through the presence of a large free surface and is in practice applicable for liquids saturated with bubbles. For the more common condition of a liquid in a flow, or in a rotating machine, liquid vaporization can only occur through the presence of “microscopic free surfaces” or microbubbles, also called “cavitation nuclei”. These are micron and sub-micron sized bubbles in suspension in the liquid or trapped in particles crevices. Several techniques have been used to measure these nuclei distributions both in the ocean and in laboratory cavitation channels. These include Coulter counter, holography, light scattering methods, cavitation susceptibility meters, and acoustic bubble spectrometers[16-22]. Therefore, any fundamental analysis of cavitation inception has to start from the observation that, any

real liquid contains nuclei which when subjected to variations in the local ambient pressure will respond dynamically by oscillating and eventually growing explosively (i.e., cavitating). Cavitation inception cannot be defined accurately independent of this liquid bubble population (sometimes characterized by liquid “strength”[23]) and independent of the means of cavitation detection. The cavitation inception is in fact a complex dynamic interaction between the nuclei and their surrounding pressure and velocity fields, interaction that can be different between the small and large scale of the same configuration. In addition, the experimental means to detect and decide on “calling” cavitation inception (practical threshold used by the experimentalists) will affect the results and could be different between laboratory experiments and full scale tests. 2.2 Inception in quasi uniform flow: spherical model When the pressure variations to which a bubble is subjected are not slow compared to the bubble response time, the nuclei cannot instantaneously adapt, inertia effects become important, and one needs to consider bubble dynamics effects. This is the case for nuclei travelling through rotating machinery. The nuclei /bubbles then act as resonators excited by the flow field temporal and spatial variations. In the case of a vortical flow field the strong spatial pressure gradients (in addition to the temporal) strongly couple with the actual bubble motion (i.e., position vs. time) to result in a driving force that depends on the resonator reaction. This makes such a case, considered in the next section, much more complex than what occurs for a travelling bubble about a foil where, relatively speaking, the position of the bubble is less coupled to its dynamics. The flow field pressure fluctuations have various time scales: e.g., relatively long for travelling cavitation bubbles over a blade or captured in a vortical region flow, or very short for cavitation in turbulent strongly sheared flow regions. The amplitudes of these fluctuations and the relationship between the various characteristic times determine the potential for cavitation inception. For cavitation inception in uniform flow fields it is common to use spherical bubble models to examine the dynamics. This is justified as long as the local velocity and pressure gradients around the bubble are weak compared to the radial gradients, and as long as proper account is made of the bubble time dependent location in the flow field i.e., account for any slip velocity relative to the liquid flow. Spherical models have been extensively studied following the original works of Rayleigh[1] and Plesset[10]. For instance, if we limit the phenomena to be modelled to inertia, small compressibility of the liquid, compressibility of the bubble content, we obtain the Gilmore[24]

318

differential equation for the bubble radius R (t ) . We modified this equation to account for a slip velocity between the bubble and the host liquid, and for the non-uniform pressure field along the bubble surface[25]. (This is further discussed later.) The resulting Surface-Averaged Pressure (SAP) equation applied to Gilmore”s equation[25,26] becomes:

§ R ·  3 § R ·  1 § R R d · RR 1 + 1   ¨ ¸ ¨ ¸ R = ¨1+ + ¸< 2 © 3c ¹ U © c c dt ¹ © c¹

2

(1)

4

where c is the sound speed, P is the liquid viscosity, u is the liquid convection velocity and ub is the bubble travel velocity. Equation (5) degenerates to the classical equation for negligible Rayleigh-Plesset[10] compressibility and slip velocity effects. If in addition, gas diffusion and heat transfer effects are neglected and a polytropic law of gas compression is assumed, the resulting modified SAP equation becomes: 3k º 3 2 1 ª § R0 ·  RR + R = « pv + pg 0 ¨ ¸  Pencounter »  U ¬« 2 © R¹ »¼

1 § 2J 4 P R ·  u  ub + ¨ ¸+ 4 U© R R ¹

the gas diffusion problem and the assumption that the gas is an ideal gas[28]. The bubble trajectory is obtained using the following motion equation[29]

dub 3 3 = ’P + CD  u  ub u  ub + dt U 4 CL  u  ub u u  ub +

§ 2V R ·  4P ¸ + ¨ pv + pg  pencounter  R R¹ ©

 u  ub

leads to a more realistic bubble dynamics. In general, the gas pressure, p g , is obtained from the solution of

2

(2)

where k is the polytropic compression law constant. In the Surface-Averaged Pressure (SAP) bubble dynamics equation, we have accounted for a slip velocity between the bubble and the host liquid, and for a non-uniform pressure field along the bubble surface. In this SAP method the definition of Pencounter as the average of the liquid pressures over the bubble surface results in a major improvement over the classical spherical bubble model which uses the pressure at the bubble center in its absence[25-27]. This has serious implications for strongly non-uniform flows as discussed in the next section. For instance, a bubble does not always continuously grow once it is captured by a vortex. Instead, it is subjected to an increase in the average pressure once it grows and this

3  u  ub R R

(3)

where the drag coefficient C D is given by an empirical equation such as that of Haberman and Morton[30]:

CD =

Reb =

24 1.38 (1+ 0.197 Reb0.63 + 2.6 u 104 Reb ), Reb

2 U R u  ub

P

(4)

2.3 Inception in vortical flow fields: Axisymmetry Spherical bubble models, as briefly described above, can be efficient tools for studying cavitation inception, scaling, bubble entrainment, and cavitation noise. They can become more powerful if they are provided with further “intelligence” based on more precise non-spherical models which account for bubble behavior near boundaries, in pressure gradients, and in high shear regions, resulting in bubble deformation, elongation, splitting, coalescence, and non-spherical sound generation. One such refinement consists of considering the case of bubbles captured on a vortex axis. The bubble then elongates along the axis and may split into two or more sub-bubbles, and/or form jets on the axis. In order to investigate this behaviour our boundary element method axisymmetric bubble dynamics code 2DYNAFS©[31-38] was exercised and was able to simulate bubble dynamics through re-entrant jet formation, jet break through, and bubble splitting. The code can handle as input vortex flow fields obtained from CFD viscous computations or from experimental measurements. 2DYNAFS© simulations of bubble dynamics on a vortex axis under a large number of conditions lead to the followings conclusions illustrated in Fig.1[32-38]: (1) If the bubble is captured by the vortex far upstream from the minimum pressure, it remains spherical while oscillating at its natural frequency.

319

(2) When the bubble reaches the axis just upstream of the minimum pressure, it develops an axial jet on its downstream side which shoots through the bubble moving in the upstream direction. Even at this stage, the spherical model provides a very good approximation because the bubble is more or less spherical until a thin jet develops on the axis. (3) The bubble behaviour becomes highly non-spherical once it passes the minimum pressure location. It elongates significantly and can reach a length to radius ratio that can exceed 10. The bubble then splits into two or more daughter bubbles emitting a strong pressure spike followed later by other strong pressure signals when daughter bubbles collapse. Two axial jets originating from the split and a strong pressure signal during the formation of the jets are observed.

conclusion that has been preliminarily confirmed experimentally[32,36,38]. Three types of tests were conducted: spark generated bubbles, laser generated bubbles, and electrolysis bubbles injected in vortex lines. Figure 2 shows high speed photography and acoustic signals of bubble splitting between two rigid walls. A small but distinct pressure spike is formed at splitting followed but a more significant spike during the collapse of the sub-bubbles. The second set of experiments was conducted in a vortex tube, where bubbles generated by electrolysis were injected and observed once captured by the vortex line. Figure 3 shows the elongated bubble dynamics and the corresponding pressure signals[35]. The third set of experiments was conducted at the University of Michigan[36,38] using laser induced bubbles (in the vortex and far upstream) and the flow field of a tip vortex behind a foil. (Shown in Fig.4) Comparisons between the observations and the 2DYNAFS© simulations showed very good correspondence as shown in Fig.5. Detailed correspondence was harder to establish because of difficulty of measurements of all experimental microscale conditions.

Fig.1 Illustration of the acoustic pressure emitted by a bubble in a vortex field as a function of the cavitation number. Note that the bubble behaviour near and above the cavitation inception is quasi-spherical[32]

Fig.3 High speed photos of an electrolysis bubble captured in a line vortex, and resulting acoustic signal indicating peak signals at splitting and collapse[32]

Fig.2 High speed photos of spark-generated bubble collapsing between two solid walls, and resulting acoustic signal. Peak signals at splitting and sub-bubbles collapses[32] Fig.4

2.4 Experimental verification This behaviour supports the hypothesis that the noise at the inception of the vortex cavitation may originate from bubble splitting and/or from the jets formed after the splitting. This is an important

Bubble

behaviour

Rc = 4.51mm ,

in

a

vortex

* =0.2123 m2 /s ,

flow

field:

V =1.72 ,

U f =10m/s , R0 = 750Pm . 3-D view of the bubble just before splitting predicted by 2DYNAFS© and observed at University of Michigan using laser-induced bubbles at the center of the vortex[36]

320

2.6

Fig.5 Bubble behaviour in a vortex flow field[26]

2.5 Bubble splitting criteria The numerical simulations resulted in the following observations on the splitting of bubbles captured in a vortex (see Fig.6), which we use in our predictive models[31-35]: (1) An explosively growing bubble splits into two sub-bubbles after it reaches its maximum volume, (equivalent radius, Rmax ) and then drops to 0.95 Rmax . (2) The two resulting sub-bubbles have the following equivalent radii: 0.90 Rmax and 0.52 Rmax . (3) The locations of the two sub-bubbles after the splitting are at –0.95 Rmax and 4.18 Rmax on the vortex axis. (4) The pressure generated by the subsequent formation of reentrant jets in the sub-bubbles can be approximated by a function of V [36]. Since the noise associated with the jet formation appears to be much higher than the pressure signal from the collapse of a spherical bubble, it is desirable to include the splitting and the associated jet noise in simulations with multiple bubble nuclei.

Fig.6 Sketch of the successive phases of a bubble behavior in a vortex flow field

Inception in vortical flow fields: Fully non-spherical In order to study the full 3-D interaction between a bubble and a complex flow field, two methods were developed. The first, using our boundary element code, 3DYNAFS©, enables study of full bubble deformations during capture but neglects the effects that the bubble may have on the underlying flow field. The second method accounts for full 2-way bubble/flow field interaction, and considers viscous interaction. This model is embedded in an unsteady ReynoldsAveraged Navier-Stokes (RANS) code, DF-UNCLE 1, with appropriate free surface boundary conditions and a moving Chimera grid scheme[26]. This full 2-way interaction non-spherical bubble dynamics model has been successfully validated in simple cases by comparing the results with base case results obtained from the Rayleigh-Plesset equation and 3DYNAFS for bubble dynamics in an infinite medium both with and without gravity.

Fig.7

Comparison of the bubble radius versus time for the spherical models (the conventional and the SAP model) and the 3-D 1-way and 2-way UnRANS computations[26]

As an illustration, Fig.7 shows results of a bubble interacting with the tip vortex of an elliptical foil. The bubble elongates once it is captured, and depending on the cavitation number, either forms a re-entrant jet directed upstream or splits into two sub-bubbles. When 2-way interaction is taken into account, further smoothing of the bubble surface is exerted by viscosity resulting in a more distorted but overall more rounded bubble. Figure 8 illustrates the various stages of the interaction between a bubble and a tip vortex flow.

1

DF_UNCLE, new dubbed 3DYNAFS-VIS is based on UNCLE developed by Mississippi State University that we extended to 3-D cavity and free surface dynamics.

321

2.7 Validation of the SAP model Here we compare the SAP spherical model and the 3-D interaction non-spherical models to predict tip vortex cavitation inception for a tip vortex flow generated by a finite-span elliptic hydrofoil[26]. The flow field was obtained by a RANS computation, which provided the velocity and pressure fields for all compared models: the classical spherical model, the SAP model, the one-way interaction model where the bubble deformed and evolved in the vortex field but did not modify it, and finally the fully coupled 3-D model in which unsteady viscous computations included modification of the flow field by the presence of the bubbles.

considered to be distributed randomly in a fictitious supply volume feeding the inlet surface or release area of the computational domain. The fictitious volume size is determined by the sought physical duration of the simulation and the characteristic velocity in the release area as illustrated in Fig.9. The liquid considered has a known nuclei size density distribution function, n( R ) , which can be obtained from experimental measurements[17-22] and can be expressed as a discrete distribution of M selected nuclei sizes. Thus, the total void fraction, D , in the liquid can be obtained by M

D = ¦ Ni i =1

4SRi3 3

(5)

where N i is the discrete number of nuclei of radius Ri used in the computations. The position and thus timing of nuclei released in the flow field are obtained using random distribution functions, always ensuring that the local and overall void fraction satisfy the nuclei size distribution function.

Fig.8

Comparison of the maximum bubble radius vs. cavitation number between the spherical models (the conventional and the SAP model) and the 3-D 1-way and 2-way UnRANS computations[26]

Comparisons between the various models of the resulting bubble dynamics history, and of the cavitation inception values obtained from many tested conditions, reveal the following conclusions, illustrated in Fig.7 and Fig.8: (1) The bubble volume variations obtained from the full 2-way interaction model deviate significantly from the classical spherical model due to the interaction between the bubble and the vortex flow field. (2) Differences between the 1-way and 2-way interaction models exist but are not major. (3) Using the SAP scheme significantly improves the prediction of bubble volume variations and cavitation inception. SAP appears to offer a very good approximation of the full interaction model. 2.8 Modelling of a real liquid nuclei field We have developed the following general approach to address bubble flow interactions in practical conditions. In order to simulate the water conditions with a known distribution of nuclei of various sizes, as illustrated in Fig.9, the nuclei are

Fig.9

Illustration of the fictitious volume feeding the inlet to the nuclei tracking computational domain (or release area) and example of resulting nuclei size distribution satisfying a given nuclei density distribution function

2.9 Effect of a propeller on a nuclei field In order to simulate the effect of a propeller on a nuclei flow field, the methods described above were used. Gas diffusion into the bubbles was also

322

included, since it is very important in this case. For the nuclei field, the fictitious volume size was determined as the product of the release area, the sought physical duration of the simulation, and the characteristic inlet velocity as illustrated in Fig.10. For the results shown below we have selected characteristic values: size range of 10 to 200 μm and void fraction, D = 3.23 u 105 . Figure 11(a) shows the selected nuclei size number density and Fig.11(b) shows the numbers of discrete bubble sizes and band widths selected.

necessary and the effects on the results were analyzed. Figure 12 shows the time variation of the total number of bubbles within the computational domain for different number of repeats. The results clearly indicate that a steady state is achieved after t > 0.03ms with a minimum of 10 repeats effectuated. Otherwise, the two-phase medium is in a transient phase. Repeating the number of bubbles more than 10 times increases the period of quasi-steady state in the computational domain.

Fig.10

Fig.12

Effect of the bubble release duration on the total number of bubbles within the computational domain. Computations were conducted for 't = 3.7 u 10-3 s and results were then repeated to cover much longer release times

Fig.13

Front and side view of the bubbles in the computational domain at a given time during the computation. Bubble sizes are to scale in Section b, They are enhanced by a factor of 5 in Sections a and c

Illustration of the fictitious volume feeding the inlet to the nuclei tracking computational domain (or release area)

Fig.11 (a) The nuclei size distribution selected to resemble the field measurement by Medwin[19] and (b) the resultant number of nuclei released

In order to minimize CPU time, a time of release, 't , during which actual bubble computations are conducted, was selected and corresponded to a meaningful bubble population distribution in the fictitious bubble supply volume. The resulting bubble behavior was then assumed to be repeatable for the next 't . This was repeated as many times as

Figure 13 illustrates the sizes and locations of bubbles that became “visible”. Both front and side

323

views are shown at a given instant for ı = 1.75 . It is noted that for the side view the bubbles in Section b were illustrated with the actual sizes, while they were amplified 5 times to make them more visible in print in Sections a and c. Figure 13 illustrates the marked increase in “visible” bubbles downstream of the propeller, traveling bubble cavitation, tip vortex cavitation and a resemblance of sheet cavitation. In order to illustrate the results in a simplified manner, a time dependent void fraction by section, D ( x, t ) , is defined by integrating bubble volumes in between x and x + dx : N ( x ,t ) 4

¦

i =1

D ( x, t ) =



3

SRi3



1 S D 2  d 2 'x 4

(6) Fig.15 Comparison of the initial nuclei size distribution and that resulting at x = 0.2m

where d is the shaft diameter, and N ( x, t ) is the number of bubble within 'x at a given time, t . The time-averaged void fraction of D ( x, t ) at t1 , over a time period dt can also be defined as:

D=

1 t1+dt ³ D ( x, t )dt dt t

x = 0.2m . It is seen that a large number of bubbles of 400 μm radius are now present in the field, while all bubbles upstream were smaller than 200 μm. The effect of the cavitation number, the initial nuclei size distribution, and the initial dissolved gas concentration, on the bubble nuclei entrainment and dynamics were shown below.

(7)

1

Figure 14 shows the distribution of D along x . It is seen that the void fraction increases significantly as a result of cavitation on the blades (first hump) and in the tip vortices (second hump). The void fraction decreases as the encountered pressure by the bubble increases. However, the downstream void fraction is 1.4 times that of upstream.

Cavitation number: The time-averaged void fraction and downstream nuclei size distribution obtained at three different cavitation numbers, V = 1.5, 1.65 and 1.75, are compared in Fig.16 and Fig.17. In the three cases the void faction increases significantly near the propeller then plateaus. As the cavitation number decreases, the void fractions within the cavitation areas increase. The relative importance of the tip vortices diminishes relative to the propeller blade cavitation. With smaller values of V , the bubbles retain a much larger size downstream of the propeller due to enhanced net gain of gas mass diffusing into the bubbles. The downstream void fraction D at x = 0.2m is 50 times the upstream value for V = 1.5, 2.35 times for V = 1.65 and 1.4 times for V = 1.75 . The largest downstream bubble radii observed are 800 μm for V = 1.5 , 600 μm for V = 1.65 and 400 μm for V = 1.75 , as compared to 200 μm at release.

Fig.14 The time-averaged void fraction, D , distribution along the axial direction for D = 3.23 u 105 and V = 1.75

The effect of propeller cavitation on the number of “observable” bubbles is much more important than on D . Figure 15 shows the bubble size distribution at

Fig.16

Comparison of the void fraction variations along x for different cavitation numbers

Nuclei distribution:

324

The effect of the initial nuclei size distribution on the results is shown in Fig.18. A bubble size distribution ranging from 10-200 μm is compared to that of a uniform distribution of 100μm bubbles with the same initial void fraction. Higher bubble void fractions are seen near the propeller blade region for the initial uniform nuclei size case. However, this does not result in any significant difference in the void fraction downstream. Figure 19 compares the resulting downstream nuclei size distributions at x = 0.2m . It is seen that the largest bubble radius is 400 μm for the non-uniform nuclei distribution, while it is 300 μm for the uniform case.

in further increase in the void fraction downstream. However, although the dissolved gas concentration was doubled, the void fraction increased only 11%.

Fig.20 Comparison of the void fraction variation along the streamwise direction for two initial gas concentrations

3. Modeling of bubbly flows 3.1 Bubbly flow and propeller/hydrofoil In order to address problems where bubbles significantly modify the flow, we have developed a two-way coupled method between the viscous solver: 3DYNAFS-VIS© and multi-bubble dynamics solver, DF_MULTI_SAP©. The coupling between the two-phase continuum flow model and the bubble discrete dynamics/tracking model is an essential part of this method and can be described as follows: (1) The bubbles in the flow field are influenced by the local densities, velocities, pressures, and their gradients in the medium. The dynamics of the individual bubbles and their tracking are based on these local flow variables. The local properties of the mixture continuum are determined by the bubble presence. The void fractions, and accordingly the mixture local densities, are determined by the bubble population and size, and evolve as the bubbles move and oscillate. The flow field is adjusted according to the modified mixture density distribution while satisfying mass and momentum conservation. (2) The void faction is obtained by subdividing the computational domain in D -cell where D is computed and measuring the volumes of all bubbles in such cells. An D -cell is large enough for the statistics to be meaningful, i.e., it should contain a large enough number of bubbles, and small enough, compared to the flow field local characteristic length, such that bubbles of equal size in the cell will behave practically the same way. This approach was applied to a NACA4412 foil at zero angle of attack with various void fractions. Figure 21 shows the void fraction distribution for a hydrofoil operating in a bubbly flow with 5 % free stream void fraction. It is seen that the void fraction level decreases near the leading edge due to screening effects. High concentrations of bubbles are observed near the suction side surface, while the concentration P

Fig.17 Bubble size distribution at different cavitation numbers

x = 0.2m

for three

Fig.18 Comparison of the void fraction variation along the streamwise direction for two initial nuclei size distributions

Fig.19 Bubble size distribution at x = 0.2m for two initial nuclei size distributions

Dissolved gas concentration: The effect of the dissolved gas concentration on the bubble dynamics is illustrated in Fig.20 by doubling the dissolved gas concentration. This results

325

is very low in the wake region.

Fig.21 Void fraction distribution for a hydrofoil operating in a bubbly flow

Figure 22 shows the effect of the void fraction on the pressure coefficients on the foil surface. As D increases, the pressure on the suction side increases, which results in loss of lift compared to the single phase flow. The drag has contributions from both skin friction and pressure drag. Skin friction decreases when bubbles are present near the hydrofoil surface, the drag due to the pressure can increase. To illustrate this, we plot the pressure coefficients with respect to the vertical coordinate in Fig.23. The BC and CD portions of the hydrofoil surface contributes to drag, while the AB and AD portions contribute to thrust. The total pressure drag is related to the area between these curves, which increases as the void fraction increases.

Fig.22 Pressure coefficient vs. x curves on the surface of the hydrofoil for various void fractions

The changes of lift and drag from the single phase flow are shown in Fig.24. It is seen that the percentage of lift reduction is greater than the percentage of void fraction increase. The total drag resulting from both pressure drag and skin friction is shown to increase as the void fraction is increased because the increase of the drag due to the pressure is greater that the reduction of drag due to the skin friction.

Fig.23 Pressure coefficient vs. z curves on the surface of the hydrofoil for various void fractions

Fig.24 Hydrofoil performance, the lift and drag changes as functions of void fraction

3.2 Bubble augmented propulsion Several publications[39-46] and patents[47-49] claim that bubble injection in a propulsor or a waterjet can significantly augment the thrust even at high speeds. Experimental work has been conducted where fluid enters the “ramjet” and is compressed by passing through a diffuser (ram effect). High pressure gas is then injected into the fluid via mixing ports and constitute the energy source for the ramjet. The multiphase mixture is then accelerated by passing through a converging nozzle. Various prototypes have been developed and tested in the past and have shown a net thrust due to the injection of bubbles. However, the overall propulsion efficiency was typically less than that anticipated from mathematical models. These results require further modelling and controlled experimentation as we are presently undertaking at DYNAFLOW. First approximation of the multiphase flow field inside the air augmented nozzle were obtained by neglecting all non-uniformities in the cross-section and considering only a quasi-one dimensional flow field[50-54]. The fluid medium is considered as the mixture of water and bubbles which satisfies the

326

following continuity and momentum equations:

wU m + ’< U m um = 0 wt

Um

(8)

Dum 2 ª º = ’pm + ’ « 2 P mG ij  Pm  ’
where the subscript m represents the mixture medium. Here the G ij in the stress tensor term is the Kronecker delta, the mixture density, U m , and the mixture viscosity, P m , are functions of D :

Um = 1  D Ul + DU g

(10)

Pm = 1  D Pl + DP g

(11)

where the subscript l represents liquid, and the subscript g represents the gas in the bubble. The flow field has a variable density because the void fraction, D , varies in space and in time. This makes the overall flow field problem similar to a compressible flow problem. One approach we are developing is to use a compressible single phase flow solver and replace the equation of state which connects the density to the pressure with a solution of the bubble phase using the Lagrangian bubble tracking model described earlier, MULTISAP©, described earlier and which computes bubble size and location and thus D ( x , t ) . The 2-way coupling between the two is as follows: (1) The bubbles in the flow field are influenced by the local density, velocity, pressure, and the pressure gradient of the mixture medium. The dynamics of individual bubbles are based on these local flow variables as described in Eqs.(8) through (11). This is achieved using our Discrete Bubble Model (DBM). (2) The mixture flow field is influenced by the presence of the bubbles. The void fraction, and accordingly the mixture density, is modified by the migration and dynamics of the bubbles. The flow field is adjusted according to the modified mixture density distribution in a way that the continuity and momentum are conserved through Eqs.(8) and (9). The 2-way interaction described above is very strong for this problem since D can vary from zero near the inlet to as high as 70% at the nozzle exit.

Fig.25 The geometry of the nozzle inside wall as modeled. Looking from the downstream of the nozzle exit

3.3 Example study The tow pool experiment in Ref.[40] was chosen as a test case because there are measured data available. The nozzle is shown in Fig.25. The water comes through the inlet, flows through a diffuser, then an expansion with a bubble injector, and finally through the converging exit nozzle. The tests were performed at 8 m/s, and thrust, air flow, and pressure were measured at four locations inside the nozzle. The flow was modelled both in 3DYNAFS-VIS and Fluent. Boundary conditions were set so that the inlet has 8 m/s velocity and the exit has the pressure of 1 atm.

Fig.26

The axial velocity distribution on the center plane of the nozzle

Fig.27

The pressure distribution on the center plane of the nozzle

327

Figure 26 shows the velocity distribution on the center plane of the nozzle. As expected, the velocity inside the nozzle decreases as the cross sectional area increases up to the bubble injector, and then increases toward the exit as the cross section becomes gradually narrower. The corresponding pressure as shown in Fig.27 is the highest just downstream of the bubble injector. The inlet pressure, which is predicted to be 84 kPa, matches with the experimental observation as reported in Ref.[40]. 3.4 Steady state 1-D approach If one is interested in the steady state solution for a relatively simple axisymmetric nozzle, a steady 1-D approach can be used where the bubbles are assumed to behave the same regardless of their locations in the cross sections of the nozzle. One then obtains the void fraction from tracking bubbles through the axial pressure field of the nozzle and conducting an iteration procedure for the 2-way interaction between the bubbles and the flow. The procedure is as follows: (1) Track the bubble in the flow field. At the first step, the liquid flow is used as an initial guess. Later in the iteration, the most recently updated flow field is used. (2) From the tracking, the bubble radius is found as a function of axial coordinate along the nozzle. (3) Deduce the medium density from the void fraction using Eq.(10). (4) Solve the flow field using the updated medium density. (5) Re-compute bubble motion and size as in step 1 and iterate until convergence. The iteration described in the previous section converges well for the simulated nozzle. In Fig.28, the convergence of the bubble radius is shown indicating good convergence as early as in the 3rd iteration. Figure 29 shows a similar fast convergence of the corresponding flow field solution, i.e., the pressure distribution along the centerline of the nozzle.

Fig.28 The convergence of the tracked bubble radius

Fig.29

The convergence of the pressure along the centerline of the nozzle

Fig.30 The bubble radius as it flows inside the nozzle

Fig.31

The void fraction distribution inside the nozzle. The initial void fraction at the injection is 0.5 (the dashed line)

Fig.32

The converged pressure distribution inside the nozzle. The pressure along the centerline (top), and the pressure distribution on the center plane (bottom)

328

A converged behavior of a 5 mm radius bubble injected at a 2 atm pressure and 1.5 m downstream of the inlet is shown in Fig.30. Since the pressure inside the nozzle at the injection location is about 1 atm, the bubble grows after the injection then oscillates as it flows down and reaches an equilibrium radius. Figure 31 shows the converged void fraction distribution inside the nozzle. It is zero from the inlet to the injection location. Then it suddenly rises at the injection point and overshoots due to the bubble oscillations. Near the outlet the void fraction increases slightly because the bubble grows due to the pressure drop near the outlet. The pressure distributions inside the nozzle are shown in Fig.32and Fig.33. The axial velocity increases toward the outlet, and reaches as high as 15 m/s which is much higher than the exit velocity of the water only flow (Fig.26).

The thrust predicted from the steady 1-D model is shown in Fig.34. As the void fraction increases, there is a large gain in the momentum component of the thrust, which is countered by a slight loss in the pressure component. As a result, the total thrust increases as the void fraction increases. At the initial void fraction of 50%, the total thrust gain above the water only flow is about 1100 N.

Fig.34 Predicted thrust for various initial void fractions

Fig.33 The converged axial velocity distribution inside the nozzle. The axial velocity distribution on the center plane (top), and the axial velocity along the centerline (bottom)

The thrust of the nozzle can be computed by integrating the pressure and the momentum flux over a surface of a control volume that contains the nozzle.

T = ³³  pm + U mum2 dA

(12)

A

where u m is the axial component of the mixture velocity. For the inside of the nozzle the thrust is defined as follows:

T   po Ao  pi Ai +  U o Ao umo 2  Ui Ai umi 2

(13)

where the subscripts i and o represent the inlet and the outlet, respectively. The component of the thrust described in the first parenthesis in Eq.(13) is the contribution from the pressure difference between the inlet and the outlet. The second parenthesis represents the thrust due to the momentum change between the inlet and the outlet.

Fig.35

Comparison of the pressure field for the quasi-steady and Rayleigh-Plesset models for the quasi 1-D problem

3.5 Results from quasi 1-D models Both quasi-steady and Rayleigh-Plesset based models were applied to the same configuration. The operating conditions were determined such that the exit pressure was atmospheric and the inlet velocity was held constant at 8 m/s. The bubbles were injected near the end of the injection section of the nozzle. For the cases presented here, the injection was implemented as a discontinuity in the void fraction. The pressure was kept constant across it while the velocity was suddenly increased in order to preserve the mass flow rate. The pressures for a void fraction injection of 20% and 50% are compared in Fig.35. It is important to note that for a 20% void fraction injection, the two models provided nearly identical answers. However, at 50% void fraction injection, the pressure field exhibited a marked difference between the two models. It is important to note that the quasi-steady approach predicts larger bubbles at the

329

nozzle exit. This is expected since the inertial effects modelled by the Rayleigh-Plesset equations will retard bubble growth.

Fig.36 Comparison of the total thrusts obtained by the DBM+Fluent approach and two 1-D models

3.6 Comparison with two-way approach The predictions using the DBM and Fluent solver are compared to two other 1-D models described earlier. The total thrust is compared in Fig.36. They match well for low void fractions, but DBM deviates from the others for higher void fractions. The pressure distribution along the nozzle is compared in Fig.37. The pressure distribution obtained by DBM+Fluent approach and by the 1-D quasi-steady model matches well for D = 20% , but is different for =50%. In these computations, the outlet pressure is fixed at 1 atm. In the case of D = 50% , the difference of the pressure gradient is prominent in the exit nozzle. Figure 38 shows the void fraction along the nozzle. It is zero from the inlet to the injection (at 1.5 m). After the injection, increases mildly for the injection D = 20% , while it increases more rapidly for the D = 50% . The large difference between DBM+Fluent approach and 1-D quasi-steady model is observed near the outlet. This is due to the bubble dynamics that includes the slip of the bubble relative to the fluid.

Fig.38 Comparison of the void fraction along the length of the nozzle predicted by DBM+Fluent approach and 1-D quasi-steady model. The injection void fractions of 20% and 50% are shown

3.7 3-D unsteady coupling For full 3-D interactions, the void fraction is defined in the 3-D space using the D -cells described earlier. The bubbles and their volumes in each D -cell are counted to compute the void fraction. Improvement in the distribution of the void fraction by introducing a smoothing over the neighbouring bins was necessary to improve convergence of the solutions. The DF_UNCLETM and DF_MULTI_SAP© coupled approach can run a simulation in unsteady 2-way coupled way. The procedure for unsteady coupled run is as follows: (1) Initialize flow field, e.g., from water steady state solution. (2) Track the bubbles in the flow field following injection. (3) Based on the current bubble size and location, compute the void fraction D  x, y, z, t , using the P

D -cells. (4) Deduce the medium density U m  x, y, z , t

from D  x, y, z, t . (5) Solve the flow field using the updated medium density, and proceed to the next time step. (6) Repeat from step 2 to step 5 until the desired simulation time is reached. Figure 39 shows an example simulation using a full unsteady two-way coupling between DF_UNCLETM and DF_MULTI_SAP©. Three different size bubbles were injected on a bottom port near the inlet. It can be observed that the larger bubbles rise faster due to the buoyancy than smaller bubbles and the largest bubbles contribute predominantly to the medium density. Figure 40 shows an instantaneous snapshot of bubbles, medium velocity vectors, and density contours for a distributed bubble size injection at the widest section of the propulsor. P

Fig.37 Comparison of the pressure distribution along the length of the nozzle predicted by DBM+Fluent approach and 1-D quasi-steady model. The injection void fractions of 20% and 50% are shown

330

This model was finally coupled with a viscous and compressible flow solver to model the two phase flow in the case of bubble injection for thrust augmentation. Work is on-going to transform this into a generalized multi-scale two-phase flow model.

Fig.39 DF_Uncle© and DF_Multi_SAP© coupled simulation for a bubbles mixture injected from a bottom port

Acknowledgements The author acknowledges the sustained support of the Office of Naval Research, Dr. Kim Ki-Han monitor, and the contributions of many DYNAFLOW colleagues, most particularly Dr. Hsiao Chao-Tsung and Dr. Choi Jin-Keun who are responsible for most of the numerical simulations presented here. References [1] [2] [3]

Fig.40

Snapshot of computed bubble distribution and mixture medium density distribution on the center plain of the bubble augmented propulsion test setup. Only the portion of the nozzle between the injection and exit is shown

These solutions are being compared against an experimental fundamental study that we are conducting in a Plexiglas quasi-2-D cross section of the propulsor.

[4] [5] [6]

[7] [8]

4. Conclusions We have presented in this communication our efforts to model bubble flow interactions under various conditions and for several families of engineering applications. These included modelling of cavitation inception and development, estimation of the effect of air bubbles on the drag and lift of a lifting surface, and the modelling of purposely injected bubbles on the thrust of a propulsor or nozzle. The models developed were based on studies of single bubble dynamics and deformation in complex flow fields and used the acquired knowledge into a simplified Rayleigh-Plesset like model using bubble surface averaged quantities to compute the bubble behaviour in the flow field. This SAP model enables modelling in reasonable computation times of a real nuclei field distribution and is presently available as a User Defined Function (UDF) in Fluent.

[9] [10] [11]

[12]

[13] [14]

BESANT W. H. Hydrostatics and hydrodynamics[M]. London: Cambridge University Press, 1859, 158, RAYLEIGH Lord, On the pressure developed in a liquid during collapse of a spherical cavity[J]. Phil. Mag., 1917, 34: 94-98. KNAPP R. T., DAILY J. W. and HAMMITT F. G. Cavitation[M]. Yew York: McGraw Hill, 1970. HAMMITT F. G. Cavitation and multiphase flow phenomena[M]. Yew York: McGraw-Hill,1980. YOUNG F. R. Cavitation[M]. Yew York: McGraw Hill, 1989. FRANC J. P., AVELLAN F. and BELHADJI B. et al. La cavitation. m l canismes physiques et aspects industrielles, collection grenoble sciences[M]. Presses Universitaires de Grenoble, 1995(in France). BRENNEN C. E. Cavitation and bubble dynamics, Oxford engineering sciences series 44[M]. Oxford: Oxford University Press, 1995. ISAY W. H. Kavitation, schiffahrts-verlag[M]. Hamburg: Hansa C. Shroedter and co., 1981(in Germany). LEIGHTON T. G. The acoustic bubble[M]. Acad. Press, 1994. PLESSET M. S. “Bubble dynamics” in cavitation in real liquids[M]. Editor Robert Davies, Elsevier Publishing Company, 1964, 1-17. HSIAO C. T., PAULEY L. L. Numerical computation of the tip vortex flow generated by a marine propeller[J]. ASME Journal of Fluids Engineering, 1999, 121(3): 638-645. CHAHINE G. L. Nuclei effects on cavitation inception and noise[C]. 25th Symposium on Naval Hydrodynamics. St. John’s, Newfoundland and Labrador, Canada, 2004. CHOI J. K., CHAHINE G. L., Discrete bubble dynamics modeling[J]. Fluent News, ANSYS, Inc., 2006, 15(3). CHEN B., STERN F. Computational fluid dynamics of four-quadrant marine-propulsor flow[C]. ASME Symposium on Advances in Numerical Modelling of Aerodynamics and Hydrodynamics in Turbomachinery. Washington, D. C., USA, 1998.

331

[15]

[16]

[17] [18] [19]

[20]

[21]

[22]

[23] [24] [25]

[26]

[27]

[28]

[29]

GORSKI J. J. Present state of numerical ship hydrodynamics and validation experiments[J]. Journal of Offshore Mechanics and Arctic Engineering, 124(2): 74-80. MACINTYRE F. On reconciling optical and acoustical bubble spectra in the mixed layer, in Oceanic whitecaps[M]. New York: Monahan and Macniocaill, Reidell, 1986,75-94. OLDENZIEL D. M. A new instrument in cavitation research: The cavitation susceptibility meter[J]. Journal of Fluids Engineering, 1982,104: 136-142. BILLET M. L. Cavitation nuclei measurements – A review[C]. ASME Cavitation and Multiphase Flow Forum, 1985, FED-23: 31-38. BREITZ N., MEDWIN H. Instrumentation for in situ acoustical measurements of bubble spectra under breaking waves[J]. J. Acoust. Soc. Am., 1989, 86: 739-743. DURAISWAMI R., PRABHUKUMAR S. and CHAHINE G. L. Bubble counting using an inverse acoustic scattering method[J]. The Journal of the Acoustical Society of America, 1998, 104(5): 2699-2717. CHAHINE G. L., KALUMUCK K. M. and CHENG L. Y. et al. Validation of bubble distribution measurements of the ABS acoustic bubble spectrometer with high speed video photography[C]. 4th International Symposium on Cavitation. California, Pasadena, USA: Institute of Technology, 2001. CHAHINE G. L., KALUMUCK K. M. development of a near real-time instrument for nuclear measurement: The ABS acoustic bubble spectrometer[C]. Proceedings of the 4th ASME-JSME Joint Fluids Engineering Conference. Honolulu, Hawai, 2003 ARNDT R. E., KELLER A. P. Water quality effects on cavitation inception in a trailing vortex[J]. ASME Journal of Fluid Engineering, 1992, 114: 430-438. GILMORE F. R. The growth and collapse of a spherical bubble in a viscous compressible liquid[R]. California Institute of Technology, Hydro. Lab. Rep., 26-4, 1952. HSIAO C. T., CHAHINE G. L. and LIU H. L. Scaling Effects on prediction of cavitation inception in a line vortex flow[J]. Journal of Fluid Engineering, 2003, 125: 53-60. HSIAO C. T., CHAHINE G. L. Prediction of vortex cavitation inception using coupled spherical and non-spherical models and Navier-Stokes computations[J]. Journal of Marine Science and Technology, 2004, 8(3): 99-108. HSIAO C. T., CHAHINE G. L. Scaling of tip vortex cavitation inception noise with a statistic bubble dynamics model accounting for nuclei size distribution[J]. ASME Journal of Fluid Engineering, 2005, 127(1): DOI:10.1115/1.1852476. CHAHINE G. L., KALUMUCK K. M. The influence of gas diffusion on the growth of a bubble cloud[J]. ASME Cavitation and Multiphase Flow Forum, Cineinnati, Ohio, 1987, 50: 17-21. JOHNSON V. E., HSIEH T. The influence of the trajectories of gas nuclei on cavitation inception[C]. Sixth Symposium on Naval Hydrodynamics. 1966, 163-179.

[30] HABERMAN W. L., MORTON R. K. An experimental investigation of the drag and shape of air bubbles rising in various liquids[R]. Report 802, DTMB, 1953. [31] CHOI J. K., CHAHINE G. L. Non-spherical bubble behavior in vortex flow fields[J]. Computational Mechanics, 2003, 32(4-6): 281-290. [32] CHOI J. K., CHAHINE G. L. Noise due to extreme bubble deformation near inception of tip vortex cavitation[C]. FEDSM”03, International Symposium on Cavitation Inception, 4th ASME/JSME Joint Fluids Eng. Conference. Honolulu, Hawaii, 2003. [33] CHOI J. K., CHAHINE G. L. A numerical study on the bubble noise and the tip vortex cavitation inception[C]. Eighth International Conference on Numerical Ship Hydrodynamics. Busan, Korea, 2003. [34] CHOI J. K., CHAHINE G. L. and HSIAO C. T. Characteristics of bubble splitting in a tip vortex flow[C]. Fifth International Symposium on Cavitation, CAV2003. Osaka, Japan, 2003. [35] CHOI J. K., CHAHINE G. L. Noise due to extreme bubble deformation near inception of tip vortex cavitation[J]. Physics of Fluids, 2004, 16(7): 2411-2418. [36] REBOW M., CHOI J. and CHOI J. K. et al. Experimental validation of BEM code analysis of bubble splitting in a tip vortex flow[C]. 11th International Symposium on Flow Visualization. Notre Dame, Indiana, 2004. [37] CHOI J. K., HSIAO C. T. and CHAHINE G. L. Tip vortex cavitation inception study using the Surface Averaged Pressure (SAP) model combined with a bubble splitting model[C]. 25th Symposium on Naval Hydrodynamics. St. John’s, Newfoundland, Canada, 2004. [38] CHOI J., HSIAO C. T. and CHAHINE G. L. et al. Growth, oscillation and collapse of vortex cavitation bubbles[J]. Journal of Fluids Mechanics, 2009, 624: 255-279. [39] ALBAGLI D., GANY A. High speed bubbly nozzle flow with heat, mass, and momentum interactions[J]. International Journal of Heat and Mass Transfer, 2003, 46: 1993-2003. [40] MOR M., GANY A. Analysis of two-phase homogeneous bubbly flows including friction and mass addition[J]. Journal of Fluids Engineering Transaction of the ASME, 2004,126: 102-109. [41] MOTTARD E. J., SHOEMAKER C. J. Preliminary investigation of an underwater ramjet powered by compressed air[R]. NASA Technical Note D-991, 1961. [42] WITTE J. H. Predicted performance of large water ramjets[C]. AIAA 2nd Advanced Marine Vehicles and Propulsion Meeting. 1969, AIAA Paper No. 69-406. [43] PIERSON J. D. An Application of hydropneumatic propulsion to hydrofoil craft[J]. Journal of Aircraft, 1965, 2: 250-254. [44] KOREN O. Water tank testing of a laboratory two-phase marine ramjet engine[D]. Master Thesis, Technion Institute of Technology, 2005. [45] VALENCI S. Parametric sea trials of marine ramjet engine performance[D]. Master Thesis, Technion Institute of Technology, 2006.

332

[46] TANGREN R. F., DODGE C. H. and SEIFERT H. S. Compressibility effects in two-phase flow[J]. Journal of Applied Physics, 1949, 20(7): 637-645. [47] SCHELL Jr. et al. The hydro-pneumatic ram-jet[Z]. US Patent 3,171,379, 1965. [48] VARSHAY H., GANY A. Underwater two phase ramjet engine[Z]. US Patent 5,598,700, 1994. [49] VARSHAY H., GANY A. Underwater two phase ramjet engine[Z]. US Patent 5,692,371, 1996. [50] MUIR T. F., EICHHORN R. Compressible flow of an air-water mixture through a vertical two-dimensional converging-diverging nozzle[C]. Proc. Heat Transfer Fluid Mechanics Institute. Stanford University Press, 1963, 183-204. [51] ISHII R., UMEDA Y. and MURATA S. et al. Bubbly

flows through a converging-diverging nozzle[J]. Physics of Fluids, 1993, 5: 1630-1643. [52] KAMEDA M., MATSUMOTO Y. Structure of shock waves in a liquid containing gas bubbles[C]. IUTAM Symposium on Waves in Liquid/Gas and Liquid/Vapor Two Phase Systems. 1995, 117-126. [53] WANG Y. C., BRENNEN C. E. One-dimensional bubbly cavitating flows through a converging-diverging nozzle[J]. Journal of Fluids Engineering, 1998, 120: 166-170. [54] PRESTON A., COLONIUS T. and BRENNEN C. E. A numerical investigation of unsteady bubbly cavitating nozzle flows[C]. Proceedings of the ASME Fluid Engineering Division Summer Meeting. Boston, 2000.