Analysis of Interaction and Flow Pattern of Multiple Bubbles in Shear-Thinning Viscoelastic Fluids

Abstract

A numerical study was conducted on the interaction of bubbles with different diameters and arrangements in shear-thinning viscoelastic fluids using OpenFOAM. The Volume of Fluid (VOF) method combined with the surface tension model was used to track the gas–liquid interface, and the rheological properties of the fluid were characterized with the Giesekus model. The numerical results are corresponded with the previous references, verifying the correctness of the simulation method. The influences of the initial bubble diameter, horizontal spacing, and arrangement on the motion state of three parallel bubbles were studied in detail. The flow pattern of the bubble rising was analyzed and compared with the motion state of parallel unequal double bubbles. As the diameter of the bubbles increases, the interaction among three equal size bubbles is changed from coalescence to detachment. Changing the diameter of one of the bubbles will significantly affect the movement of the larger diameter bubble, which is due to the enhancement in kinetic energy. The final state of some arrangement ways is consistent with the phenomenon of unequal double bubbles. The shear thinning effect, the velocity difference between bubbles, and the flow field around bubbles are considered the main reasons that decide the interaction between bubbles.

1. Introduction

The movement of bubbles in non-Newtonian fluid and the interaction between bubbles widely exist in the chemical industry, biology, environmental engineering, and other industrial processes [1,2]. In recent years, the motion characteristics of a single bubble in non-Newtonian fluid have been studied by many researchers [3,4,5,6,7]. However, for multiple bubbles, the motion laws of a single bubble cannot be directly applied. Hence the motion behavior of multiple bubbles is more instructive to the optimization of the actual process. Compared with other types of non-Newtonian fluids, the motion behavior of multiple bubbles in viscoelastic fluids seems very complicated [8,9,10,11,12,13,14], and some existing problems need to be solved urgently.

For years many researchers investigated the interactions between the bubbles. Zana and Lau et al. [15,16] found that the coalescence behavior of bubbles had a significant impact on the heat transfer, mass transfer efficiency, and chemical reaction rate between gas–liquid phases. Grienberger and Hofman [17] mentioned that the phenomenon of coalescence largely depended on the collision frequency between the bubbles. However, the collision frequency between bubbles/droplets depended on the turbulent random fluctuations, lift, shear stress, vortices generated by the rising velocity, etc. The total collision frequency was regarded as the summation of the collision frequencies caused by all factors. After the collision of bubbles, coalescence may occur. The process can be concluded that the trailing bubbles start to accelerate and stretch longitudinally, and then come into contact with the preceding bubbles to form a liquid film, which drains, thins, and ultimately ruptures [18,19,20]. For viscoelastic fluids, Lin and Lin [21,22] believed that the upward flow of leading bubbles caused by shear thinning and viscoelastic effects generated drag and thrust, leading to the coalescence of the two bubbles. Komasawa [23] investigated that the wake generated by rising bubbles in viscoelastic fluid and Newtonian fluid was very different. There were three kinds of wake generated by rising bubbles in Newtonian fluid, namely laminar wake, transitional wake, and turbulent wake. For viscoelastic fluids, a wake was generated by the rise of bubbles in the opposite direction of motion, known as a negative wake [24], which hindered bubble coalescence [25]. Imaizumi et al. [26] carried out that the generation of the negative wake was related to the viscoelasticity of the fluid, and the effect of elastic force was opposite to inertial force [27].

In recent years, the study on bubble motion in viscoelastic fluids has increasingly focused on precision, efficiency, and complexity, with more researchers preferring numerical simulation. Warnez and Johnsen [28], Gaudron et al. [29], and Johnsen and Mancia [30] in a Keller-Miksis framework, a Kelvin–Voigt viscoelastic model with full thermal effects being considered, in which the elastic term is represented by a Neo-Hookean model to account for the finite strains. Numerical results show that the relaxation time increases bubble growth and permits rebounds driven purely by residual stresses in the surroundings. Different regimes of oscillations occur depending on the relaxation time. Many different simulation methods were developed for gas–liquid interface trackings, such as Level Set Method (LSFEM), Frontal Tracking Method (FTM), Lattice Boltzmann Method (LBM), and Fluid Volume Method (VOF). Annland et al. [31] summarized the advantages and disadvantages of various methods. For viscoelastic fluids, Pillapakkam et al. [6,32] and Ohta et al. [33] used the Oldroyd-B model and FENE-CR model to characterize their rheological properties and studied bubble motion in viscoelastic fluids. Wagner et al. [34] used LBM (Lattice Boltzmann Method) to study the two-phase flow of bubbles rising in viscoelastic fluids, and obtained the sharp points at the tail of the bubbles. Pillapakkam and Singh [35] used the two-dimensional level-set method (LSFEM) and introduced dimensionless capillary number (Ca) and Deborah number (De) to characterize bubble motion in Oldroyd-B fluid. It was found that for bubbles above the critical dimensionless number, the local viscoelastic and viscous stresses caused the upper end of the bubble to be smooth and the lower end to be sharp. Hua and Lou [36] improved the traditional Frontal Tracking Method (FTM) algorithm by solving the Navier–Stokes equation to analyze the motion of bubbles. Among many interface tracking methods, the Volume of Fluid (VOF) method is relatively simple, and can accurately track the deformation of bubbles, thus making it widely used. By the Volume of Fluid (VOF) method, Liu et al. [37] studied the motion and interaction of three equally spaced parallel bubbles in a shear-thinning fluid. The dimensionless critical horizontal spacing of bubble coalescence under different physical properties was carefully developed. Zhu et al. [38] discussed the fluid concentration distribution during the rise of single bubbles and double bubbles in shear-thinning fluid. They mentioned that the concentration distribution had a great influence on the bubble deformation rate. In recent years, studies about the interaction between multiple bubbles of non-Newtonian fluid were reported widely, and most of them were shear thinning fluids. At present, although some scholars have conducted a lot of research on the rising process of bubbles in viscoelastic fluids, more in-depth research is still needed on the interaction mechanism of the rising process of multiple bubbles.

In the present study, the Volume of Fluid (VOF) method is used to track the bubble interface, and the Giesekus model is carried out to characterize the rheological properties of viscoelastic fluids. Based on the OpenFOAM CFD system, the aggregation and interaction of two and three bubbles in different states rising in the static shear-thinning viscoelastic fluid are studied. The effects of the bubble diameter, bubble spacing, and different arrangement methods on bubble rise and interaction are investigated. The interaction mechanism between bubble groups is discussed and analyzed for further exploration.

2. Mathematical Models and Numerical Methods

Figure 1 shows an initial calculation model, the spherical bubbles (A₁, A₂) with diameters d₁, d₂, and d₃ are placed in the center of the bottom of the calculation domain B, and the spacing is set as H. To eliminate the influence of the wall on the motion of bubbles [33], the calculation domain is set to a rectangle with P × Q (15 × 15 cm). The upper boundary is a zero gradient pressure outlet, the bottom and sides are regarded as walls, and the boundary condition is no slip [39,40]. At the initial moment, bubbles begin to rise freely driven by the buoyancy.

2.1. Governing Equation

Viscoelastic fluids (liquid phase) and bubbles (gas phase) are continuous and immiscible, and the system is isothermal and incompressible. The density of viscoelastic fluids ρ and relaxation time λ are considered constant. The conservation equation of the system consists of a continuity equation and a momentum equation, shown as follows:

$$\nabla \mathrm{u}=0$$

(1)

$$\rho \left(F\right)\left[\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}\times \nabla \right)\mathbf{u}\right]=-\nabla p+\nabla \times \left[2\mu \left(F\right)\mathbf{D}\right]+F_S+\rho \left(F\right)g$$

(2)

where u is the velocity vector and p is the pressure, μ(F) is the dynamic viscosity coefficient, F_S is the source term caused by surface tension, represents density, and D is the stress tensor. Local average density ρ(F) and dynamic viscosity coefficient μ(F) use the volume fraction function F for local distribution estimation.

$$\rho \left(F\right)={\rho }_l\left(F\right)+{\rho }_g\left(1-F\right)$$

(3)

$$\mu \left(F\right)={\mu }_l\left(F\right)+{\mu }_{\mathrm{g}}\left(1-F\right)$$

(4)

The subscripts g and l represent the gas and liquid phases, respectively. The stress tensor is represented as follows:

$${\mathbf{D}}_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_j}\right)$$

(5)

The volume fraction function F is defined as follows:

$$F=\left\{\begin{array}{c}1\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}\hspace{0.33em}\mathrm{t}\mathrm{h}e\mathrm{\hspace{0.33em}}\mathrm{b}\mathrm{u}\mathrm{b}\mathrm{b}\mathrm{l}\mathrm{e}\\ 0<F<1\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathrm{a}\mathrm{t}\hspace{0.33em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.33em}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\hspace{0.33em}\\ 0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathrm{o}\mathrm{u}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}\hspace{0.33em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.33em}\mathrm{b}\mathrm{u}\mathrm{b}\mathrm{b}\mathrm{l}\mathrm{e}\end{array}\right.\left.\begin{array}{l}\\ \\ \end{array}\right\}$$

(6)

$$\frac{\partial F}{\partial t}+\nabla \times \left(\mathbf{u}F\right)=0$$

(7)

2.2. Source Term of Momentum Equation Caused by Surface Tension

Surface tension is considered that it has a significant impact on the interface because the smaller curvature of a bubble generates greater additional pressure. The continuous surface force model proposed by Brackbill et al. [41] used the source term Fs to describe the influence of the surface tension:

$$F_S={\displaystyle \sum _{i<j}{\sigma }_{ij}\frac{F_i{\rho }_i{K}_j\nabla F_j+F_j{\rho }_j{K}_i\nabla F_i}{\frac{1}{2}({\rho }_i+{\rho }_j)}}$$

(8)

While there are only liquid and gas phases, K_i = K_j; ∇F_i = −∇F_j, which can be represented as

$$F_S={\sigma }_{ij}\frac{\rho K_i\nabla F_i}{\left(1/2\left({\rho }_i+{\rho }_j\right)\right)}$$

(9)

K is the interface curvature, represented by the divergence of the unit normal vector of the interface:

$$K=-\left(\nabla \times \stackrel{⌢}{n}\right)=\frac{1}{\left|\mathbf{n}\right|}\left[\left(\frac{\mathbf{n}}{\left|\mathbf{n}\right|}\times \nabla \right)\left|\mathbf{n}\right|-\left(\nabla \times \mathbf{n}\right)\right]$$

(10)

While $\stackrel{⌢}{n}=\mathbf{n}/\left|\mathbf{n}\right|$, $n=\nabla \tilde{F}$, $\tilde{F}$ is expressed as a stationary continuous volume fraction function.

2.3. Constitutive Equation of Continuous Phase

The Giesekus model [42] can simultaneously represent the shear thinning characteristics and elasticity of fluids, which is widely used in the simulation of gas–liquid two-phase flow [43,44]. However, the selection and determination of some of the parameters of the Giesekus model may rely on experimental data or experience, leading to limitations in the application of the model. The differential equation is defined as:

$${\tau }_{P_K}+{\lambda }_K{\stackrel{\nabla }{\tau }}_{P_K}+{\alpha }_{}\frac{{\lambda }_{}}{{\eta }_{P_K}}\left({\tau }_{P_K}\cdot {\tau }_{P_K}\right)=2{\eta }_{P_K}D$$

(11)

${\tau }_P$ is the symmetry tensor obtained by the summation of polymer contributions under each relaxation mode, lying on the viscoelastic constitutive equation:

$${\tau }_P={\displaystyle \sum _{K=1}^n{\tau }_{P_K}}$$

(12)

${\stackrel{\nabla }{\tau }}_{P_K}$, ${\lambda }_K$, and ${\eta }_{P_K}$ are applied to the upper convection derivative, relaxation time, and viscosity at zero shear rates of stress tensors, respectively. Nonlinear parameters in the Giesekus model α is a dimensionless “mobility factor” of the relaxation mode, which is related to the various anisotropic forces on macromolecular polymers that prevent flow.

The expression for the derivative of upward convection is expressed as

$${\stackrel{\nabla }{\tau }}_{P_K}=\frac{D}{Dt}{\tau }_{P_K}-\left[\nabla {\mathbf{u}}^T\cdot {\tau }_{P_K}\right]-\left[{\tau }_{P_K}\cdot \nabla \mathbf{u}\right]$$

(13)

For symmetric tensors:

$$\stackrel{\nabla }{{\tau }_{P_K}}=\frac{D}{Dt}{\tau }_{P_K}-\left[{\tau }_{P_K}\cdot \nabla \mathbf{u}\right]-{\left[{\tau }_{P_K}\cdot \nabla \mathbf{u}\right]}^T$$

(14)

The key derivative of the additional stress tensor is represented as follows:

$$\frac{D}{Dt}{\tau }_{P_K}=\frac{\partial }{\partial t}{\tau }_{P_K}+\mathbf{u}\cdot \nabla {\tau }_{P_K}$$

(15)

2.4. Numerical Methods

OpenFOAM is used for numerical calculation and Finite Volume Method is developed for solution. The fluid volume method (VOF) is carried out to capture the interface between gas–liquid two-phase flow [45]. Mesh refinement technology can ensure the numerical simulation of viscoelastic two-phase flow is more accurate. The SIMPLEC algorithm is used to couple pressure velocity [46], ensuring the correct and stable flow of viscoelastic fluids [39]. The Giesekus model is performed to calculate highly elastic two-phase flow [47]. Using the logarithmic conformational tensor method, the polymer stress tensor ${\tau }_p$ is related to the conformational tensor A [48]:

$${\tau }_p=\frac{{\eta }_p}{\lambda }(\mathbf{A}-\mathbf{I})$$

(16)

Here η_p is the viscosity of the viscoelastic fluid, λ is the coefficient of fluid relaxation time, and I is the unit tensor [25]. A new tensor $\mathrm{\Theta }$ is created to represent the logarithm of conformational tensors:

$$\mathrm{\Theta }=\ln (\mathbf{A})=R\ln (\mathrm{\Lambda })R^T$$

(17)

Since A is positive definite and symmetric and it has eigenvalues, which is easy to diagonalize $\mathbf{A}=\mathbf{R}\mathsf{\Lambda }{\mathbf{R}}^T$. R is a rotation matrix formed by the eigenvector of A [49].

3. Results and Discussion

3.1. Validation

Table 1. Physical and rheological properties of gas–liquid two-phase flow.

	ω	ρ	η	σ	λ	α
liquid phase	0.8	1000.9 kg/m3	0.001 Ns/m2	0.076 N/m	0.207 s	0.6
gas phase		1 kg/m3	1.7 × 10−5 Ns/m2

shows the physical and rheological properties of the continuous phase, the experimental data is cited from Pilz and Brenn [50]. A 0.8 wt% Polyacrylamide(PAA) solution is used. Due to the shear thinning effect of the solution, the Giesekus model is studied caused of the rheological properties of the continuous phase. Migration factor α is achieved by fitting the formula between the relaxation time λ used by Pilz and Brenn [50] and the zero shear viscosity value. The variation of shear viscosity with the shear rate for the Giesekus model can be obtained according to Equations (18) and (19), as shown in Figure 2. The error is similar to that in the work of Yuan et al. [25] and the effect on the results can be neglected.

$$\eta ={\eta }_s+{\eta }_p\frac{4(1-\alpha )}{\sqrt{f+1}[\sqrt{f+1}+\sqrt{2}(1-2\alpha )]}$$

(18)

$$f=\sqrt{1+16\alpha (1-\alpha ){\lambda }^2{\dot{\gamma }}^2}$$

(19)

Three grid sizes of 0.1 mm, 0.2 mm, and 0.3 mm are selected for the independence test. Taking the rise of bubbles with a diameter of 6 mm as an example, Figure 3 shows the bubble velocity calculated by three kinds of grid computing. It can be seen from Figure 3 that the velocity is fluctuating. When the bubble movement tends to be stable, the difference between grids of 0.1 mm and 0.2 mm can be ignored. However, the error of grid computing with a size of 0.3 mm is relatively large. Taking into account the calculation duration and the accuracy of simulation results, finally, a grid size of 0.2 mm is selected.

As we know, the most reliable way to verify the accuracy of the simulation method is to compare the experimental data with the results obtained from the simulation. However, experimental studies on the motion of parallel double and multiple bubbles in viscoelastic fluids have not been reported in the open literature. Therefore, this study presents a comparison between the experimental results of a stable rising single bubble in viscoelastic fluids with the simulation results. Figure 4 shows the functional relationship between the terminal velocity of a single bubble and the volume of a bubble in a 0.8 wt% PAA non-Newtonian liquid. It is observed that the numerical simulation results are in good agreement with the experimental results of Pilz and Brenn [50], both of which exhibit velocity jumps [50,51].

3.2. Parallel Unequal Diameter Double Bubble Upflow Pattern

Double bubbles freely rising experiment is developed in a viscoelastic solution of 0.8 wt% PAA. Large bubble diameter is set as d₁ = 6 mm, for small bubbles of d₂ = 4 mm, 5 mm, and edge spacing of H = 1 mm. Figure 5 shows the upward flow pattern of double bubbles with unequal diameters, and the bubbles coalesce. From the figure, it can be seen that as the diameter of the small bubble increases, under the combined action of the wake effect, shear thinning, and buoyancy, the velocity of the small bubble increases, thus gradually catching up with the large bubble. The small bubble is then squeezing the liquid below it, pushing the large bubble up and causing its velocity to enhance obviously. Finally, the liquid film ruptures and coalesces into a single bubble. The simulation results are almost consistent with the conclusions obtained by Bracbill et al. [41]. The process of online bubble coalescence is as follows: two bubbles first collide to produce a liquid film, and then continue to compress, causing the liquid film to drain and dry. Finally, the liquid film ruptures and the two bubbles coalesce. Some scholars also observed this phenomenon [18,19,20].

Figure 6 presents the interaction process between two bubbles with diameters of d₁ = 6 mm and d₂ = 4 mm respectively, as well as the changes in the surrounding flow field. It can be seen that at t = 0.01 s, due to the buoyancy and surface tension leading to the rise of bubbles, the bubble velocity is slower and the shape is close to spherical. The two bubbles start to interact with each other, and the flow field is intertwined, with vortices appearing on both sides. While at t = 0.06 s, the vortex on the left side of the large bubble disappears, and the flow fields of the two bubbles show the same performance as well as that of the individual large bubble. The small bubble is pushed towards the large bubble, and the wake effect of the large bubble causes the small bubble to stretch laterally. When t = 0.2 s, based on the flow field results in Figure 6, it can be seen clearly that the small bubble is significantly affected by the large bubble. The wake effect of the large bubble causes the small bubble to accelerate and elongate in the direction of the large bubble, thus gradually moving toward the center of the wake. When t = 0.34 s, the two bubbles and their flow field are symmetrical in the vertical direction, which displays the same phenomenon of the rise of linear unequal double bubbles. At this time, the small bubbles accelerate to catch up with the large bubbles in the vertical direction due to the shear thinning effect and wake effect of the fluid, thus pushing the liquid above to uniformly act on the lower surface of the large bubbles. The shape of the large bubbles is flat and gradually turns into a spherical cap shape [52]; for t = 0.36 s, the liquid between the two bubbles is squeezed out, and then forming a very thin liquid film. The upward force of the small bubble is greater than the horizontal force. It can be observed that the velocity at both ends of the large bubbles increases sharply, and the bubbles further stretch horizontally with greater concave curvature, leading to the rupture of the liquid film. At this moment, the flow field distribution is very uniform, and merging into the large bubble.

We can see the rising process of bubbles of d₂ = 5 mm is similar to that of the bubble of d₂ = 4 mm. The difference is that a small bubble gets coalescence before reaching the bottom of the large bubble. It can be explained that the bubble with a diameter of 5 mm has a larger volume, which is easily captured by the rising wake of the large bubble. In addition, the shear thinning effect accelerates the small bubble to catch up with the large bubble. As can be seen from the flow field diagram in Figure 7, the streamline at the bottom left region of the large bubble is denser than that of the bottom right region. Since the small bubble pushes the liquid upwards, the viscoelastic stress and thrust of the bubble act intensively on the bottom left side, leading to the early coalescence of the two bubbles. According to the experiments of Li and Gupta [52] and Stewart [53], it is concluded that the coalescence between bubbles occurs while one bubble is captured by the wake of another bubble.

3.3. Three Bubbles Upflow Pattern with Parallel Equal Diameter

Figure 8 shows the motion trajectory of bubbles with diameters of 4 mm, 5 mm, and 6 mm and spacing of 1 mm. In Figure 8a, three bubbles come into contact and coalesce with each other simultaneously. At the initial moment, when the bubble of 4 mm rises, under the action of shear thinning, the viscosity around the middle bubble is decreased, as shown in Figure 9. The motion speed increases and the velocity gradient difference formed between the bubbles cause the bubbles on both sides to be attracted until direct contact. When t = 0.275 s, the bubbles at both ends move toward the middle bubble and contact with each other. From Figure 10a, it can be seen that the vortices of the two below bubbles act on the upper sides of the two bubbles, respectively, converging on both sides below the middle bubble. The middle bubble interacts with the two bubbles to form two steam bridges [34]. The curvature of the gas bridge is small, and under the action of surface tension, the merged bubbles can still stretch longitudinally, which is attributed to differences in bubble size and rheological property.

As shown in Figure 8b, the left and right bubbles first coalesce below the middle bubble, and then come into contact with the lower end of the middle bubble to coalesce. This phenomenon can be explained that as the volume of the bubble increases, the initial rising speed of the middle bubble increases, and the upward force of the bubbles at both ends of the middle bubble enlarges. Later, due to the shear thinning effect, the lower bubble appears inclined motion, and the two bubbles come into contact, causing the liquid film to become thinner and eventually rupture, finally leading to coalescence. Interestingly, with the coalescence of three bubbles, a satellite bubble appears at the lower end. That means that when the three bubbles coalesce, a portion of the liquid remains in the bubble after completing a series of processes such as the liquid film contact, thinning, and rupture.

As seen in Figure 10b, the flow field performance around the formation of the satellite bubble is studied. Due to the high mass and density of the liquid, which is affected by gravity with a low rising speed, the merged bubbles generate eddies in the gas valley. The eddies push the liquid downward and compress the bubbles at both ends, ultimately causing the bubbles to rupture and form satellite bubbles. Figure 8c shows that the three bubbles do not merge after rising, and the rising speed of the middle bubble is slowly increased, ultimately exceeding the speed of the bubbles on both sides. While the diameter of the bubble increases to 8 mm, due to the small initial spacing between the bubbles, the surrounding flow fields overlap with each other during the rising process, and the interaction between the bubbles is significant. After the bubble rises for a period of time, it can be seen from Figure 10a that there is a tendency for vortices to appear between the middle bubble and the right bubble. The vortex of the right bubble also appears clearly, indicating that the horizontal force of the bubble gradually decreases. From the velocity cloud map in Figure 11, it can also be observed that as the bubble rises, the vertical velocity is not changed obviously, as well as the vertical distance. This indicates that the interaction force between bubbles has a significant impact on the horizontal distance; however, it has little effect on the vertical distance. This phenomenon is consistent with the results of Liu [37].

In order to better explore the deformation and interaction of bubbles, Figure 12 shows the first normal stress difference contour ($N_1=\left|{\tau }_{\mathrm{xx}}-{\tau }_{yy}\right|$) at different times under three configurations. The non-zero region in the figure can clearly display the viscoelastic stress of the liquid. It can be observed that during the rising process of bubbles, viscoelastic stress is mainly distributed on the surface of the bubbles and resonates, and the interaction between bubbles affects the stress field. As the diameter of the bubble increases, while at time 0, the aggregation stress does not occur, and the bubble is accelerated under the action of buoyancy. With forming the aggregation stress, the stress between the bubbles increases, and the deformation of the bubbles on both sides is greater. There is a certain degree of stagnation stress, which hinders the rise of the middle bubble. This phenomenon is similar to the result obtained by Vahabi et al. [54]. Considering the viscoelastic stress between bubbles, it can be observed that the enhancement of the viscoelastic stress in the gap between bubbles would be more likely to form the coalescence. In Figure 12c, while at t = 0.1 s, the viscoelastic stress appears at the lower end of the bubbles on both sides, and the lower end of the middle bubble is relatively high, showing a trend of coalescence. However, under the acceleration of buoyancy and shear thinning, the bubbles ultimately fail to coalesce, indicating that viscoelastic stress can remarkably affect the interaction between bubbles.

The upward motion trajectory of three bubbles with a spacing of H = 3 mm is shown in Figure 13. It can be seen that when the diameters of the three bubbles are less than 4 mm, the bubbles at both ends coalesce firstly. This can be explained that as the spacing increases, the middle bubble is pushed by the liquid above the bubbles at both ends. Under the combined effect of buoyancy and the driving force provided by the bubbles at both ends, the velocity increases and rises first. The shear influence of the wake of the middle bubble causes the viscosity of the liquid below it to decrease, and we observe the bubbles at both ends approach each other and coalesce clearly. After the bubbles on both sides coalesce, they are subjected to the action of the wake of the upper bubbles and stretched longitudinally. Due to the increase in volume and velocity, they eventually collide with the previous intermediate bubbles and then recoalescence. The rising phenomenon of three bubbles with other sizes is almost consistent with that of the spacing of 1 mm.

3.4. Upflow Pattern with Parallel Unequal Diameters and Complex Arrangement of Three Bubbles

In industrial processes, the interactions between bubbles of different diameters are considered a common phenomenon. For the application, we change the diameter of the bubble to observe the interaction between bubbles. Figure 14 and Figure 8a show the rising process of three bubbles and the variation of the bubble shape with the different diameters of the right bubble. During the rising process, there is a force between three bubbles, causing them to deform and approach each other, which results in the middle bubble coming into contact with the right bubble and quickly merging. Unlike the coalescence of medium-diameter bubble in Figure 8a, the middle bubble does not rise first because the bubble diameter increases and so the volume increases. Under the action of buoyancy, the initial rising speed increases, and then it is faster than the other two bubbles. Under the shear thinning effect and the circulation flow of large bubbles, the rising speed of the middle bubble is higher than that of the left bubble, which is consistent with the principle of the middle bubble rising firstly in equal diameter of three bubbles. Due to the velocity difference between bubbles, the middle bubble firstly contacts with the large bubble to form a gas bridge, and the two bubbles start to form a large bubble. Similarly, due to the velocity difference caused by the shear thinning effect, a small bubble catches up with a large bubble. At t = 0.48, the two bubbles come into contact to form an air bridge, which is consistent with the results in Figure 7a. Finally, under the action of surface tension, viscoelastic stress, and inertial force, a whole large bubble is formed. When the bubbles on the right increase to 6 mm, the initial rising stage is basically consistent with that of Figure 14a, but in the middle and later stages, the small bubbles and merged bubbles do not coalesce. This is because after the merger process, the volume rapidly increases and the rising speed increases obviously. Although the bubbles below are accelerated under the action of shear thinning, a wake opposite to the rising direction appears below the merged bubbles, thus hindering the small bubbles from catching up with the large bubbles. This factor causes the three bubbles to not merge into one large bubble.

Figure 15 shows the motion state of the three bubbles rising freely when the bubble size remains unchanged and the spacing between the bubbles is enlarged. The phenomenon in Figure 14a is consistent with that in Figure 13b, and the principle of interaction between bubbles is basically identical. When the diameter of the right bubble increases to 6 mm, the motion state of the bubble will be changed clearly. Figure 16 shows the cloud diagram of the velocity of bubbles and surrounding liquid during the rising process. When t = 0.1 s, the velocity of the right bubble is much higher than that of other bubbles, and the surrounding fluid flow drives the middle bubble to move toward the direction of the large bubble. As the bubbles continue to rise, the movement of the large bubbles above gradually stabilizes, and the shape of the bubbles is an inverted teardrop. When t = 0.7 s, the interaction between bubbles gradually decreases, and the rise of bubbles has a significant impact on the flow of surrounding fluid, and bubbles cannot coalesce. Here, we predict that the two bubbles below will merge at some point. This is because after the stable rise of a large bubble, a negative wake will be generated at the lower end of the bubble, which will hinder the rise of the middle bubble. When the two bubbles are close together, under the action of shear thinning, the lower bubble accelerates to catch up with the middle bubble. Under the influence of inertia, viscoelastic stress, and other factors, the two bubbles will converge at a certain time.

As discussed, the size of the intermediate bubble strongly affects the upward trajectory of the bubble, as shown in Figure 17a,b. Intermediate bubbles with d = 4 mm and d = 5 mm rise in viscoelastic fluids, and the bubble shape and upward trajectory are basically consistent. Taking the 5 mm middle bubble as an example shown in Figure 18, while at 0.15 s, the viscosity of the liquid above the middle bubble decreases. However, under the circulation of bubbles on both sides and the squeezing of the liquid between the bubbles, the bubble fails to accelerate its movement, causing the bottom of the bubble to flatten, and the bubbles on both sides are subjected to reaction forces and tilt upwards. As the bubbles continue to rise, the liquid between the two sides of the bubble and the middle bubble decreases in apparent viscosity due to shear thinning. The rise of the middle bubble is accelerated, but due to the volume differences, wake effects, and repulsion between bubbles, the middle bubble cannot catch up with the two sides of the bubble. From the perspective of shape, the interaction force between bubbles decreases, and the shape of bubbles on both sides gradually tends to be a single bubble. There is a shear thinning effect with the middle bubble produced by the bubbles on both sides, and the rising speed becomes faster, leading to the appearance of sharp points.

As the spacing between bubbles increases to 3 mm, the interaction between bubbles is weakened, and the horizontal and vertical spacing between bubbles are both enlarged. The overall phenomenon is almost consistent with that of a spacing of 1 mm. The motion trajectory is shown in Figure 13c and Figure 19.

4. Conclusions

In this study, the Volume of Fluid (VOF) method and the Giesekus model were carried out to study the motion behavior of multiple bubbles in shear-thinning viscoelastic fluids under different conditions. The SIMPLEC method was applied to couple pressure and velocity smoothly and accurately. Mesh refinement technology was used at special behavior locations. By varying the velocity, shape, and viscoelastic stress obtained from the rise of unequal double bubbles, in-depth research was conducted on the rising of three parallel bubbles under different conditions. The mechanism of the separation, coalescence, and deformation between bubbles were discussed and analyzed through a complete flow state, bubble state, flow field structure, bubble deformation, and stress diagram.

In the rise of parallel and unequal bubbles, with the increasing diameter, the two bubbles attract each other caused by the pressure difference and bubble circulation flow. The small bubble will be pushed by the upward wake of the large bubble. The rising of the large bubble leads to shear thinning and viscoelastic stress in the surrounding fluid, ultimately resulting in a bubble coalescence.

During the motion of equal diameter of three bubbles, as the diameter of the bubbles increases, the viscoelastic stress, velocity difference, pressure difference, and circulating flow force of the bubbles are weakened, thus transitioning from bubble aggregation to bubble detachment. Afterward, due to the repulsion between bubbles, the horizontal spacing gradually increases and ultimately remains stable. Varying in bubble spacing is almost consistent with the above situation.

By varying the diameter of bubbles at a certain position, it can be found that the bubbles coalesce or not largely depends on the initial spacing between bubbles and the initial size of bubbles. When bubbles rise, the surrounding flow field affects the shape, rising speed, and wake structure of the bubbles, which have a significant impact on bubble coalescence. While a bubble rises, the side where the large bubble is located has a greater impact on the interaction between bubbles. The size of the bubble determines the shape and flow field structure of the bubble and influences the interaction between the bubbles.