Numerical Study of Formation of a Series of Bubbles from an Orifice by Applying Dynamic Contact Angle Model

*e dynamic contact angle model is applied in the formation process of a series of bubbles from Period-I regime to Period-II regime by using the VOF method on a 2D axisymmetric domain. In the first process of the current research, the dynamic contact angle model is validated by comparing the numerical results to the experimental data. Good agreement in terms of bubble shape and bubble detachment time is observed from a lower flow rate Q� 150.8 cm/min (Re� 54.77, Period-I regime) to a higher flow rateQ� 603.2 cm/min (Re� 219.07, Period-III regime). *e comparison between the dynamic contact angle model and the static contact angle model is also performed. It is observed that the static contact angle model can obtain similar results as the dynamic contact angle model only for smaller gas flow rates (Q≤ 150.8 cm/min and Re≤ 54.77)). For higher gas flow rates, the static contact angle model cannot produce good results as the dynamic contact angle model and has larger relative errors in terms of bubble detachment time and bubble shape.


Introduction
e process of bubble formation is widely applied in various fields, such as chemistry, biology, and industrial manufacturing [1][2][3], depending on different operating parameters, such as gas flow rates, orifice diameters, contact angles, and liquid properties. Different regimes of bubble formation are classified as Period-1, Period-2, Period-3, Period-4, and chaotic bubbling regimes [4]. Period-1 is within the regime for smaller gas flow rates, bubbles with constant detachment time and volume are formed, and no coalescence is observed between neighboring bubbles after detaching from the orifice. When the gas flow rates further increase, the odd-even phenomenon with different detachment times of the bubble formation is found, and bubble coalescence happens after detaching from the orifice [5]. In this regime, the definition of Period-N is further classified by different coalescence phenomena: Period-2 for pairing and double coalescence, Period-3 for triple coalescence, and Period-4 for quadruple coalescence. And finally, when the gas flow rates are large enough, the phenomena in respect to the bubble detachment and coalescence become more complicated. In order to understand these typical physical phenomena, several experimental studies have been performed: Snabre and Magnifotcham [6], Yuan et al. [7], and Zhang and Shoji [8] researched the bubble formation with different flow rates and the force acting on the bubble and found that the higher flow rate leads to bubble coalescence near the orifice and found that the higher flow rate leads to bubble coalescence near the orifice and the smaller bubble detachment time. Di Bari and Robinson [9] studied the process of bubble formation and the influence of contact angle. ey found that the bubble formation process contains three distinct periods: nucleation stage, growth stage, and detachment stage. Davidson and Amick [10] and Albadawi et al. [11] investigated the whole process of bubble formation based on the experiment and observed when the flow rate further increases, the large coalesced bubbles themselves undergo the second coalescence as they rise. e bubbles are erratic, and the liquid becomes filled with small circulating bubbles. In addition, various influence factors such as density and viscosity were studied for a series of dynamic bubble formation process [8,12,13]. Zhang and Shoji [8] studied a series of bubble rising process after formation and interaction of rising bubble, and the main three types of bubble interaction are (a) bubble interference, (b) bubble collision, and (c) bubble coalescence. Based on the various flow fields, the bubble formation is influenced by the wake of the previous bubble and the surrounding liquid field, and the movement of the previous bubble leads to different periodic behaviours of a series of bubble formation. Dou et al. [14], Loubière et al. [15], Dietrich et al. [16], Nahra and Kamotani [17], and Simmons et al. [18] researched the bubble formation and growth stage from a submerged orifice with a liquid cross-flow, and the formation of the bubbles is strongly influenced by the surrounding liquid. In this situation, the drag force is the main factor for bubble detachment instead of buoyancy force. In addition, the influence of liquid properties such as viscosity and density is also intensively studied by numerical simulation and compared with experiment [19][20][21][22]. e detachment bubble time is larger with a higher viscosity and surface tension. Other than the numerical study, several simulation methods are generated for simulating different researching aspects, such as immerse boundary (IB) method [23], level set (LS) method [24], front tracking (FT) method [25], the lattice Boltzmann (LB) method [26,27], and the volume-of-fluid (VOF) method [28].
As discussed before, the surface tension force plays an important role in bubble formation. In order to specifically study the influence of surface tension and contact angle, several experimental and numerical studies were carried out to understand the physical background of the surface tension force. Gerlach et al. [29][30][31] found that the wettability has strong effects on the bubble formation process and investigated the behaviour of series bubble formation under different inlet flow rates. Chen et al. [32], Cano-Lozano et al. [33], and Salman et al. [34] studied the bubble formation and rising in a small tube; the larger surface tension force stabilised the bubble longer and formed a Taylor bubble. Brackbill et al. [35] used the continuum surface force (CSF) model by considering the surface tension force as a volume force for simulation. Siddiqui et al. [36], Gnyloskurenko et al. [37], Drelich et al. [38], and Yang et al. [39] researched different orifice materials with different contact angles and observed that the bubble size increases when the surface tension coefficient is higher, especially in the case of larger contact angle. Ahmed et al. [40] and Das et al. [41] studied the wettability of bubble formation from submerged orifice through experiment, and the experimental result shows that the bubble grows with a spherical shape and smaller diameter when the contact angle increases. Corchero et al. [42] and Yuan and Lee [43] used a simple model for different static contact angles, and they found numerical results agree on the experiment at small flow rates. However, some experimental studies discovered that the contact angle is not a constant value during the whole process of bubble formation. Sun et al. [44], Tadmor [45], and Chibowski and Terpilowski [46] studied the dynamic contact angle and obtained the correlation between advancing, receding, and equilibrium contact angles by experiment. Pasandidehfard et al. [47] and Nichita et al. [48] developed the dynamic contact angle model for liquid droplet impact on the solid surface, and the predictions of droplet diameter are more accurate than that of static contact angle model during spreading and at equilibrium. Mukherjee and Kandlikar [49] and Ardron et al. [50] investigated dynamic contact angle for nucleate boiling bubble based on the equilibrium thermodynamics theory. For the bubble formation processes, Mirsandi et al. [51] and Chen et al. [52,53] used two dynamic contact angle models to simulate the single bubble formation process at lower flow rates. e numerical result is improved, and the bubble detachment characteristics in terms of bubble shape show a good agreement with experiment. Göhl et al. [54] summarized different dynamic contact angle models and simulated the single nonaxisymmetric droplet and found that Shikhmurvaev's model predicts the spreading dynamics best, especially on hydrophilic surfaces. is model is also dependent on the contact line velocity through the Ca number and the equilibrium contact angle.
From the summary of the literature study above, the previous research mainly focused on the motion of liquid droplet and single bubble detachment from the orifice. On the contrary, it is also found that when the gas flow rates become larger, bubble formation becomes complicated and different periodicities will be shown with different contact angles [5,15]. So it is necessary to dig deeper into the bubble formation process with more physical contact angle model. In this paper, the dynamic contact angle model is utilized to simulate the formation process of a series of bubbles, and numerical simulations are performed with different inlet gas flow rates from Period-I to Period-III regimes. Furthermore, in order to examine the accuracy and improvement of this model, the results are also compared with both experimental results and results from the static contact angle model, which is widely used in the previous simulations, and the improvement of the dynamic contact angle is also analysed and explained. e structure of the thesis is as follows: theory and numerical models are presented in Section 2; Section 3 contains the results and discussion, including numerical validation, mesh study, and results explanation. e final conclusion is made in the last chapter.

Mathematical Model.
e mass and momentum conservation equations are used to describe the motion of the fluid, which is presented in the mass conservation equation as shown in the following equations: Here, f is the additional force term per unit volume and composed of surface tension force in the bubble formation. e stress tensor σ can be decomposed into the vertical pressure p and viscous shear stress τ: σ � −pI + τ, in which I

2
Advances in Materials Science and Engineering is the unit tensor. In Newtonian fluids, the viscous stress can be expressed as where µ is the dynamic viscosity and e is the rate of strain tensor based on velocity gradient. Inserting equation (3) into (2), the equation of momentum conservation can be written as follows: In the liquid-gas system, the numerical technique of capturing the physics of the interface between liquid and gas is one of the vital problems to obtain meaning results. In the previous research, several numerical methods are developed to simulate the interface, including volume-of-fluid (VOF) method, level set method, immersed boundary method, front tracking method, and lattice Boltzmann method [31]. e 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 reinitialization [25,55,56]. In this method, it utilizes an additional variable ∅ (defined as volume fraction) to indicate different phases: if ∅ � 1, it means liquid; if ∅ � 0, it means gas; and the interface is defined when ∅ is between 0 and 1.
It can be noticeable that the sum of volume fraction of different phases must be 1 for the whole domain. By using the volume fraction, the mixture properties such as density ρ and viscosity μ are given by the following equations, respectively [55]: where the subscripts f and g denote the liquid and gas, respectively. e volume fraction ∅ equation follows the general transport equation [55]: In the VOF method, equations (1), (2), and (7) are solved. For the liquid-gas system, the external force f in equation (2) on the interface should be taken into consideration. Based on the research of Brackbill et al. [35], this force is expressed in the following equation: where κ is the curvature defined as κ � ∇ · n, n is the unit outer-normal vector defined as n � (∇∅/|∇∅|), ρ is the density of the mixture calculated in equation (5), and ρ g and ρ f are density of gas and liquid, respectively. For time-dependent multiphase simulations, the pressure implicit with split operator (PISO) [57] algorithm is used to calculate the pressure and velocity field, and then, the volume fraction is updated with the new velocity field. e basic procedure of the VOF method is shown as follows [28]: (1) e intermediate velocity field at the new time step is calculated based on a guessed pressure field by using equation (2). (2) e guessed pressure field and the intermediate velocity field are further modified by using equation (2). In the PISO algorithm, the correction step will be performed twice to reduce error and enhance mass conservation [28]. (3) Finally, the volume fraction ∅ is updated by solving equation (7) and a new interface is reconstructed based on the PLIC algorithm [58].
All simulations are performed by using the software Ansys Fluent. e basic settings are as follows: the implicit Euler method is used for time discretization; the PRESTO algorithm is used for the pressure discretization, and the 2nd order upwind scheme is chosen for the convection and diffusion terms in the momentum equation. e detailed description of the PLIC algorithm is shown in [58,59], and from the previous study [5], the sharpness of the discontinuous volume fraction extends over 3 cells. More description can be found in Appendix A.

Dynamic Contact Angle Model. Experimental studies
show that the contact angle does not stay constant during the process of bubble expansion, shrinkage, and detachment. Tadmor [45] and Chibowski and Terpilowski [46] study the dynamic contact angle from the experiment and found the relationship between static contact angle and dynamic contact angle, which is shown as follows: where Here, θ A is the contact angle during the bubble expansion phase, which is called as advancing contact angle, and θ R is the contact angle during the bubble neck shrinking phase when the bubble is about to detach from the orifice and is called the receding contact angle. Please notice that these two angles are defined in the air phase (secondary phase) in the process of bubble formation, while in the research regarding droplet movement, these two angles are defined in the liquid phase. ese two angles define the upper and lower bounds of the contact angle during the whole bubble formation process: in the dynamic contact angle model, the apparent contact angle θ is restricted within the range between receding contact angle θ R and advancing Advances in Materials Science and Engineering [45,46,51]. θ c is the equilibrium contact angle, while r A and r R are the coefficients for advancing and receding contact angle, respectively (demonstrated in Figure 1).
Compared to the static contact angle model, this dynamic contact angle model utilized two different angles to define the direction of the surface tension force for different phases during the whole process of bubble formation. Chen et al. [52,53] developed contact line velocity-dependent model for modelling the transition between the advancing and receding contact angles: when the instantaneous contact line velocity is between the given velocity −u ct, max and u ct, max , the contact angle is changed between the advancing contact and receding contact angle; during the bubble expansion phase, if the contact line velocity is larger than the u ct, max , the contact angle no longer grows and maintains the advancing contact angle; the same treatment is also applied during the bubble neck shrinking phase when the contact velocity is lower than the −u ct, max , and the contact angle remains the receding contact angle [51,52]. Equation (11) summarized this concept for the contact line velocity-dependent model: In our study, the VOF method and the contact line velocity-dependent model are used to simulate the bubble formation process. e basic solving procedure is briefly summarized as follows: (1) e contact line position is found through interface capturing by search the location with the value of volume fraction between 0 and 1. (2) e velocity of the contact line is obtained, and the contact angle is calculated based on equation (11). (3) e contact angle value is calculated and returned for the PISO algorithm. (4) Repeat Steps 1-3 for the next time step.

Dimensionless Number.
ere are several dimensionless numbers used to characterize specific multiphase problems. Reynolds number (Re), Bond number (Bo), capillary number (Ca), and Weber number (We) are utilized in the present study. e Reynolds number is defined as the ratio between the inertial force and the viscous force: where ρ g is the density of gas, Q and d a is flow rate and the diameter of orifice, u is the inlet velocity, and μ g is the dynamic viscosity of the gas.
e Bond number (Bo) is defined as the ratio between the gravitational force and surface tension force: Ca is the capillary number for the rate between the viscous force and the surface tension force, defined as We is the Weber number for the rate between the inertial force and the capillary force, defined as e model of the dynamic contact angle is mainly based on the contact line velocity and Ca, where the bubble formation is mainly based on the Re, Bo, and We.

Computational Domain and Boundary Condition. A 2D
axisymmetric domain is used to simulate the process of the bubble formation from the orifice, as shown in Figure 2. From the previous research of Gerlach et al. [30,31] and Vafaei et al. [60], this assumption is physically applicable for gas injection from a single orifice. e geometry setting up is proportionally based on the inlet orifice diameter d a : the sidewall is 20d a away from the orifice; the height of the computational domain is 50d a . Both are chosen in order to minimize the effects of the boundary condition and visualize enough bubbles within the computational domain. e gas is injected from the gas inlet at different flow rates for different cases. e boundary condition at the inlet is a fully developed parabolic flow profile and set by the userdefined function (UDF), where the maximum velocity is two times of the average velocity u max � 2u 0 ; the pressure outlet condition is defined on the outlet boundary; the side and bottom wall boundaries are set as nonslip walls, and since the bottom wall is in contact with the bubble, the dynamic contact angle is imposed here.
Based on the previous study [5], the mesh with 260000 number of cells (416 × 625) is used and the time step is chosen as 5 × 10 − 6 s. A more detailed description of the mesh and time step study can be found in Appendix B.

Numerical Study and Comparison of the Dynamic Contact
Angle Model. Based on the previous study, the dynamic contact angle model is used to simulate the formation process of a series of bubbles and the numerical results are compared to the experiments [31]. Simulations with four different flow rates from a lower flow rate Q � 150.8 cm 3 /min (Re � 54.77, Period-I regime) to a higher flow rate Q � 603.2 cm 3 /min (Re � 219.07, Period-II regime) are investigated. e equilibrium contact angle is 100°given in [31]. rough the measurement of advancing contact angle from the experimental picture of flow rate Q � 301.6 cm 3 / min and Q � 603.2 cm 3 /min, the average of advancing contact angle θ A � 130°is used (shown in Figure 3). Based on equation (10), the value of the receding contact angle θ R is calculated as 84°. e general input data are listed in Table 1; based on the previous study for maximum contact line velocity u ct, max [32], the value u ct, max � 0.02 m/s is used for the simulation. Figure 4(a) shows the snapshots of bubble formation with a lowest gas flow rate Q � 150.8 cm 3 /min (Re � 54.77). Since this gas flow rate is within the Period-1 regime, the bubble rises individually after formation without further splitting or coalescence, and the detachment time maintains constant. e numerical results by using the dynamic contact angle model are shown in Figure 4(b). Comparing to the experimental results, the detachment time of the dynamic contact angle model numerical result maintains 51.5 ms with a relative error of 0.96% and the bubble shape is also similar to the experimental results. e numerical simulation with static contact angle model is also performed, and the results are shown in Figure 4(c); the bubble detachment time by using a static contact angle model is 51 ms with 1.92% relative error. It can be seen from the results that, for a lower gas flow rate, both dynamic contact angle model

Advances in Materials Science and Engineering
and static contact angle model can produce the comparable results in terms of bubble detachment time and bubble shape, while the improvement of applying dynamic contact angle model is not obvious in this case. Figure 5 shows the bubble profile comparison between the dynamic and static contact angle. It can be seen that, for the case with a lower gas flow rate, the differences in bubble shapes between these two models are not so obvious during the whole process of bubble formation. is is probably due to the smaller value of the contact line velocity: when the gas flow rate is small, the contact line velocity is not high enough such that the upper and lower bounds of the contact angles are not reached yet. Furthermore, this conclusion further confirms the results of the previous numerical study by using the static contact angle model [14,[29][30][31]: good agreements are also observed between the experimental data and numerical results for a smaller gas flow rate.
When the flow rate increases and Period-2 regime is reached, the bubble detachment time is not a constant value. Figure 6(a) shows the experimental snapshots with gas flow rate Q � 301.6 cm 3 /min (Re � 109.54), and two different detachment times are shown during the bubble formation. With this gas flow rate, the coalescence of the rising bubbles also appears after the formation process and the influence of the flow field from the previous bubble becomes more significant. By checking the numerical results with the static contact angle model and the dynamic contact angle model, both models can still obtain two different detachment times. However, the difference in these two models becomes larger: the detachment time from dynamic contact angle is closer to experimental results with a relative error of 4.9%, while the detachment time of the static contact angle model has a larger relative error around 12.2%. Furthermore, when comparing to the differences in the two detachment times in the experimental data, the dynamic contact angle model has the same value (2 ms), while the value obtained by the static contact angle model is larger (4 ms). Figure 7 shows the bubble profile comparison between the two contact angle models during the bubble formation process. It can be seen that, by using the dynamic contact angle model, the pinch-off of the bubble neck is slightly thinner and the overall bubble shape is slightly higher in the advancing contact angle period because of the larger advancing contact angle. Moreover, these differences remain during the receding periods and finally lead to a smaller bubble detachment time. e reason for the different profiles in the early stage during the formation process is due to the larger gas flow rate, which leads to a larger contact line velocity. Figure 8 shows the comparison of the surrounding liquid flow field of the corresponding snapshots in Figures 6(b) and 6(c) for the dynamic contact angle model and static contact angle model, respectively. From this figure, the coalescence of the bubbles shown in the second snapshots of Figures 8(a) and 8(b) is slightly different between these two models, and this causes slightly different gaps as well as the flow fields. From the results shown in Figures 7 and 8, two factors cause the differences in the detachment time: bubble profile in the beginning and the surrounding liquid field. e reason lies in the fact that there are the two different apparent contact angles in the dynamic contact angle model, and these two angles are used in different stages during the bubble formation so that this model can obtain closer results to the experiments than the static contact angle model.
When the gas flow rate further increases to Q � 452.4 cm 3 /min (Re � 164.30), the bubble detachment time decreases further and the bubble coalescence happens closer to the orifice. As is shown in Figure 9, both contact angle models show good agreement in terms of bubble shape and two different detachment times can also be obtained. However, compared to the experimental value of bubble detachment time (40 ms and 35 ms), the numerical results of dynamic contact angle model have the relative error of the bubble detachment time 2.5% and 2.9%, while the bubble detachment time relative error of a static contact angle model is 7.5% and 11.4%. Based on these results, it concludes that when the gas flow rate increases further, when comparing to the static contact angle model, the dynamic contact   Figure 10 shows the comparison of profile between dynamic and static contact angle models. Similar to the previous case (Re � 109.54), the profile from the dynamic contact angle model is higher especially for the second formed bubble with a smaller detachment time, and the difference between two contact angle models is already shown up during the advancing contact angle period. In addition, due to the further increment of Re, the profile difference is more obvious. Figure 11 shows the flow field comparison between two contact angle models. e bubble shape shown in the fourth snapshots between Figures 11(a) and 11(b) is already slightly different: the bubble shape from the dynamic contact angle model is slightly higher and thinner, and the gap between the previous bubble and the 52ms 52 ms Advances in Materials Science and Engineering 7 current bubble is slightly smaller. Due to this smaller bubble gap, the influence of the former bubble becomes slightly stronger and the surrounding fluid compresses the bubble from both sides of the pinch-off neck even harder. As a result, these differences cause different bubble detachment times between these two models. When the flow rate increases to Q � 603.2 cm 3 /min (Re � 219.07), the Period-III is reached. Within this regime, the bubble detachment time decreases further and the time difference of two detachment time becomes even larger. As a result, the coalescence between the former rising bubble and the current bubble is closer to the orifice and the formation of the current bubble is strongly influenced by the previous one, which leads to a long flat shape. In this case, both models can still produce two different bubble detachment times yet with different values. Comparing to the experimental results in terms of bubble detachment time (40 ms and 33 ms), the numerical bubble detachment time of dynamic contact angle model has relative error 2.5% and 3.0%, while the corresponding relative error of static contact angle model is around 7.5% and 9.1%. e comparison of the bubble profile is shown in Figure  12. Similar to the previous conclusion, the difference in bubble profile between two contact angle models becomes even larger during the advancing contact angle period, and the bubble shape obtained by using the dynamic contact angle model is higher due to a larger contact line velocity and the stronger influence of surrounding flow field when Re further increases.
is phenomenon obeys the same tendency from the results with smaller gas flow rates, and the same conclusion can be applied. Figure 13 shows the comparison of flow field between two contact angle models. From the flow field shown in Figure 13(b), due to the splitting of the previous rising bubble, the difference in the flow field between two contact angle models becomes even larger and finally leads to different detachment times. e reason for the bubble splitting is further explained in Figure 14.
e first snapshot in Figure 14 shows different contact angles generating different bubble shapes for the first bubble, and a relatively longer bubble is formed by using the static contact angle model.
is bubble leads to splitting of the previous rising bubble in the later stage of the formation. On the contrary, the bubble shape obtained from the dynamic contact angle model is not too high to split the previous bubble, which agrees with the experimental snapshots. Table 2 summarizes the detachment time for two contact angle models and the corresponding relative errors when comparing with experimental data. Based on the study of all four gas flow rates, one can safely come to the conclusion that the detachment time of the dynamic contact angle 8 Advances in Materials Science and Engineering model has a smaller relative error and shows obvious improvement at the higher flow rates when compared to the static contact angle model, especially in Period-II and Period-III. e advantage of dynamic contact angle model is utilizing two different values to adjust the shape of the bubbles for different phases during the formation process, and this model is closer to the actual physics at a higher Re [51]. For a smaller Re of single bubble formation within the regime of Period-I, both dynamic and static contact angle models can deliver similar results and good agreements to the experimental results are observed because of smaller contact line velocity and similar profile. When Re increases

Conclusion and Outlook
In this study, the dynamic contact angle model for the formation of a series of bubbles is studied. e dynamic contact angle model is compared to the experimental results and static contact angle model from lower to higher Reynold's numbers covering Period-I regime to Period-III regime. e improvement of the dynamic contact angle model is also studied and explained by means of bubble shape and surrounding liquid fields. Based on the current study, the following conclusions can be drawn. For the lower Re in the Period-I regime with individual bubble formation, both models deliver good results in terms of detachment time and bubble shape comparing the experimental results, and this is because the apparent contact angle, as well as the fluid field, is still similar during the formation process and the fluid field.
When flow rate increases to Period-II regime, the bubble detachment time is not a constant value anymore and the differences between the dynamic contact angle model and the static contact angle model become larger with the increment in the gas flow rates. Results by using the dynamic contact angle model are closer to the experimental data in terms of bubble shape and detachment time, while the results with a static contact angle model have a larger relative error. is is because of the variation in the bubble profile in the advancing contact angle period and the bubble coalescence behaviour. e influence of the surrounding flow field in combination with the effect of contact angle also causes different bubble detachment times and bubble shapes. In the regime of Period-III, the static contact angle model is not comparable to experimental in terms of bubble shape, while the dynamic contact angle model still has good agreements.
For the further study, the bubble coalescence problem can be investigated in Period-IV and even chaotic regimes together with the dynamic contact angle model, and a 3D domain can be used for such high flow rates to solve the axis-symmetric problem. Appendix e content of the appendix is from the published paper [3] from the same author. In order to avoid plagiarism, the content will not be repeated again. However, in order to make it convenient for the reviewers and the readers, the contents are shown in Appendix.
A. PLIC Algorithm and the Sharpness of the Interface e calculation and reconstruction of the volume fraction is based on PLIC (1) 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. e calculation details are shown in (1,2). e advantage of the PLIC algorithm is the higher accuracy and sharpness. Figure 16 demonstrates the sharpness after smoothing over the discontinuous volume fraction: the interface extends over only 3 cells.

B. Mesh Independence Study and Time
Step Study e mesh independence study is based on the previous study of  Table 3 shows the quantitative comparison of the numerical results with different mesh sizes, and Figure 17 shows the snapshots at detachment times with different mesh sizes. Hence, the numerical results are based on the mesh size 0.06 mm, and the total number of mesh cells is 260000, about 44% less than if the mesh size was taken to be 0.05 mm. Figure 18 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.

Data Availability
All the data used are included in the paper and can be used for further verification.

Conflicts of Interest
e authors declare that there are no conflicts of interest.