Simulation Analysis of Gas Bubble Formation and Escape in Non- Newtonian Drilling Fluids

In this study, the formation and escape movements of a bubble injected in non-Newtonian drilling fluid through a pore were numerically simulated using a volume of fluid method. The pattern of a single bubble and the pressure and velocity fields of the surrounding liquid phase during the bubble formation were analyzed and compared with experimental results; based on the comparison, the formation and escape properties of the bubble were further studied. In particular, the effects of static shear force, consistency coefficient, and flow behavior index on the growth and escape time of the bubble were analyzed. The results show that, owing to the effect of velocity on the viscosity of a non-Newtonian drilling fluid, the escape time and volume of the bubble increase with an increase in static shear force, consistency coefficient, and flow behavior index. Among the three parameters, the flow behavior index has the greatest effect. This is because the shear disturbance of a bubble to its surrounding fluid during its growth and escape, caused by the shear thinning of a yield-power-law fluid, reduces the fluid viscosity. The shear thinning decreases, and the resistance to the bubble increases as the flow behavior index approaches 1, leading to larger bubble formation times and separation volumes. An empirical formula for predicting the equivalent radius of bubbles considering the liquid yield stress, inertial force, viscous force, and surface tension is established. The average error of predicting the equivalent radius of detached bubble is 0.80%, which can provide a reference for the better study of bubble migration and flow pattern in non-Newtonian fluid.


Introduction
The contact between the gas and liquid phases is an important phenomenon in petroleum engineering [1][2][3][4][5][6][7]. The gas-liquid contact usually occurs in two forms: either the liquid enters the gas in the form of droplets, or the gas is bubbled into the liquid. In engineering applications, such as absorption, distillation, flotation, and gas cut during well drilling, all the gas-liquid interactions occur in the latter form. The hydraulic characteristics and flow regimes of the gas-liquid system are affected by the physical and chemical properties of the liquid phase (such as viscosity, surface tension, and density) and by the characteristic parameters of the dispersed phase (such as bubble size and bubble rising velocity and among others) [8][9][10][11][12][13]. The characteristics of a bubble in the liquid phase, the different flow regimes, and the distribution of the flow field of the liquid under the effect of bubble movement are helpful for the indepth understanding of the transfer process during gas-liquid contact and have been extensively studied. For instance, Raymond and Rosant [14] conducted experimental and numerical studies on the behavior of spherical and ellipsoidal bubbles. Bhaga and Weber [1] studied experimentally the different shapes of bubbles rising in an immobile fluid, obtaining the shapes of bubble trails and the general characteristics of a velocity streamline. Zenit and Magnaudet [15] and Cano-Lozano et al. [16] conducted experimental studies on the starting point of path instability. To date, great progress has been made in the study of the bubble rising behavior in viscous liquids [10,[16][17][18][19]. Numerical simulations of bubble formation and escape have also been seldom conducted. Albadawi et al. [20] studied the growth and escape of an isolated bubble injected through a single pore using four different interface-capture methods. Boubendir et al. [21] studied the role of the surface tension in the growth and separation mechanism of bubbles in the gas-water flow in micropipes. Gumulya et al. [12] simulated numerically the movement of bubbles in a viscous Newtonian fluid and calculated the dynamic shape, resistance, and average bubble rising velocity for different bubble sizes in a viscous fluid and the corresponding flow fields. Additionally, they analyzed the expansion and formation of the closed/open laminar flow and trail after bubble rising using the 3D volume of fluid (VOF) method. Islam and Bergstrom [22] discretized the Reynolds-averaged Navier-Stokes (RANS) equations of a two-fluid model using the finite volume method and analyzed the effect of turbulence on the liquid flow field using the transport equation for the turbulence kinetic energy and dissipation rate.
The aforementioned studies focused primarily on the pressure, velocity, disturbance, and bubble collision and coalescence during the two-phase gas and liquid movement and their effect on bubble formation and migration. Most fluids used in petroleum engineering and especially in drilling, such as drilling fluids and mud slurry, are yield-power-law fluids. The Herschel-Bulkley (HB) model is widely used to characterize the dynamic shear force of a fluid and to represent the shear thinning of a non-Newtonian fluid. However, few studies have been conducted on the migration of yieldpower-law drilling fluids, especially those addressing gas. The characteristics of bubble movement in a non-Newtonian drilling fluid include the bubble shape, movement trajectory, rising velocity, rupture, and the coalescence between bubbles [23][24][25]. The initial gas flow into the liquid phase is the basis for the movement analysis, including the bubble formation velocity, size, internal and external flow, pressure field, and initial bubble shape. Numerical analyses can substitute a considerable amount of experimental work in determining the effect of different parameters [26][27][28][29]. Through numerical simulation, this study analyzes the bubble formation and escape from the wall surface when the gas phase enters the yield-power-law drilling fluid phase from a pore. Furthermore, this study reveals the effect of the fluid-dynamic shear force, consistency coefficient, power-law index, and other physical parameters on the initial shape of a bubble, providing a reference for further studies on the migration and flow regimes of bubbles in non-Newtonian drilling fluids.

Model Development
A momentum equation describes the properties of momentum conservation in a flow system. For the formation and migration of bubbles in the liquid phase, the effect of surface tension is to be considered. Then, the momentum equation is where u ! is the flow velocity, p is the flow field pressure, μ is the dynamic viscosity, g is the gravitational acceleration, F sf ! is the surface tension, and ρ is the average density for volume fraction, ρ = α 2 ρ 2 + ð1-α 2 Þρ 1 . When the flow is turbulent, the k − ε turbulence model is used to approximate the RANS equations, using Here, G k is the kinetic energy of the turbulent flow generated by the average velocity gradient, G k = −ρu i ′ u j ′ ð∂u j /∂x i Þ, G b is the kinetic energy of the turbulent flow generated by buoyancy, G b = βg i ðμ t /Pr t Þð∂T/∂x i Þ, and β is the coefficient of thermal expansion, β = −1/ρð∂ρ/∂TÞ p , Pr t = 0:85. μ t is the viscosity of the turbulent flow, determined using the equation μ t = ρC μ ðk 2 /εÞ. C μ = 0:0845, C 1ε = 1:42, C 2ε = 1:68: 2.2. Interface Tracking Method. The VOF method [30] is used to track and reconstruct an interface. In this method, the volume fraction α q is defined as 1 at a point of fluid q and 0 at a point outside it. Thus, the average value of α q in a cell represents the fraction of fluid in that cell. If the average value of α q is 1, the grid is filled with the fluid; if it is 0, there is no fluid in the cell; if it is between 0 and 1, the grid is a free-surface one.
Because it tracks the volume of fluid in the grid rather than the movement of fluid particles, the VOF method is very practical, requiring low computation power and providing high accuracy. The volume fraction equations used in this method are [31,32] ∂α q ∂t and for the two-phase flow, Modified Hershel − Bulkley model : where τ 0 and τ are the dynamic and static shear force in Pa, respectively, k is the consistency coefficient in Pa·s n , γ · is the second invariant of the rate-of-strain tensor, and 5 is the fluidity index.
For the two-dimensional flow, the explicit form of the shear rate is given by [34][35][36] where ðx, yÞ are the dimensional coordinates, and ðu, vÞ are the velocity components along ðx, yÞ.
The Reynolds number of Hershel-Bulkley fluid was derived as follows: The greater deviation of n from 1, the stronger non-Newtonian property of fluid. When the value of n is less than 1, the apparent viscosity decreases with the increase of the shear rate. When the value of n is greater than 1, the apparent viscosity increases with the increase of the shear rate. Both the Bingham model and power-law model are considered in the Hershel-Bulkley model.
In the simulation process, there is a critical shear stress. For the modified Hershel-Bulkley model, set τ = τ 0 , ð γ · Þ c is given by the following equation.  γ · Þ c is the critical shear rate, s -1 . When τ > τ 0 , the modified Hershel-Bulkley model can be turned into the Hershel-Bulkley model. When τ < τ 0 , it reveals very low shear rate, the fluid is affected by yield stress and produces strain, and the modified Hershel-Bulkley model is used. Set τ 0 = 2 Pa, k = 5 Pa · s n , and n = 0:5, N = 10000, and the relationship between shear force and shear rate is obtained. The relationship between shear force and shear rate is shown in Figure 1.

Meshing
Partition. The numerical discrete method was used to solve the finite volume method, and a structured grid was used to generate the space grid. Before the numerical calculation, we compared the triangular mesh with the quadrilateral mesh. It is found that quadrilateral mesh is better than triangular mesh in terms of computation time and convergence speed. In this work, quadrilateral mesh is used. In terms of the number of grids, too sparse grid will lead to inaccurate calculation results and even cause the calculation results not to converge in some cases, and too dense grid will not only increase the amount of calculation greatly, but also have higher requirements for computer hardware. We selected 3 kinds of quadrilateral mesh numbers (200000, 246000, 260000). In order to make the results closer to the reality, the number of grids should be large enough that with the increase of the number of grids, the results will not change significantly. Through simulation, the number of grids selected in this paper is 246000. The mesh generation is shown in Figure 2. According to the grid quality inspection function, the worst grid skewness is less than the standard value of 0.97; so, the grid quality is better, and all the grids are evenly distributed in the whole fluid calculation domain. Therefore, the calculation with this grid can meet the engineering accuracy requirements.

Solving Method and Boundary Condition.
For the numerical simulations, we used a single-precision finite volume method. For discretization, an implicit scheme was applied to time, and a second-order upwind scheme was applied to the convective term. The pressure-velocity field coupling was obtained by solving the corresponding equation using the pressure-implicit with splitting of operators (PISO) algorithm in the unsteady compressible or incompressible fluid flow field. This was proposed to solve the asynchronous correction problem of the momentum and mass continuity equations in the semi-implicit method for pressure-linked equations (SIMPLE) series algorithm. The main idea is to add a velocity correction step after the pressure one in the SIMPLE algorithms; so, that the iterative equation satisfies the mass conservation equation explicitly and the momentum conservation equation implicitly. The time step is 1 × 10 −3 s, velocity convergence criteria is 1 × 10 −5 m/s, and pressure convergence criteria is 1 × 10 −5 Pa. respectively. The rectangular simulation domain is first meshed into different amounts of the structured grid (see grid independence analysis) by using a commercial software pack-age named GAMBIT (Ansys Inc., US). Thereafter, the CFD code FLUENT (Ansys Inc., US) is applied to calculate equations.
2.6. Validation. Using the numerical method established, the bubble formation and escape from the wall surface in water were first simulated and then compared with the experimental results of Albadawi et al. [20]. Figure 1 shows a comparison of the simulated and experimental results on bubble formation at different times. The experimental conditions and computational parameters used are listed in Tables 1  and 2. Figures 3(a)-(e) show the bubble shapes when t/t det = 0 , 0:4, 0:6, 0:8, and 1, respectively. Gas-liquid interfaces can be seen in different working conditions. In the figure, the left column displays the bubble shapes and sizes obtained experimentally by Albadawi et al., and the right column displays those obtained from the numerical simulation. As seen, the simulated results are in close agreement with the experimental ones. Therefore, it is concluded that the bubble formation and separation can be accurately reproduced through computational simulation, and therefore, numerical analysis of the processes can be safely conducted.

Results and Discussion
Using the above method, the process of gas entering a non-Newtonian fluid at the velocity of 0.1 m/s and    3.1. Bubble Formation Process and Shape Analysis. Figure 4 shows a simulation of the process starting from the bubble formation to its escape for a pore size of 1 mm and pore velocity of 0.1 m/s. At this point, the liquid phase parameters were τ 0 = 10 Pa, k = 2 Pa · s n , and n = 0:4. The volume fractions of gas are represented in different colors. As observed, the bubble went through two stages in its formation process: an expansion stage and a stretching-escape one. First, the bubble expanded radially, though its lower part was still in contact with the pore. Then, with the continuous growth of the bubble, the lower side rose and formed a thin neck, by which the bubble contacted with the pore. When the buoyancy of the bubble was greater than the resistance, the bubble escaped from the pore and moved upwards. Figure 5 shows a cloud diagram of the pressure distribution for the formation and escape of a bubble. As seen, the pressure was not constant during the formation. At a defined gas injection velocity, a bubble emerged gradually, and the pressure at the initial time was the highest. At an internal pressure of 4920 Pa and an external hydrostatic fluid column pressure of 4870 Pa, the bubble kept expanding and growing. Similar to a balloon being blown, as the bubble grew, the pressure decreased gradually. After the bubble grew to a certain size and escaped from the wall surface, the internal pressure became slightly higher than the external hydrostatic fluid column pressure, nearly equal to the ambient pressure, and equivalent to the hydrostatic fluid column pressure at the bottom of the escape position, which was 4830 Pa. The hydrostatic fluid column pressure at the top of the bubble was 4800 Pa. Figure 6 shows the velocity distribution in the gas and liquid phases during the formation, growth, and escape of a bubble. When the bubble formed, the velocity became distributed laterally and vertically inside the bubble; simultaneously, the bubble expanded in the lateral and vertical directions. After growing to a certain size, under the joint actions of the buoyancy, viscous force, and surface tension, the bubble began to grow vertically. Then, the internal vertical velocity became dominant. Its magnitude is related to the bubble shape. The gas flow rate at the necking increased to 5-10 times of the inlet injection velocity, and the pressure decreased, which led to the further development of the necking until the bubble detached. After the bubble separation, the internal velocity decreases gradually. With the growth of the bubble, the external liquid phase fluid will be disturbed, and the eddy current will appear gradually outside the bubble, but the eddy speed is not large. In the simulation conditions, the maximum transverse velocity of liquid is 0.06 m/s.

Effect of the Yield-Power-Law Properties on Bubble
Formation and Escape. In addition to the density difference between gas and liquid and the surface tension of liquid, the non-Newtonian fluid properties such as liquid static shear force, consistency coefficient, and rheological index also have significant effects on bubble shape formation, separation time, the volume, etc. ( [37,38]). However, the current study-only qualitative analysis has been carried out, and the quantitative relationship between the characteristic parameters of bubble separation at the orifice and the yield power-law characteristics of non-Newtonian fluids has not been established. In this section, the influence of different static shear force, consistency coefficient, and rheological index on bubble formation and separation is simulated, which provides a basis for establishing the quantitative relationship between the bubble volume and non-Newtonian fluid properties.

Static Shear
Force τ 0 . The method developed in this study was used to simulate the formation and escape processes of a bubble by changing the static shear τ 0 while maintaining the gas physical properties, liquid density, flow behavior index, and consistency coefficient constant of the non-Newtonian fluid. The results obtained with a τ 0 value of 5 Pa and 25 Pa are shown in Figures 7 and 8.
As shown in Figure 7, when τ 0 increased from 5 Pa to 25 Pa while maintaining the other parameters constant (k = 2 Pa · s n , n = 0:4), the bubbles initially had a uniform growth process involving formation, radial growth, stretching, and neck reducing breaking. However, they displayed different escape times and volumes. In fact, with the bub-ble growth, the greater the static shear force of the liquid, the longer the time for the bubble to escape, and the larger the bubble volume when it escaped. When τ 0 increased from 5 Pa to 25 Pa, the bubble escape time increased from 0.30 s to 0.335 s, and the escape volume increased from 23.55 mm 3 to 26.30 mm 3 . At different static shear forces, the velocities concentrated mainly in the bubble and liquid phase near the bubble wall surface. With the growth of the bubble, the bubble wall surface caused disturbance to the liquid phase, and a small vortex appeared. In the formation and escape of a bubble, the velocity reached its maximum just before the bubble escaped, and it decreased with the increase of the static shear force of the liquid phase, as shown in Figure 8. The main cause is that, when the static shear force in the non-Newtonian fluid is large, the bubble has a large lateral size. At a constant gas flow rate, the bubble has a greater size than its reduced neck, and the velocity distribution is relatively more uniform. The statistical results on bubble escape time, maximum velocity, and volume in liquids with different static shear forces are shown in Table 3.

Consistency Coefficient k.
Using the same method, the effect of the consistency coefficient k of a non-Newtonian yield-power-law fluid on the formation and escape processes of a bubble was analyzed. Figure 9 shows the formation and escape processes of bubbles when k = 1 Pa · s n and 5 Pa·s n . In the first part of the formation, the bubbles had similar shapes and sizes in the two fluids. In the last part of the process, the bubbles in the fluid with the smallest consistency coefficient had a faster stretching-breaking velocity than those in the fluid with a larger consistency coefficient. For example, at 0.25 s, the bubbles in the former had a reduced neck, whereas those in the latter displayed nearly constant sizes in the upper and lower parts, as shown in Figure 9(b). At 0.29 s, the bubbles in the former broke and escaped with a volume of 22.76 mm 3 , whereas it took 0.36 s for the bubbles in the latter to escape, with a volume of 32.97 mm 3 . The main cause is that, when increasing the consistency coefficient, the viscosity of the liquid phase and the viscous force that the bubble has to overcome to escape increases. Therefore, it took a longer time for the bubbles to escape, and the bubbles grew larger. The statistical results on   Table 4.

Flow Behavior Index n.
For a non-Newtonian fluid, the flow behavior index k has a greater effect on the formation and escape of a bubble than k and τ 0 . Figure 10 shows the formation and escape processes of bubbles when n = 0:2 and 1.4. The commonly used drilling fluid has shear dilution property (n < 1). If n > 1, the circulation of drilling fluid will be difficult. As observed, the formation and escape processes of a bubble in the liquid having n = 0:2 were faster than those of a bubble in the liquid having n = 1:0. In fact, the bubble escape time was 0.28 s in the former and 0.89 s in the latter. In the latter, the lower part of the bubble had a long tail that connected it with the next bubble formed under the action of a highly viscous liquid phase, as shown in Figure 10(b). A yield-power-law fluid is a pseudoplastic fluid, in which a bubble grows fast. The shear thinning of a pseudoplastic fluid can reduce both the viscosity of the fluid around a bubble and the viscous resistance acting on the bubble. Therefore, it was more likely for the bubble to move in the tangent direction at the early stage of growth, increasing the probability of bubble elongation. That is, the smaller the value of n, the greater the shear thinning, and the greater the probability of bubble elongation. When n = 1, the bubble had a large roundness, and it was not likely to escape from the wall surface. The statistical results of the bubble escape time, maximum velocity, and volume in liquids with different flow behavior indexes are shown in Table 5.

Analysis of the Results.
According to the results of Tsamopoulos and Dimakopouloset [39] and Dimakopoulos et al. [40], at a certain bubble volume, the gas-liquid interface is in equilibrium under the actions of the pressure inside a bubble, hydrostatic pressure of the fluid column outside it, surface tension, yield stress, viscous force, inertial force, etc. That is, when the gas enters from a pore, the bubble shape depends on the properties of the fluid. In a yield-stress fluid, the suspended and migrating bubbles generally have spherical or inverted teardrop shapes, and the distribution of the fluid stress and velocity fields around the bubbles differs. To date, it is generally accepted that the main forces inducing bubble escape are the rising velocity of the gas, the buoyancy of the bubble, and the convection of the liquid caused by the formation, growth, and upward movement of the bubble [41,42]. In contrast, the main forces that prevent bubble escape are the surface tension, viscous resistance, and reactive force of the liquid during the formation and growth of the bubble. When the resultant upward force on the bubble is greater than the resultant downward one, the bubble escapes [43][44][45].
In a non-Newtonian yield-power-law fluid, the vertical expansion velocity of a bubble is greater than its radial expansion velocity, and the bubble has an elliptical shape connected to the nozzle; that is, the bubble shape is easier elongated than in a Newtonian fluid [19]. The reason for the different shapes is that the rheological properties of the two fluids are different. During the formation and expansion of a bubble, the disturbance to the liquid flow and the rheological properties of the fluid result in residual stress accumulation, which leads to a decrease in the viscosity of the fluid inside the bubble tail during stress dissipation, reducing the effect on the viscous force of the bubble. The viscosity distribution of a yield-power-law fluid during bubble formation and escape is shown in Figure 11. According to the results,     10 Geofluids increasing τ 0 , k, and n all induced a viscosity increase. Additionally, the smaller the value of n, the greater the shear thinning and the lower the viscosity at constant velocity. The effect of viscosity can better illustrate the stretching and escape processes of a bubble during its growth. An increase in viscosity will make the drag force increase. When the resistance to a bubble increases, both the bubble formation time and separation volume increase. In general, the greater the viscous force, the longer the bubble tail at the escape time.

Dimensionless Analysis.
The buoyancy is the driving force while the yield stress and viscous force are the resistance during the growth and rise of the bubble. In order to analyze the variation of the bubble size with fluid properties, dimensionless analysis of numerical simulation results is carried out in this section. The dimensionless number Y g is introduced to represent the ratio of yield stress to inertial force (buoyancy), the dimensionless number B n represents the ratio of yield stress to viscous force, and the dimensionless number B o represents the ratio of inertial force to surface tension [39,40], which are defined in equations (8)- (10).
The dimensionless number is defined as the ratio between the bubble volume radius and the orifice radius.
where R b is the equivalent radius of the bubble, R 0 is the radius of the orifice, and V b is the volume of the detached bubble.
The dimensionless parameters of the detached bubbles are obtained under different yield stress, consistency coefficient, and fluidity index, as shown in Table 6. The variation of dimensionless radius of bubbles with various parameters is shown in Figures 12-14. The dimensionless equivalent radius of detached bubbles has no obvious change with Y g , B n , and B o .
Both the yield stress and the viscous force are the resistance of bubbles falling off. In addition, the surface tension affects the resistance by affecting the shape of bubbles, which is also an important reason for the volume of bubbles falling off. In this work, the influence factors of three forces on the equivalent radius of bubbles falling off are considered, and the empirical relationship as shown in Eq. (14) is established.

Geofluids
The results are shown in Figure 15. Compared with the results of numerical simulation, the average error of the model is 0.80%, and the maximum error is 2.42%.

Conclusions
A numerical model has been developed for the formation and escape of a bubble introduced into a non-Newtonian yieldpower-law fluid through a single channel. The velocity, pressure, and viscosity distribution around a bubble during its formation and escape have been obtained. The effect of the physical parameters of the fluid (dynamic shear force, consistency coefficient, and power-law index) on the formation and escape behavior of a bubble has been analyzed. The results showed that, because of the different viscosities, the escape time and volume of a bubble are larger for increasing static shear force, consistency coefficient, and fluidity index. Among the parameters considered, the flow behavior index has the greatest effect, mainly because the shear disturbance of a bubble to its surrounding fluid during its growth and escape, caused by the shear thinning of yield-power-law fluids, reduces the fluid viscosity. The closer the value of n to 1, the lower the shear thinning and the greater the resistance to a bubble, causing larger bubble formation times and separation volumes. The buoyancy is the driving force in the process of bubble growth and rising, and the yield stress and viscous force are resistance. Through dimensionless analysis, an empirical formula for predicting the equivalent radius of bubbles considering the yield stress, inertial force, viscous force, and surface tension is established. The average error of pre-dicting the equivalent radius of detached bubble is 0.80%, which can provide a reference for the better study of bubble migration and flow pattern in non-Newtonian fluid.

Data Availability
All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
The author(s) declare(s) that they have no conflicts of interest.  12 Geofluids