Numerical study of formation of a series of bubbles at a submerged orifice

Applied Mathematical Modelling 73 (2019) 668–694

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Numerical study of formation of a series of bubbles at a submerged orifice Xiyu Chen, Nan Chen∗, Antonio Delgado Institute of Fluid Mechanics, Friedrich-Alexander University, Erlangen, Germany

a r t i c l e

i n f o

Article history: Received 31 August 2018 Revised 26 March 2019 Accepted 4 April 2019 Available online 15 April 2019 Keywords: Bubble Contact angle Detachment time Aspect ratio Period doubling Cascade to chaos

a b s t r a c t Bubble formation from a submerged orifice is widely applied in bio-process and chemical reaction systems. In this study, the effects of different orifice diameters and contact angles in Period-I and Period-II regimes are studied systematically on a 2D axisymmetric domain. Simulation results are presented from the formation of the first bubble and explained by means of the surrounding fluid field, bubble interaction, and bubble aspect ratio. The orifice diameter is varied from 0.6 mm to 3 mm. The numerical results show that the detachment time of all bubbles remains constant (in time) for smaller orifice diameters (da ≤ 1.5 mm), while the detachment time of the first bubble is different from the rest of the bubbles for larger orifice diameters (da ≥ 2 mm), which is due to the different surrounding flow field. Contact angles from 60° to 165° are considered for the gas flow rates in the regime of bubble pairing, and it is observed that the bubble detachment time decreases when the contact angle increases, and it converges to a constant value when the contact angle is larger than 135°. In addition, the transition from period doubling to deterministic chaos (in which there is a variable number of bubbles within each period) is observed. A new scenario of inserting a submerged tube upward into the liquid is considered and compared to the previous cases. It is observed that when the tube is vertically inserted into the liquid, the bubble detachment time is even smaller because of higher influence from the surrounding liquid field, leading to a different phenomenon from the non-inserted tube cases. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Bubble formation from an orifice is of great interest because it is widely applied in various fields, such as chemistry, biology and industrial manufacture etc. The process of bubble formation from an orifice is fairly complicated due to the influence of the forces from the wall and the surrounding liquid: these forces cause varied detachment time and bubble shapes. Different bubble detachment situations and shapes influence the chemical reaction speed for gas as catalyst in chemical reaction and oxygen for biological cell growth processes [1–3], so it is necessary to dig this process in depth and understand the background physical reason. In order to understand the physical phenomena of bubble formation from a submerged orifice, considerable research efforts have been devoted by means of experimental study. Davidson [4] built up a theory of bubble formation based on ∗

Corresponding author. E-mail addresses: [email protected] (X. Chen), [email protected] (N. Chen), [email protected] (A. Delgado).

https://doi.org/10.1016/j.apm.2019.04.016 0307-904X/© 2019 Elsevier Inc. All rights reserved.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

669

Nomenclature Symbol da d

ρf ρg μf μg Q

σ θ

A vg v ∅ P

Meaning (Units) orifice diameter (m) bubble diameter (m) density of fluid (kg/m3 ) density of gas (kg/m3 ) viscosity of fluid, Pas viscosity of gas (Pas) gas flow rate (m3 /s) surface tension coefficient (N/m) contact angle (degree) area (m2 ) gas velocity (m/s) expansion velocity (m/s) volume fraction (-) pressure in the liquid phase (Pa)

the motion of a bubble in a viscous liquid. The theory assumed the bubble has spherical shape and found that the liquid viscosity has a major effect on bubble size, dominating the effects of gas flow rate, orifice diameters, liquid density and surface tension. Davidson and Amick [5], Snabre and Magnifotcham [6] and Albadawi et al. [7] studied different procedures of formation of bubble through a set of experiments and did some theoretical analysis based on the result of the obtained experimental data. They found that the bubble coalescence happens close to the orifice and observed that when the flow rate further increases, the large coalesced bubbles themselves undergo a second coalescence as they rise. Das et al. [8] introduced a new methodology for measuring the frequency of the formation and presented a simple analytical model for prediction of the bubble formation frequency. The bubble frequency increases when the flow rate increases. Zhu et al. [9] and Qu et al. [10] studied the dynamic bubbling behaviors in a bubble column with a micro-orifice of 54 μm in diameter by visual experiments and found a bifurcation structure of bubble radius and aspect ratio when increasing the distance from the orifice and presented a bubbling regime map to define the bubble formation behavior based on Bond number and We number. Xie et al. [11] and Zhang and Shoji [12] investigated the dynamic characterization of bubble formation and found three main types of bubble interaction: a) bubble interference, b) bubble collision, c) bubble coalescence. Other than the development of the theoretical models, a set of intensive experimental studies on the influences of different factors such as viscosity, surface tension, gas flow rates had already been carried out. McCann and Princes [13] subdivided the dynamic regime as single bubbling, pairing, double bubbling and double pairing, they observed that the individual bubble formation occurs at the low flow rate and can be described by the previous theory, when the flow rate increases, the interaction of bubble is stronger which cause the double bubbling and double pairing behavior. Irons and Guthrie [14] performed an experimental study to determine how several variables affect the size of gas bubbles formed at nozzles in liquid pig iron and reported that the bubble formation changes from uniform size at low flow rate to "doublets" and "pairs" as the chamber volumes vary. Mittoni et al. [15] presented a short analysis including the variation of chamber volume, injection nozzle diameter, liquid viscosity, and gas flow rate at room temperature. They found that the classic period doubling sequence leading to irregular behavior is observed only over certain system parameter ranges. Other than the experimental research, numerical study has become more and more popular due to cheaper and more flexible computer processing, and it has been widely used to extend the range of research into bubble formation beyond experimental studies. Oguz and Prosperetti [16] studied the bubble formation from an orifice by considering the Rayleigh– Plesset model and analyzed the results with respect to the mathematical formula. Lesage et al. [17,18] utilized the capillary equation to simulate the bubble formation from boiling liquid. Due to the limitation of the mathematical model, they mainly focus on the bubble formation under small gas flow rate in the direction more towards the theoretical model. Valencia et al. [19] used the VOF method to study the bubble formation process at the low inlet velocity for a fixed orifice diameter, and found the first two bubbles have different formation behavior due to the influence of the surrounding flow field. Gerlach et al. [1,20,21] and Buwa et al. [22] compared the accuracy of several numerical methods and used Coupled Level Set and Volume of Fluid method (CLSVOF) for bubble simulation and performed a systematic numerical study on the bubble formation from an orifice including the bubble interactions and wettability. They mainly focused on the phenomena of bubble pairing and doubling bubbles with low and moderate gas flow rates and the transition between single periodic and double periodic bubble formation and found that the strong vortex generated by larger leading bubble influences the growth of the following bubbles for the larger flow rates within the bubble paring and doubling bubble regime. Charkraborty et al. [23, 24] also performed a systematic numerical study of single bubble rising and bubble formation for a wide range of Bond number and Weber number by using the CLSVOF method. In addition, they also further investigated the influence of a co-flowing stream [25] and high density liquid with high Reynolds numbers [26]. Ma et al. [27] researched the bubble formation process with different viscosities and coalescence behaviors when bubble is rising and found that the bubble coalescence becomes more

670

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

frequent with increasing orifice gas flow rates because the two neighboring bubbles are closer near the orifice. Loubière et al. [28] investigated the forces acting on the bubbles during formation from a submerged orifice during different growth stages, and found that the drag force is the main factor for bubble detachment instead of buoyancy force. In addition, a considerable amount of research papers have investigated the effects of liquid viscosity, gas viscosity, surface tension coefficient, gravity, gas flow rates, orifice diameters including comparing the numerical results and experimental data for 2D and 3D simulations, for example, Dietrich et al. [29], Nahra et al. [30] and Simmons et al. [31], Zhang and Tan [32], Quan and Hua [33], Georgoulas et al. [34] and Zhang et al. [35]. The previous numerical studies focus mainly on the influence of liquid properties such as density, viscosity and surface tension as well as the boundary conditions such as the gas flow rates and orifice diameters. However, the influence of contact angles has not been paid enough attention. From the previous study on this field, Chen et al. [36] researched the movement of the macroscopic contact line during the bubble formation process of bubble. Vafaei et al. [37] studied the dynamics of bubble growth and bubble detachment from a substrate nozzle immersed in water experimentally by using the model Young–Laplace equation based on an axisymmetric bubble shape assumption for small gas flow rate. They have found that the contact angle of the bubble does not remain constant during the process of bubble formation, but it varies due to the movement of the contact line. Gnyloskurenko et al. [38, 39] presented experimental results of the surface phenomena effect on bubble formation from a single orifice (1 mm diameter) submerged in water with air blowing at an extremely small flow rate (2 cm3 /min) with different contact angle range 68°–110°, and observed that the bubble volume increases by more than half as the equilibrium contact angle increases. Chigarev and Chigareva [40] and Kandlikar and Steinke [41] investigated the influence of contact angle in terms of bubble shape during the boiling process. The bubble is higher when the contact angle increases, and the contact angle is larger for a very smooth surface. Gerlach et al. [20] has also simulated a set of cases with different contact angle within the regime of single bubble formation (Period-I) and found out that the bubble detachment time decreases when the contact angle becomes smaller. From the literature review above, the influence of contact angle needs to be investigated further, especially in the regime with bubble interaction and that is the original motivation of the current study. In this study, we extended previous work and performed our simulations with higher gas flow rate in the bubble pairing and doubling bubble (Period-II) regime. The influence of contact angle on periodicity of bubble formation is the main focus. Furthermore, a different scenario of a vertically inserted tube inside submerged liquid is simulated and compared to the non-inserted tube case. All simulation cases are performed with commercial software ANSYS Fluent and the results of the bubble formation from the very beginning are shown so that the complete process of bubble formation from an orifice can be studied. The structure of this paper is as follows: theory and numerical models are presented in Section 2, and Section 3 contains the results and discussion, including numerical validation, mesh study and results explanation. Conclusion is finally made in the last section. 2. Mathematical model and VOF method 2.1. Mathematical model In the present study, the conservation laws of mass and momentum for incompressible Newtonian liquid are used as follows:

∇ · u = 0

(1)

   ∂ (ρ u ) u  ) = −∇ p + ∇ μ ∇ u  + (∇ u  )T + Fs + ρ g + ∇ · (ρ u ∂t

(2)

For the surface tension force term Fs , the CSF (continuum surface force) model [42] is used. In the CSF model, the surface tension is assumed as constant along the surface and only the force normal to the interface is considered, so the force at the interface is described as a volume force in the momentum equation by the following equation [42]:

FS = σ

ρκ∇ ∅ f   0 . 5 ρg + ρ f

(3)

where κ is the vector curvature: κ = ∇ · nˆ , nˆ = n/|n|, n = ∇ · ∅f . The 2nd order derivative in the formula for curvature κ is not directly calculated, while the surface tension force Fs in the momentum Eq. (2) is re-formulated and calculated with the divergence of the outer normal vector nˆ . The calculation detail can be found in [42]. 2.2. VOF method There are several numerical methods to capture the advection of the interface, such as: Marker and Cell [43], Level Set [44], or Front Tracking method [45]. The VOF method constitutes a powerful and efficient technique for complex free surface

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

671

Fig. 1. The sharpness of the interface reconstruction; red region: bubble; blue region: liquid; transition region: interface; bold black line: actual interface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

problems. This method uses a color function ∅f to account for the fractional amount of fluid at a certain position (x, y, z) and time t. The volume fraction is described by the following equation:



∅f = 0 0 < ∅f < 1 ∅f = 1

gas phase interface liquid phase

(4)

If ∅f is zero, the content of the volume at that cell is gas, and if ∅f is unity, the content of the volume at that cell is liquid. Values of ∅f between one and zero represents a mixture of gas and liquid and generally take place at the gas-liquid interface. The color function ∅f follows the usual continuity relations of mass and momentum for a conserved quantity. The advantage of the VOF method compared to the other methods is that it guarantees mass conservation and gives accurate results for substantial topology changes in the interface without the requirement of re-initialization [46]. In the VOF model, the phases are described through solving the volume fraction continuity equation. As for the liquid phase, it is expressed by the following equation:

∂∅f  · ∇∅ f = 0 +u ∂t

(5)

The volume fraction of gas and liquid phases follow the equation below:

∅ f + ∅g = 1

(6)

The calculation of density and viscosity of the fluid mixture in a computational cell is shown by the following equations.

  ρ (x, t ) = ∅ f (x, t )ρ f + 1 − ∅ f (x, t ) ρg

(7)

  μ(x, t ) = ∅ f (x, t )μ f + 1 − ∅ f (x, t ) μg

(8)

2.3. Discretization scheme and solving algorithm In VOF method, Eqs. (1), (2) and (5) are solved based on the finite volume method. The implicit Euler method is used for time discretization, the pressure discretization method is pressure staggering option (PRESTO), and the momentum discretization is 2nd order upwind scheme [47]. The calculation and reconstruction of the volume fraction is based on PLIC [48] algorithm: the new field of volume fraction is calculated at new time-level and then the new interface is captured by calculating the new outer normal of each cell. The calculation details are shown in [48] and [49]. The advantage of the PLIC algorithm is the higher accuracy and sharpness. Fig. 1 demonstrates the sharpness after smoothing over the discontinuous volume fraction: the interface extends over only 3 cells. For time dependent multiphase simulations, the Pressure Implicit with Split Operator (PISO) algorithm [50] is used to calculate the pressure and velocity field with under relaxation factor to enhance stability, and then the volume fraction is updated with the new velocity field. The basic procedure of the algorithm is shown as follows: (1) Based on the initial conditions or previous time-level values of pressure field to calculate the intermediate velocity fields at new time level by using Eq. (2).

672

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 2. Computational domain; (a) domain of the original bubble formation case; (b) domain of the submerged vertically inserted tube case.

(2) To satisfy mass conservation, the pressure field is adjusted first by using the Eq. (1), then the intermediate velocity at new time level is also corrected through the new pressure field. The corrector step is performed one more time in the PISO algorithm by considering the neighboring corrector for higher accuracy. (3) Finally, the volume fraction is updated to give the new fluid configuration and the new interface is captured through PLIC algorithm. 2.4. Computational domain and boundary conditions Fig. 2 shows the domain of bubble formation in this study; a 2D axisymmetric domain is used to simulate the process of the bubble formation from orifice. Based on the previous research of Gerlach et al. [20, 21] and Vafaei et al. [37], this assumption is physically applicable. The inlet orifice diameter is da , the side wall is 20da away from the orifice in order to minimize wall effects; the height of the computational domain is 50da to mimic the influence of an outlet boundary and provide sufficient vertical distance to observe the bubble interaction behavior after formation. The bubble is injected from the gas inlet at different flow rates. The boundary condition at the inlet is fully developed parabolic flow profile, set by the user defined function (UDF), where the maximum velocity is two times of the average velocity umax = 2u0 , the outlet boundary is connected to the atmosphere, so the pressure outlet condition is imposed here. The side and bottom wall boundaries are set as non-slip walls. Since the bottom wall is in contact with the bubble in Fig. 2(a), the contact angle is imposed here. For the computational domain with the vertically inserted tube case shown in Fig. 2(b), the inserted portion with length 8.3da is set as non-slip wall. Initially, the computational domain is filled with liquid, and gas is injected at the bottom. A time step size of 5 × 10−6 s is used for the simulation cases. Results of different time step sizes and mesh independence study are shown in Section 3.2. It should also be mentioned that the definition of the contact angle in this paper is defined inside the bubble; it may be different from other papers such as Gnyloskurenko et al. [38, 39]. In this paper, the wetting condition is for contact angle θ ≤ 90°, while the non-wetting condition is for θ > 90°. 2.5. Dimensionless number Several dimensionless numbers are used to describe the flow characteristics, such as Reynolds number (Re), and Bond number (Bo). The Reynolds number is the ratio of inertial forces to viscous force shown in Eq. (9):

Re =

ρg u 0 d a ρg Q d a ρg Q d a 4 ρg Q = = = π da 2 μg μg A μ g π da μg 4

(9)

where the ρ g is the density of gas, Q and da are flow rate and the diameter of orifice, respectively; u0 is the average inlet velocity and μg is the dynamic viscosity of gas.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

673

Table 1 General properties for the validation of single bubble formation at higher gas flow rate. Property

Unit

Value

μf ρf μg ρg

Pa·s kg/m3 Pa·s kg/m3 cm3 /min m degrees N/m

1.003 × 10−3 998.2 1.7894 × 10−5 1.225 407.4 3 × 10−3 90 72.75 × 10−3

Q da

θ σ

Bond number (Bo) is the ratio of gravitational forces to surface tension forces, and the equation for Bo is given as follows:

  ρ f − ρg gda 2 Bo = σ

(10)

The symbols ρ f and σ are density of liquid and surface tension coefficient, respectively. 3. Results and discussion 3.1. Numerical validation Three different scenarios are used for validation of the model before using it to investigate other cases, namely: • • •

Single bubble formation with higher gas flow rate Formation of a series of bubbles with lower gas flow rate Single bubble formation with submerged tube vertically inserted into liquid

3.1.1. Validation of single bubble formation with higher gas flow rate This validation case refers to the experimental study of Chen et al. [36] and is based on air injected into water at different gas flow rates. Table 1 shows the property data of the gas-liquid system. At the beginning (t = 0), the gas injection is turned on for both experiments and numerical simulations, and a series of bubbles will then be formed. The detachment time of each formed bubble is calculated based on the time interval between the local starting point and the pinch-off of the neck from the orifice as the bubble detaches. Fig. 3(a) and (b) shows a snapshot comparison between experiment and simulation for detachment of single bubble. The experimental data shows that the bubble detachment time is 50.5 ms, and the simulation result shows that the detachment time of the 2nd bubble is 51 ms: it is close to the experimental result and the relative error is only around 1%. Compared to the previous numerical results shown by the green curve in Fig. 3(a), our numerical result (Fig. 3(a), blue curve) is closer to the experimental snapshot in terms of bubble height. However, since the simulation is based on a 2D axisymmetric domain while the bubble in the experiment is not perfectly symmetric, there exists some deviations regarding the profile of the bubble, but overall the comparison is acceptable. Fig. 3(c) shows the bubble profile at several snapshots during the formation process as predicted previously by Mirsandi et al. [51] by using the 3D Local Front Reconstruction Method (LFRM) with relatively lower grid resolution (mesh size of 0.25 mm), and comparison with the experimental photo and the current simulation results. In order to compare the bubble profile at specific time, the last snapshot is taken at 50.5 ms, which is consistent with the experimental data, and it is noticeable that the neck of the bubble is thicker than in Fig. 3(a) and (b). It can be seen that 2D axisymmetric simulations could still give comparable results to the experiments and the current numerical results are closer to the experiments in terms of the bubble height and a narrower neck due to higher grid resolution used compared to the work of Mirsandi et al. Fig. 4 shows the development of apparent contact angle of the formed bubble shown in Fig. 3 during the formation process and the present result is also compared to the experimental data as well as the numerical result from the previous study. It can be seen that the present result stays within the experimental range and is also comparable to the previous numerical result but with gradual decrease of the contact angle during the receding phase. It is probably another reason for the different bubble profiles shown in Fig. 3(a) and (b). 3.1.2. Validation of formation of a series of bubbles with lower gas flow rate In this case, formation of a series of bubbles is simulated, and the numerical results are compared to the experimental study from Zhang and Shoji [12] and the previous numerical study from Gerlach et al. [20] in terms of bubble shape and bubble detachment time. The general input properties are given by the Table 2, and in this case, liquid and gas are also water and air. All results are presented in Fig. 5. It can be seen that the numerical results shown in Fig. 5(b) are quantitatively comparable to the experimental snapshot in Fig. 5(a) in terms of bubble detachment time. Moreover, the bubble shape and bubble location shown in the red curves

674

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 3. Comparison for detachment of single bubble (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): (a) comparison of bubble profile at detachment time with experiments (green curve: simulation results from Chen et al. [36]; blue curve: present simulation profile of the 2nd formed bubble at 51 ms; dash-dot curve: indicator of axisymmetric line); (b) present simulation snapshot of the 2nd formed bubble at 51 ms; (c) comparison of several snapshots of bubble profile at specific times, (dash curve: the simulation results from Mirsandi et al. [51]) the blue curve in the last frame shows the prediction of the present simulation at 50.5 ms in order to be consistent with the experimental snapshot. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2 General properties in the validation case of formation of a series of bubbles. Property

Unit

Value

μf ρf μg ρg

Pa·s kg/m3 Pa·s kg/m3 cm3 /min m degrees N/m

1.003 × 10−3 998.2 1.7894 × 10−5 1.225 100 2 × 10−3 125 72.75 × 10−3

Q da

θ σ

in Fig. 5(a) have also good agreement to the experiments. Fig. 5(c) presents the results from Gerlach et al. [20] by using the CLSVOF method to perform the same simulation. It is noticeable that the numerical results from Gerlach start from the first bubble formed, but the snapshots from the experiments are taken from at least the 3rd bubble formed. Further study in Section 3.3.1 shows that the first bubble formed has slightly different detachment time from the rest of the bubbles with

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

675

Fig. 4. Comparison of the development of apparent contact angle (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): red dash curve: numerical result from Chen et al. [36]; blue curve: present result; black dot: experimental data from Chen et al. [36]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3 Properties of reference fluid in the insert tube validation case. Property

Unit

Value

μf ρf μg ρg

Pa·s kg/m3 Pa·s kg/m3 cm3 /min m N/m

43 × 10−3 1250 1.7894 × 10−5 1.225 30 0.6 × 10−3 60 × 10−3

Q da

σ

orifice diameter da = 2 mm, so it is not precise enough to perform the comparison from the beginning of the numerical simulation. This is probably the reason that our numerical results are closer to the experiment in this aspect.

3.1.3. Validation of bubble formation with vertically inserted tube submerged into liquid The previous validation cases are based on the orifice located at the bottom. When the gas is injected through a vertically inserted submerged tube, the phenomena of bubble formation could be different due to different surrounding flow fields. For the validation of this case, we use an experiment from Snabre and Magnifotcham [6] with the air injected into a Newtonian water-glycerol solution. The general properties are represented in Table 3. Fig. 6 shows the snapshot from the simulation when the first bubble is just detached from the inserted tube. The experimental data shows that the detachment bubble volume is 20.97 mm3 . In the simulation case, the area of the simulated bubble is calculated by using numerical integration at the detachment time 42 ms, with the volume fraction integrated over cells with volume fraction larger than 0.5. Hence, the equivalent diameter can be calculated from the area obtained. And the corresponding volume is further obtained by calculating from the equivalent diameter. The final numerical result of the detachment volume is 20.19 mm3 with the relative error around 3.8%, so it satisfies the quantitative agreement very well.

676

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 5. Formation process comparison for a series of bubbles: (a) experimental results (Zhang et al. [12]) and numerical profile (red curve); (b) present simulation results (Q = 100 cm3 /min, Re = 72.63, da = 2 mm); (c) previous numerical results (Gerlach et al. [20]). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3.2. Time step size and mesh independence study In order to reduce the computational cost and increase the efficiency of the simulation, it is necessary to perform a mesh independence study. The case used here refers to the experimental data [36] described in Section 3.1.1. Six mesh sizes 0.05 mm (500 × 750), 0.06 mm (416 × 625), 0.075 mm (333 × 500), 0.1 mm (250 × 375), 0.15 mm (166 × 250) and 0.25 mm (100 × 150) are used in the mesh independence study. Table 4 shows the quantitative comparison of the numerical results with different mesh sizes and Fig. 7 shows the snapshots at detachment time with different mesh sizes. For the coarsest mesh, the bubble detachment time and bubble shape are not comparable to the experimental results due to the lack of resolution. The numerical results become better when the mesh is refined, and good agreement can be obtained when mesh size is 0.06 mm, with the relative error of detachment time being 1.0%. When the mesh is further refined, there is no

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

677

Fig. 6. Numerical validation for the case of bubble formation from vertically inserted submerged tube, snapshot is taken at 42.0 ms when the first bubble is detached (μf = 43 mPa·s, Q = 30 cm3 /min, Re = 72.63, da = 0.6 mm). Table 4 Bubble detachment time with different mesh size, orifice diameter da = 3 mm. Mesh size

0.05 mm

Numerical value(ms) Mesh size Numerical value(ms) Reference value(ms)

Value 50.9 0.1 mm Value 52.0

0.06 mm Error 0.7% Error 3.0%

Value 51.0 0.15 mm Value 52.2 50.5

0.075 mm Error 1.0% Error 3.4%

Value 51.8 0.25 mm Value 48.6

Error 2.5% Error 3.8%

Fig. 7. Simulation snapshots at detachment time for different mesh sizes.

significant change in terms of bubble detachment time and shape. Hence, the numerical results are based on the mesh size 0.06 mm and the total number of mesh cells is 260,0 0 0, about 44% less than if the mesh size was taken to be 0.05 mm. Fig. 8 shows the comparison of the bubble profile at detachment time with different time step t sizes with the mesh size 0.06 mm. It can be seen that the difference in bubble profile is negligible when t ≤ 5 × 10−6 s since the blue and green curves are almost overlapped. So finally the time step 5 × 10−6 s is selected for further simulations.

678

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 8. Simulation results at detachment time with different time steps: (a–c) the bubble snapshots at detachment time with different time step sizes; (d) the comparison of the bubble profile with different time step sizes. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 9. Detachment time and volume for different orifice diameter from 0.6 mm to 3 mm (average of inlet gas velocity vg = 0.53 m/s is used for all cases; the volumetric flow rate varies from 9 cm3 /min to 225 cm3 /min): (a) detachment time for a series of bubbles; (b) bubble detachment volume for a series of bubbles.

3.3. Numerical results 3.3.1. Effects of different orifice diameters Extending the validation case in Section 3.1.2, the bubble formation process from different orifice diameters (from 0.6 mm to 3 mm), keeping the inlet gas velocity constant, is studied. The change of the surrounding flow field and development of aspect ratio with different orifice diameters are also investigated. The material properties are the same as the data given in Table 2 except for specifically mentioned da . According to the theoretical study of Oguz and Prosperetti [16], the ratio to the critical gas flow rate is within the range [0.19, 1.22] and the ratio to the Fritz volume is within the range [1.45, 2.35]. Simulation results from the first bubble are shown together with the following bubbles until repeatable behavior of bubble formation is accomplished. Fig. 9(a) and 9(b) show respectively the detachment time and detachment bubble volume for different orifice diameters from the 1st bubble until the 5th bubble. From Fig. 9(a), it can be observed that the detachment time decreases when the orifice diameter increases. This is because gas momentum increases quadratically compared to the linear increase in surface tension force [6], so it is easier for the bubble to detach from the orifice when the orifice diameter increases. Moreover, starting from the smallest orifice diameter da = 0.6 mm, at fixed inlet velocity, the detachment time maintains a constant value for the formation of a series of bubbles up to da = 1.5 mm. A difference in detachment time between the first bubble and the following bubbles arises with larger orifice diameter (da = 2 mm and 3 mm). An obvious difference between the first bubble and subsequent bubbles can be observed in the case of da = 3 mm in Fig. 9(b) when comparing the bubble volume. Similar phenomena are also observed in the experimental study from [9, 10]. In order to further understand these phenomena, the surrounding flow field and aspect ratio development during bubble formation process are presented in the following figures. Fig. 10 shows the flow field of the case with orifice diameter da = 3 mm. It is clearly shown that the flow field surrounding the first bubble (Fig. 10(a)) is different from the field around the second bubble (Fig. 10(c)) since there is no vortex generated from a previous bubble for the first one. Due to this effect, the detachment time of the second bubble is 10% larger than the first one. It can also be seen that the flow field during formation process of the 3rd bubbles (Fig. 10(e)) is similar to the second one (Fig. 10(c)), so the detachment time of 3rd bubble is the same as the 2nd one. This phenomenon continues for the rest of the bubbles, so that the detachment time afterwards maintains a constant value. Fig. 11 presents the aspect ratio (height/width) development during the bubble formation process. It can be observed that the process of bubble formation consists of three phases: in the first phase, the aspect ratio increases rapidly because of the domination of gas momentum force and it pushes the bubble growth preferentially in the vertical direction; in the second phase, the bubble expands in both horizontal and vertical directions and the increment of aspect ratio becomes smaller, (this is probably due to the influence of the previous bubble, which delays the bubble being further growth in the vertical direction). In the third phase, the aspect ratio starts to grow again until detachment. Bari and Robinson [52] also proved

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

679

Fig. 10. Flow field of bubble formation for the orifice diameter 3 mm (Q = 225 cm3 /min, Re = 108.84): a, b) first bubble; c, d) second bubble; e, f) third bubble.

Fig. 11. Aspect ratio of bubble formation for orifice diameter da = 3 mm (Q = 225 cm3 /min, Re = 108.84).

these three-phase phenomena by experimental observation. From Figs. 10 and 11, the formation process of the second and the third bubbles is almost the same due to the similar surrounding flow field; considering the formation of the first bubble, since there is no influence of being compressed, the aspect ratio of the first bubble has shorter time in Phase II compared to the next two bubbles, leading to a smaller detachment time. Fig. 12 shows the flow field of orifice diameter da =1 mm. The small vortex generated by the first bubble is still observed during the formation process of the following bubble (Fig. 12(b)). However, since the bubble volume is small and the gap between the neighboring bubbles becomes larger and larger (Fig. 12(c) and (d)), the influence of the wake from the first bubble becomes smaller, so that the change of detachment time of second bubble is only 2.4% different from the first one. A similar phenomenon is also shown during the formation of the third bubble (Fig. 12(e) and (f)), resulting in the same detachment time as well as detachment volume and this phenomenon continues for the rest of the bubbles. Fig. 13 shows a comparison between aspect ratio history for the first and second bubbles, and it is fairly understandable that both bubbles have similar curves because of the similar flow field. Based on the discussion above, one can come to the conclusion that for small gas flow rate in the Period-I regime, the bubble detachment time remains constant except for the first one formed from larger orifices (da ≥ 2 mm). This is because when the orifice is larger, the detachment time becomes smaller while the bubble volume becomes larger. Both factors lead to a narrower gap between the neighboring bubbles as well as a larger surrounding vortex (Figs. 10(c) and 12(c)). As a result,

680

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 12. Flow field of bubble formation for the orifice diameter da = 1 mm (Q = 25 cm3 /min, Re = 36.28): a, b) first bubble; c, d) second bubble; e, f) third bubble.

Fig. 13. Aspect ratio of bubble formation for orifice diameter da = 1 mm (Q = 25 cm3 /min, Re = 36.28).

when the influence from the first bubble cannot be neglected (da ≥ 2 mm), the formation of the second bubble is changed with the extra force from the wake of the first one. This difference only happens between the first and subsequent bubbles; the detachment time remains the same from the second one because of the same surrounding flow field and development of aspect ratio.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

681

Fig. 14. Detachment time for different contact angles (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm).

3.3.2. Effects of different contact angles There are several factors that influence the bubble formation, and lots of research has been performed investigating the effects of liquid properties such as viscosity, density, surface tension coefficient and boundary conditions such as orifice diameter and gas flow rate. However, the effects of contact angle have been given less attention yet it is an important factor for bubble detachment. Gnyloskurenko et al. [38, 39] presented experimental results of bubble formation from a single orifice at an extremely small flow rate (2 cm3 /min) with contact angle range 70°–112° (with contact angle as defined in this paper). They found that the bubble volume increases by more than half as the equilibrium contact angle decreases. Gerlach et al. [20] studied the influence of contact angle for a single bubble at low flow rate (Q ≤ 200 cm3 /min). Buwa et al. [22] researched the formation of bubbles with higher gas flow rate numerically with single contact angle. In the current study, we focus on the formation of a series of bubbles at higher flow rate (v = 0.96 m/s, Q = 407.4 cm3 /min) for larger orifice diameter (da = 3 mm) within the contact angle range from 60° to 165°. According to the theoretical study of Oguz and Prosperetti [16], the ratio to the critical gas flow rate is 2.21 and the ratio to the Fritz volume is within the range [3.30, 7.89]. For this gas flow rate, the phenomenon of bubble formation is fairly complicated due to bubble interaction, and the influence of contact angle also becomes more significant. Based on the validation case shown in Section 3.1.1, a series of simulations with different contact angles are performed and the results of detachment time are shown in Fig. 14. It can be clearly seen that the detachment time increases when the contact angle decreases, which was also found by Gnyloskurenko et al. [38]. This is because when the surface changes from wetting to non-wetting, surface tension force between bubble and bottom wall decreases and the bubble can detach from the orifice easier. Moreover, the detachment time for the first bubble is also smaller than for the others for the reason explained in Section 3.3.1. It can also be observed that the detachment time of each bubble is almost independent of contact angle when the contact angle is larger than 135°. Due to higher gas flow rate in the current study, different contact angles show different periodicities: the periodicity become more complex as the contact angle increases, with the period of the repeating sequence involving a greater number of bubble formation events. Different periodicities can be divided into three contact angle ranges: for smaller contact angles from 60° to 90°, one period contains two bubbles; for medium contact angle 105° and 120°, one period contains four bubbles; for large contact angle θ > 135°, the period continues increasing and the periodicity is not found. The reason why different contact angles show different periodicities is also due to the surrounding flow field and bubble coalescence, which will be discussed in detail for each range in the following section. Fig. 15 shows the flow field for the contact angle 90°. Similar to the previous case, the generation of the second bubble is influenced by the vortex from the first one. However, for higher gas flow rate, the gap between bubbles is variable. As is shown in Fig. 15, the gap between the 1st and the 2nd bubbles (Fig. 15(c)) is slightly larger than the gap between the 2nd and 3rd ones (Fig. 15(e)). Due to this phenomenon, the detachment time is different. The same process repeats again as the 4th bubble has the same surrounding flow field as the 2nd one (Fig. 15(c), (d) and (g), (h)), and then the flow field around the 3rd and 5th bubbles is also the same (Fig. 15(e), (f) and (i), (j)). As a result, for contact angle θ = 90°, one period contains two bubbles with different detachment time. Fig. 16 shows the development of aspect ratio of the case θ = 90°. It can be observed that, due to the gap and the vortex from the previous bubble, the detachment time is larger with a relatively smaller aspect ratio. As a result, the development of aspect ratio contains three phases for the bubbles with longer detachment time (3rd and 5th bubbles); while for the shorter-detachment-time bubbles (2nd and 4th bubbles), the development of aspect ratio contains only two phases: without the 2nd one in which the growth of aspect ratio is fairly small. Comparing the bubble shape shown in Fig. 17(b) between the 2nd and 3rd bubbles, the height of the 2nd bubble is slightly larger than the 3rd one while the expansion in the horizontal

682

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 15. Flow field of bubble formation for the contact angle 90° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): a, b) first bubble; c, d) second bubble; e, f) third bubble; g, h) fourth bubble; i, j) fifth bubble.

Fig. 16. Aspect ratio of bubble formation for contact angle θ = 90° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm).

direction is similar. This phenomenon indicates that due to smaller gap (Fig. 15(e)), the bubble gets compressed and the growth in the vertical direction is delayed. From the discussion and explanation above, due to the flow field and development of aspect ratio, for small contact angle θ ≤ 90° with higher gas flow rate, a bubble pairing phenomenon is observed with the odd numbered bubbles (3rd, 5th) chasing the even numbered bubbles (2nd, 4th) and the similar pairing phenomenon repeats for every two bubbles which shows a two-bubble period. Fig. 18 shows the bubble detachment time for contact angles θ = 105°and 120°. Although the bubble detachment time varies between odd and even bubbles, it has two different longer (3rd and 5th bubbles) and two different shorter (4th and 6th bubbles) detachment times, so one period contains four bubbles for contact angle θ = 105° and 120°. Fig. 19 shows the flow field for θ = 120°: it can be seen that for bubbles with longer detachment time (3rd, 5th, 7th and 9th bubbles), the surrounding flow field is different since the gap between the 4th and 5th bubbles (Fig. 19(e)) is larger than the gap between the 2nd and 3rd ones (Fig. 19(a)), so the bubble detachment time will also be different. A similar phenomenon can also be observed by comparing the gap and flow field of the 8th and 9th bubbles (Fig. 19(m)) to the ones between the 6th and 7th

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

683

Fig. 17. Snapshot of the bubble shape at labeled time (a, b, c and d correspond to the time in Fig. 16).

Fig. 18. Periodicity of detachment time for contact angle 105° and 120° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm).

bubbles (Fig. 19(i)). On the other hand, the gap and flow field are similar between the formation of the 3rd and 7th bubbles (Fig. 19(a) (b) and (i) (j)), leading to the same detachment time, so is the formation of the 5th and 9th bubbles (Fig. 19(e) (f) and (m) (n)). Similar phenomena can also be observed for the bubbles with shorter detachment time (4th, 6th, 8th and 10th bubbles), when comparing Fig. 19(d) and (l) to Fig. 19(h) and (p), respectively. So finally, a four-bubble period is found with contact angle θ = 120°. Further explanation is given by considering the aspect ratio curve shown in Fig. 20: the aspect ratio of the 5th bubble is lower than the 3rd bubble in Phase-II, so the detachment time of the 5th bubble is larger; on the other hand, the aspect ratio curves for 4th and 6th bubble are similar, so the relative difference of detachment time is only 2.1%. These phenomena could also be reviewed by comparing the flow field shown in Fig. 19(c) (d) to Fig. 19(g) (h), it can be observed that the difference in terms of previous bubble shape and flow field are similar. Since the detachment time is further decreased with larger contact angles, the development of aspect ratio contains only two phases. Similar behavior is shown for both longer and shorter detachment time bubbles in the first phase; in the second phase, for longer detachment time bubbles the aspect ratio is smaller while for shorter detachment time bubbles the aspect ratio is larger. This phenomenon is further verified in the snapshots of bubble shape in Fig. 21. From the comparison of Fig. 21(c) and (d) between the 4th and 5th bubble, it is noticeable that the 5th bubble has larger size in the horizontal direction, leading to a smaller aspect ratio. Fig. 22 shows the bubble detachment time for contact angle θ = 150°. The bubble detachment time also shows similar behavior as the case θ = 120°: with larger detachment time for odd numbered bubbles and smaller detachment time for even numbered bubbles. However, this odd-even fluctuation is irregular and the periodicity in this case is not clearly shown (at least at this stage of the simulation). Fig. 23 shows the flow field for contact angle θ = 150°: the flow field of the even numbered bubbles with smaller detachment time is shown in the first two columns, while the flow field of the odd numbered bubbles with larger detachment time is shown in the third and fourth column, respectively. Snapshots of Fig. 23(b), (f) and (j) indicates that all the previous bubbles shown in this group have different shapes and the flow fields generated from these leading bubbles are

684

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 19. Flow field of bubble formation for the contact angle 120° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): from (a) to (h) is one period including the 3rd to 6th bubble (each bubble contains two snapshots); from (i) to (p) is the other period including the 7th to 10th bubble (each bubble formation contains two snapshots).

Fig. 20. Aspect ratio of bubble formation for the contact angle 120° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm).

slightly different, too. As a result, the formation of the corresponding following bubbles is influenced by different wakes from the previous ones, and the detachment time is finally different from each other. Similar phenomena can also be found in Fig. 23(c), (g), (k) and (d), (h), (l) groups with different leading bubble shapes and flow fields. From Fig. 24 showing the curves of aspect ratio development, the aspect ratio development of these bubbles is also slightly different from each other. Although all curves follow the same tendency: the higher the aspect ratio, the smaller the bubble detachment time; there exist special curves with individual characteristics, especially for the 5th and 6th bubbles.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

685

Fig. 21. Snapshot of the bubble shape at labeled time (a, b, c and d correspond to the time in Fig. 20).

Fig. 22. Detachment time with contact angle θ = 150° (Q = 407.4 cm3 /min, Re = 197.28, da =3 mm).

For the 5th bubble, it has smaller aspect ratio in Phase II; for the 6th bubble, it has larger aspect ratio in Phase I, this detail is further explained in Fig. 25. From the comparison of the snapshot at instants ‘a’ and ‘b’, it can be observed that the 6th bubble has slightly larger height than the 8th bubble due to the influence from the 5th bubble, so the detachment time of the 6th bubble is smaller. On the other hand, the comparison of the snapshot at instants ‘c’ and ‘d’ between the 5th and 7th bubble shows that the 5th bubble is compressed more in the vertical direction, leading to a smaller aspect ratio in Phase-II, so the corresponding detachment time of the 5th bubble is larger. Since the flow field and the aspect ratio of all these bubbles are slightly different, the detachment time is not repeatable up to the 15th bubble, so no clear periodicity is observed in this case. Mittoni et al. [15] has also demonstrated this transition to chaos for bubbling experimentally. And here we explore further from our numerical results by studying the velocity fields, aspect ratio as well as bubble interaction. 3.3.3. Bubble formation from inserted vertical submerged tube When the orifice wall is extended into the liquid by means of a vertical tube, the influence of the flow field becomes stronger from the bottom side of the bubble, while the influence from the bottom wall becomes weaker, so new detachment behavior could be observed, and Snabrea and Magnifotcham [6] has done some experimental research for such cases. Based on the validation case from Section 3.1.3, we perform our numerical analysis and focus on the different flow physics that occurs for bubble formation at the vertical tube. According to the theoretical study of Oguz and Prosperetti [16], the ratio to the critical gas flow rate is 1.22 and 2.21 for lower and higher gas flow rate, respectively; and the ratio to the Fritz volume is within the range [2.22, 3.73].

686

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 23. Flow field of bubble formation for the contact angle θ = 150° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): a, b) fourth bubble; c, d) fifth bubble; e, f) sixth bubble; g, h) seventh bubble; i, j) eighth bubble; k, l) ninth bubble.

Fig. 24. Aspect ratio of bubble formation for the contact angle 150° (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm).

Fig. 26 shows the detachment time for the inserted tube case. Due to stronger surrounding flow field from the bottom side of the bubble and weaker surface tension force between bubble and tube wall, even for the lower gas flow rate case, there is a small perturbation in detachment time with a four-bubble periodicity, compared to the constant detachment time for bubbles without tube insertion. When the inlet flow rate increases further, the phenomenon of a period of four bubbles is preserved, however, the fluctuation of detachment time between odd and even numbered bubbles becomes larger because of the larger influence of the flow field. Fig. 27 shows the detachment time comparison between inserted tube and non-inserted tube for lower gas flow rate (Q = 225 cm3 /min, Re = 108.85). Generally speaking, when the tube is vertically inserted into the liquid, the detachment time is smaller than in the non-inserted tube case. In addition, it shows small fluctuation because of the extra force generated from the flow field underneath the bubble.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

687

Fig. 25. Snapshot of the bubble shape at labeled time (a, b, c and d correspond to the time in Fig. 24).

Fig. 26. Detachment time for inserted tube with different flow rates (Q = 225 cm3 /min, Re = 108.84, da = 3 mm) and (Q = 407.4 cm3 /min, Re = 197.16, da = 3 mm).

Fig. 28 shows the flow field for the inserted tube case with lower flow rate (Q = 225 cm3 /min, Re = 108.85). It can be seen that the gap of the longer detachment time bubbles (the 5th, 7th and 9th bubbles in Fig. 28(d), (h), (l)) is slightly larger than the gap of the shorter detachment time bubbles (the 4th, 6th and 8th bubbles in Fig. 28(b), (f), (j)), so small deviations of detachment time among different bubbles are shown. Although this difference is fairly small (around 3%), the fluctuation of detachment time is still clearly observed in Fig. 27 especially when compared to the constant bubble detachment time in the non-inserted tube case. Due to this phenomenon, the bubble detachment time is considered to be variable and a period of four bubbles is found. Fig. 29 shows the development of the aspect ratio in the inserted tube case. Based on the conclusion from the previous discussion, since both bubbles generated from the inserted tube have slightly larger aspect ratio during the formation process, their detachment time will be smaller than in the non-inserted tube case. In addition,

688

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 27. Comparison of detachment time (blue curve: submerged vertical tube is inserted into the liquid; orange curve: non-inserted tube; Q = 225 cm3 / min, Re = 108.84, da = 3 mm). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

because the bubble detachment time becomes smaller, the aspect ratio development contains only two phases with different growing speeds instead of three phases. Fig. 30 shows the comparison of detachment time between inserted tube and non-inserted tube (contact angle θ = 105°, 120° and 150°) cases for higher gas flow rate (Q = 407.4 cm3 /min, Re = 197.28). The periodicity of the inserted tube case contains four bubbles, including two different larger detachment time and two different smaller detachment time, which is similar to the contact angle case θ = 120°. However, compared to the smallest bubble detachment time with non-inserted tube case shown in the curve θ = 150°, the bubble detachment time of the inserted tube case is even smaller because of the influence of flow field from the bottom of the bubble. Fig. 31 shows the flow field of the inserted tube case with higher gas flow rate (Q = 407.4 cm3 /min, Re = 197.28). It can be seen that the flow field and the leading bubble shape are quite different from the non-inserted case (Fig. 19). For the bubbles with smaller detachment time (4th and 6th bubbles), the gap between neighboring bubbles is even smaller (Fig. 31(d) and (h) compared to Fig. 19(d) and (h)), so the influence from the wake of previous bubble is even stronger and the bubble is formed with a higher aspect ratio (shown in Fig. 32). For the bubbles with larger detachment time (the 5th and 7th bubbles), the bubble gap is fairly large so the influence from the flow field is not as strong as the even-numbered bubbles; however, because there is no wall contact, the detachment time is still smaller than the corresponding bubbles in the non-inserted tube case. Fig. 32 presents the aspect ratio of the bubbles. Similar to the phenomenon shown in Fig. 20, it is quite obvious to observe the aspect ratio difference between longer and short detachment time, just with a bigger difference. From the discussion above, one can come to the conclusion that when the tube is vertically inserted into the liquid, the influence of the flow field especially from the bottom side will be changed, so is the development of bubble aspect ratio. These two factors lead to different bubble detachment time from the non-inserted tube case. Generally speaking, due to stronger influence from the flow field from the bottom side of the bubble and larger aspect ratio during the formation process, the bubble detachment time is smaller than in the non-inserted tube case, and the fluctuation of detachment time between odd and even numbered bubble becomes larger.

4. Conclusion and outlook In this study, formation of a series of bubbles is investigated and the results are shown from the first formed bubble. The influence of orifice diameter, contact angle, gas flow rate is studied, which covers the range of [0.19, 2.21] regarding the ratio between the imposed gas flow rate to the critical gas flow rate and the range of [1.45, 7.89] regarding the ratio between the bubble volume to the Fritz volume, together with further comparison to the cases with vertical inserted submerged tube. The following conclusions can be drawn: When the inlet velocity is fixed, the detachment time increases when the orifice diameter decreases because gas momentum decreases quadratically compared to the linearly-decreasing surface tension force. When the orifice diameter is small (da ≤ 1.5 mm), the bubble detachment time remains constant; for larger orifice diameters (da ≥ 2 mm), the formation process of the first bubble is different from the following ones with smaller detachment time. For the bubbles with larger detachment time, the development of the aspect ratio (height/length) contains three phases: growth-stable-growth; while for the bubbles with smaller detachment time, it contains two phases with the 2nd phase being shorten. Generally speaking, during the development of aspect ratio, larger aspect ratio leads to smaller detachment time.

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

689

Fig. 28. Flow field of bubble formation for the inserted tube(Q = 225 cm3 /min, Re = 108.84, da = 3 mm): a, b) fourth bubble; c, d) fifth bubble; e, f) sixth bubble; g, h) seventh bubble; i, j) eighth bubble; k, l) ninth bubble.

When the gas flow rate and the orifice diameter are fixed, as the contact angle is increased, the bubble detachment time decreases when contact angle θ < 135°, but the detachment time becomes independent of contact angle when the large contact angle θ ≥ 135°. In addition, different contact angles lead to different periodicities of bubble formation: • • •

60° ≤ θ ≤ 90°, two-bubble period 105° ≤ θ ≤ 120°, four-bubble period θ ≥ 135°, no obvious period observed

With the increment of the apparent contact angle, the transition of the period doubling to the deterministic chaos with different number of bubbles within one period is observed in the bubbling dynamics. This phenomena is known as cascade to chaos. Different periodicities are mainly due to the different flow fields and aspect ratios during the process of bubble formation, both of which lead to different detachment time of each bubble. When the tube is vertically inserted into the liquid, because of the stronger influence of the flow field from the bottom side of the bubble, the detachment time is generally smaller, and the fluctuation of detachment time between odd and even numbered bubbles starts to appear for the smaller gas flow rate in Period-I and becomes larger for the higher gas flow

690

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 29. Comparison of aspect ratio of bubble formation for the inserted tube case (6th and 7th bubble) and non-inserted tube case (3rd bubble) (Q = 225 cm3 /min, Re = 108.84, da = 3 mm).

Fig. 30. Comparison of detachment time between non-inserted tube (contact angle θ = 105°, 120° and 150°) and the inserted tube cases (Q = 407.4 cm3 / min, Re = 197.28, da = 3 mm).

rate in Period-II. This can also be explained by the larger differences observed the surrounding flow field as well as the development of aspect ratio. Due to the limits of computational capacity, a 2D axisymmetric model is applied here and the numerical simulation only focuses on the gas flow rate within Period-I and II regime. This 2D domain is not able to capture the whole trajectory of rising bubbles, especially the wobbling bubbles. Further research could focus on further increment of the gas flow rate in Period-III and IV regime, and the dynamic contact angle model could also be implemented and applied for numerical study of bubble formation. Moreover, a 3D model with enough grid resolution can be adopted for real physical representation when the computational capacity fulfils the requirement. Appendix: Ratio to the Fritz volume and critical gas flow rate RF is calculated based on the force balance between the buoyancy force and surface tension force, derived by Fritz [53]



RF = VF =

3σ da 4ρ f g

 13

4 π R3F 3

Oguz and Prosperetti [16] stated that there is a critical gas flow rate, when the gas flow rate Q < Qcr , the detachment 6

volume V ≈ VF , while for the bubbles growth at a faster rate than Qcr , the detachment volume is proportional to Q 5

16 16  σ d  56 a Qcr = π 2ρ f 3g2

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

691

Fig. 31. Flow field of bubble formation for the inserted tube (Q = 407.4 cm3 /min, Re = 197.28, da = 3 mm): from (a) to (h) is one period including the 3rd to 6th bubble (each bubble contains two snapshots); from (i) to (p) is the other period including the 7th to 10th bubble (each bubble formation contains two snapshots).

692

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

Fig. 32. Comparison of aspect ratio of bubble formation for the inserted tube (Q = 407.4 cm3 /min, Re = 197.28, da =3 mm).

Validation Case 1 Validation Case 2 Validation Case 3

R (mm)

RF (mm)

R/R F

V/VF

Q (cm3 /min)

Qcr (cm3 /min)

Q/Qcr

4.34 2.64 1.71

2.56 2.23 1.30

1.70 1.18 1.32

4.90 1.64 2.28

407.4 100 30

184.01 131.25 33.98

2.21 0.76 0.88

Case1: Effect of different orifice diameters Orifice diameter(mm)

0.6

1

1.5

2

3

R(mm) RF (mm) R/RF V/VF Qcr (cm3 /min) Q/Qcr

1.85 1.50 1.24 1.89 48.12 0.19

2.04 1.77 1.15 1.52 73.66 0.34

2.30 2.03 1.13 1.45 103.27 0.54

2.70 2.23 1.21 1.77 131.25 0.76

3.40 2.56 1.33 2.35 184.01 1.22

Case2: Effect of different contact angles Contact angle(°) Rmin (mm) Rmax (mm) RF (mm) (R/RF )min (R/RF )max (V/VF )min (V/VF )max Contact angle(°) Rmin (mm) Rmax (mm) RF (mm) (R/RF )min (R/RF )max (V/VF )min (V/VF )max Q Qcr (cm3 /min) Q/Qcr

60 5.01 5.09 1.96 1.99 7.53 7.89 120 3.91 4.13 1.53 1.62 3.59 4.23

75 90 4.67 4.32 4.81 4.54 2.56 1.83 1.69 1.88 1.78 6.09 4.84 6.65 5.60 135 150 3.84 3.81 4.07 4.05 2.56 1.50 1.49 1.59 1.58 3.39 3.32 4.03 3.98 407.4 184.01 2.21

Case3: Vertically inserted tube Flow rate (cm3 /min) Rmin (mm) Rmax (mm) RF (mm) (R/RF )min (R/RF )max (V/VF )min (V/VF )max Qcr (cm3 /min) Q/Qcr

225 3.34 3.38

407.4 3.39 3.97 2.56 1.30 1.32 1.32 1.55 2.22 2.33 2.30 3.73 184.01 1.22 2.21

105 4.05 4.32 1.58 1.69 3.98 4.81 165 3.81 4.03 1.49 1.58 3.30 3.92

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

693

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]

D. Gerlach, G. Biswas, F. Durst, V. Kolobaric, Quasi-static bubble formation on submerged orifices, Int. J. Heat Mass Transf. 48 (2005) 425–438. J. Davidson, B. Schueler, Bubble formation at an orifice in an inviscid liquid, Trans. Inst. Chem. Eng. 38 (1960) 335–342. J. Davidson, B. Schueler, Bubble formation at an orifice in a viscous liquid, Trans. Inst. Chem. Eng. 38 (1960) 144–154. J.F. Davidson, A.M.I. Mech, B.O.G. Schuler, Bubble Formation at an Orifice in a viscous Liquid, Chem. Eng. Research and Design 75 Supplement (1997) 105–115. L. Davidson, E. Amick, Formation of gas bubbles at horizontal orifices, A.I.Ch.E. J. 2 (3) (1956) 337–342. P. Snabrea, F. Magnifotcham, Formation and rise of a bubble stream in a viscous liquid, Eur. Phys. J. B 4 (1998) 369–377. A. Albadawi, D.B. Donoghue, A.J. Robinson, D.B. Murray, Y.M.C. Delaure, On the analysis of bubble growth and detachment at low Capillary and Bond numbers using Volume of Fluid and Level Set methods, Chem. Eng. Sci. 90 (2013) 77–91. A.K. Das, P.K. Das, P. Saha, Formation of bubbles at submerged orifices – Experimental investigation and theoretical prediction, Exper. Therm. Fluid Sci. 35 (2011) 618–627. X. Zhu, J. Xie, Q. Liao, R. Chen, H. Wang, Dynamic bubbling behaviors on a micro-orifice submerged in stagnant Liquid, Int. J. Heat Mass Transf. 68 (2014) 324–331. C. Qu, Y. Yu, J. Zhang, Experimental study of bubbling regimes on submerged micro-orifices, Int. J. Heat Mass Transf. 111 (2017) 17–28. J. Xie, X. Zhu, Q. Liao, H. Wang, Y.D. Ding, Dynamics of bubble formation and detachment from an immersed micro-orifice on a plate, Int. J. Heat Mass Transf. 55 (2012) 3205–3213. L. Zhang, M. Shoji, Aperiodic bubble formation from a submerged orifice, Chem. Eng. Sci. 56 (2001) 5371–5381. D.J. Mc-Cann, R.G.H. Princes, Regimes of bubbling at a submerged orifice, Chem. Eng. Sci. 26 (1971) 1505–1515. G.A. Irons, R.I.L. Guthrie, Bubble Formation at Nozzles in Pig Iron, Metall. Trans. B 9B (1978) 101–110. L.J. Mittoni, M.P. Schwarz, R.D. La Nauze, Deterministic chaos in the gas inlet pressure of gas–liquid bubbling systems, Phys. Fluids 7 (4) (1995) April. H.N. Oguz, A. Prosperetti, Dynamics of bubble growth and detachment from a needle, J. Fluid Mech. 257 (1993) 111–145. F.J. Lesage, F. Marois, Experimental and numerical analysis of quasi-static bubble size and shape characteristics at detachment, Int. J. Heat Mass Transf. 64 (2013) 53–69. F.J. Lesage, J.S. Cotton, A.J. Robinson, Analysis of quasi-static vapour bubble shape during growth and departure, Phys. fluids 25 (2013) 067103. A. Valencia, M. Cordova, J. Ortega, Numerical simulation of gas bubbles formation at a submerged orifice in a liquid, Int. Commun. Heat Mass Transf. 29 (6) (2002) 821–830. D. Gerlach, N. Alleborn, V. Buwa, F. Durst, Numerical simulation of periodic bubble formation at a submerged orifice with constant gas flow rate, Chem. Eng. Sci. 62 (2007) 2109–2125. D. Gerlach, G. Tomar, G. Biswas, F. Durst, Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows, Int. J. Heat Mass Transf. 49 (2006) 740–754. V.V. Buwa, D. Gerlach, F. Durst, E. Schlücker, Numerical simulations of bubble formation on submerged orifices: period-1 and period-2 bubbling regimes, Chem. Eng. Sci. 62 (2007) 7119–7132. I. Chakraborty, B. Ray, G. Biswas, F. Durst, A. Sharma, P.S. Ghoshdastidar, Computational investigation on bubble detachment from submerged orifice in quiescent liquid under normal and reduced gravity, Phys. fluids 21 (2009) 062103. I. Chakraborty, G. Biswas, P.S. Ghoshdastidar, A coupled level-set and volume-of-fluid method for the buoyant rise of gas bubbles in liquids, Int. J. Heat Mass Transf. 58 (2013) 240–259. I. Chakraborty, G. Biswas, P.S. Ghoshdastidar, Bubble generation in quiescent and co-flowing liquids, Int. J. Heat Mass Transf. 54 (2011) 4673–4688. I. Chakraborty, G. Biswas, S. Polepalle, P.S. Ghoshdastidar, Bubble formation and dynamics in a quiescent high-density liquid, Am. Inst. Chem. Eng. 61 (11) (2015) 3996–4012. D. Ma, M. Liu, Y. Zu, C. Tang, Two-dimensional volume of fluid simulation studies on single bubble formation and dynamics in bubble columns, Chem. Eng. Sci. 72 (2012) 61–77. K. Loubière, V. Castaignède, G. Hébrard, M. Roustan, Bubble formation at a flexible orifice with liquid cross-flow, Chem. Eng. Process. 43 (2004) 717–725. N. Dietrich, N. Mayoufi, S. Poncin, N. Midoux, H.Z. Li, Bubble formation at an orifice: a multiscale investigation, Chem. Eng. Sci. 92 (2013) 118–125. H.K. Nahra, Y. Kamotani, Bubble formation from wall orifice in liquid cross-flow under low gravity, Chem. Eng. Sci. 55 (20 0 0) 4653–4665. J.A. Simmons, J.E. Sprittles, Y.D. Shikhmurzaev, The formation of a bubble from a submerged orifice, Eur. J. Mech. B Fluids 53 (2015) 24–36. W. Zhang, R.B.H. Tan, A model for bubble formation and weeping at a submerged orifice, Chem. Eng. Sci. 55 (20 0 0) 6243–6250. S. Quan, J. Hua, Numerical studies of bubble necking in viscous liquids, Phys. Rev. E 77 (2008) 066303. A. Georgoulas, P. Koukouvinis, M. Gavaises, M. Marengo, Numerical investigation of quasi-static bubble growth and detachment from submerged orifices in isothermal liquid pools: the effect of varying fluid properties and gravity levels, Int. J. Multiph. Flow 74 (2015) 59–78. Y. Zhang, M. Liu, Y. Xu, C. Tang, Three-dimensional volume of fluid simulations on bubble formation and dynamics in bubble columns, Chem. Eng. Sci. 73 (2012) 55–78. Y. Chen, S. Liu, R. Kulenovic, R. Mertz, Experimental study on macroscopic contact line behaviors during bubble formation on submerged orifice and comparison with numerical simulations, Phys. fluids 25 (2013) 092105. S. Vafaei, P. Angeli, D. Wen, Bubble growth rate from stainless steel substrate and needle nozzles, Colloids Surf. A Phys. Chem. Eng. Asp. 384 (2011) 240–247. S.V. Gnyloskurenko, A.V. Byakova, O.I. Raychenko, T. Nakamura, Influence of wetting conditions on bubble formation at orifice in an inviscid liquid. Transformation of bubble shape and size, Colloids Surf A Phys. Chem. Eng. Asp. 218 (2003) 73–87. A. Byakova, S.V. Gnyloskurenko, T. Nakamura, O. Raychenko, Influence of wetting conditions on bubble formation at orifice in an inviscid liquid: mechanism of bubble evolution, Colloids Surf A Physicochem Eng Asp. 229 (1) (2003) 19–32. N. Chigarev, T. Chigareva, Variation of the contact angle in the quasistatic growth of a vapor bubble on a horizontal surface in a boiling liquid, J. Eng. Phys. Thermophys. 50 (4) (1986) 398–401. S.G. Kandlikar, M.E. Steinke, Contact angles and interface behavior during rapid evaporation of liquid on a heated surface, Int. J. Heat Mass Transf. 45 (18) (2002) 3771–3780. J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modelling surface tension, J. Comput. Phys. 100 (1992) 335–354. F. Harlow, J. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182–2189. M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed—algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49. M. van Sint Annaland, N.G. Deen, J.A.M. Kuipers, Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method, Chem. Eng. Sci. 60 (2005) 2999–3011. ANSYS Fluent User’s Guide, v18.0, 2017. J. Li, Piecewise Linear Interface Calculation [Calcul d’Interface Affine par Morceaux], Comptes Rendus - Academie des Sciences, Serie II: Mecanique, Physique, Chimie, Astronomie, vol. 320, 1995/01/01, pp. 391–396.

694

X. Chen, N. Chen and A. Delgado / Applied Mathematical Modelling 73 (2019) 668–694

[49] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski, Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows, J. Comp. Phys. 152 (1999) 423–456. [50] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1986) 40–65 (ISSN 0021-9991). [51] H. Mirsandi, A.H. Rajkotwala, M.W. Baltussen, E.A.J.F. Peters, J.A.M. Kuipers, Numerical simulation of bubble formation with a moving contact line using Local Front Reconstruction Method, Chem. Eng. Sci. 187 (2018) 415–431. [52] S. Di Bari, A.J. Robinson, Experimental study of gas injected bubble growth from submerged orifices, Exper. Therm. Fluid Sci. 44 (2013) 124–137. [53] W. Fritz, Berechnung des maximale Volume von Dampfblasen, Phys. Z. 36 (1935) 379–384.