Single bubble rising behaviors in Newtonian and non‐Newtonian fluids with validation of empirical correlations: A computational fluid dynamics study

The bubble rising (BR) dynamic is a common phenomenon in numerous processes of industries. Here, a single air BR behavior is studied using computational fluid dynamics (CFD) modelling in a Newtonian fluid (NF) and non‐Newtonian fluid (nNF). The volume of fluid formulation with the continuum surface force equation is used to track the air bubble in a NF, while the viscosity of the nNF fluid is estimated by using the power‐law equation. The bubble terminal velocity and its shape deformation, as well as the influence of different dimensionless numbers on BR characterization are investigated. The bubble rise in NF, the bubble terminal velocity increases with decreasing Morton number, and bubble moves up in a zigzag way with shape oscillation for the case of low Morton number of the NF. The bubble rise in nNF, the bubble terminal velocity increases with the increases in the rheological index, and bubble size as well as its shape transforms from an ellipsoidal to ellipsoidal cap types. It is found that the drag coefficient is high in a low rheological index compared with the high rheological index. The CFD results are compared with experimental results and empirical correlations reported in the literature. Good agreements are found between the CFD and the literature data for both fluids.

an orifice, free BR characteristics, and its shape deformation) in viscous, viscoelastic, and shear-thinning non-Newtonian fluid (nNF). [7][8][9][10][11] The bubble terminal velocity and its shape deformation in different glycerin solutions were investigated in a wide range of Reynolds numbers and Weber numbers. 8 A correlation was developed by using the Weber number, Morton number, and Reynolds number. A good agreement was found between the experimental and numerical results. Another study of the unsteady motion of a single BR in high viscous quiescent liquid (different concentrations of sucrose solution) was carried out. 9 An empirical correlation was also developed to estimate the drag coefficient by utilizing the Reynolds number, Archimedes number, and Acceleration number. However, BR behaviors were investigated in a viscoelastic fluid. 10 It was found that an intense shear layer formed at the tip of the bubble's rear pole and helps to change the bubble shape. Similarly, a study was performed to investigate single BR characteristics in a viscoelastic polymer solution (Praestol 2500). 12 It was found that the polymer molecules were moved around the bubble from the top to the bottom by aligning its coordinate system with the motion of BR. As a result, the polymer molecules in the bubble lower hemisphere contributed to the BR velocity. Furthermore, the rising behavior of single bubbles was investigated in water and xanthan solution (XGS). 11 The bubbles rose in straight way with tiny oscillations that were convoyed by a steady bubble shape in water and a marginally low rising velocity in XGS. The effect of channel depth, bubble size, superimposed velocity, and rheology of liquid was examined on single BR dynamics in water and XGS. 2 The velocity magnitude, vorticity, strain, and shear stress increased with increasing bubble size and superimposed liquid velocity. Also, the shear stresses acting in XGS were almost two times higher than that of water.
Nowadays, researchers have utilized an advanced computational technique for a single BR study in Newtonian fluid (NF) [13][14][15][16] and nNF. [17][18][19] There are three types of computational model namely the Eulerian-Eulerian, the Eulerian-Lagrangian, and the direct numerical simulation model. The Eulerian-Eulerian model was used when the gas bubble treats as a pseudo-continuum phase so no need to provide flow specifics for distinct bubbles. 20 The Eulerian-Lagrangian model was used for a nondeformable bubble analysis. The direct numerical simulation model is used to analyze the details feature of a single BR dynamic in liquid. In this model, the motion of the bubble-liquid interface is determined by using the front capturing (Eulerian framework system) approach or front tracking (Lagrangian grid system) approach. However, the computational power and time in the front tracking approach are high in comparison to the computational power and time in the front capturing approach, which tracks the bubble-liquid interface by a scalar-indicator or volume of fluid (VOF) fraction that is calculated by solving an advection equation.
Researchers have used the VOF formulation to investigate single or multiple bubbles rise characterization and bubble formation through an orifice at atmospheric pressure [21][22][23][24] or at elevated pressure conditions in viscous liquid. 25 Also, the VOF model was utilized to investigate the BR analysis in membrane bioreactor processes, 26,27 ionic liquids, 27,28 and shear-thinning nNF. 18,19,29,30 For example, the effect of bubble size, column height, and fluid properties on a BR dynamics was investigated in viscous liquid. 31 The increasing viscosity and surface tension of liquid decreased and increased the bubble oscillation, respectively. The dynamics of a sphere BR in a viscous liquid have been investigated for different Eotvos number, Eo. 32 The authors developed five distinct regimes with bubble shape based on the Galilei number (Ga) and Eo.
The VOF with the level set formulation is also used of single BR in carboxymethyl cellulose (CMC) solution for a wide range of Reynolds number and Eo. 33 The bubble deformed from sphere to oblate shapes with increasing Reynolds number and Eo and a symmetrical flow field was also observed in the bubble upper and lower hemisphere. Another study was carried out by using the VOF model to understand the BR dynamics in different nNF, while the Carreau-Yasuda model was used to calculate the liquid viscosity. 18 The bubble deformation increased with increasing Ga and Eo or with reducing the rheological index. The bubble terminal velocity increased with increasing the rheological index and Ga, but it decreased with increasing Eo. Similarly, numerical and experimental studies of a single air BR in CMC solution were investigated and the Carreau-Yasuda model was used for viscosity calculation. 19 It was found that decreasing rheological index and increasing inelastic time constant promoted the bubble path instabilities. Another study found that the viscosity around the bubble decreased significantly as the rheological index of nNF became small so that the bubble rose easily through the shear-thinning fluid. 29 In literature, most of the single BR studies were analyzed by using dimensionless numbers namely Reynolds number, drag coefficient, Weber number, Eotvos number, Galilei number, and dimensionless liquid properties. Very few studies focused on the BR characterization by other dimensionless numbers such as Velocity number, Flow number, Archimedes number, and Capillary number. 9 These dimensionless numbers are also useful for single BR analysis. In this study, single BR characterization is investigated by using those dimensionless numbers in NF and shear-thinning nNF. The VOF formulation with the continuum surface force (CSF) equation is used to track the bubble-liquid interface in the NF, while the power-law equation is used to calculate the viscosity of the shear-thinning nNF. The bubble terminal velocity and its deformation, and the influence of dimensionless number on BR behaviors are examined and compared with experimental data and empirical correlations reported in the literature.

VOF method
The VOF tracks the interface on fixed meshes and is used to track the bubble interface by using the volume fraction function. The transport equation is described as follows: where − → W and L are the liquid velocity and the liquid volume fraction, respectively. The fully liquid and gas phase are defined by using the value of L = 1 and 0, respectively. The interface boundary between the bubble and the liquid is known as 0 < L < 1. Here, the piecewise linear interface construction algorithm 34 is used to solve Equation (1) since it can provide a high precision to reconstruct the bubble-liquid interface. The local-averaged density ( ) and viscosity ( ) of liquid in the interface are calculated as follows 18 : where the subscripts L and G denote the liquid and the gas phase, respectively.

Continuity and momentum equations
The continuity and the momentum equations for incompressible flow are written as follows: where P, t, and g represent the liquid pressure, reference time, and gravitational acceleration, respectively. − → F is the volume force that is computed by using the CSF method developed by Brackbill et al. 35 The equation is expressed as follows: where is the liquid surface tension. k represents the curvature of the bubble-liquid interface and can be described as follows:

Constitutive equation of nNF
In an nNF, the power-law equation is used to calculate the liquid viscosity and is written as follows 17 : where K and f are the consistency coefficient and the power-law index or rheological index, respectively. The local shear rate (̇) can be expressed as follows:̇= where − → is the strain rate tensor and can be described as follows: In Equation (5), the viscosity of the nNF can be estimated by combining Equations (8) to (10). For nondimensional analysis of a single rising bubble, the Morton number, Reynolds number, Capillary number, Eotvos number, Weber number, drag coefficient, and Archimedes number are described as follows: For NF, For nNF, 36 where U T is the bubble terminal velocity. d eq is the bubble equivalent diameter which is calculated as follows 37 : d eq = (d h × d 2 w ) 1∕3 ; here, d h and d w are the shortest distance and longest distance of the bubble, respectively.

Physical model
Literature review shows that several researchers investigated the BR behaviors in a full two-dimensional (2D) computational domain. 28,38,39 Here, a 2D domain is selected since the present study is mainly focused on the BR characterization in a high Mo case. Figure 1 shows a schematic view of the 2D computational domain, which is used to investigate single air BR characterization in NF and shear-thinning nNF. In the initial condition, a sphere bubble with a diameter (d b ) is placed at 10 + 0.

Numerical procedure
The governing equations were solved by using Ansys-Fluent 14.0 41 that is based on the finite volume method. The QUICK upwind and PRESTO schemes were used to solve the momentum equation and the liquid pressure, respectively. The pressure-velocity coupling is solved by using the PISO scheme which can able to give a fast convergence rate with min-

Mesh
A uniform structured with different sizes of the grid such as 0.25, 0.20, and 0.15 mm (which generates total elements number of 6.93 × 10 4 , 1.08 × 10 5, and 1.93 × 10 5 , respectively) are examined for the mesh dependency study. Figure 2 shows the 2.5 mm BR distance with time for the different mesh sizes. For all cases, the CFD results are almost the same up to t = 0.10 second and after that, there is some difference in results until 0.5 second. Very less difference can be seen between the results in the size of the mesh 0.20 and 0.15 mm. While 0.25 mm mesh size produces a low estimation of the BR distance. Considering all, mesh size of 0.20 mm is selected.

Terminal velocity
The bubble terminal velocity is one of the key parameters for BR characterization. Figure 3 shows the instantaneous BR velocity with simulation time for different Morton numbers, Mo. For all cases, the BR velocity increases until 0.2 second and then there is no variation on its velocity up to 0.5 second. It means that the instantaneous rise velocities are in a steady condition. Values of 0.196, 0.162 and 0.063 m/s are found as the terminal velocity for low, medium and high Mo, respectively. The experimental data of 2.5 mm bubble in low Mo (Case 1) is also included. 9 Up to 2.36% relative error (([Experimental-CFD] ÷ Experimental) × 100%) is found by comparing between the CFD and the experimental data. 9 This low error represents a good benchmark for the present CFD methods. The CFD results are also compared with the empirical correlation developed by Rodrigue,42 to predict single BR velocity based on velocity number (V N ) and flow number (F N ). The correlation was developed by using a large range of physical properties and hydrodynamic regimes, and can be described as follows: The velocity number (V N ) and flow number (F N ) are calculated as follows:

F I G U R E 3 The instantaneous rising velocity of a bubble with time for different
Morton number (Mo). Experimental data (d b = 2.5 mm) is taken from Zhang et al. 9 Figure 4 shows a comparison between the CFD results and the correlation (Equation (19)) results. The velocity number increases with increasing flow number and the relative error is less than 10% which is also a reasonable agreement for the present CFD method. Figure 5 shows the BR trajectory for different Morton number, Mo. In a low Mo, the bubble moves up and follows a zigzag path from the left to right. During zigzag motion, the bubble changes its axis (or shape oscillation) from 150 • to 45 • until the top of the domain ( Figure 5B). In this case, a strong vortex of liquid jet acts behind the bubble that rises faster with shape oscillation. The bubble shape oscillation was also observed by other researchers. 7,11 In medium Mo ( Figure 5A), the bubble moves up from the left to right with less oscillation. However, there is no trajectory of the bubble for the case of high Mo because the surface tension force maintains the shape with minimum surface energy. As a result, the bubble aspect ratio is high, and a low terminal velocity is observed. In addition, the bubble shapes are examined with a well-known shape regime chart which was developed by using a vast amount of experimental data and theoretical analysis. 43 Figure 6 shows the bubble shape regime chart with Reynolds number (Re) vs Eotvos number (Eo) for a wide range of Mo. The computed bubble shapes agree well with the shape regime chart for all cases.  43 and the contour shows the computed bubble shape at t = 0.4 second. Note: s, oe, oed, oec, sc, and sk denote the spherical, oblate ellipsoidal, oblate ellipsoidal-disk, oblate ellipsoid-cap, spherical-cap, and skirted type shape of a bubble, respectively Also, the bubble shape can be examined with the Weber number (We) which is the ratio of dynamic pressure to the surface tension. When We < 1, the bubble maintains a spherical shape that alters to an oblate ellipsoidal or oblate ellipsoidal disk shape due to the velocity pressure at We > 1. 44 For example, We of 0.123 and 1.51 represent a spherical and oblate ellipsoidal shape for low and high Mo, respectively (refer to Figure 6). The calculation of the We is compared with available empirical correlation developed by Raymond & Rosant. 8 The authors developed the correlation with the range of Mo = O (10 −8 )-O (7) and can be expressed as follows: Figure 7 shows the Weber number (We) as a function of Reynolds number (Re) for different Mo. It is found that We increases with increasing Re. Results from CFD match well with correlation Equation (22) and less than 5% relative error is obtained. Similarly, Cai et al. 45 investigated a single BR behavior in high viscous liquid and a wide range of Mo = O (10 −11 )-O (163). The authors developed a relationship between the bubble aspect ratio (the ratio of shortest distance to longest distance of the bubble, d h /d w ) and dimensionless numbers such as Re, Eo, and Mo. The equation is described as follows: Figure 8 shows the bubble aspect ratio vs the Reynolds number for different Mo. The bubble aspect ratio decreases with increasing Re. At low Re < 3, the bubble contains its original shape which becomes an ellipsoidal shape due to the surface tension force and the buoyancy when Re > 3. The computed results are well fit with correlation (Equation (23)) and up to 5% relative error is found.

Drag coefficient
The drag force is another important parameter for BR characterization. The computed results are examined with Tomiyama et al. correlation 46 which is expressed as follows:  Figure 9 illustrates the drag coefficient vs Reynolds number. The drag coefficient decreases with increasing Re and good agreement can be seen between the computed results and the correlation results (Equation (24)). Also, Zhang et al 9 experimental data are included and up to 8% relative error is found by comparing the correlation and the experimental data. Similarly, Zhang et al 9 developed a correlation to estimate the drag coefficient based on bubble Archimedes number (Ar) and Re. The authors reported that if the bubble velocity reached into a steady condition which represented the Acceleration number is equal to zero, that is same for this present work. The correlation for the steady condition is expressed as follows: ) .
(25) Figure 10 shows the drag coefficient vs Archimedes number. The drag coefficient decreases with increasing Archimedes number. An overprediction (up to 5%) can be seen by comparing with Zhang et al.'s estimations.

Terminal velocity
In this section, single BR dynamics are examined in the shear-thinning nNF. Three distinct sizes of bubble, d b = 10, 12 and 17 mm are considered in XGS of 0.025% and 0.10% concentration solution, corresponding to rheological index (f ) of 0.8248 and 0.5481, respectively. Figure 11 shows the bubble terminal velocity vs bubble diameter for f = 0.8248 and 0.5481. The  Figure 12 illustrates a BR trajectory with shape deformation (like as dimpled ellipsoidal cap shape) which happens at the beginning of the rising stage at t = 0.05 second. The changes in shape occur due to the unbalanced pressure difference, buoyancy force, and drag force. However, the bubble deformation and swing increase with increasing bubble size and the large bubble experiences more resistance by the nNF as well as a large deformation due to the effect of strong wake at bubble lower surface. In Figure 12A (d b = 10 mm), the bubble shapes are close to a spherical-cap which changes to wobbling shape when d b is increased to 12 mm ( Figure 12B) at t = 0.2-0.5 second. In Figure 12C (d b = 17 mm), the dimpled ellipsoidal cap shape (t = 0.05 second) changes to skirted type shape (t = 0.1 second) and then the edge of the bubble breakup into small bubbles and form ellipsoidal cap shape bubble (t = 0.2-0.5 second).

Drag coefficient
In a non-Newtonian power-law fluid, the Reynolds number based on the aspherical bubble with a vertical symmetry axis can be calculated as follows 17 : Figure 13 shows the drag coefficient vs Reynolds number for two different rheological indexes of 0.8248 and 0.5481. The drag coefficient decreases with increasing Reynolds number. The lower drag coefficient can be seen for a high rheological index (f = 0.8248) than that of the low rheological index, f = 0.5481. The CFD results are also compared with the available correlation: For Up to 12% relative error is found by comparing between the CFD results and these correlation results.

CONCLUSION
In this study, single air BR behaviors in an NF and nNF are investigated by using the VOF with the continuum surface force formulations. The power-law equation is used to estimate the viscosity of the nNF. The bubble terminal velocity increases with decreasing Morton number, and bubble moves up in a zigzag way with shape oscillation for the case of low Morton number of NF. In the nNF, the bubble terminal velocity increases with the increases of the rheological index and the bubble size as well as its shape changes from an ellipsoidal to ellipsoidal-cap. Also, the drag coefficient is high in a low rheological index compared to the high rheological index. Finally, the computed results show good agreements with experimental results and empirical correlations reported in the literature.

CONFLICT OF INTEREST
The authors declare no conflict of interest.