Computational simulation of a single Taylor bubble rising in a vertical column with stagnant liquid

This work presents a computational simulation of a single Taylor bubble rising in a vertical column of stagnant liquid. The computational simulation was based on the Navier-Stokes equations for isothermal, incompressible, and laminar flow, solved by using the open source software OpenFOAM. The two fluids were assumed immiscible. The governing equations were discretized by the volume-of-fluid (VOF) method and solved using the Gauss iteration method. Parametric mesh was used to improve the modeling of curvilinear geometry. Numerical results were obtained for the rise velocities and shapes of the bubbles which are in excellent agreement with experimental data and correlations from literature.


INTRODUCTION
As a common gas-liquid two-phase flow pattern, slug flow is encountered in a variety of industrial applications, such as nuclear reactor cooling systems, evaporators, boilers, condensers, among others. When accompanied by fluctuations in flow rate and pipe temperature, high pipe wall temperatures in slug flow may result in "dryout" and cause damages in the nuclear power generating systems. Slug flow is characterized by long bullet-shaped bubbles, also called Taylor bubbles, which occupy nearly the entire cross-section of the pipe, and liquid slugs between successive bubbles. The liquid moves around the Taylor bubbles as a thin film with thickness δ and expands at the rear of the bubble, inducing a liquid wake ( Figure 1).

Figure 1: Schematic of a Taylor bubble with length Lb rising in a vertical column with the liquid film around the bubble.
Dumitrescu performed the first relevant study concerning individual bubbles rising in stagnant liquids and derived a Taylor bubble profile for air-water systems from potential flow theory, with a correlation to estimate the bubble velocity U0 [1]. From the experimental point of view, the pioneering work is attributed to Davies and Taylor, who carried out experiments using both water and nitrobenzene, and studied bubble velocity and profile, among other parameters [2]. The absolute velocity of a Taylor bubble in slug flow is usually expressed by the characteristic rise velocity U0 and a velocity component due to the movement of liquid phase, as follows [3]: where Ub is the velocity of the bubble, c is a constant, UL is the velocity of the liquid and U0 is the velocity of the bubble rising in stagnant liquid. For low viscosity liquids, the expression 0 = 0.35√ is in agreement with the works previously cited [1,2], where g is the gravity acceleration and D is the inner tube diameter. The shape profile derived by Dumitrescu was divided into two regions [1,4]: For Z/D < 0.25 (nose region), (2) For Z/D > 0.25 (film region), where Z is the axial distance from the tip of the bubble nose and r is the radial distance from the tube axis.
Commercial CFD codes and computers have been advanced significantly in past decades.
However, computational simulation of multiphase flows is still a considerably difficult task [5]. For a gas-liquid two-phase flow, the main difficulty is the existence of deformable interfaces.
Lagrangian and Eulerian strategies have been employed in develop a wide range of methods for advecting a sharp interface [6]. Currently, most CFD codes use variants of the volume-of-fluid (VOF) method for the interface advection step in their interfacial flow solvers. The VOF method can track the motion of gas-liquid interface using a transport equation for volume fraction occupied by each phase [7], which is suitable for simulations involving two immiscible fluids, and can be used to accurately predict the profile of the interface between the fluids [8].
A numerical method for the advection of a surface, such as the interface between two incompressible fluids, across a computational mesh has been devised by Roenby et al. [9]. The the InterFlow solver and the IsoAdvector method. Numerical results were obtained for the rise velocities and shapes of the Taylor bubbles, which were compared favorably with experimental data [10] and correlations from literature.

MATHEMATICAL FORMULATION
We consider a single Taylor bubble rising in a vertical column with an inner diameter D filled with stagnant liquid (Ub = U0), as defined in Eq. (1) and illustrated in Figure 1.

Governing Equations
The mathematical model is formed by the Navier-Stokes equations for the conservation of mass and momentum, respectively, applied in both gas and liquid phase: where, ρ is the fluid density, u is the velocity vector, p is the pressure, T is the stress tensor, g is the gravity vector, and fσ is the surface tension. A laminar flow and a Newtonian incompressible fluid is assumed to solve the mass and momentum conservation equations.
The VOF method employs an additional advection equation for the transport of the volume fraction α: where α = 1 corresponds to a control volume entirely occupied by the liquid and α = 0 corresponds to a control volume containing only air. The value of α is averaged in each of the mesh cells. The interface between the phases is found in cells where 0 < α < 1. The second term in Eq.(6) is referred to as the advection term [11].
In the VOF method, it is difficult to properly represent the value transitions at the interface between two phases, which can be overcome by the use of additional surface capture methods. In recent versions of OpenFOAM, the newly developed method (IsoAdvector) was introduced together with a solver named Interflow, which is an improvement of the InterFOAM solver.

The Interflow Solver
The Interflow solver uses a scheme named MULES for improving the surface sharpness.
MULES is a numerical scheme where the advection term in Eq. (6) is modified to compress the surface [11]. The scheme is obtained by firstly rewriting the advection equation in integral form: where Ωi represents each cell, ∂Ωi is the cell boundary and n is the cell boundary normal vector.
Equation (7) is then discretized, using a time-stepping scheme for the first term and writing the second term as a sum over each face of the cell: where Fu and Fc are the advective fluxes and λM is a delimiter taking the value 1 at the surface and 0 elsewhere.

The IsoAdvector Method
The IsoAdvector method uses the concept of isosurfaces to calculate more accurate face fluxes for the cells containing the interface [9,11]. The value for the phase fraction in cell i at time t, αi(t), is calculated from a function H(x; t) describing the continuous phase fraction field: where Vi is the volume of cell i and Ωi represents each cell.
Knowing the phase fraction in each cell at time t, it is necessary to calculate the phase fractions at the next time step using the following equation, where the flux of α over each cell face is integrated in time and added together: where Bi is the list of all faces Fj belonging to cell i, sij is used to orient the flux to going out from the cell and τ is the variable of integration used in the time step. dS is the differential area vector pointing out of the volume. sij is either +1 or -1 to ensure that the product sij dS is always in the direction out form the cell boundary even when the orientation of face j makes dS point into the cell.

Mesh Refinement and Boundary Conditions
It has been verified in a previous work that computational simulations of Taylor bubbles in stagnant water column using the InterFlow solver with parametric meshes showed to be better than those using nonparametric meshes [12]. When a nonparametric mesh is used to represent a cylindrical geometry, the central region of the cylinder assumes a square profile, which proved to be inappropriate for the bubble geometry.  The initial velocity of the phases was taken to be zero and the pressure was just the hydrostatic pressure since the simulation was performed for a closed tube.

RESULTS AND DISCUSSION
The numerical results were compared with correlations from literature and experimental data previously obtained from a vertical column consisting of acrylic tubes with 2.0 m in length and inner diameter of 0.024 m sealed at the ends. In this work, a Taylor bubble with length Lb was formed from different air pocket lengths inside tubes partially filled with water-glycerin mixtures, including pure distilled water. The rise velocities and shapes of the bubbles were determined by using a pulse-echo ultrasonic technique, with uncertainties estimated in 2% for Ub (relative uncertainty) and between 100μm and 280μm for δ on the experimental bubble shapes, respectively [10].
The rise velocities of the bubbles U0 were determined numerically after the full development of their flows, recovering the bubble axial positions at different time steps of the simulations and calculating the angular coefficients of the best linear fits of these points. Conservatively, the relative uncertainty on these determinations were estimated as smaller than 4%. An excellent linear fit was observed for all the simulated cases, which confirmed that the bubble motion was fully developed, as shown in Figure 3.   Table 1.   Figure 5a indicates that the simulated velocity of the bubble tended to improve the agreement with the theoretical and experimental data as the number of cells in the height of the cylinder increased. Together with the simulated bubble shapes, Figure 5b also shows the correlation proposed to predict the shape of a Taylor bubble rising in stagnant water (Eqs. (2) and (3)) and the experimental bubble shape [1,4,10]. It can be observed that the bubble shape also improves the agreement with the theoretical and experimental shapes as the number of cells in the height of the cylinder increased. The results showed in Figures 5a and 5b are physically consistent since the velocity of a Taylor bubble rising in stagnant liquid is closely related with its shape [1,2]. Thus, a better simulation of the bubble shape resulted in better simulated results for bubble velocity.   The results presented in Figures 5-7 indicate that the IsoAdvector method of the Interflow solver proved to be a powerful tool for computational simulation of a Taylor bubble moving inside cylindrical tubes and for the analysis of the hydrodynamic characteristics of slug flow in tubes.
The method was shown to determine accurately the interface between the fluids, resulting in an excellent prediction of the velocity and shape of the bubble. However, these figures also show that the mesh construction for a good performance simulation must be done with criteria. For the simulation proposed in the present work, the best results were obtained for mesh constructions with a minimum number of 6 cells at the inner square, 30 cells between the square and the circle and 160 cells in cylinder height. Naturally, a definition of the best mesh construction for a given simulation must take into account the balance between the accuracy of results and the computational costs.
It is important to note that previous simulations of single Taylor bubbles rising in vertical tubes performed by using the InterFOAM solver, without the IsoAdvector method, did not produce satisfactory results. The bubble tending to fragment as it moved inside the tube and this behavior was not compatible with the reality or with experimental observations. The use of the IsoAdvector method seems to be essential in solving this problem.
Thus, will be possible in future works to perform computational simulations of the movement of Taylor bubbles inside inclined cylindrical tubes by using the Interflow solver of OpenFOAM. For vertical tubes it is expected to determine other bubble flow parameters besides velocity and shape of the bubble and compare them with experimental measurements [10,13] and with data available in literature. Since IsoAdvector method was able to accurately determine the interface between the fluids for the vertical case, the next step will test it for the simulation of elongated bubbles inside vertical and inclined slug flows.

CONCLUSION
Computational simulations of single Taylor bubbles rising in a vertical column with stagnant liquid were performed using the InterFlow solver of the OpenFOAM software with parametric meshes. The numerical results were compared with experimental data and correlations available in the literature. Some conclusions of this work are listed below: -The IsoAdvector method of the Interflow solver was able to determine accurately the interface between the fluids, resulting in an excellent prediction of the rise velocity and shape of the bubble.
-An accurate simulation depends critically on the appropriate mesh construction. For the simulation reported in the present work, the best results were obtained for mesh constructions with a minimum number of 6 cells at the inner square, 30 cells between the square and the circle and 160 cells in cylinder height.
-The obtained results indicate that the IsoAdvector method of the Interflow solver proved to be a powerful tool for computational simulation of a Taylor bubble moving inside cylindrical tubes and for the analysis of the hydrodynamic characteristics of slug flow in tubes.