Tomomi Uchiyama^{1*} and Yusuke Kishimoto^{2}  
^{1}EcoTopia Science Institute, Nagoya University, Furocho, Chikusaku, Nagoya 4648603, Japan  
^{2}Fujitsu Limited, 152 HigashiShimbashi, Minatoku, Tokyo 1057123, Japan  
Corresponding Author :  Tomomi Uchiyama Eco Topia Science Institute Nagoya University, Furocho Chikusaku, Nagoya 4648603, Japan Tel: +81527895187 Fax: +81527895187 Email: [email protected] 

Received August 29, 2014; Accepted September 26, 2014; Published September 30, 2014  
Citation: Uchiyama T, Kishimoto Y (2014) Numerical Simulation of Airwater Bubbly Jet Issuing from a Square Nozzle by a Vortex in Cell Method. J Chem Eng Process Technol 5:207. doi:10.4172/21577048.1000207  
Copyright: © 2014 Uchiyama T, et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.  
Related article at Pubmed Scholar Google 
Visit for more related articles at Journal of Chemical Engineering & Process Technology
A water jet entrained with small air bubbles is simulated. The Vortex in Cell method for bubbly flow proposed by the authors in a prior study is used for the simulation. The method discretizes the vorticity field into vortex elements and computes the time evolution of the flow by tracing the convection of each vortex element using the Lagrangian approach. The bubbly jet issues vertically upward from a nozzle of square crosssection. The Reynolds number based on the water velocity is 5000. The bubble diameter is 0.2 mm, and the bubble volumetric flow rate ratio at the nozzle is 0.0025. The bubbles increase the water turbulence intensity around the jet centerline. This makes the water momentum diffusion in the lateral direction larger at the initial region of the jet. The bubble effect lessens with increasing axial distance from the nozzle. These simulated results are favorably compared with existing measurements, demonstrating the validity of the current simulation method for bubbly jet. In the developed region, the water velocity decay is relaxed and the halfwidth reduces. This is because the water is accelerated by the bubbles, which have higher axial velocity due to the buoyant force.
Keywords  
Gasliquid bubbly flow; CFD; Jet; Vortex in cell method; Bubble behavior  
Notations  
B: side length of nozzle crosssection [m]; b: halfwidth for singlephase water jet [m]; d: bubble diameter [m]; Eo: Eötvös number ; g: gravitational constant [m/s^{2}]; p: pressure [Pa]; Q: second invariant of the water velocity gradient tensor [1/s^{2}]; Re: Reynolds number=U_{0}B / ν_{l} []; Reb: bubble Reynolds number= t: time [s]; u: velocity [m/s]; U_{0}: water velocity at nozzle exit [m/s]; ν: volume of single bubble=πd^{3}/6 [m^{3}]; W: redistribution function of vorticity []; x, y, z: orthogonal coordinates [m]; α: volume fraction []; β: density ratio=ρ_{g} / ρ_{l} []; Δt: time increment [s]; Δx, Δy, Δz: grid width in direction of x, y, z [m]; ν: kinematic viscosity [m2/s]; ρ: density [kg/m^{3}]; σ: surface tension [kg/s^{2}]; ϕ: scalar potential [m2/s]; ψ: vector potential [m^{2}/s]; ω: vorticity of liquidphase=[1/s]; c: jet centerline; g: gasphase; l: liquidphase; rms: root mean square; x, y, z: components in directions of x, y, z  
Introduction  
Vortex in Cell (VIC) method is one of the vortex methods simulating incompressible flows [1]. It discretizes the vorticity field into vortex elements and computes the time evolution of the flow by tracing the convection of each vortex element using the Lagrangian approach. The Lagrangian calculation markedly reduces numerical diffusion and also improves numerical stability. Thus, the VIC method is eminently suitable for direct numerical simulation (DNS) of turbulent flows. Cottet and Poncet [2] applied the VIC method to the wake simulation of a circular cylinder mounted in a uniform flow field and successfully captured the stream wise vortices occurring behind the cylinder. Cocle et al. [3] analyzed the behavior of two vortex system near a solid wall and made clear the interaction between the two counterrotating vortices and the eddies induced in the vicinity of the wall. Chatelain et al. [4] simulated the interaction between trailing vortices and visualized the unsteady phenomena caused by disturbances. These studies demonstrate the applicability of the DNS by the VIC method to free turbulent flows which are dominated by the convection of largescale eddies.  
The authors [5] have previously proposed two improvements of the VIC method. First, a staggeredgrid discretization method guarantees consistency among the discretized equations and prevents numerical oscillations in the solution. Second, a correction method for the vorticity enables the computation of vorticity fields satisfying the solenoidal condition. The improved VIC method was applied to the DNS of a turbulent channel flow [5]. The DNS highlighted the successful capture of organized flow structures, such as streaks and streamwise vortices in the near wall region, demonstrating that the VIC method is applicable to DNS of wall turbulent flows. Uchiyama [6] presented a simulation method for bubbly flow by using the VIC method. The behavior of vortex element and the bubble motion are simultaneously analyzed by the Lagrangian method. Uchiyama [7] also proposed a VIC method for incompressible gas flows laden with small solid particles. This method was favorably validated in simulations of small freefalling solid particles in unbounded air [7]. Uchiyama and Shimada [8] used the VIC method to simulate the interaction between a vortex pair and small solid particles near a horizontal wall in the air. The simulation made clear the agitation of particles by the vortex pair approaching the wall, the production of vorticity fields by the particles, and the particleinduced changes in the strength and behavior of the vortex pair.  
Liquid jets entraining small gas bubbles are utilized in equipment, such as chemical reactors, heat exchangers and aerators. A number of experimental researches were conducted to investigate the relationship between the largescale eddies and the bubble distribution [9] as well as the bubble effect on the liquid velocity field [1014]. But there are no numerical simulations except for a few studies. Sun and Faeth [10,11] simulated a bubbly jet. But they assumed an axisymmetric flow to employ a kε model [15] developed for gasparticle twophase flow. Thus, the validity of the simulation is unclear. Xing et al. [16] reported a DNS for a cavitating jet. As they applied a homogeneous model which assumes the gas and the liquid have the same velocity, the interactions between the two phases were not fully analyzed. Consequently, a reliable simulation method for bubby jets is not presented.  
In this study, the author’s VIC method for bubbly flow [6] is applied to simulate an airwater bubbly jet issuing vertically upward from a nozzle into quiescent water. The Reynolds number based on the water velocity is 5000, and the bubbles having a diameter 0.2 mm are issued from the nozzle with the same velocity as the water. The bubble volumetric flow rate ratio at the nozzle is 0.0025. The simulation clarifies the effects of the bubbles on the water flow, such as the changes in the vorticity field and the relaxation of the velocity decay at the jet initial region. It also highlights that such simulated results are in good agreement with the measured ones of Kumar et al. [12], demonstrating the validity and applicability of the VIC method to bubbly jet simulations.  
Basic Equations and Numerical Method  
Assumptions  
The following assumptions are employed for the simulation:  
(a) The mixture is a gasliquid bubbly flow entraining small bubbles.  
(b) Both phases are incompressible and no phase changes occur.  
(c) The mass and momentum of the gasphase are very small and negligible compared with those of the liquidphase.  
(d) The bubbles maintain their spherical shape, and neither fragmentation nor coalescence occurs.  
Governing equations for bubbly flow  
The mass conservation equation for the liquidphase and that for the gasphase are independently derived. Taking the summation of them and rearranging the resultant equation with the assumptions (a)(c), the mass conservation equation for the twophase mixtures is obtained [6]:  
(1)  
where α_{l} is the liquid volume fraction satisfying the following relation with the gas volume fraction α_{g}:  
(2)  
The momentum conservation equation for the twophase mixtures, which is also derived from the equations for each phase by the same manner as the mass conservation equation, is expressed as [6]:  
(3)  
where  
It is postulated that the virtual mass force, the drag force, the gravitational force, and the lift force act on the bubble. In this case, the equation of motion for the bubble is expressed by the following equation [17,18] with the assumption (d):  
(4)  
where d is the bubble diameter, and β is the density ratio defined as ρ_{g} ρl. C_{V}, C_{D} and C_{L} are the virtual mass coefficient, the drag coefficient, and the lift coefficient respectively.  
Vorticity equation and orthogonal decomposition of liquid velocity  
When taking the curl of Eq. (3), the vorticity equation for the bubbly flow is derived:  
(5)  
where ω is the vorticity of the liquidphase.  
(6)  
According to the Helmholtz theorem, any vector field can be represented as the summation of the gradient of a scalar potential ϕ and the curl of a vector potential ψ. The liquid velocity u_{l} is thus expressed as  
(7)  
The velocity calculated from Eq. (7) remains unaltered when any gradient of a scalar function is added to ψ. To remove this arbitrariness, the following solenoidal condition is imposed on ψ:  
(8)  
When substituting Eq. (7) into Eq. (1), the following equation is obtained:  
(9)  
Taking the curl of Eq. (7) and substituting Eq. (8) into the resultant equation, the following vector Poisson equation for ψ is derived:  
(10)  
Simulation Based on VIC Method  
Discretization of vorticity field into vortex elements  
Once φ and ψ have been computed from Eqs. (9) and (10) respectively, the velocity u_{l} is calculated from Eq. (7). The vorticity ω in Eq. (10) is estimated from Eq. (5). The VIC method discretizes the vorticity field into vortex elements and calculates the distribution of ω by tracing the convection of each vortex element.  
It is postulated that the position vector and vorticity for the vortex element p are x_{p}=(x_{p}, y_{p}, z_{p}) and ω_{p}, respectively. The Lagrangian form of the vorticity equation, Eq. (5), is written as follows:  
(11)  
(12)  
When the position and vorticity of a vortex element are known at time t, the values at t+Δt are computed from Eqs. (11) and (12). In the VIC method, the flow field is divided into computational grid cells to define ψ, φ and ω on the grids. If ω is defined at a position x_{k}=(x_{k}, y_{k}, z_{k}),the vorticity ? is assigned to x_{k}, or a vortex element with vorticity ? is redistributed onto x_{k}.  
(13)  
where N_{ν} is the number of vortex elements, and Δx, Δy and Δz are the grid widths. For the redistribution function W, the following equation is employed [19].  
(14)  
Calculation of gas volume fraction  
The abovementioned grid cells are also used to calculate the gas volume fraction αg. The bubble diameter is prescribed to be smaller than the grid width on the basis of the assumption (a). When a bubble with volume ν_{r} exists at a position x_{g}=(x_{g}, y_{g}, z_{g}), the α_{g} value at the grid point of x_{q}=(x_{q}, y_{q}, z_{q}) is calculated by the following equation:  
(15)  
where N_{b} is the number of bubbles. For the function Wα, the following equation, which is the redistribution function of vortex element [19], is employed.  
(16)  
Discretization with staggered grid and correction of vorticity  
This simulation employs a staggered grid to solve Eqs. (9) and (10) so as to ensure consistency between the discretized equations and to prevent numerical oscillations of the solution. Figure 1 shows the grid. The scalar potential φ and the gas volume fraction α_{g} are defined at the center of a grid cell. The liquid velocity ul is defined at the sides, while the vorticity ω and the vector potential ψ are defined on the edges.  
The vorticity field is discretized with vortex elements, and it is expressed by the superposition of the vorticity distribution around each vortex element as found from Eq. (13). Thus, the vorticity field obtained from Eq. (13), denoted by , does not always satisfy the solenoidal condition. The curl of the velocity calculated by using results in the vorticity satisfying the solenoidal condition. This correction method for the vorticity is used in this simulation.  
Numerical procedure  
Given the flow at time t, the flow at t+Δt is simulated by the following procedure:  
1. Calculate the bubble motion from Eq. (4).  
2. Calculate α_{g} from Eq. (15), and compute α_{l} from Eq. (2).  
3. Calculate the time variation of ω at every grid point from Eq. (12).  
4. Calculate the convection of each vortex element from Eq. (11).  
5. Calculate ω from Eq. (13).  
6. Calculate ψ from Eq. (10).  
7. Calculate φ from Eq. (9).  
8. Calculate u_{l} from Eq. (7).  
9. Correct the vorticity, or calculate the corrected vorticity from the curl of u_{l}.  
Simulation Conditions  
bubbly jet issuing from a nozzle having a square crosssection is simulated. The flow direction is vertically upward. The nozzle and the coordinate system are shown in Figure 2. The yz plane is horizontal, and the xaxis is taken in the vertically upward direction along the jet centerline.  
Table 1 lists the simulation conditions, Airwater bubbly mixtures issues from the nozzle. The side length of the nozzle exit B is 10 mm. The computational domain consists of a hexahedral region of 45B×30B×30B. It is divided into 150×300×300 grid cells. A simulation using more computational cells 180×360×360 was preliminarily performed. But  
the bubble distribution and the water velocity field were less affected by the increase in the number of cells. Therefore, the cell size is found to be appropriate. The water has a uniform velocity U_{0} at the nozzle exit, and the Reynolds number Re (=U_{0}B/ν_{l}) is 5000. The time increment Δ_{t} is set at U0Δ_{t}/B=0.015.  
The bubble diameter d is 0.2 mm. The bubble volumetric flow rate ratio at the nozzle exit is 0.0025. The bubbles are released from the nozzle exit at a time interval Δt with the same velocity as the water. The number of released bubbles is specified by the volumetric flow rate ratio. The releasing positions are determined by using random numbers. The relation of C_{L}=C_{V}=0.5 is used in Eq. (4) [18], and the drag coefficient C_{D} is given as:  
(17)  
where Re_{b} is the bubble Reynolds number defined as Equation (17) combines the theoretical formula of Moore [20] and the Stokesian drag coefficient. It was also used for a bubble plume simulation by Murai and Matsumoto [21].  
A nonslip condition is imposed on the boundary (x/B=0), on which the nozzle exit is mounted. The velocity gradient is set at zero on the other boundaries [22]. The Poisson equations, Eqs. (9) and (10), are solved by the SOR method.  
Results and Discussions  
Singlephase water jet  
First, the present method is applied to simulate the singlephase water jet. When the developed jet occurs, the isosurface for the second invariant of the velocity gradient tensor Q distributes as shown in Figure 3. The vortical structure, in which vortex tubes having various scales entangle, is grasped.  
Figure 4 shows the vorticity component parallel to the jet axis at the same time point as Figure 3, where the isosurfaces of ω_{x}/(U_{0}/B)=± 0.002 are plotted. Pairs of positive and negative vortex tubes exist and entangle, suggesting the appearance of highly threedimensional vertical flow.  
The mean velocity on the jet centerline u_{k} changes as a function of the axial distance from the nozzle as shown in Figure 5. The velocity decays markedly at x/B ≥ 2. Mi et al. [23] measured the velocity distribution of a jet of Re=15000 at the fullydeveloped region. The present result agrees nearly with the measurement. The velocity decay at the fullydeveloped region of the jet satisfies a relation of [24]. In this simulation, such relation is obtained.  
The change in the halfwidth b is shown in Figure 6. The present result is slightly larger than the measured one by Tsuchiya et al. [25] and the experimental formula of Quinn and Militzer [26]. The change in the axial direction is favorably computed.  
Figure 7 shows the lateral distribution of the timeaveraged axial velocity u_{lx}, where the centerline velocity u_{lc} and the halfwidth b are used for the nondimensional expression. Though the effect of the potential core appears on the section at x/B=3, the similarity of velocity distribution is analyzed on the sections at x/B ≥ 6.  
The axial component of the turbulent intensity on the jet centerline varies as shown in Figure 8. It takes the maximum value at x/ B=2.4. The present result agrees well with the measured one of Mi et al. [23] at the fullydeveloped region.  
From the abovementioned results, one can confirm the validity of the present method for singlephase jet simulations.  
Bubble distribution and vortical structure  
Figure 9 shows the bubble distribution of the fullydeveloped bubbly jet. The bubbles disperse markedly in the lateral (horizontal) direction at x/B ≥ 4.3. This is attributable to the water momentum diffusion in the direction. The dispersion becomes larger with increasing axial distance from the nozzle.  
According to Clift et al. [27], the shape of a bubble rising in a viscous fluid under the influence of gravity is specified by the bubble Reynolds number, ,and the Eötvös number, , where σ is the surface tension. In this simulation, R_{eb} is less than 16.1 and Eo=0.0527. The bubble shape is specified to be spherical [27], demonstrating the validity of the assumption (d).  
The isosurface of Q for the bubbly jet at the same time point as Figure 9 is presented in Figure 10. When compared with the result for the singlephase water jet shown in Figure 3, the vortical region spreads more in the horizontal direction at x/B ≤ 7. One can grasp the increase in the jet width due to the bubbles. But the width of the vortical region reduces at x/B ≥ 30. The bubbles hardly affect the scale of the vortex tubes.  
Figure 11 shows the isosurfaces of the axial component of the vorticity at the same time point as Figure 10. The isosurfaces are divided into finer ones when compared with the result for the singlephase water jet (Figure 4). The threedimensional features of the vortical structure become remarkable owing to the bubbles distributing nonuniformly in the jet.  
Bubble motion and water velocity field  
The mean velocities of water and bubble on the jet centerline change as function of the axial distance from the nozzle as shown in Figure 12. The water velocity of the bubbly jet is slightly lower than that of the singlephase jet at x/B ≤ 9.6 while it is higher downstream of the region. The bubble velocity is slightly higher than the water, because the buoyant force acts on them. The lesson of the water velocity of the bubbly jet at x/B ≤ 9.6 is caused by the increase in the water momentum diffusion in the horizontal direction due to the bubbles as explained later. The relaxation of the water velocity decay at x/B>9.6 is attributable to the fact that the water is accelerated by the bubbles flowing faster than the water. The increase in the water momentum diffusion due to the bubbles at jet initial region was also clarified experimentally by Kumar et al. [12].  
A numerical simulation of an air jet laden with small solid particles (diameter 100 μm, density 2590 kg/m^{3}) also made clear such air velocity reduction at the jet initial region [28]. It seems that a dispersed phase has the similar disturbances on the continuous phase in jet flows.  
Figure 13 shows the axial evolution of the halfwidth of the water velocity. At x/B ≤ 7.5, the halfwidth of the bubbly jet is larger than that of the singlephase jet. This is caused by the abovementioned increase in the water momentum diffusion in the lateral direction due to the bubbles. At x/B ≥ 7.5, it is smaller because the velocity decay is relaxed as shown in Figure 12.  
The lateral distribution of the water and bubble velocities is presented in Figure 14, where the water velocity at the nozzle exit U_{0} and the halfwidth for the singlephase jet b are used for the nondimensional expression. On the sections at x/B=3 and 12, the bubble velocity is higher than the water. On the section at x/B=3, the water velocity of the bubbly jet is slightly lower and higher than the singlephase jet at z/b ≤ 0.86 and z/B ≥ 0.86, respectively. This is because the water momentum diffusion in the lateral direction increases due to the bubbles. On the section at x/B=1, the water velocity of the bubbly jet is higher than the singlephase jet. This is attributable to the acceleration effect by the bubbles. The velocity difference caused by the bubbles is lower than that on the section at x/B=3.  
The axial component of the water turbulence intensity distributes on the jet centerline as shown in Figure 15. The turbulence intensity of the bubbly jet is larger than that of the singlephase jet. This result is parallel with the measured one of Kumar [12].  
Figure 16 shows the lateral distribution of the axial component of the water turbulence intensity. The intensity is heightened by the bubbles. The marked increase is observed on the section at x/B=3, where the water momentum diffusion is promoted by the bubbles.  
In Figures 14, 15 and 16, it is found that the bubble effects on the water flow lessen with increasing axial distance from the nozzle. This is because the bubbles disperse in the lateral (horizontal) direction, and accordingly the bubble volume fraction around the jet centerline reduces. Such decrease in the bubble effect was also reported by Kumar [12].  
Conclusions  
An airwater bubbly jet, issuing vertically upward from a nozzle of a squarecross section, is simulated. The Vortex in Cell method for bubbly flow proposed by the authors in a prior study is used for the simulation. The Reynolds number of the water is 5000, the bubble diameter is 0.2 mm, and the bubble volumetric flow rate ratio at the nozzle exit is 0.0025. The results are summarized as follows:  
(1) The bubbles heighten the water turbulence intensity and cause the increase in the water momentum diffusion in the lateral direction at the initial region of the jet. This simulated result is in good agreement with the existing measured one, demonstrating the validity of the simulation method.  
(2) At the developed region, the water velocity decay is relaxed and the spread of the jet in the lateral direction decreases. This is because the bubbles have the higher velocity than the water.  
(3) The threedimensional features of the water vortical structure become more remarkable owing to the bubbles.  
References  

Table 1 
Figure 1  Figure 2  Figure 3  Figure 4  Figure 5  Figure 6 
Figure 7  Figure 8  Figure 9  Figure 10  Figure 11  Figure 12 
Figure 13  Figure 14  Figure 15  Figure 16 