Simulation analysis of bubble coalescence behavior characteristics in 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 significantly impacts bubbles’ size, shape, and movement. These parameters are essential for interphase mass, heat transfer processes, and equipment performance. This paper presents a simulation study of the bubble coalescence behavior of 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 Newtonian fluids. The results show that the relative distance of bubbles is the crucial factor for the occurrence of bubble coalescence in Newtonian fluids, and different distributions show different laws on the motion of bubbles.


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, investigating bubble coalescing law is integral to studying bubble dynamics behavior and chemical equipment process enhancement.
For the study of bubble coalescence behavior, Chen et al. [2] used the idea of microlayer modeling to simulate the single bubble growth and detachment process numerically.They analyzed the effect of the gas-liquid contact angle on the device's mass and heat transfer process.Huang et al. [3] discussed the deformation characteristics and laws of single bubbles experimentally in the rise of fluids with different viscosities through a self-designed bubble generation device.Ye et al. [4] numerically solved the kinetic regulation of the double bubble coalescence process.They used the level set method combined with MATLAB to simulate the coalescence process of double bubbles numerically.Wen [5] conducted an experimental study on the factors influencing the coalescence of multiple bubbles in pure water systems utilizing an in-house designed experimental setup, verified by numerical simulations using a multiple grid algorithm.Maxim K [6] empirically studied double bubble collision and fragmentation phenomena.
From the available literature, it is known that there are few studies related to coalescence behavior of bubbles in 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 Newtonian fluids under different arrangements, and to explore the motion law of bubbles in this process, which provides some reference for related fields.

Physical model
The coalescing behavior of bubbles affects the equipment's heat and mass transfer index.According to the conclusion of the literature [7], the bubble diameter is the main factor affecting its coalescence, the double bubble diameter in this paper are set to 4mm.In addition, in the simulation analysis of this paper, the initial bubble is located at a distance of 5mm from the lower edge of the geometric model.According to the literature [8], when the distance between the bubble's center and the wall is greater than five times the bubble's radius, 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 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, centered 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. Adjusting these two parameters makes the position relationship between the double bubbles appear vertical, parallel, and random.
2.2 Control Equations 2.2.1 Continuity Equation.ρm is the average fluid density; um is the mean fluid velocity vector, whose value can be obtained from Equation (3); xm is the average fluid parameter; αi is the volume fraction of phase i; xi is the fluid parameter of phase i.

(
) .p is the pressure on the fluid micro-element; μm is the average hydrodynamic viscosity; F is the surface tension source term.

Phase field interface Equation.
Numerical simulation methods commonly used for gas-liquid twophase flow include the volume of fluid, level set, and phase field methods.The volume of fluid method requires strict quality of the Mesh, and the horizontal set method has limited accuracy and is prone to mass non-conservation [9].Compared with the first two methods, 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 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.The phase field interface control equation used in this paper is the Cahn-Hilliard's non-linear diffusion equation containing the fourth-order diffusion term, which can achieve continuous variation of physical quantities such as density and viscosity in this region.
In Equation ( 6), φ 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 water are -1 and 1.The gas-liquid interface is where the phase field variable φ=0 can determine the bubble's position at different moments.ψ denotes the chemical potential also known as the phase field auxiliary variable and is calculated from Equation (7).Equation (7) in the ε 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 generally accepted as γ=ε 2 ; λ denotes the mixing energy density, whose value satisfies Equation (8).In Equation (8), σ is the surface tension coefficient, and the surface tension is calculated from Equation (4) to find.

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
In the geometric model shown in Figure 1, the initial and relative positions of the double bubbles can be adjusted accordingly depending on the study target.The Newtonian fluid used in this paper is the most common form of air bubble flow in production water, where the gas-liquid two-phase's physical parameters are shown in Table 1.Since the inflow of 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 condition is the gas phase inlet velocity and is set to a constant flow rate of 0.03m/s.The upper part of the geometric model is the outlet boundary condition, and the pressure value is set to 1 standard atmosphere.The Mesh's size and quality directly determine the simulation results' computational accuracy.Use the Mesh refine function in COMSOL to encrypt the local mesh at the junction of the two phases, use the boundary layer mesh for the model's boundary, and use the refine mesh option in other areas of the model.After using the above settings, the mesh distribution is shown in Figure 3, with 27,000 calculated mesh cells and an average mesh mass of 0.92.  is slightly slower than that of the centered distribution, and the bubbles show "zigzag" ascending path after coalescence.This is because the double bubble in this distribution condition, in the process of rising, due to the Newtonian fluid internal surface without shear conditions generated by the vortex volume in the respective wake area will interact with the wall surface to cause repulsive lift, resulting in the bubble rising process away from the wall phenomenon.Therefore the rise rate of the bubbles under this condition also exhibits an oscillatory form as shown in Figure 6, which also agrees with the literature [11] findings.

Coalescence 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=16mm can be achieved by the coalescence of vertical double bubbles in Newtonian fluid at a bubble diameter of 4mm. Figure 7 shows the rising process of the centered double bubble with vertical spacing Δy=18mm, it can be found that under this condition, the interaction between two bubbles is weakened, and the two bubbles rise separately, and chasing effect between bubbles in the rising process is not apparent.The phenomenon of coalescence does not appear.

Coalescence behavior of parallel double bubbles in Newtonian fluid
The coalescence behavior law of parallel bubbles is more complicated, and the horizontal spacing Δx is the main factor for the coalescence of parallel bubbles or not.In addition, the flow after coalescence of bubbles in different distribution cases also has various phenomena.
t=0.04S t=0.07S t=0.1S t=0.17S t=0.23S t=0.36S(a) Δx=3.5mm double bubble-centered distribution.The horizontal double bubbles are an ideal form of symmetry, so only the left-leaning and centered distributions are presented in Figure 8 for discussion.By adjusting Δx, it is found that the coalescence of horizontal double bubbles can be completed when 3mm ≤ Δx ≤ 5mm, and Figure 8 shows the motion of flat double bubbles in different distribution cases when Δx=3.5mm.From Figure 8(a), it can be seen that the two bubbles gradually come together under the force provided by the individual vortices until the completion of the coalescence.A narrow local vortex region is formed above the two bubbles, which drives coalescence [12].The new bubbles formed after aggregation continue to move in a vertically upward direction.The horizontal double bubbles with the leftward distribution shown in Figure 8(b) have the same movement pattern as the case of centered allotment before coalescence.Still, it takes a longer time to complete the coalescence.In this case, the side of the bubble near the wall shows a flattened shape due to the wall effect, and the lower surface has been clearly concave.In addition, its motion also offers a non-linear motion characteristic.9, when Δx=6mm, the two bubbles will move in the opposite direction at the same rate, resulting in the inability to complete the coalescence.As Δx increases, the difference between the rise velocity of the double bubble and the single bubble becomes smaller and smaller, as shown in Figure 10, which means that the interaction force between the double bubble becomes smaller and smaller under this condition..

Coalescence behavior of random bubbles in Newtonian fluids
When the double bubbles are randomly distributed, the coalescence behavior is more difficult to have regularly; simultaneously, 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 bubbles with diameter 4mm, Δx=2mm, Δy=6mm for simulation analysis.As shown in Figure 11, When t=0.04S, the double bubbles each rise vertically, and their shapes change from the initial round to ellipsoidal.When t=0.07~0.1S, the velocity of the lower bubble gradually increases and approaches the upper bubble under the effect of the wake of the upper bubble, and the rising routes of the two bubbles no longer maintain the vertical direction at this stage.When t=0.17S, the two bubbles began to merge, and after the coalescence, the bubbles show "zigzag" rising motion law.Random bubble coalescence process.Doubles bubble with diameter 4mm, Δx=10mm, Δ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 in this position, the two bubbles each move in a vertically upward direction and do not affect each other.When t=0.17~0.36S,due to the vortex effect of the surrounding flow field and the effect of the trailing area of the bubble below by the bubble above, the two bubbles eventually rise in a spiral motion.

Conclusion
The coalescence behavior of double bubbles in Newtonian fluids with different distributions and relative positions is discussed using COMSOL, and the results of the study show that： (1) The initial spacing in Newtonian fluid determines the possibility of bubble aggregation.In the Newtonian fluid of this paper, when the bubble diameter is 4 mm, the maximum distance between vertical double bubbles to achieve coalescence is 16 mm; (2) The spacing conditions for achieving the coalescence of parallel bubbles are more demanding.In the Newtonian fluid of this paper, when the bubble diameter is 4 mm, once the horizontal spacing exceeds 5 mm, the two bubbles cannot be aggregated; (3) When random bubbles coalesce, the horizontal spacing has a more significant effect on the coalescence of bubbles than the vertical spacing; (4) When the bubbles are distributed to the left or the right, the wall effect will cause the bubbles to show a "zigzag" ascending path after aggregation.

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.According to the literature[8], when the distance between the bubble's center and the wall is greater than five times the bubble's radius, 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 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, centered and right-sided

Figure 5 .
Figure 5. Coalescence of Δy=9mm vertical double bubbles with different distributions.isslightly slower than that of the centered distribution, and the bubbles show "zigzag" ascending path after coalescence.This is because the double bubble in this distribution condition, in the process of rising, due to the Newtonian fluid internal surface without shear conditions generated by the vortex volume in the respective wake area will interact with the wall surface to cause repulsive lift, resulting in the bubble rising process away from the wall phenomenon.Therefore the rise rate of the bubbles under this condition also exhibits an oscillatory form as shown in Figure6, which also agrees with the

Figure 8 .
Figure 8. Coalescence of Δx=3.5mm horizontal double bubbles with different distributionsAs shown in Figure9, when Δx=6mm, the two bubbles will move in the opposite direction at the same rate, resulting in the inability to complete the coalescence.As Δx increases, the difference between the rise velocity of the double bubble and the single bubble becomes smaller and smaller, as shown in Figure10, which means that the interaction force between the double bubble becomes smaller and smaller under this condition..

Figure 11 .
Figure 11.Random bubble coalescence process.Doubles bubble with diameter 4mm, Δx=10mm, Δ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 in this position, the two bubbles each move in a vertically upward direction and do not affect each other.When t=0.17~0.36S,due to the vortex effect of the surrounding flow field and the effect of the trailing area of the bubble below by the bubble above, the two bubbles eventually rise in a spiral motion.
distribution.In this case, due to the influence of the vertical wall, the catching speed of the lower bubbles