Bubble Formation in Yield-Stress Fluids Using a Coupled Level-Set and Volume-of-Fluid Method

The bubble formation under the influence of orifice diameter submerged in yield-stress fluids was studied numerically using a coupled level-set (LS) and the volume-of-fluid (VOF) method and was in contrast with that seen in Newtonian fluids. The bubble formation process had a good consistency by virtue of comparing numerical simulation and experiment. The process of bubble formation could be divided into two parts, bubble expansion stage and stretch stage. The influence of orifice diameter and liquid rheological characteristics (consistency coefficient, flow index, and yield stress) on the bubble formation parameters (the bubble formation time, detachment volume, and aspect ratio) was investigated. The results revealed that the bubble detachment volume increases with the increase of orifice diameter, consistency coefficient, flow index, and yield stress. In different kinds of fluids, the formation time and detachment volume of bubbles in the shear-thinning fluid were the lowest, followed by the Newtonian fluid, and finally, in yield-stress fluids.


■ INTRODUCTION
The phenomenon of the bubble formation in fluids is widely used in a number of industrial processes, such as chemical, environmental, biological, and mineral processes. Virtually, the knowledge of the bubble motion characteristics is of the enormous importance to optimize the above-mentioned processes. The typical motion characteristics of bubbles are principally comprehensive of formation, rising, coalescence, deformation, and fracture. Especially, the formation characteristics of bubbles have an important and high impact on the shape, volume, and velocity, which further affects the specific surface area and residence time of bubbles in the liquid phase, and finally, the mass and heat transfer as well as the reaction rate between gas and liquid phases are significantly influenced. 1 Consequently, a deep-and-meaningful comprehension of bubble formation characteristics will be vitally favorable for completing the optimization design of gas−liquid two-phase related industrial processes. In view of this, a number of scholars have been interested in the formation characteristics of bubbles in fluids. In the early year, Davidson and Schuler probed into the bubble formation in the low-viscosity fluid, and a mathematical model predicting the spherical bubble formation volume was put forward based on the analysis of forces. 2 However, the Davidson model is a one-stage model. Subsequently, the bubble formation process was divided into two stages (expansion and stretching) by Kumar and Kuloor, and the two-stage formation volume prediction model was proposed on the basis of the Davidson model. 3 More information on the bubble formation models can be achieved from the review of Tsuge and Hibino. 4 As a matter of fact, there are multifarious factors affecting the bubble formation, such as fluid viscosity, 5 surface tension, 6 density of gas and liquid phase, 7 orifice form, 8 and more. Nevertheless, those above-mentioned investigations prevailingly involved Newtonian fluids. Until now, most of the fluids in the industrial processes belong to non-Newtonian fluids. Compared to Newtonian fluids, very few understand the bubble formation characteristics in non-Newtonian fluids owing to its intricate rheological characteristics. Given that, some achievements in the bubble formation characteristics in non-Newtonian fluids have also been incorporated, which can be acquired from the work of Chhabra. 9 In addition, limited research conducted thus far indicates that they are especially dominant on the basis of experiments.
Numerical simulation has turned into a momentous method for exploring the microcosmic flow structure around the bubbles due to the limitations of experimental conditions, such as the high temperature and pressure as well as opacity systems. Nevertheless, the crux of simulating bubble motion characteristics is to track the gas−liquid interface accurately. So, a number of research studies by scholars have explored the way in which pervasive tracing methods of gas−liquid interface can simulate the bubble motion characteristics in fluids, such as the front-tracing (FT) method, the lattice-Boltzmann (LB) method, the level-set (LS) method, and the volume-of-fluid (VOF) method. Deen et al. carried out the numerical simulation of complex multiphase flow by the front-tracing method. 10 Yu et al. used the lattice-Boltzmann method to explore the interaction between bubbles. 11 Zhang et al. simulated the rising process of bubbles in shear-thinning fluids using the level-set method. 12 Cho and Son developed the motion of droplets and bubbles based on the level-set method. 13 Rabha and Buwa simulated the formation of a single bubble and multiple bubbles in shear-thinning fluids by virtue of the volume-of-fluid method. 14 Tomiyama and Sou applied the volume-of-fluid method to the numerical analysis of a single bubble. 15 The above investigation methods have their own characteristics, such as the front-tracking method is accurate and stable but prone to mass loss. The lattice-Boltzmann method is accurate but difficult to implement. The level-set method is simpler and more effortless to operate but the accuracy is lower. The volume-of-fluid method has a high calculation accuracy and can capture complex interface deformation, but the interface automatically merges and ruptures and numerical divergence may occur. Hence, some researchers attempted to couple two interface tracking methods to simulate the motion characteristics of bubbles in the liquid phase. 16−18 Recently, Fan et al. examined the formation volume and shape changes of bubbles in shearthinning non-Newtonian fluids through coupling the volumeof-fluid and level-set method (CLSVOF). 19 The yield-stress fluids are a kind of vitally crucial non-Newtonian fluids, which behave as solids when the applied stress is smaller than the yield stress and will flow when the applied stress is greater than the yield stress. A number of fluids, foams, slurries, and emulsions, belong to yield-stress fluids in industrial processes. In recent years, many researchers have studied the characteristics of bubble motion in yield-stress fluids. These studies have focused on the critical conditions for determining the onset of bubble motion 20,21 and the heat and mass transfer between gas and liquid, 22,23 while the study of bubble formation characteristics in yield-stress fluids is relatively rare.
In this study, a novel numerical simulation approach, CLSVOF, is given priority, which is able to make the simulation interface clear and the flow field changes well. However, there are no reports on the bubble formation characteristics of yield-stress fluids via coupling the level-set method and the volume-of-fluid method. In the first place, the numerical simulation results and experimental results are compared; subsequently, the influence of orifice diameter as well as rheological characteristics of fluids on the bubble formation time, volume, and aspect ratio is discussed. Last, but not least, the results are compared with Newtonian fluids and other types of non-Newtonian fluids, which is effortless to acquire a series of conclusions that play a dominant role in the bubble formation process and contribute to the reinforcement of chemical process facilities.

■ RESULTS AND DISCUSSION
Influence of Orifice Diameter. A legible graph describing the influence of the orifice diameter (examples 1, 2, 6) on the formation time and detachment volume of the bubble is shown in Figure 1a, where it can be seen that the bubble formation time decreases with the increase of the orifice diameter. It can be attributed to the forces acting on the bubble, such as . There, the adhesive force of the orifice impedes the bubble detachment from the orifice and, as a result, the bubble formation time decreases with the increase of the orifice diameter. It also can be observed from Figure 1a that the detachment volume of the bubble increases with the increase in the orifice diameter due to the setting of the gas inlet velocity as a constant. Thus, the gas flow increases with the increase of the orifice diameter under the same gas inlet velocity. In addition, larger bubble volume leads to the larger buoyancy force, which also decreases the formation time. This is consistent with the simulation results of Terasaka and Tsuge. 21 Figure 1b indicates that the orifice diameter d has a conspicuous bearing on the aspect ratio E during the bubble formation process. Here, the aspect ratio E is defined as h/w, where h is the vertical diameter of the bubble and w is the horizontal diameter of the bubble. It can be concluded from Figure 1b that the aspect ratio of the bubble increases with time under three different orifice diameters. In the initial stage of bubble formation, the aspect ratio increases at a high rate, which can be attributed to the pressure accumulation at the orifice outlet due to the successive gas flow. Once the pressure exceeds the critical value, the gas runs out of the orifice and the bubble formation process begins. On account of the inertial effect, the gas has a larger velocity in the vertical direction. There, the inertial effect reduces and disappears rapidly. Moreover, it can also be observed from Figure 1b that the

ACS Omega
http://pubs.acs.org/journal/acsodf Article growth rate of the aspect ratio increases when the orifice diameter increases. This is because an increase of the orifice diameter causes an increase of the gas flow and leads to a larger aspect ratio growth rate. Influence of Fluid Consistency Coefficient. The influence of the fluid consistency coefficient k (examples 3, 6, 9) on the formation time and detachment volume of bubbles in yield-stress fluids is demonstrated in Figure 2a. Obviously, an increase of the consistency coefficient k results in an increase of the bubble formation time and detachment volume. This increase is relevant to the difficulty for the bubbles to expand in the liquid phase and detach from the orifice, which increases with the increase of the liquid viscosity and then with the consistency coefficient. Homologous conclusions were maintained in the study of Fan et al. 19 The effect of consistency coefficient on the aspect ratio of the bubble in the formation process is shown in Figure 2b, which indicates that an increase of the consistency coefficient in the liquid phase causes a decrease of the aspect ratio growth rate. This decrease is related to the difficulty for the bubble to stretch in the perpendicular direction owing to the increase in the viscosity that has higher consistency coefficient in the liquid phase.
Influence of Fluid Flow Index. Figure 3a Figure 4, which indicates that the viscosity distribution (nondimensionalized by the consistency coefficient k) around the bubble at the end of the formation. The color bar reveals that the figures at the left indicate the degree of reducing viscosity around the bubbles. It can be observed that the viscosity around the bubbles decreases and the reduced region also decreases with the deceases of the flow index. The reason being that the detachment volume and formation time of the bubble increase with the increase of the flow index in Figure 3a. The influence of the flow index on the aspect ratio during the formation of the bubble is shown in Figure 3b. Distinctly, the aspect ratio growth rate decreases with the increase of the flow index in the yield-stress fluids. This is because the fluid with a smaller flow index has lower apparent viscosity, thereby reducing the viscosity resistance and increasing the growth rate of the aspect ratio of the bubble.
Influence of Yield Stress and Comparison with Newtonian Fluids. Figure 5a illustrates the bubble formation   This is ascribable to the fact that the shear-thinning fluid has the shear-thinning property and the surrounding fluid has the least resistance on the bubble, and so the formation time and detachment volume are lowest. In the case of the Newtonian fluid, the bubble is not affected by the yield stress, and so the formation time and detachment volume are shorter than that in yield-stress fluids. In yield-stress fluids with different yield stresses, the bubble formation time and detachment volume increase with the increase of yield stress. The influence of the fluids' rheological characteristics on the aspect ratio variation is seen in Figure 5b. Visibly, the growth rate of the aspect ratio in the shear-thinning fluid is largest, followed by the Newtonian fluid, and finally, in the yield-stress fluids. With regard to yield-stress fluids, the growth rate of the aspect ratio decreases with the increase of yield stress. This is because the shear-thinning fluid has the lowest viscosity around the bubbles and the resistance on bubble formation is less than other fluids, which leads to the fastest growth rate of the aspect ratio. Moreover, the growth rate of the aspect ratio in the Newtonian fluid is faster than that in yield-stress fluids due to the effects of the yield stress. For yield-stress fluids, the growth rate of the aspect ratio decreases with increasing yield stress.

■ CONCLUSIONS
A CLSVOF method is proposed and validated for the simulation of bubble formation in yield-stress fluids, shearthinning fluid, and Newtonian fluid. The simulation results are compared with the experimental results and it is found that the simulation results are in agreement with the experimental results. The influence of the orifice diameter, consistency coefficient, flow index, and yield stress on the bubble formation parameters (the bubble formation time, detachment volume, and aspect ratio) is investigated. The conclusions can be drawn as follows. The formation of the bubble can be divided into two stages: expansion dominated by surface tension and stretching dominated by buoyancy. For yield-stress fluids, the detachment volume increases with the increase of orifice diameter, consistency coefficient, flow index, and yield stress; the bubble formation time increases with the decrease of the orifice diameter but increases, respectively, with the consistency coefficient, flow index, and yield stress; and the growth rate of the bubble aspect ratio increases with the decrease of the consistency coefficient, flow index, and yield stress and decreases with the increase of the orifice diameter. In different kinds of fluids, the formation time and detachment volume of the bubble in the shear-thinning fluid are the lowest, followed by the Newtonian fluid, and finally, in yield-stress fluids. However, the growth rate of the aspect ratio in the shear-thinning fluid is the largest, followed by the Newtonian fluid, and finally, in yield-stress fluids. (1) where u⃗ is the velocity vector, p is the static pressure, F ⃗ represents the volumetric surface tension, and the density ρ(H) and the dynamic viscosity μ(H) of the mixed fluid are calculated with the following formula l g (4) where subscripts g and l represent the gas phase and the liquid phase, respectively. In the VOF model, the volumetric surface tension can be calculated by the continuum surface force model (CFS) as where σ is the surface tension coefficient of a fluid, κ(α) is the local surface curvature, and H(α) is the Heaviside function, defined as The level-set function α is defined as a function at a distance from the interface. When α > 0, it represents the liquid phase and, when α < 0, it represents the gas phase. The interface is given by α = 0. The local surface curvature κ can be obtained by Since the free surface is expressed via the CLSVOF method, the solution of the equation is replaced by the VOF function F and the LS equation H(α), given by eqs 8 and 9.
The local shear rate γ̇can be calculated as follows γ̇= · ⃗ ⃗ D D 2 ( : ) (10) The matrix form of the shear rate can be expressed as Because of the equation of continuity ∇·u⃗ = 0 and combined with type (12), D ⃗ 2 can be calculated as follows The shear rate γ̇is Numerical Simulation Method. The geometric model was established by SolidWorks, which was a rectangular domain with a width of 0.102 m and a height of 0.15 m. The orifice was located in the middle of the wide edge of the rectangular region and the diameters were 1.0, 1.5, and 2.0 mm, respectively. The air was set to the first phase and the density was 1.225 kg m −3 . Meanwhile, the liquid was set to different n, k, and τ 0 values through the Herschel−Bulkley (H−B) model in the second phase, respectively. The inlet was set as the inlet velocity boundary condition and it remained consistent. Subsequently, the top outlet was set as a pressure outlet boundary condition whose value was 0.1 MPa, and the remaining solid wall without a sliding boundary was provided. Finally, as the simulation belonged to the transient flow question, the calculation process changed at any time. Based on the pressure solver, the calculation model was CLSVOF multiphase laminar flow model. The pressure implicit with the splitting of the operator (PISO) was selected as the pressure− velocity correction method and the pressure interpolation format was selected. The convection term discrete equation in the equation was set as a first-order upwind equation. The convergence criterion of each scale residual was 1.0 × 10 −3 and the time step was 1.0 × 10 −3 s in the simulation.
Experimental Material. Carbopol 940 is a high molecular polymer that forms a stable, transparent, and thixotropic colloid when dispersed in water and neutralized. The yield stress of Carbopol gel can be tuned by controlling the polymer concentration and the pH value. When the bubble is moving at a low velocity, the surrounding fluids are at low shear rate conditions. In this case, the elastic behavior of the Carbopol gel

ACS Omega
http://pubs.acs.org/journal/acsodf Article can be neglected and such material can be approximated as an ideal yield-stress fluid.
Test of Grid Independence. The grid independence verification was carried out under the same conditions that included the orifice diameter d = 2 mm, fluid consistency coefficient k = 5 Pa·s n , flow index n = 0.7, and yield stress τ 0 = 3 Pa (example 6). The results of the grid independence verification are shown in Figure 6. Four mesh sizes of 0.5, 0.3, 0.2, and 0.1 mm were tested. The tests revealed that the results of 0.1 and 0.2 mm mesh sizes were better than others and the 0.1 mm mesh size needed too much simulated time. Thus, a mesh size of 0.2 mm was finally adopted.

■ TEST OF MODEL
The simulated results of example 6 in Table 1 were compared with the experimental results under slightly different conditions (d = 2 mm, k = 5.1 Pa·s n , n = 0.68, τ 0 = 2.9 Pa). The experimental device and the method were described in our previous work. 25 Notably, the simulation results were in good agreement with the experimental results in Figure 7, which proved the accuracy of the simulation method. It could be seen that the bubble formation was principally divided into two stages, expansion and stretching, which was similar to the simulation results. 26,27 The bubble underwent both horizontal and vertical expansion in the expansion stage between 0 and 1.2 s, which displayed a spherical shape due to the predominance of the surface tension. Nevertheless, longitudinal expansion of the bubble was much larger than the horizontal expansion in the stretching stage between 1.2 and 3.0 s. As the bubble was subjected to drag and buoyancy at this point, it caused the bubble to be stretched into a teardrop shape until it left the orifice.  ACS Omega http://pubs.acs.org/journal/acsodf Article