Simulation analysis of bubble coalescence behavior characteristics in non-Newtonian fluids based on the phase field method

Gas-liquid two-phase flow is widely used in various fields of chemical production because of its stability. The bubble coalescence behavior has a significant impact on the size, shape and movement of bubbles. These parameters are important parameters for affecting the interphase mass and heat transfer processes and equipment performance. This paper presents a simulation study of the bubble coalescence behavior of non-Newtonian fluids using the phase field method through simulation software. The relative location of bubbles were divided into three types: vertical, parallel and random; and also the distribution of bubbles can be divided into three cases: symmetric, left-skewed and right-skewed. This paper focuses on the influence of the arrangements of the bubbles on their coalescence behavior in non-Newtonian fluids. The results show that in non-Newtonian fluids, the relative distance of the bubbles and the flow characteristics of the fluid are both key factors in whether bubble coalescence can occur, and they exhibit different patterns of action when the bubbles are in different arrangements. Also bubbles show different patterns of motion in different distributions.


Introduction
Gas-liquid two-phase flow is widely used in chemical equipment such as bubble towers, fluidized beds and bioreactors because of the small specific surface area and good stability of bubbles, which can provide high rates of mass and heat transfer and chemical reaction efficiency [1].The deformation and velocity of the bubbles in the production equipment are related to their coalescence behavior, and these factors directly affect the mass and heat transfer performance and process index of the above equipment.Therefore, the investigation of bubble coalescing law is an important part of the study of bubble dynamics behavior and chemical equipment process enhancement [2].
Jiang et al [3] carried out experimental and numerical simulations on the deformation law of a single bubble of a non-Newtonian fluid, and the study revealed the influence of its aspect ratio on the deformation and viscosity during the rise of the bubble under this condition.Chen et al [4] used the idea of microlayer modeling to numerically simulate the single bubble growth and detachment process.Huang et al [5] discussed experimentally the deformation characteristics and laws of single bubbles in the rise of fluids with different viscosities by means of a self-designed bubble generation device.Ye et al [6] solved the kinetic law of the double bubble coalescence process numerically and used the level set method.Wen [7] carried out an experimental study on the factors influencing the coalescence of multiple bubbles in pure water systems by means of an in-house designed experimental setup.Maxim K [8] conducted an experimental study of double bubble collision and fragmentation phenomena.Zhao et al [9] used the univariate method to simulate and analyze the deformation influence laws of different size parameters on the process of double bubble coalescence.
From the existing literature, it is known that there are fewer studies related to the coalescence behavior of bubbles in non-Newtonian fluids under different arrangements.The coalescing behavior of bubbles can affect processes such as distillation, wastewater treatment, absorption, etc.This paper uses the phase field method combined with numerical computation software to analyze the coalescence behavior of double bubbles in non-Newtonian fluids under different arrangements by changing the model and boundary conditions, and to explore the motion law of bubbles in this process, which provides some reference for related fields.[10], the bubble diameter is the main factor affecting its coalescence, the double bubble diameter in this paper are set to 4mm.The Reynolds number of the bubble is less than 200 and its motion shows two-dimensional characteristics [11], in this paper, the Reynolds number is about 25, so the geometrical model in this paper can be simplified to a two-dimensional model.According to the literature [12], when the distance between the center of the bubble and the wall is greater than five times the radius of the bubble, the effect of the wall on the flow can be ignored and the width of the geometric model is set accordingly.The geometric model of this paper selects a rectangular area of 80×60mm, which is filled with stationary liquid.The distribution of bubbles is based on the vertical symmetry axis of the geometric model, and three distributions are considered: left-sided distribution, center-symmetric distribution and right-sided distribution, as shown in Figure 1.The horizontal spacing Δx and the vertical spacing Δy are used to describe the position relationship between the bubbles, as shown in Figure 2. By adjusting these two parameters, the position relationship between the double bubbles can be made to appear vertical, parallel and random.2.1.2Non-Newtonian fluid properties.The non-Newtonian fluid in this paper is a solution of air bubbles in transformer oil, a form of flow that has been widely discussed in studies on the insulation and safety properties of transformers [13], and whose viscosity has the following mathematical relationship with the shear strain rate: n is the flow characteristic index.When n and 1 exhibit different magnitude-value relationships, the viscosity and shear velocity of the fluid exhibit different laws of action.The values of the physical properties parameters of the non-Newtonian power law fluids in this paper are shown in Table 1.For non-Newtonian fluids, the smaller the flow index n, the stronger the non-Newtonianity; n=0.6 in this example belongs to pseudoplastic fluids or non-Newtonian fluids with weak shear thinning..Numerical simulation methods commonly used for gas-liquid twophase flow include the volume of fluid method, the level set method and the phase field method.The volume of fluid method requires strict quality of the grid, and the horizontal set method has limited accuracy and is prone to mass non-conservation [14].The phase field method is a mathematical model for solving boundary problems by introducing the phase field variable φ and considering the combined effect of ordering potential and thermodynamic driving force to establish the phase field equations, thus obtaining the instantaneous state of the study system in time and space under the assumption that the interface is diffuse.So the phase field method can better satisfy the conservation of mass and deal with the abrupt changes of physical quantities of smooth cross sections.The phase field interface control equation used in this paper is the Cahn-Hilliard's nonlinear diffusion equation containing the fourthorder diffusion term.
In Equation (7), φ denotes the normalized phase field variable, which varies between -1 and 1; ρ1 and ρ2 are the densities of the gas and liquid phases; C is the corresponding two-phase concentration.This paper's corresponding φ for air and transformer oil are -1 and 1.The gas-liquid interface is where the phase field variable φ=0 can determine the bubble's position at different moments [15].ψ denotes the chemical potential and is calculated from Equation (8).ε for the two-phase interface thickness control parameter usually taken ε=hc/2, where the hc for the mesh feature size.γ denotes the migration regulation parameter, which is accepted as γ=ε 2 ; λ denotes the mixing energy density, whose value satisfies Equation (9).2.2.4 Intrinsic Equation.It forms a closed set of equations to characterize the flow of the fluid; in this paper, the intrinsic equation of the power-law fluid is:

Basic assumptions
The following assumptions are made for the simulation model developed in this paper: (1) The liquid phase is incompressible fluid, the gas phase is compressible fluid; (2) The bubble motion is subject to gravity, viscous force and surface tension; (3) Phase change occurs only at the gas-liquid contact surface.

Initial and boundary conditions
The main physical parameters used in the simulation analysis of gas-liquid two-phase flow in this paper are shown in Table 2.The effect of air bubbles on the insulating properties of transformer oil is mainly based on the bubble breakdown theory, which mainly describes the phenomenon of breakdown caused by air bubbles in a static liquid dielectric, so the simulation analysis in this paper is calculated according to the following boundary conditions.Since the inflow of the liquid phase is not considered, the lower part of the geometric model and the left and right sides are no-slip boundary conditions.Flow properties are set to laminar flow.The inlet boundary conditions are gas-phase inlet velocity set to 0.06m/s, 1 standard atmospheric pressure at room temperature.The upper part of the geometric model is the outlet boundary condition, and the pressure value is set to 1 standard atmosphere.Surface tension(N•m -1 ) 0.048 In order to avoid the influence of mesh scale on the results and at the same time minimize the solution volume as much as possible, four meshing schemes, Mesh1, Mesh2, Mesh3, and Mesh4, are established to check the influence of mesh scale on the computational results, corresponding to the minimum mesh sizes of 0.45mm, 0.3mm, 0.2mm, and 0.15mm, respectively.In verifying the grid independence, the bubble diameter was 4 mm, the inlet velocity was 0.06m/s, and other physical parameters were set according to the data in Tables 1 and 2. The variation of single bubble rise velocity with time for the four meshing schemes is shown in Figure 3.As can be seen from Figure 3, the calculations derived from Mesh3 and Mesh4 are basically the same, which shows that continuing to refine the mesh size does not have much impact on the improvement of the calculation accuracy, and in order to take into account the calculation accuracy and the amount of calculations, the Mesh3 meshing scheme with a minimum mesh size of 0.2mm is selected for the calculation.left-leaning, centrosymmetric and right-leaning distributions, the law of the influence of location and distribution factors on the coalescence behavior is derived.In the simulation analysis of this paper, the double bubble diameter is set to 4mm, and the initial bubble is located at a distance of 5mm from the lower edge of the geometric model.As can be seen from Figure 4, under the parametric conditions of this paper, the single bubble rise in a non-Newtonian fluid exhibits an axisymmetric ellipsoidal spherical appearance, which matches well with the experimental results in the literature [16], and also proves the reliability of the model and the algorithm used in this paper.According to Figure 5, it can be obtained that a single bubble in a non-Newtonian fluid produces a symmetric vortex structure on both sides of it, and there is a vertically upward component of the velocity in this part, thus driving the bubble upwards.In addition, the bubble squeezes the non-Newtonian fluid, and a shear thinning effect occurs, forming a region of reduced viscosity around the bubble, a symmetric viscosity distribution on both sides below the bubble.

Coalescence behavior of vertical double bubbles in non-Newtonian fluids
The coalescence behavior of vertical double bubbles in different arrangements is discussed using the vertical spacing Δy=9mm as an example.

Coalescence behavior of vertical double bubbles with different distributions.
As shown in Figure 6(b), the rise in the vertical double bubble can be roughly described in three stages.When t=0.06~0.12S, the upper bubble is stretched along the horizontal side and the lower bubble is stretched along the vertical direction.At the same time, the two bubbles rise slowly, which is due to the high viscosity exhibited by the pseudoplastic fluid in its stationary state.When t=0.18~0.36S, the two bubbles close to each other, the lower bubble by the role of the upper bubble wake, its acceleration becomes larger and after a period of time to catch up with the upper bubble.When t=0.36~0.48S, the two bubbles complete the aggregation and form a long-tailed new bubble to continue to rise.It can be seen from Figure 7 that in the case of a bubble diameter of 4 mm, the rising velocity of the vertical double bubble is larger than that of the single bubble, and the velocity of the lower bubble is larger than that of the upper bubble, which is due to the nature of the non-Newtonian fluid in the wake region of the upper bubble in terms of shear-thinning, which results in accelerated movement of the lower bubble into the region.It is also worth noting that the shape of the new bubble after aggregation is significantly different from that of the single bubble represented in Figure 4.    distribution.It can also be seen that the longitudinal size of the bubbles after aggregation is significantly larger than that of the centered distribution in both cases due to wall effects, which do not have much effect on the skirt size of the bubbles for the weakly shear-thinned fluid in this algorithm [17].In addition, the bubbles show a "zigzag" rising path after aggregation, which is because the vortex structure around the bubbles under this distribution condition will interact with the wall to cause repulsive lifting force, resulting in the phenomenon that the bubbles are far away from the wall in the rising process.3.1.2Coalescence behavior of vertical double bubbles with different spacing.Simulation and comparison of the rising process of several centered double bubbles with different vertical spacing, the maximum vertical spacing Δy=14mm that can be achieved by the coalescence of vertical double bubbles in a non-Newtonian fluid at a bubble diameter of 4mm.This value is significantly smaller than in a Newtonian fluid under the same conditions, which is due to the high viscosity characteristic of the pseudoplastic fluid at rest or at low flow rates in this example.Figure 8 shows that when Δy≥16mm, the interaction between the two bubbles is weakened, and the catching up effect between the bubbles is not obvious, and the two bubbles rise individually, and the phenomenon of coalescence does not occur.

Coalescence behavior of parallel double bubbles in non-Newtonian fluid
The horizontal spacing Δx is the main factor for the coalescence of parallel bubbles or not.The horizontal double bubble is an ideal form of symmetry, so only the left-leaning and centered distributions are presented in Figure 9 for discussion.Figure 9 shows the motion of horizontal double bubbles in different distribution cases when Δx=5mm.As shown in Figure 9(a), When t=0.06~0.12S, a shear-thinning region of intersection is formed between the bubbles, which generates a viscosity reduction gradient and creates conditions for the completion of coalescence.When t=0.12~0.24S, the two bubbles collide and squeeze each other, the connection surface below between the bubbles gradually flattens, and the fluid above produces a viscosity reduction gradient, which also promotes the realization of bubble coalescence.When t=0.24~0.48S, the bubbles after completing the coalescence continue to rise with an elliptical appearance and exhibit the same shape as the single bubble rising process shown in Figure 4.By comparing the shapes of the bubble after aggregation, this also shows that the parallel arrangement of bubbles has less effect on bubble deformation.As shown in Figure 10, the rising speed of the two bubbles in the above case exhibits the same trend and eventually stabilizes at about 0.066 m/s, while the rising speed of the single bubble is slightly larger, and the reason for this difference lies in the fact that the vortex structure between the bubbles will have a hindering effect on the bubbles approaching each other horizontally.By comparing Figure 9(a) and (b), it is found that the aggregation patterns of parallel double bubbles are very similar in the two different distribution cases, but the rising speed of bubbles after aggregation is significantly slower in the left-skewed distribution case than in the centered case.This is because for non-Newtonian fluids with weak shear thinning effect, the pressure difference between the top and bottom sides of the bubbles gradually increases with the enhancement of the wall effect, which leads to a decrease in the velocity of the bubbles after aggregation [18].In addition, for the case of leftward distribution, the bubbles show a clear trend of zigzag upward movement, which shows that non-Newtonian fluids, the parallel arrangement of bubbles is more susceptible to the wall effect.Figure 11 shows the aggregation process of parallel double bubbles at Δx=18mm, under this condition the double bubbles failed to complete the coalescence, this is because as Δx increases to a certain extent, the local viscosity reduction characteristics of non-Newtonian fluids and the flow field structure of the inter-bubble interactions are significantly reduced, making it difficult for the aggregation to take place.After analysis, it is found that the bubble diameter of 4mm parallel double bubbles in the non-Newtonian fluid can be completed in the aggregation of horizontal spacing conditions for 3mm≤Δx≤16mm, more than this range, the interaction between the bubbles from mutual attraction to mutual repulsion.This indicates that in non-Newtonian fluids, although the local viscosity reduction property between parallel double bubbles gradually decays as the bubble spacing increases [19], it is also this property that makes the horizontal spacing at which parallel double bubbles in non-Newtonian fluids can be realized to aggregate much larger than that of the Newtonian fluids in the same case.bubble rising process.regularity; at the same time, the symmetry of the simulation analysis cannot be guaranteed, so this part only discusses the coalescence behavior of random double bubbles from the relative positions of the double bubbles.Double bubble with diameter 4mm, Δx=8mm, Δy=6mm for simulation analysis.As shown in Figure 12, When t=0.04~0.12S,each of the double bubbles rises vertically, and the shape of the bubbles remains spherical.When t=0.12~0.24S, the wake region of the upper bubble starts to act on the lower bubble, which produces a local viscosity reduction effect, under which the lower bubble's velocity gradually increases, while the rising routes of the two bubbles no longer keep the vertical direction.When t=0.24~0.48S, the two bubbles occur aggregation, the formation of new bubbles of surface oscillation to produce two directions of the lift, one so that the bubble to produce lateral motion, the other to maintain the bubble circular motion, so that it shows a "zigzag" upward movement tendencies.Random bubble coalescence process.Doubles bubble with diameter 4mm, Δx=16mm, Δy=8mm for simulation analysis, at this time, the double bubble did not happen to merge, its motion process is shown in Figure 12.As can be seen from the diagram, at the beginning of the movement of the bubble, the two bubbles each move in a vertically upward direction and do not affect each other.The two bubbles then move in opposite directions under the influence of the vortex structure of the fluid, and eventually the two bubbles rise in a spiral motion.

Conclusion
The coalescence behavior of double bubbles in non-Newtonian fluids with different arrangements is discussed using COMSOL, and the results of the study show that： (1) Under the same positional conditions, the stronger the shear rarefaction of the non-Newtonian fluid, the more likely the bubbles to merge; the initial spacing of the bubbles is also an important influence on the aggregation of bubbles; When vertical bubbles merge, the wake field of the upper bubbles drags the lower bubbles to form a catch-up effect to make the coalescence phenomenon occur.In this paper, when the bubble diameter is 4mm, the maximum distance between vertical bubbles is 14mm, and the maximum distance between horizontal bubbles is 16mm.Horizontal spacing has a greater effect on bubble coalescence than vertical spacing in random bubble coalescence.
(2) When the bubbles are distributed to the left or right, the shear rarefaction of the non-Newtonian fluid enhances the wall effect, making the bubble movement show a "zigzag" upward path; (3) In non-Newtonian fluids, parallel double bubbles rise similarly to single bubbles, the forms of vertical and random arrangements have a greater effect on the deformation of the bubbles.

Figure 1 .
Figure 1.Schematic diagram of the distribution of double bubbles.

Figure 2 .
Figure 2. Schematic diagram of the relative positions of the double bubbles.

Figure 3 .
Figure 3. Rising velocity of bubbles under four meshing schemes.

Figure 4 .
Figure 4. Deformation of a single bubble of an power-law fluid.

Figure 5 .
Figure 5. Velocity field and viscosity field of a single bubble of a power-law fluid.As can be seen from Figure4, under the parametric conditions of this paper, the single bubble rise in a non-Newtonian fluid exhibits an axisymmetric ellipsoidal spherical appearance, which matches well with the experimental results in the literature[16], and also proves the reliability of the model and the algorithm used in this paper.According to Figure5, it can be obtained that a single bubble in a non-Newtonian fluid produces a symmetric vortex structure on both sides of it, and there is a vertically upward component of the velocity in this part, thus driving the bubble upwards.In addition, the bubble squeezes the non-Newtonian fluid, and a shear thinning effect occurs, forming a region of reduced viscosity around the bubble, a symmetric viscosity distribution on both sides below the bubble.

For
both the left-skewed and right-skewed distributions, the double bubbles are 10 mm from the vertical symmetry axis of the geometric model.The bubble motion patterns shown in Figure 6(a) and (c) for the left-leaning and right-leaning distributions are generally the same as those for the centered t=0.06S t=0.12S t=0.18S t=0.24S t=0.36S t=0.48S(a) Δy=9mm double bubble leftward distribution.

Figure 6 .
Figure 6.Coalescence of Δy=9mm vertical double bubbles with different distributions.distribution.It can also be seen that the longitudinal size of the bubbles after aggregation is significantly larger than that of the centered distribution in both cases due to wall effects, which do not have much effect on the skirt size of the bubbles for the weakly shear-thinned fluid in this algorithm[17].In addition, the bubbles show a "zigzag" rising path after aggregation, which is because the vortex structure around

Figure 7 .
Figure 7. Rising velocity of vertical double bubbles as a function of time.

Figure 9 .
Figure 9. Coalescence of Δx=5 mm horizontal double bubbles with different distributions.Figure11shows the aggregation process of parallel double bubbles at Δx=18mm, under this condition the double bubbles failed to complete the coalescence, this is because as Δx increases to a certain extent, the local viscosity reduction characteristics of non-Newtonian fluids and the flow field structure of the inter-bubble interactions are significantly reduced, making it difficult for the aggregation to take place.After analysis, it is found that the bubble diameter of 4mm parallel double bubbles in the non-Newtonian fluid can be completed in the aggregation of horizontal spacing conditions for 3mm≤Δx≤16mm, more than this range, the interaction between the bubbles from mutual attraction to mutual repulsion.This indicates that in non-Newtonian fluids, although the local viscosity reduction property between parallel double bubbles gradually decays as the bubble spacing increases[19], it is also this property that makes the horizontal spacing at which parallel double bubbles in non-Newtonian fluids can be realized to aggregate much larger than that of the Newtonian fluids in the same case.

3. 3 Figure 10 .
Figure 10.Parallel Double Bubble Rising Velocity as a function of Time.

Figure 11 .
Figure 11.Centered distribution Δx=18mmbubble rising process.regularity; at the same time, the symmetry of the simulation analysis cannot be guaranteed, so this part only discusses the coalescence behavior of random double bubbles from the relative positions of the double bubbles.Double bubble with diameter 4mm, Δx=8mm, Δy=6mm for simulation analysis.As shown in Figure12, When t=0.04~0.12S,each of the double bubbles rises vertically, and the shape of the bubbles remains spherical.When t=0.12~0.24S, the wake region of the upper bubble starts to act on the lower bubble, which produces a local viscosity reduction effect, under which the lower bubble's velocity gradually increases, while the rising routes of the two bubbles no longer keep the vertical direction.When t=0.24~0.48S, the two bubbles occur aggregation, the formation of new bubbles of surface oscillation to produce two directions of the lift, one so that the bubble to produce lateral motion, the other to maintain the bubble circular motion, so that it shows a "zigzag" upward movement tendencies.

Figure 12 .
Figure 12.Random bubble coalescence process.Doubles bubble with diameter 4mm, Δx=16mm, Δy=8mm for simulation analysis, at this time, the double bubble did not happen to merge, its motion process is shown in Figure12.As can be seen from the diagram, at the beginning of the movement of the bubble, the two bubbles each move in a vertically upward direction and do not affect each other.The two bubbles then move in opposite directions under the influence of the vortex structure of the fluid, and eventually the two bubbles rise in a spiral motion.

Figure 13 .
Figure 13.Random bubble non-aggregation process.(1)Under the same positional conditions, the stronger the shear rarefaction of the non-Newtonian fluid, the more likely the bubbles to merge; the initial spacing of the bubbles is also an important influence on the aggregation of bubbles; When vertical bubbles merge, the wake field of the upper bubbles drags the lower bubbles to form a catch-up effect to make the coalescence phenomenon occur.In this paper, when the bubble diameter is 4mm, the maximum distance between vertical bubbles is 14mm, and the maximum distance between horizontal bubbles is 16mm.Horizontal spacing has a greater effect on bubble coalescence than vertical spacing in random bubble coalescence.(2)When the bubbles are distributed to the left or right, the shear rarefaction of the non-Newtonian fluid enhances the wall effect, making the bubble movement show a "zigzag" upward path;(3) In non-Newtonian fluids, parallel double bubbles rise similarly to single bubbles, the forms of vertical and random arrangements have a greater effect on the deformation of the bubbles.

Table 1 .
Values of flow characteristic parameters for non-Newtonian fluids.