Using a Dynamic and Constant Mesh in Numerical Simulation of the Free-Rising Bubble

: A two-phase bubbly ﬂow is often found in the process industry. For the efﬁcient operation of such devices, it is important to know the details of the ﬂow. The paper presents a numerical simulation of the rising bubble in a stagnant liquid column. The interFOAM solver from the open source Computational Fluid Dynamics (CFD) toolbox OpenFOAM was used to obtain the necessary data. The constant and dynamic computational grids were used in the numerical simulation. The results of the calculation were compared with the measured values. As expected, by using the dynamic mesh, the bubble trajectory was closer to the experimental results due to the more detailed description of the gas–liquid interface.


Introduction
Buoyancy-driven bubbly flow represents very demanding two-phase time-dependent three-dimensional flow, which is found in many natural and technological processes. Good knowledge of such flow is extremely important for the chemical and pharmaceutical industries where mixing vessels [1], bubble columns [2,3], and gas-lift reactors [4] can be found. In order to control the chemical reactions in these devices, the knowledge of the contact area between liquid and gas is crucial. The interfacial area depends on the flow regime or bubble size distribution. Another very important application of bubbly flow is in the field of heat transfer. Free-rising bubbles can cause the mixing of water layers with a temperature gradient, which is the so-called destratification [5], or increase heat transfer in heat exchanging devices by breaking the thermal boundary layer along the wall [6].
The complexity of the bubble flow represents a major challenge for theoretical, experimental, and numerical research. Sometimes it is very difficult to provide control over the experimental process and thus, obtaining accurate data due to the interference between the measuring probe and fluid flow can be challenging. Therefore, numerical simulation is increasingly used as a tool for the analysis of the bubbly flow [7][8][9].
The velocity fields of the two-phase flow can be significantly modified by the interfacial interactions, meaning that they must be taken into account through adequate closure laws. The determination of the interface depends on the considered problem. There are two approaches to describe the development of the interface. One is the interface tracking method [10] (Lagrangian), which is based on the use of markers on the interface. Another one is the interface capturing (Eulerian) method, which uses fixed mesh to solve the transport equation of the scalar quantity. Level-Set [11,12], Volume of Fluid (VOF) [13,14] and Euler-Euler [15,16] methods are representatives of this approach. When the bubble is substantially larger than the cell size of the computational grid, the VOF method is most commonly used. The drawback of this method is that we obtain a smearing interface when using the interpolation schemes of the lower order and the numerical oscillations when using higher order interpolation schemes. A possible solution involves the refinement of the computational mesh, but this requires a significant increase in computer resources (memory, computational time). This can be avoided by using the "Adaptive Mesh Refinement" (AMR) technique [17], in which the computational grid is refined only in the range of rapid changes in physical quantities (e.g., cells near the gas-liquid interface) and therefore, the computational mesh is dynamically changing as the refined mesh follows the movement of the bubble.
We have presented the numerical simulation of free-rising air bubble in a stagnant water column using constant (CM) and dynamic (DM) computational mesh in the following sections. The results are compared with experimental data obtained in the Laboratory for Fluid Dynamics and Thermodynamics during 1980-1992. The paper is organized as follows. In Section 2, the numerical model is described. The numerical simulation is presented in Section 3. In Section 4, the results of the simulation are presented and discussed. Finally, the conclusions are drawn in Section 5.

Governing Equations
To determine the gas-liquid interface, we used the VOF method. The surface tension force was taken into account using the Continuum Surface Force (CSF) model [18]. Due to a low Reynolds number, the conditions were considered to be laminar. The time-dependent laminar incompressible two-phase flow can be described with the solution of the conservation equations: where α is the liquid volume fraction, u i is the velocity vector, p is the pressure, σ ij is the stress tensor due to the molecular viscosity, g i is the gravity acceleration vector, σ is the surface tension coefficient and κ is the gas-liquid interface curvature. Equation (1) and Equation (2) are connected through material properties. The mixture density ρ and mixture dynamic viscosity µ are determined by the following equations: where ρ l and ρ g are the liquid and gas density, whereas µ l and µ g are liquid and gas dynamic viscosity. The system of Equations (1)-(4) was solved within the OpenFOAM [19,20] toolkit, which has proven to be an effective tool for solving such problems [21]. Table 1 shows the material properties for air-water two-phase system, which was considered in our case.

Discrete Model
The computational domain is represented by a hexahedron with dimensions L x = 40 mm, L y = 80 mm and L z = 290 mm ( Figure 1). The domain is discretized with the structured computational grid. The number of cells in x, y and z coordinate directions is N x = 100, N y = 200 and N z = 725. The total number of cells is 14.5 million and the size of individual cell is 0.4 mm. This mesh resolution is a consequence of the VOF two-phase flow simulation practice that there should be at least 15 cells per bubble diameter [5,22]. This recommendation applies in the case of constant computational mesh ( Figure 2a). In the case of the dynamic mesh, the cell size varies in the part of computational domain according to the prescribed condition (e.g., 0.001 ≤ α ≤ 0.999). This means that cell is near gas-liquid interface ( Figure 2b). The number of cells varied between 14.79 and 16.09 million.
As an initial condition, we put a spherical bubble with a diameter of 5.8 mm at the bottom of the domain. The bubble center is located at coordinates x 0 = 20 mm, y 0 = 40 mm and z 0 = 10 mm. The free surface occurs at the height H = 280 mm. consequence of the VOF two-phase flow simulation practice that there should be at least 15 cells per bubble diameter [5,22]. This recommendation applies in the case of constant computational mesh ( Figure 2a). In the case of the dynamic mesh, the cell size varies in the part of computational domain according to the prescribed condition (e.g., 0.001 ≤ α ≤ 0.999). This means that cell is near gas-liquid interface ( Figure 2b). The number of cells varied between 14.79 and 16.09 million.
As an initial condition, we put a spherical bubble with a diameter of 5.

Numerical Simulation
The combination of the Pressure Implicit with Splitting of Operator (PISO) and Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), which is the so-called PIMPLE algorithm, was used for velocity-pressure coupling [19,23]. The volume fraction, momentum and turbulence equations were solved with smoothSolver using symGaussSeidel smoother. For the pressure equation, we used Geometric Agglomerated Algebraic Multigrid (GAMG) solver with the Diagonal Incomplete-Cholesky (DIC) preconditioner for symmetric matrices.
The discretization schemes were: • Euler explicit for the time derivative,  consequence of the VOF two-phase flow simulation practice that there should be at least 15 cells per bubble diameter [5,22]. This recommendation applies in the case of constant computational mesh (Figure 2a). In the case of the dynamic mesh, the cell size varies in the part of computational domain according to the prescribed condition (e.g., 0.001 ≤ α ≤ 0.999). This means that cell is near gas-liquid interface (Figure 2b). The number of cells varied between 14.79 and 16.09 million.
As an initial condition, we put a spherical bubble with a diameter of 5.

Numerical Simulation
The combination of the Pressure Implicit with Splitting of Operator (PISO) and Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), which is the so-called PIMPLE algorithm, was used for velocity-pressure coupling [19,23]. The volume fraction, momentum and turbulence equations were solved with smoothSolver using symGaussSeidel smoother. For the pressure equation, we used Geometric Agglomerated Algebraic Multigrid (GAMG) solver with the Diagonal Incomplete-Cholesky (DIC) preconditioner for symmetric matrices.
The discretization schemes were: • Euler explicit for the time derivative, • Gauss linear for the gradient operator,

Numerical Simulation
The combination of the Pressure Implicit with Splitting of Operator (PISO) and Semi-Implicit Method for Pressure-Linked Equations (SIMPLE), which is the so-called PIMPLE algorithm, was used for velocity-pressure coupling [19,23]. The volume fraction, momentum and turbulence equations were solved with smoothSolver using symGaussSeidel smoother. For the pressure equation, we used Geometric Agglomerated Algebraic Multigrid (GAMG) solver with the Diagonal Incomplete-Cholesky (DIC) preconditioner for symmetric matrices.
The discretization schemes were: • Euler explicit for the time derivative, • Gauss linear for the gradient operator, • Gauss upwind for the divergence operator • Gauss linear corrected for the Laplacian operator.
We used the adjustable time step. This means that ∆t was changed so that the Courant number did not exceed the prescribed value. In our simulation, we used the following maximum parameter values max(Co) = maxAlpha(Co) = 0.5 and consequently, the time step was varied between 1. The difference is very large, but it is necessary to take two things into account: 1.
The division of computational cells increases the number of equations considerably.

2.
The ∆t decreases due to the small cells and consequently increases the number of time steps. Figure 3 shows the trajectory of the bubble center or precisely the deviation s x = x c − x c and s y = y c − y c from the average value in x and y direction, respectively. • Gauss linear corrected for the Laplacian operator. We used the adjustable time step. This means that Δt was changed so that the Courant number did not exceed the prescribed value. In our simulation, we used the following maximum parameter values max(Co) = maxAlpha(Co) = 0.5 and consequently, the time step was varied between 1. The difference is very large, but it is necessary to take two things into account: 1. The division of computational cells increases the number of equations considerably. 2. The Δt decreases due to the small cells and consequently increases the number of time steps. Initially, the stationary spherical bubble begins to move upwards due to the buoyancy force. As the bubble velocity rises, surface waves appear as a result of viscous friction and surface tension forces, which destroy symmetry and cause its shape to change from spherical to ellipsoidal. The bubble begins to wobble and travel along the spiral path, leaving behind the characteristic decaying vorticity structures [24], which is the so-called 'horseshoe' vortex (Figure 4), which further deforms  Initially, the stationary spherical bubble begins to move upwards due to the buoyancy force. As the bubble velocity rises, surface waves appear as a result of viscous friction and surface tension forces, which destroy symmetry and cause its shape to change from spherical to ellipsoidal. The bubble begins to wobble and travel along the spiral path, leaving behind the characteristic decaying vorticity structures [24], which is the so-called 'horseshoe' vortex (Figure 4), which further deforms its shape. The bubble mean lateral displacement is defined by:

Results and Discussion
where sj is the bubble center deviation from initial position: where tj is the time and N is the number of saved time steps. The bubble mean lateral displacement is s ∞ = 2.26 mm in the case of a constant mesh and s ∞ = 1.42 in the case of a dynamic mesh. The comparison of these two values with the experiment [25] is shown in Figure 5. As we can see, the dynamic mesh value is closer to the measured trend-line than the constant mesh result. The bubble terminal rise speed is defined by: For a constant mesh, the bubble terminal velocity is v ∞ = 22.56 cm/s and for a dynamic mesh, this is v ∞ = 22.59 cm/s. Both values match well with the experimental data ( Figure 6).
By using the Fast Fourier Transformation (FFT) for the time series of the bubble center coordinates, we calculate the mean oscillation frequency in the bubble wobbling motion. The result is f ∞ = 5.21 Hz for a constant mesh and f ∞ = 4.55 Hz for a dynamic mesh. The latter value is closer to the measured trend line (Figure 7).
The non-dimensional Weber number We = ρ σ ⁄ represents the ratio of inertial force and surface tension force. The dependence of the Weber number on the bubble equivalent sphere diameter is shown in Figure 8. For the bubble with de = 5.8 mm, the calculated value for a constant mesh is We = 4.05 and We = 4.07 for a dynamic mesh. Both values are slightly lower than the measured The bubble mean lateral displacement is defined by: where s j is the bubble center deviation from initial position: where t j is the time and N is the number of saved time steps. The bubble mean lateral displacement is s ∞ = 2.26 mm in the case of a constant mesh and s ∞ = 1.42 in the case of a dynamic mesh. The comparison of these two values with the experiment [25] is shown in Figure 5. As we can see, the dynamic mesh value is closer to the measured trend-line than the constant mesh result. The bubble terminal rise speed is defined by: For a constant mesh, the bubble terminal velocity is v ∞ = 22.56 cm/s and for a dynamic mesh, this is v ∞ = 22.59 cm/s. Both values match well with the experimental data ( Figure 6). By using the Fast Fourier Transformation (FFT) for the time series of the bubble center coordinates, we calculate the mean oscillation frequency in the bubble wobbling motion. The result is f ∞ = 5.21 Hz for a constant mesh and f ∞ = 4.55 Hz for a dynamic mesh. The latter value is closer to the measured trend line (Figure 7).
The non-dimensional Weber number We = ρ l v 2 ∞ d e /σ represents the ratio of inertial force and surface tension force. The dependence of the Weber number on the bubble equivalent sphere diameter is shown in Figure 8. For the bubble with d e = 5.8 mm, the calculated value for a constant mesh is We = 4.05 and We = 4.07 for a dynamic mesh. Both values are slightly lower than the measured ones. This means that the contribution of the inertial force is lower than the actual one in the numerical simulation. . Figure 6. Bubble terminal rise speed. . Figure 6. Bubble terminal rise speed.    The spherical bubble equivalent diameter = 6 ⁄ = 5.8 mm, where V B is the volume of the bubble, is calculated using the Gauss's theorem for closed iso-surface α = 0.5, which represents the gas-liquid interface. The free-rising bubble mode can be determined with the help of Grace's diagram [26], which connects the bubble Eötvos (Eo), Reynolds (Re) and Morton (Mo) dimensionless numbers: According to the Grace's diagram, the bubble from the present case is in the wobbling mode. Figure 9 shows the bubble shape at the time t = 1 s. It is in the region of the wobbling mode (time changing) shape with dynamic mesh, while it is closer to ellipsoidal mode shape with constant mesh. The spherical bubble equivalent diameter d e = 3 √ 6V B /π = 5.8 mm, where V B is the volume of the bubble, is calculated using the Gauss's theorem for closed iso-surface α = 0.5, which represents the gas-liquid interface. The free-rising bubble mode can be determined with the help of Grace's diagram [26], which connects the bubble Eötvos (Eo), Reynolds (Re) and Morton (Mo) dimensionless numbers: According to the Grace's diagram, the bubble from the present case is in the wobbling mode. Figure 9 shows the bubble shape at the time t = 1 s. It is in the region of the wobbling mode (time changing) shape with dynamic mesh, while it is closer to ellipsoidal mode shape with constant mesh.
According to the Grace's diagram, the bubble from the present case is in the wobbling mode. Figure 9 shows the bubble shape at the time t = 1 s. It is in the region of the wobbling mode (time changing) shape with dynamic mesh, while it is closer to ellipsoidal mode shape with constant mesh.

Conclusions
The numerical simulation of the free-rising air bubble in a stationary water column was performed by using the open source CFD toolbox OpenFOAM. The constant and dynamic computational meshes were used and the results of the calculation were compared with the experimental data. The interFoam solver was used for the case with time-constant mesh and the interDymFoam for the case with time-varying mesh.

Conclusions
The numerical simulation of the free-rising air bubble in a stationary water column was performed by using the open source CFD toolbox OpenFOAM. The constant and dynamic computational meshes were used and the results of the calculation were compared with the experimental data. The interFoam solver was used for the case with time-constant mesh and the interDymFoam for the case with time-varying mesh.
A bigger difference between the calculations is evident in the bubble trajectory. In the case of a constant mesh, the deviation amplitude of the helical path is much smaller than in the case of a dynamic mesh due to a poorer description of the gas-liquid interface.
The bubble terminal velocity and the average frequency of oscillation of the bubble movement fits well with the experimental regression curve in both cases, but the matching is slightly better with the dynamic mesh. The Weber numbers partially deviate from the measured value. The wobbling bubble flow regime, which is obtained by the dynamic mesh, closely matches the Grace's diagram, while the bubble shape is closer to the ellipsoidal flow regime in the case of the constant mesh.
From the above-mentioned points, we can conclude that the VOF method is suitable for describing and analyzing the two-phase flow with large (according to the size of the computational cells) bubbles. A more detailed description of the gas-liquid interface is obtained by using a dynamic computational mesh, but we must be aware that this means a significant increase in computing time.