Vortex Simulation of the Bubbly Flow around a Hydrofoil

This study is concerned with the two-dimensional simulation for an air-water bubbly flow around a hydrofoil. The vortex method, proposed by the authors for gas-liquid two-phase free turbulent flow in a prior paper, is applied for the simulation. The liquid vorticity field is discrerized by vortex elements, and the behavior of vortex element and the bubble motion are simultaneously computed by the Lagrangian approach. The effect of bubble motion on the liquid flow is taken into account through the change in the strength of vortex element. The bubbly flow around a hydrofoil of NACA4412 with a chord length 100 mm is simulated. The Reynolds number is 2.5 × 10, the bubble diameter is 1 mm, and the volumetric flow ratio of bubble to whole fluid is 0.048. It is confirmed that the simulated distributions of air volume fraction and pressure agree well with the trend of the measurement and that the effect of angle of attack on the flow is favorably analyzed. These results demonstrate that the vortex method is applicable to the bubbly flow analysis around a hydrofoil.


INTRODUCTION
Turbo pumps are originally designed for handling singlephase liquid flows. But they are frequently operated under gas-liquid two-phase flow conditions because of their diverse use. As the pump performance for two-phase flow generally deteriorates, it is practically important to grasp and improve the performance. A number of researchers [1][2][3] have investigated the performance of pumps at various specific speeds, and they made clear the relation between the performance and the flow in impeller. Some numerical simulations on two-phase flows in impeller have been carried out so as to predict the effect of the pump design parameters, such as the impeller geometry and the rotating speed, on the pump performance [4,5]. With the intention of clarifying the twophase flow in impeller from a basic standpoint, the bubbly flow around a hydrofoil has also been studied. Ohashi et al. [6] and Ichikawa et al. [7] found from their measurement that the gas volume fraction is distributed quite unevenly around a hydrofoil. Nishikawa et al. [8] performed the numerical simulation for a bubbly flow around a hydrofoil to examine the slip motion of bubble relative to the water flow. Uchiyama [9,10] analyzed the bubbly flow around a hydrofoil undergoing heaving and pitching motions to predict the unsteady characteristics of the hydrofoil.
The authors [11] proposed a vortex method for gasliquid bubbly flow in a prior paper. The liquid vorticity field is discretized by vortex elements, and the behavior of vortex element and the bubble motion are simultaneously computed by the Lagrangian approach. The method can simulate directly the development of vortical structure, such as the formation and deformation of vortices. Therefore, it is usefully employed for the simulation of free turbulent flows, in which the large scale eddies play a dominant role. The authors applied the vortex method to simulate a bubbly flow around rectangular cylinder [11], a bubble-laden plane mixing layer [12], and a bubble-driven plume [13]. The simulations demonstrated that the interaction between the largescale eddies and the bubbles is successfully analyzed and that the simulated two-phase flow field is favorably compared with the corresponding experimental results. Since the two-phase flow around a hydrofoil is chiefly governed by the liquid vortices and the bubble motion, the authors' vortex method seems to be effectively employed for the simulation.
In this study, the vortex method is applied to the twodimensional simulation for an air-water bubbly flow around a hydrofoil. The flow was experimentally investigated by Ichikawa et al. [7]. A hydrofoil of NACA4412 is mounted in a channel. The chord length is 100 mm, and the Reynolds number is 2.5 × 10 5 . The angle of attack is set at 0 • and 10 • .

International Journal of Rotating Machinery
The bubble diameter is 1 mm, while the volumetric flow ratio of bubble to whole fluid is 0.048. It is confirmed that the simulated distributions of air volume fraction and pressure agree well with the trend of the measurement and that the effect of angle of attack on the flow is favorably analyzed. These demonstrates that the vortex method is applicable to the bubbly flow analysis around a hydrofoil.

Assumptions
The following assumptions are employed for the simulation.
(1) The mixture is a bubbly flow entraining small bubbles.
(2) Both phases are incompressible and no phase changes occur. (3) The mass and momentum of the gas phase are very small and negligible compared with those of the liquid phase. (4) The bubbles maintain their spherical shape, and neither fragmentation nor coalescence occurs.

Conservation equations for bubbly flow
The conservation equations for the mass and momentum of the bubbly flow are expressed by the following equations under the assumptions (1)- (3): where the third term on the right-hand side of (2) expresses the buoyancy effect, and the volume fraction satisfies the following relation:

Equation of motion of bubble
It is postulated that the virtual mass force, the drag force, the gravitational force, and the lift force act on the bubble. The equation of motion for the bubble is expressed by the following equation under the assumption (4): where d is the bubble diameter, u r = u g − u l and β = ρ g /ρ l . C V , C D , and C L are the virtual mass coefficient, the drag coefficient and the lift coefficient, respectively. C D is given as [14], where Re b = d|u r |/v l .

Decomposition of liquid velocity field
According to the Helmholtz' theorem, any vector field can be represented as the summation of the gradient for a scalar potential φ and the rotation for a vector potential ψ. Thus, the liquid velocity u l is written as The velocity calculated from (6) remains unaltered even when any gradient of scalar function is added to ψ. To remove this arbitrariness, a solenoidal condition is imposed on ψ: When substituting (6) into (1) and rewriting the resultant equation by using a relation ∇·(∇ × ψ) = 0, the following equation is obtained: Taking the rotation of (6) and substituting (7) into the resultant equation, the vector Poisson equation is derived: where ω is the vorticity for the liquid phase.

Lagrange analysis for vortex element and bubble
When taking the rotation of (2), the vorticity equation for the bubbly flow is derived. For the two-dimensional calculation, it is expressed by the following equation: The vorticity field is discretized by vortex elements. It is postulated that the vortex element has a viscous core and that the vorticity distribution around the element is represented with the Gaussian curve. When the vortex element γ at x γ has a circulation Γ γ and a core radius σ γ , the vorticity at x induced by the element is expressed as The convection of the vortex element γ is estimated by the Lagrangian calculation of the following equation: Equations (8) and (9) are solved by a finite element method. The two-dimensional computational domain in the x-y plane is resolved into quadrilateral elements. An element is shown in Figure 1. The scalar potential φ, the z-component for vector potential ψ z , the vorticity ω, and the liquid velocity u l are defined on the nodes, while the volume fraction is defined at the center of element. When solving (9) for ψ z , the T. Uchiyama and T. Degawa ω value on the nodes is determined by taking the summation of the vorticity produced by every vortex element. The α l value in (8) is estimated from (3) after computing the gas volume fraction α g from the Lagrangian calculation of (4). When the bubbles are assumed to distribute uniformly with an interval Δz in the z-direction normal to the x-y plane, α g for a quadrilateral element with an area A is expressed as where n g is the number of bubbles in the element. It is found from (10) that the vorticity varies with the lapse of time due to the viscous diffusion and the gradient of the volume fraction. These variations are considered through the changes in σ and Γ for the vortex element, and they are computed simultaneously with the Lagrangian calculation of (12).

Change in core radius due to viscous diffusion
The viscous diffusion is simulated by applying a core spreading method, which was used for single-phase flow simulations [15],

Change in circulation due to bubble motion
When employing the Reynolds transport theorem and (10), the time rate of change in the circulation Γ along any closed curve (area element dS) is expressed as where the viscous diffusion term is omitted because it is already considered by (14). The application of (15) to a quadrilateral element (Figure 1) yields the time rate of change in Γ, ΔΓ/Δt, in the element. In the case where the number of vortex elements in a quadrilateral element is n v , the change in the circulation for each vortex element during Δt is supposed to be ΔΓ/n v . In the case where there are no vortex elements in the quadrilateral element, a vortex element with a circulation ΔΓ is generated at the center of the quadrilateral element.

Numerical procedure
When the flow field at t = t is known, the flow at t = t + Δt is estimated by the following procedures.

SIMULATION CONDITIONS
The vortex method is applied to the two-dimensional simulation for an air-water bubbly flow around a hydrofoil. The hydrofoil of NACA4412 is mounted in a channel with 200 mm × 100 mm cross-sectional area. The chord length l is 100 mm. The Reynolds number, based on the water velocity upstream of the hydrofoil u l0 and l, is 2.5 × 10 5 . The volumetric flow ratio of bubble to whole fluid α g0 is 0.048. Figure 2 shows the computational domain. The inlet and outlet boundaries, B 0 and B 3 , are located: l upstream and 3l downstream of the leading edge of the hydrofoil. The simulations are performed at the angles of attack of 0 • and 10 • .
The vorticity field is generated by the bubble motion and the velocity shear layer originating from the hydrofoil. To simulate the vorticity field due to the velocity shear layer, the vorticity layer on the hydrofoil is divided into segments with a length s, and a vortex element is released from each segment into the flow field at a time interval Δt. In this simulation, the number of segments is 50, the core radius σ 0 of the vortex element at the release is 0.5s, and the height of vorticity layer is set at 10/ √ Re. The σ 0 value for a vortex element generated from a quadrilateral element by the bubble motion is determined so that πσ 2 0 coincides with the area for the quadrilateral element.
The bubbles are released into the water flow field with a velocity u g0 from the inlet boundary B 0 at a time interval Δt. The diameter is set at the measured value 1 mm, and u g0 is estimated from a drift-flux model [16]. To reproduce the bubbly flow condition, the bubble releasing positions are  determined by using random numbers. The second-order Adams-Bashforth method is used for the Lagrangian calculation of vortex element and bubble. The time increment Δt is set at u l0 Δt/s = 0.025.
Equations (8) and (9) are solved by using the following boundary conditions: the uniform flow on the inlet boundary B 0 ; the convective outflow on the outlet boundary B 3 ; the slip on the channel walls B 1 and B 2 , and the nonslip on the hydrofoil surface B 4 . The nonslip condition is also applicable to B 1 and B 2 by representing the vorticity layer with the same method applied on the hydrofoil surface. Since the Reynolds number is so high, this study assumes that the vorticity layer scarcely affects the flow around the hydrofoil and imposes the slip condition.
The above-mentioned boundary conditions are summarized as follows: where n is the outward unit vector normal to the boundary. The Sommerfeld radiation condition is imposed for ψ z on the boundary B 3 . The pressure is calculated from the Poisson equation derived by taking the divergence of (2).

Distribution of water velocity
The distribution of water velocity around the hydrofoil is shown in Figure 3. When the angle of attack θ is 0 • , the flow parallels the hydrofoil surface and remains steady. But the flow at θ = 10 • separates from the upper surface of the rear half for the hydrofoil. The closeup of the velocity field near the separated region is shown in Figure 4. The flow exhibits an unsteady behavior.
The vortex elements introduced from the hydrofoil surface distribute as plotted in Figure 5, where the distribution at the same instant as Figure 3 is indicated. As the vortex elements generated by the bubble motion distribute in the whole region, their depiction is omitted. When θ = 0 • , the vortex elements flow along the hydrofoil surface as shown in Figure 5(a), since the water flow is in parallel with the surface. They flow almost straight in the wake of the hydro-  foil. The distribution of vortex element at θ = 10 • is shown in Figure 5(b). The flow separation from the upper surface of the hydrofoil is confirmed. The vortex elements meander markedly in the wake, and they form clusters. The clusters, corresponding to the large-scale eddies, flow downstream at almost even interval. The authors [12] analyzed such largescale eddies in a bubble-laden plane mixing layer by the vortex method. From the present simulation, one can grasp the generation and development for a velocity shear layer downstream of the hydrofoil. Figure 6 shows the distribution of air volume fraction α g at the same instant as Figures 3 and 5. When θ = 0 • , α g is low near the hydrofoil surface, as found from Figure 6(a). This is because the bubbles are prevented from approaching the leading edge due to the stagnation pressure, and accordingly they flow downstream in parallel with the hydrofoil surface. The α g value is also low in the wake. This is caused by the above-mentioned sparse distribution of bubble near the hydrofoil surface. Near the rear half of the upper-surface, α g takes its maximum value. The distribution for θ = 10 • is shown in Figure 6(b). The α g value reduces near the hydrofoil surface. The reduction near the lower surface is larger than that for θ = 0 • . Because the pressure rises in the downstream direction near the leading edge on the lower surface, accordingly such positive pressure gradient prevents the bubbles from approaching the hydrofoil surface. It should be noted that α g is higher near the rear half of the upper-surface. Since the bubbles are entrained into the large-scale eddies shed from the trailing edge, α g is also higher in the wake. It reaches its maximum value at the center of every eddy. Rotating blades are used in engineering devices handling the chemical reaction of gas-liquid two-phase mixtures. It is considered that the shed of largescale eddies and the bubble entrainment into the eddies can attain the active contact and mixing between the two phases. Murakami and Minemura [17] performed the experimental investigation on an axial-flow pump operated under air-water two-phase flow conditions. Their flow observation inside the impeller made clear that the entrained bubbles accumulate on the rear half of the suction surface for the blade. The present vortex method simulates such bubble accumulation near the rear half of the upper surface (suction surface) for a hydrofoil, as demonstrated in Figure 6. Consequently, it seems to be usefully employed for the simulation of twophase flow in rotating impeller.

Distribution of air volume fraction
Calculating the time-averaged air volume fraction α g on the six lines L 1 -L 6 indicated in Figure 7, it distributes as shown in Figure 8, where α g is plotted against the distance y from the hydrofoil surface. The result at θ = 0 • is shown by the broken lines. In the upper-surface regions (L 1 , L 2 , L 3 ), α g takes its maximum value near the surface. The maximum value increases in the downstream direction. In both sides of the hydrofoil, α g reduces with approaching the hydrofoil surface. The distributions of α g at θ = 10 • are plotted by the solid lines. Comparing with the results at θ = 0 • , α g is lower near the point C which locates in the separated region. It is also lower near the lower surface. Figure 9 shows the distribution of α g on the line L 0 passing through a point G downstream of the trailing edge. When θ = 0 • , α g reaches its minimum value at the center of wake (y /l = 0), indicating quite uneven distribution. It is found that the wake is markedly affected by the bubble distribution near the hydrofoil surface. Similarly, uneven distribution is also observed when θ = 10 • .
The above-mentioned distribution of α g agrees well with the trend of the measured result [7]. But the maximum value for α g is lower than the measurement. This may be because the present simulation ignores the bubble deformation. Figure 10 shows the distribution of pressure coefficient C p . The result at θ = 0 • is plotted in Figure 10  hydrofoil surface and causes the bubble accumulation there, accordingly α g heightens as found in Figure 6(a). The bubble accumulation does not occur near the lower-surface, because the pressure remains unaltered in the flow direction. Figure 10(b) shows the result at θ = 10 • . The pressure is locally high near the leading edge of the lower-surface. As the pressure gradient prevents the bubble from approaching the hydrofoil surface, it causes the marked decrement in α g near the lower-surface as shown in Figures 6 and 8. The time-averaged value for C p , C p distributes on the hydrofoil surface as shown in Figure 11. On the lower-surface, the distribution for the two-phase flow (α g0 = 0.048) is almost the same as that for the single-phase flow (α g0 = 0). This is because the air volume fraction is lower near the lowersurface, as found from Figure 6. On the upper-surface, however, the absolute value |C p | decreases at the two-phase flow condition. According to the experiment at Re = 4 × 10 5 [7], |C p | on the upper-surface reduces with an increment in α g0 . Therefore, the present simulation seems to be reasonable.

CONCLUSIONS
The air-water bubbly flow around a hydrofoil is simulated with the vortex method proposed by the authors in a prior study. A hydrofoil of NACA4412 with a chord length 100 mm is mounted in a channel. The water Reynolds number is 2.5 × 10 5 , the bubble diameter is 1 mm, and the volumetric flow ratio of bubble to whole fluid is 0.048. The results are summarized as follows.
(1) When the angle of attack θ is 0 • , the flow parallels the hydrofoil surface, and the air volume fraction α g takes its maximum value near the rear half of the uppersurface. When θ = 10 • , the flow separates from the upper-surface of the rear half, and large-scale eddies are shed from the trailing edge. The bubbles are en-trained into the eddies, and α g reaches its maximum value at the center of each eddy. (2) The time-averaged value for α g , α g takes its maximum value near the hydrofoil surface, and it reduces with approaching the surface. The maximum value increases in the downstream direction. The α g value at θ = 10 • is lower than that at θ = 0 • in the separated flow region and near the lower-surface. (3) The absolute value of the pressure on the uppersurface for the bubbly flow is lower than that for the single-phase flow. (4) The simulated results are in good agreement with the trend of the measurement. It is confirmed that the authors' vortex method is indeed applicable to the bubbly flow analysis around a hydrofoil.