A pair of bubbles’ rising dynamics in a xanthan gum solution: a CFD study

Abstract

The motion and interaction of a bubble pair in a non-Newtonian fluid (xanthan gum solution) were numerically simulated using a volume of fluid (VOF) method, in which the continuous surface tension model and the power-law model were adopted to represent the surface tension and rheological properties of non-Newtonian fluids, respectively. The effects of the initial horizontal bubble interval, oblique alignment and rheological properties of non-Newtonian fluids on a pair of bubbles rising side-by-side are evaluated in this study. The results indicate that for the case with a non-dimensional initial horizontal bubble interval h* = 4.0, the interaction between the bubbles shows a minimal repulsive effect. Moreover, for the oblique angle alignment a greater repulsive force between the bubbles was seen when the angle was reduced. However, oblique coalescence occurred due to the higher attraction between the bubbles at higher angles, which is independent of the flow index. It is also found that the repulsion effect as well as the variation of the bubble shape from spherical to irregular are more significant at a lower flow index (n < 0.5) due to the shear-thinning effect as well as the differences in their flow field structures.

---

1. Introduction

A bubble column reactor is a device in which gas is usually dispersed as bubbles which rise through a liquid. The overall performance of a bubble column mainly depends on the bubble flow characteristics, bubble coalescence and bubble break up phenomenon which alters the bubble size distribution and enhances the gas–liquid contact area by the action of vortices via stretching, tearing etc., resulting in a significant increase in heat and mass transfer and chemical reaction rate.^1–3 Nowadays, bubble column reactors are widely used as gas–liquid interface equipment in various industrial sectors such as chemical, petrochemical, pharmaceutical etc., in reactions such as oxidation, chlorination, alkylation, polymerization and hydrogenation, and are employed as bioreactors to produce valuable products like enzymes, proteins, antibiotics, etc.^1,2 Recently, a number of experimental and theoretical studies^4–10 have been carried out on the formation and rise characteristics of single bubbles in stagnant fluids under the buoyancy force in Newtonian fluids. The literature provides almost complete information regarding single bubble formation and rise characteristics. However, in industrial environments, swarm bubbles or multiple bubbles are more commonly encountered than the single bubble. Thus, the obtained information on the single bubble might not be broadened to multiple bubble systems. On the other hand, most of the materials encountered in both nature and industry are non-Newtonian fluids. Compared to Newtonian fluids, non-Newtonian fluids usually show many peculiar properties, for instance, shear-thinning, viscoelasticity, rod-climbing and tubeless siphoning.^11–13 Hence, the bubble rising behavior and interactions in non-Newtonian fluids are more complex and also the relevant research is quite rare. Therefore, the fundamental knowledge of bubble rising behavior and interactions in non-Newtonian fluids could be improved to more scientifically confirm the modeling of gas–liquid flow in bubble columns.^14,15

Most experimental studies have focused on single or in-line bubble rise dynamics in non-Newtonian fluids.^16–19 For example, Hassagar¹⁶ investigated the rise dynamics of two in-line bubbles in non-Newtonian fluids. It was found that a negative wake encouraged by elasticity pushes the liquid away from the bubble. Lin et al.¹⁷ also found that for two in-line bubbles, the acceleration of the trailing bubble to the leading bubble is caused by negative pressure; and the shear-thinning effect in addition to the pushing force is caused by viscoelastic effects. Hassan et al.²⁰ investigated the zigzag trajectory of bubbles with maximum amplitude for smaller bubbles in low xanthan gum concentration solutions.

Few attempts have been made to investigate experimentally the interaction between bubble pairs rising side-by-side in non-Newtonian fluids.^21–24 Sanada et al.²¹ experimentally investigated the motion of a horizontally aligned pair of rising bubbles in silicone oil. They found that the velocities of the bubbles decrease after coalescence by as much as 50% when the bubbles rising side-by-side bounce off each other. The bubble repulsion effect was observed due to a large amount of fluid in the space between the bubbles, which resulted in the bubbles’ movement in a direction away from each other. Vélez-Cordero et al.²³ experimentally investigated the interaction of two bubbles rising in shear-thinning fluids of a xanthan gum solution (0.55 < n < 1.0). They observed that the attractive motion between the bubbles was increased with the amount of shear-thinning (decreasing the flow index). The bubble pairs showed an oscillatory motion due to the reduced viscosity behind the leading bubble. Fan et al.²⁴ focused on the rise and interaction between two parallel rising bubbles by analyzing the velocity field around the bubbles using particle image velocimetry (PIV), and found that within a certain distance between two bubbles, the interaction between two neighboring bubbles changes from mutual repellence to attraction with decreasing the angle of the line which links the two bubbles’ centers to the vertical direction. Legendre et al.²² investigated the rising behavior of two parallel spherical bubbles in a viscous fluid. They reported that the hydrodynamic interactions between the bubbles would be cohesive or repulsive, primarily depending on the Reynolds number.

Recently, more and more researchers^25–35 are making use of various numerical methods such as the volume of fluid method (VOF), level set method (LS), lattice Boltzmann method (LB), and front tracking method (FT) to investigate bubble dynamics or interactions. In general, their results are in reasonable agreement with the existing experimental data. Fan et al.²⁵ successfully analyzed the dynamics of two bubbles rising side-by-side in concentrated carboxymethyl cellulose (CMC) solution by using the VOF method. The results showed a good agreement with the experimental measurements. They found that the repulsive effect between two bubbles decreases with increasing the initial center-to-center distance of the bubbles and an increase of the oblique angle between them. Similarly, Yu et al.²⁶ investigated two parallel bubbles’ rising behavior in a viscous fluid using an adapted LB method. The authors found that the repulsive behavior of two spherical bubbles occurs at lower Reynolds numbers, but cohesive behavior and finally coalescence of the bubbles occur at higher Reynolds numbers. Li Zhang et al.²⁸ studied the motion of a single bubble rising freely through a CMC sodium salt, sodium hydroxyl-ethyl cellulose (HEC) and xanthan gum (XG) solution using a level set method for tracking the bubble interface. They investigated the shear rate and viscosity distribution and shape of a bubble rising in the CMC, HEC, XG solution and compared them to those in a sodium acrylate polymer (SAP) shear-thinning solution. Liu et al.³² performed a numerical investigation on three parallel bubbles at equal intervals rising in CMC non-Newtonian fluids. The governing equations were solved using the VOF method. It was seen that the three parallel bubbles coalesced when the horizontal interval between them was less than 1 mm, otherwise the bubbles would experience a repulsive effect.

In summary, a detailed study of a pair of bubbles’ rising dynamics in a pure concentrate of xanthan gum solution (i.e., non-Newtonian) and their change in shape is not available to the best of the authors’ knowledge. There are a few works in which xanthan gum mixture solutions (xanthan–glycerin/water) are tested, which have different physical properties from the pure concentrate of xanthan gum solution, e.g. see ref. 20 and 23. The present study provides detailed information on a pair of bubbles’ rising dynamics and shows how the shape of the bubbles evolves with time in a pure concentrate of xanthan gum solution. Specifically, a pair of 6 mm bubbles’ rising dynamics with different initial configurations in such a concentrate of a non-Newtonian fluid are studied. The effects of initial horizontal bubble interval distance, non-horizontal configuration, and rheological properties of shear-thinning fluids are simulated in detail using the VOF method to study the interactions between bubble pairs and the flow field structures and pressure distribution around the bubbles. The knowledge from the present study could be useful for a multi-scale approach to predict bubble swarm behavior as well as to improve the understanding of the mechanisms of multiple bubble rising dynamics.³¹

2. Numerical methods

In this work, the continuous surface tension model and the rheological properties of non-Newtonian fluids are incorporated into the VOF method for a pair of bubbles rising in shear-thinning fluids.

2.1 Physical model

The physical model of the simulation was simplified to a 2D computational domain with dimensions of 120 mm height and 90 mm width, as shown in Fig. 1a. The bottom and two sides of the domain are assigned as no-slip wall boundary conditions. The top of the domain is assigned as a pressure outlet boundary condition. The operating pressure is set to be equal to ambient pressure, i.e. 101 325 Pa, and a gravitational force (g) of −9.81 m s⁻² is assigned along the y-axis. At the initial stage of simulation, a pair of bubbles with diameters of 6 mm is imposed at the centre and 12 mm height from the bottom of the domain. The bubbles in quiescent non-Newtonian fluids rise under the action of buoyancy and the bubbles’ rising dynamics are numerically investigated. The effect of the column wall on the bubble is negligible since the bubble size is small in comparison to the column width. In the present study, a xanthan gum solution is used as the non-Newtonian fluid. The non-Newtonian fluid data are taken from Hassan et al.²⁰ The non-Newtonian fluid has the following properties: density, ρ_l = 997.0 kg m⁻³ and surface tension, σ = 0.063 N m⁻¹. Rheological data of the consistency coefficient K = 0.095 Pa sⁿ and flow index n = 0.548 were used to fit the power-law model.

Fig. 1 (a) Physical model of the computational domain and (b) uniform structured grid.

2.2 Governing equations

2.2.1 Equations of continuity and momentum. The continuity and momentum equations for an incompressible fluid can be written as:³⁶

∇·V⃑ = 0 (1)

where ρ is the fluid density; V⃑ is the velocity vector of the fluid; P is pressure; F⃑_s is the body force; μ is the dynamic viscosity coefficient. The strain rate tensor, D⃑ is written as below:

The local averaged density ρ(F) and kinematic viscosity coefficient μ(F) are evaluated from the local distribution of the phase volume function F:

ρ(F) = ρ_l(F) + ρ_g[1 − F] (4)

μ(F) = μ_l(F) + μ_g[1 − F] (5)

where ρ_l and ρ_g, are the density of the liquid and gas, and μ_l and μ_g are the viscosity of the liquid and gas, respectively. The volume fraction (F) is defined as the fraction of the liquid inside a control volume or cell, in which F takes the value of 0 for a pure gas cell; 1 for a pure liquid cell and between 0 and 1 for an interface of gas and liquid in the cell. The volume fraction equation is defined as follows:³⁷

2.2.2 Source term of the momentum equation induced by surface tension. Surface tension has an important impact on the interface because the minor curvature of a bubble could generate major additional pressure. The continuum surface force (CSF) is used to calculate the gas–liquid interface motion, which is incorporated as a source term in the momentum equation (eqn (2)) by introducing a body force F⃑_s as described by Brackbill et al.³⁸ This body force is calculated by the following equation:where, F_l is the liquid phase fraction and k is the surface curvature of the interface, which is defined in terms of divergence of the normal vector, n̂ and is calculated by using the following equation:

2.2.3 Constitutive equation of continuous phase. In non-Newtonian fluids, the shear-thinning effect is presented by the power-law model:³⁹

μ(F) = K[small gamma, Greek, dot above]ⁿ⁻¹ (9)

where K and n are the consistency coefficient and flow index of shear-thinning fluids, respectively. The viscosity μ can be obtained from the local shear rate [small gamma, Greek, dot above] which can be written as below:

Therefore, the viscosities of non-Newtonian fluids can be obtained by the combination of eqn (3), (9) and (10).

2.3 Model validation

2.3.1 Grid analysis. A uniform structured mesh is used for the present study, as shown in Fig. 1b. Fig. 2 shows the effect of the mesh interval sizes (0.15, 0.25, 0.35 and 0.45 mm) on the simulation results of the single bubble rising velocities with an initial diameter of 6 mm in the xanthan gum solution. The density, surface tension, consistency coefficient and flow index of the power-law fluid were 997.0 kg m⁻³, 0.064 N m⁻¹, 0.095 Pa sⁿ and 0.548, respectively.²⁰ It could be found that the single bubble rising velocities at a mesh interval size of 0.25 mm are almost the same as those at 0.15 mm. Therefore, a mesh interval size of 0.25 mm was adopted throughout this study to take into account both the computational accuracy and the time consumption. Similarly, satisfactory simulation results were obtained between the current study and the one carried out by Ma et al.⁴⁰ on the influence of different factors on single bubble formation and bubble dynamics with a mesh interval size 0.25 mm.

Fig. 2 Effect of mesh interval size on the rising velocity of a 6mm diameter bubble.

2.3.2 Code validation. In order to validate the reliability of the computational method, two types of study have been considered. First, the processes of two parallel bubbles rising in carboxymethyl cellulose sodium (CMC) solution were simulated and the results agreed well with the experimental results,³² as depicted in Fig. 3a and b. The density, surface tension, rheological parameters of the consistency coefficient and flow index of the CMC aqueous solution are 1005.6 kg m⁻³, 0.06875 N m⁻¹, 0.048 Pa sⁿ and 0.922, respectively. For the second validation, single bubbles with initial diameters of 6 mm and 10 mm in the present of non-Newtonian fluids (i.e., xanthan gum solution) were simulated and compared with the experimental results reported in ref. 20. The relative errors of the Reynolds numbers were about 2.43% and 1.22% for the initial diameters of 6 mm and 10 mm, respectively. The Reynolds number is calculated as:^32,41,42where d_e = [(d_h × d_w²)^1/3]²⁰ and U_T are the equivalent diameter and the bubble terminal velocity, respectively. The equivalent diameter (d_e) was calculated using the short bubble diameter (d_h) and larger bubble diameter (d_w) of the elliptic bubbles after they reached steady conditions, or once there was no more variation of the instantaneous bubble velocity.

Fig. 3 Comparisons of the height and shape of two parallel bubbles between (a) experimental shape of bubbles from Liu et al.³² and (b) the present simulation.

The relative error of the Reynolds number was found to be about 3.76% for the two parallel bubbles. Therefore, the above comparison results indicated that the VOF computational method can be reliable for predictions in the present investigation (Table 1).

Table 1 Comparison of equivalent diameter (d_e); bubble rising velocities (U_T) and Reynolds numbers between experiments and the present simulation

	de (mm)	UT (m s^−1)	dh/dw	Re [—]	Error of Re (%)
Two parallel bubbles
Experimental^32	4.1034	0.181	0.648	20.96	3.76
Simulation	4.213	0.172	0.636	20.17
Single bubble
Experimental^20	6	0.231	[—]	77.75	2.43
Simulation	6.452	0.225	0.683	75.86
Experimental^20	10	0.242	[—]	107.32	1.22
Simulation	10.824	0.237	0.562	108.63

---

2.4 Numerical procedures

In this study, the CFD code FLUENT was used to solve the governing equations employing a pressure based solver. The pressure–velocity coupling equation was solved using the pressure implicit with the splitting of operators (PISO) algorithm. The geometric reconstruction approach was adopted to track the interface between two phases using a piecewise linear interface calculation (PLIC) method.⁴³ The discretization schemes of pressure and momentum were the pressure staggering option (PRESTO!) and second order upwind, respectively. The iteration and time steps were 1.0 × 10⁻⁶ and 0.0001 s, respectively. At the beginning, initial spherical bubbles were patched at the bottom of the computational domain containing only quiescent liquid.

2.5 Simulation cases

A total of 17 simulation cases were carried out. The first two cases were carried out for validation purposes, in which bubbles with diameters of 4 mm (case 1) and 6 mm (case 2) were set to rise from a resting condition with the initial position of 12 mm vertical height from the bottom of the computational domain. Three cases (cases 3–5) were used to investigate the effect of non-dimensional horizontal intervals (i.e., h* = 1.5, 2.0 and 4.0) between a pair of bubbles with d_o = 6 mm. The non-dimensional horizontal interval is defined as h* = X_i/d_o, where X_i is the initial distance between the centres of the bubbles and d_o is the initial bubble diameter. Cases 6–8 were used to investigate the effect of the flow index (n) of the non-Newtonian fluid between the bubble pair on the rising dynamics. In cases 9–17, the effect of the oblique angle (θ) and flow index (n) on the fluid flow were studied. A summary of all of the simulation cases is given in Table 2.

Table 2 Simulation cases

Case	ρl [kg m^−3]	K [Pa s^n]	n [—]	h* [—]	θ [°]	Purpose
1	1005.6	0.048	0.922	2.0	0	Validation of the simulation with the available literature^20,32
2	997	0.095	0.548	0.0
3–5	997	0.095	0.5	1.5	0	To investigate the effect of the horizontal interval between a pair of bubbles
				2.0
				4.0
6–8	997	0.095	0.1	2.0	0	To investigate the effect of the flow index (n) on the dynamics of a pair of rising bubbles
			0.5
			0.9
9–17	997	0.095	0.1	2.0	10	To investigate the effect of the oblique angle and flow index (n) on the dynamics of a pair of rising bubbles
			0.5		22.5
			0.9		45

---

3. Results and discussion

3.1 Effect of different initial horizontal intervals between a pair of bubbles

A bubble pair’s rising dynamics in a horizontal direction are closely related to the initial horizontal bubble interval and the physical properties of non-Newtonian fluids. The interaction of bubble pairs is investigated with three different non-dimensional initial horizontal intervals h* = 1.5, 2.0 and 4.0. Fig. 4 shows the 6 mm bubble pair’s rising trajectories with three different initial bubble intervals. The results reveal that the path, shape and velocity experience a significant variation with increasing distance between the centres of the two bubbles. The bubble shape changes from spherical to irregular for the three different initial bubble intervals. The bubble pair’s rising trajectory is asymmetric along the perpendicular line in the middle of the vertical column, as can be seen from Fig. 4. For h* = 1.5 (see Fig. 4a) and h* = 2.0 (see Fig. 4b), a stronger repulsive interaction between the two parallel bubbles is observed compared to the interval of h* = 4.0 (see Fig. 4c). Consequently, the results for bubble horizontal interval ratio, bubble aspect ratio (or the ratio of the minimum to maximum deformation of a bubble), bubble rising velocity and static pressure on the bubble as a function of time at three different initial bubble intervals are shown in Fig. 5a–d, respectively.

Fig. 4 Two parallel bubbles’ rising trajectory at different initial intervals when the flow index, n = 0.5, (a) h* = 1.5; (b) h* = 2.0; (c) h* = 4.0.

Fig. 5 (a) Bubble interval ratio; (b) bubble aspect ratio; (c) bubble rising velocity and (d) static pressure as a function of time at different initial bubble intervals.

According to Fig. 5a, the highest variation in the bubbles’ rising trajectory in the horizontal direction occurs for the low initial bubble intervals of h* = 1.5 and 2.0. This is due to the stronger effect of repulsion as a result of a large amount of vortices generated between the bubble intervals. While in the case of h* = 4.0, the bubble rising trajectory in the horizontal direction fluctuates less, indicating the weak repulsive effect. These results are in agreement with the experimental work of Vélez-Cordero et al.²³ Consequently, from Fig. 5b, for h* = 1.5, the curve of the bubble aspect ratio fluctuated intensely, but less fluctuation was observed in the h* = 2.0 curve. When the initial interval increased to h* = 4.0, the bubble shape did not vary so significantly compared to the other non-dimensional intervals. The change of bubble aspect ratio occurs due to the stronger vortex field lying in the gap between the bubble pair that keep interacting with each other and also the interacting force in the horizontal direction which varies with time. Meanwhile, it is also worth noting that the bubbles’ rising velocity depends on the distance between the bubbles, as can be observed in Fig. 6a and b. As can be seen in Fig. 5c, the bubble pair with the initial interval h* = 4.0 was observed to reach a terminal velocity of 0.16 m s⁻¹, which is higher than those with h* = 2.0 and 1.5. This means that the vertical motion of the bubbles with smaller initial intervals is weakened by their repulsive effect in the horizontal direction, and the terminal velocity of the bubble pair would increase with the increase of the initial bubble interval due to a weaker vortex field between the bubble gap, as can be observed in Fig. 6c. So, it is reasonable to conclude that bubble–bubble interaction exhibits a repulsive effect, which decreases with the increase of distance between the bubbles and the repulsive effect can be considered negligible at h* ≥ 4. The effect of initial bubble interval on the bubble static pressure is calculated as P = (H − h)ρg, where H is the height of the bubble column and h is the height to which the bubble had risen. The results are shown in Fig. 5d. It is found that the static pressure on the bubble decreases with an increase of h* and time. In comparing the value of h* = 1.5 and 4.0, on average, the drop in static pressure is about 30%.

Fig. 6 Velocity flow field around the rising bubble pairs at t = 0.05 s and n = 0.5; when (a) h* = 1.5; (b) h* = 2.0; (c) h* = 4.0.

3.2 Effect of flow index on a pair of bubbles’ rising dynamics

The flow index n reflects the deviating degree of non-Newtonian fluids from Newtonian fluids (n = 1). Most non-Newtonian fluids are pseudo plastics with 0 < n < 1. If n is less than one, the power law (eqn (9)) predicts that the effective viscosity would decrease with an increasing shear rate. The simulation was performed by setting different flow indexes; however the other properties such as air bubble density, density of the liquid, surface tension, column shape and the consistency of the fluid, were kept constant. In this study, three different flow indexes (n = 0.1, 0.5 and 0.9) were taken into consideration to see the effect of the flow index on the bubbles’ rising dynamics in non-Newtonian fluids. Fig. 7a and b display the rising trajectory of two bubbles with 6 mm diameter rising in a non-Newtonian fluid at a low flow index (n = 0.1) and high flow index (n = 0.9), respectively. As shown in Fig. 7a, when the bubbles with an initial diameter of 6 mm and a flow index of n = 0.1 rose in the column, the bubbles’ shape changed from spherical to irregular and they began to deviate from each other in the horizontal direction. Whereas for the case of n = 0.9, the bubbles’ shape changed from spherical to ellipsoidal and the rising trajectory showed a nearly straight path (see Fig. 7b).

Fig. 7 Trajectory of the 6 mm diameter bubbles at (a) low flow index, n = 0.1 and (b) high flow index, n = 0.9.

Fig. 8a and b illustrate the velocity fields around the bubbles at low flow index (n = 0.1) and high flow index (n = 0.9), respectively. For n = 0.1, at time t = 0.05 s, it is observed that a stronger circulation of vortexes between the bubbles push at the bubbles’ edges. As a result, the bubbles’ shapes change to an oblate ellipsoidal disk (oed) shape at t = 0.1 s. When the time is increased to 0.2 s, the oed shape again changes to an irregular shape, which is attributed to the stronger shear-thinning effect of non-Newtonian fluids. Note that an excess vortex is seen at t = 0.2 s. This excess vortex may be the primary reason for the deviation of the bubbles from each other in the horizontal direction (see Fig. 7a at t = 0.2–0.5 s). Whereas, weak circulations occur between the bubbles for the case with a flow index of n = 0.9, which are not able to push at the bubbles’ boundaries. Therefore, the shape of the bubbles does not change too much, as can be seen in Fig. 8b.

Fig. 8 Velocity field around the bubbles at different times when h* = 2, (a) low flow index, n = 0.1; (b) high flow index, n = 0.9.

Fig. 9a–d show the bubble horizontal interval ratio, bubble aspect ratio, bubble rising velocity and static pressure on the bubbles, respectively at three different flow indexes of n = 0.1, 0.5 and 0.9. In addition, the results of Fan et al.²⁵ for a pair of 6 mm diameter bubbles’ rising behavior when h* = 1.7, which are close to that in the present study, for a carboxymethyl cellulose (CMC) solution of n = 0.904 are also included in Fig. 9a. Comparing the results obtained by Fan et al.²⁵ with the present study, on average the results are quite close to one another; some differences are observed and this may have resulted from the difference in the properties of the non-Newtonian fluids used, whereby Fan’s ²⁵ study was based on a CMC solution and the present study is based on a xanthan gum solution. For example, in comparing n = 0.904 of Fan et al.²⁵ and n = 0.9 from the present study, the maximum difference is about 9%. However, the horizontal interval ratio between the bubbles remains unchanged until t = 0.25 s for n = 0.1–0.9. After that, the values of the interval ratio of the non-Newtonian fluid with n = 0.1 increases dramatically, indicating high repulsion effects between the bubbles. However, much lower values were obtained for the fluid with n = 0.5 and 0.9, along with a similar trend for them.

Fig. 9 (a) Bubble interval ratio and data for the CMC solution from Fan et al.²⁵ is also included; (b) bubble aspect ratio; (c) bubble velocity and (d) static pressure as a function of time at different flow index n, when h* = 2.0 and d_o = 6 mm.

Additionally, the effect of change in the flow index on the bubble aspect ratio (or bubble shape) is shown in Fig. 9b. It is observed that the curve fluctuates severely with n = 0.1, but less fluctuation was observed at high flow index n = 0.5. When the flow index increased to n = 0.9, the variation in the bubble shape was not so remarkable, compared to the case of n = 0.1. Changes of the bubble shape occurred due to the development of a stronger push of fluid jet around the bubbles in the gap between the bubbles (see Fig. 8a). Accordingly, Fig. 9c shows that the bubbles’ rising velocity increases with a decrease in n. This is may be related to the shear-thinning effects in non-Newtonian fluids, which reduced the viscosity of the fluid around the gap between the bubbles, thus the rising resistance of the bubbles was faster in both the vertical and horizontal directions.

Fig. 9d shows the effect of the flow index on the bubble static pressure with time. It has also been found that the static pressure decreases with an increase of ‘n’ and time. Comparing three different n, on average, the results are relatively close to each other; but some differences are seen from t = 0.2 s to 0.5 s and this may have resulted from the start of strong shear-thinning. For example, comparing n = 0.1 and 0.9, the maximum difference is about 16% on average from t = 0.2–0.5 s. Subsequently, Fig. 10 shows the normalized static pressure as a function of the bubble position in the column. The static pressure and bubble position are normalized by atmospheric pressure and column height, respectively. The results show that the normalized static pressure decreases with the height of the bubble in the column and the maximum difference is found to be less than 6% for the cases of different flow index (n) and different initial bubble interval (h*).

Fig. 10 Normalized static pressure as a function of the bubble position at different flow index n and different initial bubble interval.

3.3 Effect of oblique angle θ between bubble pairs at different flow indexes

In this section, bubble pairs with non-horizontal configurations are investigated. The simulation results reveal that in the rising process, the bubble path and bubble interval experience a huge variation at different oblique alignments. The oblique configuration between two bubbles is shown in Fig. 11. The contours of the bubbles’ shape with the variation of time for three initial angles θ = 10°, 22.5° and 45° and n = 0.1 and 0.9, are shown in Fig. 12a–d. By considering the identical horizontal distance of h* = 2.0, it is manifested that when the configuration angle θ increases, the interaction between the two bubbles alters from repellence to attraction. This behavior can be observed with the comparison of Fig. 12a and b. Finally, the bubbles coalesce by increasing θ to 45° for both flow index n = 0.1 and 0.9 of the non-Newtonian fluid, as displayed in Fig. 12c and d, respectively. The velocity fields around the bubbles at the beginning stage for different θ and low flow index n = 0.1 are displayed in Fig. 13a–c. From Fig. 13a, for a small angle (θ = 10°), when the bubbles begin moving, the lower surface of the bubbles is pushed up by vortexes. The upper portion of the following bubble encounters a negligible wake developed behind the lower surface of the upper bubble. As a result, the upper bubble has experienced a stronger pushed-away effect by the following bubble. Thus, both bubbles’ shapes deform significantly and their rising paths deviate gradually from the original vertical line, resulting in enlargement of the bubble distance. The pushed-away effect gradually gets weaker (or the wake effect gradually gets stronger) behind the upper bubble with increasing θ from 22.5° (see in Fig. 13b) to 45° (see in Fig. 13c).

Fig. 11 Schematic of the configuration of a bubble pair; noted: FB and UB mean following bubble and upper bubble, respectively.

Fig. 12 Bubble rising trajectory at different angles when h* = 2, (a) β = 10°, n = 0.1; (b) β = 22.5°, n = 0.1; (c) β = 45°, n = 0.1 and (d) β = 45°, n = 0.9.

Fig. 13 Velocity flow field around the rising 6 mm diameter bubbles at n = 0.1 and t = 0.05 s, when the oblique angle (a) θ = 10°; (b) θ = 22.5°; (c) θ = 45°.

The bubble horizontal interval ratio and the change of the angle between the centre line of the bubbles and the horizontal direction at three different initial angles (θ = 10°, 22.5° and 45°) with three different flow indexes (n = 0.1, 0.5 and 0.9) are illustrated in Fig. 14 and 15, respectively. According to Fig. 14a and b, at low flow indexes (i.e., n = 0.1 and n = 0.5), the fluid generates a stronger repulsive effect for the cases with θ = 10° and 22.5°, respectively. Consequently, higher values of the interval ratio are obtained compared to the highest flow index (n = 0.9). The stronger repulsive effect is generated in the horizontal direction due to the stronger shear-thinning effect. In addition, the results of Fan et al.²⁵ for an initial angle of 17° when h* = 1.7, which is close to that in our simulation of 22.5°, are also included in Fig. 14b. Comparing between Fan et al.²⁵ and the present study, some differences are seen because of the differences in the non-Newtonian fluid properties as discussed earlier. For example, in comparing n = 0.904, θ = 17° of Fan et al.²⁵ and n = 0.9, θ = 22.5° from the present simulation, the maximum difference is about 13%.

Fig. 14 6 mm diameter bubble interval ratio as a function of time at different flow indexes, n; (a) θ = 10°; (b) θ = 22.5°; data for the CMC solution from Fan et al.²⁵ is also included for θ = 17°; (c) θ = 45°.

Fig. 15 Variation of the angles as a function of time at different flow indexes, n; (a) θ = 10°; (b) θ = 22.5°; (c) θ = 45°.

When θ = 45° (see Fig. 14c), the decreasing trend of the curves represents the weaker repulsive effect for the three different flow indexes (n). Additionally, a value of the interval ratio (ΔX/X_i) of zero for the three different n represents the coalescence of the bubbles. Our results show noticeable differences with the results of Vélez-Cordero et al.,²³ who found that the bubbles have an attractive and repulsive behavior at θ = 42° for n = 0.5–0.8. Moreover, as θ increased to 61° and 78°, the bubbles coalesced for n = 0.5 and 0.7, but no coalescence can be found for n = 0.8. In contrast to these results, our study has showed bubble coalescence only at θ ≥ 45° for n = 0.1–0.9. This may be explained by the relatively low physical properties of the pure xanthan gum solutions that facilitate the attractive behaviour of bubble coalescence, while the higher physical properties of mixture solutions (i.e., xanthan gum and glycerine/water) showed both attraction/coalescence and repulsive phenomena.²³

Fig. 15a–c show the variations of the configuration angle against time for various flow indexes. For the first two graphs (θ = 10° and 22.5°), the curves related to both flow indexes of n = 0.1 and 0.5 experience a sudden drop after a time step at around t = 0.25 s, indicating that the shear-thinning effect is strong enough to change the bubble angle configuration, while smaller changes are observed for the curves with n = 0.9. However, for the initial angle of θ = 45°, the configuration angle curve increases steeply for the flow index n = 0.9 and 0.5 with no fluctuation (see in Fig. 15c). Moreover, as indicated above, a fluctuating configuration angle is seen for n = 0.1 due to its higher shear-thinning effect. It is worth mentioning that the end of each curve represents the bubble coalescence time at approximately t = 0.30 s for θ = 45° (see Fig. 15c). By comparing with the results for small angles, there always exists a repulsive effect between the bubbles in both directions as well as the variation of the angle for the low flow index of a 0.1–0.5 non-Newtonian fluid; whereas for large angles, the attraction between bubbles is stronger than the repulsive effect. In the other words, the bubble repulsion effect can be controlled by using initial oblique angle variation.

The velocity flow fields, outside and inside the pressure profile on the bubbles before the start of coalescence (t = 0.30 s) for low flow index n = 0.1 and high flow index n = 0.9 are also displayed in Fig. 16a and b, respectively. For both flow indexes, the horizontal attraction between the bubble pair increases which results in the bubbles moving closer to each other and subsequently coalescing at t = 0.30 s. It is noted that in case of n = 0.l (Fig. 16a), the upper surface of the trailing (or following) bubble moves more quickly due to the wake (or low pressure region) created by the upper bubble (refer to the outside pressure profile in Fig. 16a). While the lower surface of the trailing (or following) bubble is altered significantly due to the stronger vortex of the fluid (see the velocity flow field in Fig. 16a) for n = 0.1, which is due to the high pressure region (see the outside pressure profile in Fig. 16a) as compared to the velocity flow field for a higher flow index of n = 0.9 (see Fig. 16b). However, a noticeable difference in the values between the outside and inside bubble pressure can be found for n = 0.1, e.g. the upper and trailing bubble pressure differences are 29.8 Pa and 97 Pa, respectively. However, for n = 0.9, the upper bubble pressure difference is 1.5 Pa, which means less deformation of the bubble, while a pressure difference of 40.24 Pa represents relatively more deformation of the trailing bubble. This can be attributed to the reduction of viscosity in local regions produced by the shear-thinning effect of the wake behind the upper bubble as well as the higher pressure difference on the bubble, which results in more deformation of the bubble.

Fig. 16 Velocity flow field, outside and inside pressure profiles on the rising of 6 mm bubble before coalescence at t = 0.30 s and θ = 45°, when (a) low flow index, n = 0.1; (b) high flow index, n = 0.9.

4. Conclusion

The dynamical characteristics of a pair of bubbles rising through non-Newtonian fluids were simulated numerically by the VOF method. The conclusions of this study are drawn below:

• In a horizontal configuration, the interaction between bubble pairs in non-Newtonian fluids shows a repulsive effect which increases with decreasing initial center-to-center distance of the bubbles due to the reduction of the shear-thinning properties of the fluid. This leads to an intensification of vortexes between the bubbles as well as a strong push-away effect between them.

• For an oblique angle configuration of a pair of bubbles, it was found that there is a repulsive effect between the bubbles when θ < 45°, while an initial angle configuration of θ > 45° shows an attraction between the bubbles, resulting in collision in a wide range of flow indexes, n ∼ 0.1–0.9.

• The flow field around the rising bubbles experiences significant changes with the variation of the flow index which can bring about a noticeable difference in the rising velocities. In addition, as the flow index decreases, a dramatic oscillation occurs between the bubbles and their shape alters from spherical to irregular.

Nomenclature

D⃑	Strain rate tensor, N m^−1
Fs	Body force, N m^−3
F	Volume fraction function
g	Gravitational acceleration, m s^−2
K	Consistency coefficient, Pa s^n
n	Flow index
n̂	Normal vector
n⃑	Unit normal vector
P	Pressure, Pa
do	Initial bubble diameter, mm
de	Equivalent bubble diameter, mm
dh	Bubble height as considered for the shorter bubble diameter, mm
dw	Bubble width as considered for the larger bubble diameter, mm
H	Column height, mm
h	Bubble position in the column, mm
Re	Reynolds number
ΔX	Actual distance between bubbles, mm
Xi	Initial distance between bubbles, mm
h*	Non-dimensional horizontal interval, [—]
t	Time, s
V⃑	Velocity vector, m s^−1
UT	Terminal velocity, m s^−1

Greek symbols

θ	Angle between the centre line of the bubbles in the horizontal direction, (°)
[small gamma, Greek, dot above]	Shear rate, s^−1
τ	Shear stress, Pa
μ	Viscosity, Pa s
μ(F)	Kinematic viscosity coefficient, m^2 s^−1
ρ	Density, kg m^−3
σ	Surface tension, N m^−1
k	Interfacial curvature between gas and liquid

Subscript

g	Gas phase
l	Liquid phase

Acknowledgements

This research is financially supported by University of Malaya, Ministry of Higher Education High Impact Research (UM.C/HIR/MOHE/ENG/20), and University of Malaya Research Grant (UMRG: RG121/11AET).

References

1. N. Kantarci, F. Borak and K. O. Ulgen, Process Biochem., 2005, 40, 2263–2283.
2. N. Yang, Z. Wu, J. Chen, Y. Wang and J. Li, Chem. Eng. Sci., 2011, 66, 3212–3222.
3. M. Pourtousi, J. Sahu and P. Ganesan, Chem. Eng. Process., 2014, 75, 38–47.
4. Z. Cai, Y. Bao and Z. Gao, Chin. J. Chem. Eng., 2010, 18, 923–930.
5. V. Badam, V. Buwa and F. Durst, Can. J. Chem. Eng., 2007, 85, 257–267.
6. S. D. Bari and A. J. Robinson, Exp. Therm. Fluid Sci., 2013, 44, 124–137.
7. P. Di Marco, N. Forgione and W. Grassi, Congresso Nazionale UIT sulla Trasmissione del Calore, 2005.
8. X. Zhang, H. Dong, D. Bao, Y. Huang, X. Zhang and S. Zhang, Ind. Eng. Chem. Res., 2013, 53, 428–439.
9. H. Dong, X. Wang, L. Liu, X. Zhang and S. Zhang, Chem. Eng. Sci., 2010, 65, 3240–3248.
10. L. Parkinson, R. Sedev, D. Fornasiero and J. Ralston, J. Colloid Interface Sci., 2008, 322, 168–172.
11. W. Fan, S. Jiang, C. Zhu, Y. Ma and H. Li, Opt. Laser Technol., 2008, 40, 389–393.
12. J. Herrera-Velarde, R. Zenit, D. Chehata and B. Mena, J. Non-Newtonian Fluid Mech., 2003, 111, 199–209.
13. T. Miyahara and S. Yamanaka, J. Chem. Eng. Jpn., 1993, 26, 297–302.
14. H. A. Jakobsen, H. Lindborg and C. A. Dorao, Ind. Eng. Chem. Res., 2005, 44, 5107–5151.
15. X. Yuan and G. Yu, Chin. J. Chem. Eng., 2008, 16, 497–502.
16. O. Hassager, Nature, 1979, 279, 402–403.
17. T.-J. Lin and G.-M. Lin, Chem. Eng. J., 2009, 155, 750–756.
18. R. C. Coke, R. Nagel and R. Etcheverry, Can. J. Chem. Eng., 2008, 81, 476–482.
19. H. Z. Li, X. Frank, D. Funfschilling and Y. Mouline, Chem. Eng. Sci., 2001, 56, 6419–6425.
20. N. Hassan, M. M. K. Khan and M. Rasul, WSEAS Trans. Fluid Mech., 2008, 3, 261.
21. T. Sanada, A. Sato, M. Shirota and M. Watanabe, Chem. Eng. Sci., 2009, 64, 2659–2671.
22. D. Legendre, J. Magnaudet and G. Mougin, J. Fluid Mech., 2003, 497, 133–166.
23. J. R. Vélez-Cordero, D. Sámano, P. Yue, J. J. Feng and R. Zenit, J. Non-Newtonian Fluid Mech., 2011, 166, 118–132.
24. W. Fan, Y. Ma, X. Li and H. Li, Chin. J. Chem. Eng., 2009, 17, 904–913.
25. W. Fan and X. Yin, Chin. J. Chem. Eng., 2013, 21, 705–713.
26. Z. Yu, H. Yang and L.-S. Fan, Chem. Eng. Sci., 2011, 66, 3441–3451.
27. M. Cheng, J. Hua and J. Lou, Comput. Fluids, 2010, 39, 260–270.
28. L. Zhang, C. Yang and Z.-S. Mao, J. Non-Newtonian Fluid Mech., 2010, 165, 555–567.
29. M. Ohta, E. Iwasaki, E. Obata and Y. Yoshida, J. Non-Newtonian Fluid Mech., 2003, 116, 95–111.
30. Z. Yu and L. S. Fan, Can. J. Chem. Eng., 2008, 86, 267–275.
31. J. Liu, C. Zhu, T. Fu and Y. Ma, Ind. Eng. Chem. Res., 2014, 53, 4850–4860.
32. J. Liu, C. Zhu, T. Fu, Y. Ma and H. Li, Chem. Eng. Sci., 2013, 93, 55–66.
33. M. van Sint Annaland, W. Dijkhuizen, N. Deen and J. Kuipers, AIChE J., 2006, 52, 99–110.
34. J. P. Singh and M. M. Denn, Phys. Fluids, 2008, 20, 040901.
35. J. Tsamopoulos, Y. Dimakopoulos, N. Chatzidai, G. Karapetsas and M. Pavlidis, J. Fluid Mech., 2008, 601, 123–164.
36. G. K. Batchelor, An introduction to fluid dynamics, Cambridge University Press, 2000.
37. C. W. Hirt and B. D. Nichols, J. Comput. Phys., 1981, 39, 201–225.
38. J. Brackbill, D. B. Kothe and C. Zemach, J. Comput. Phys., 1992, 100, 335–354.
39. E. C. Bingham, Fluidity and Plasticity, McGraw-Hill, New York, 1922, p. 440.
40. D. Ma, M. Liu, Y. Zu and C. Tang, Chem. Eng. Sci., 2012, 72, 61–77.
41. K. Dewsbury, D. Karamanev and A. Margaritis, Chem. Eng. Sci., 1999, 54, 4825–4830.
42. A. Margaritis, D. W. te Bokkel and D. G. Karamanev, Biotechnol. Bioeng., 1999, 64, 257–266.
43. D. Youngs, Numerical methods for fluid dynamics, 1982, vol. 24, pp. 273–285.

---