Numerical simulation of bubble formation with a moving contact line using Local Front Reconstruction Method

The process of adiabatic bubble formation from an orifice plate occurs in various industrial applications. It is important to understand the dynamics of bubble formation and to develop numerical models to accurately predict the formation dynamics under various operating conditions. For the numerical models, an appropriate contact line boundary condition is necessary since this process may involve a moving contact line, which significantly affects the bubble departure size. In this paper, we extend the Local Front Reconstruction Method by incorporating contact angle dynamics. The predictions of the improved model are extensively verified and validated with experimental and numerical data available in the literature. The problem of three-dimensional bubble injection from an orifice into quiescent water using various volumetric flow rates is used to assess the numerical model under capillary dominant conditions and conditions where the interplay between inertial, viscous, surface tension, and buoyancy forces cause a complex interface deformation. 2018 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The formation of gas bubbles by submerged needles or orifices is of great interest for diverse applications in chemical, nuclear, and metallurgical industries, because it influences the bubble size and thereby the bubble rise velocity. To optimize the bubble formation (e.g., the formation time, bubble volume and motion of the contact line), the fundamental physics should be understood. However, some of this fundamental knowledge is still lacking (Kulkarni and Joshi, 2005).
Several studies on the formation of bubbles have already been performed using both experimental and numerical approaches. The majority of these studies focused on the case of bubble formation using a constant volumetric gas flow rate through a submerged orifice (McCann and Prince, 1971;Zhang and Shoji, 2001;Badam et al., 2007;Das et al., 2011). Experimental and theoretical studies have shown that the forces acting on the bubble during the growth are divided into two main groups according to their influence on the formation. The first group causes the bubble detachment. This group includes the buoyancy, gas injection momentum, and contact pressure forces, which emerges due to the pressure difference inside and outside of the bubble over the contact area. The second group includes viscous, surface tension, and inertial forces. This group resists the bubble detachment and tends to keep the bubble attached to the orifice (Duhar and Colin, 2006;Di Bari and Robinson, 2013).
For single bubble growth at a constant gas flow rate, three formation regimes have been identified. The first regime is the static regime, which prevails at low gas flow rates. In this regime, there is equality between buoyancy and surface tension forces. The volume of the formed bubbles is independent of flow rate; decreasing the flow rate will increase the formation time. At higher flow rates, bubble formation enters the dynamic regime. Because the bubbling mode is more complicated, this regime is divided into six bubbling modes. Contrary to the static regime, the formation time is constant while the bubble volume is affected by the flow rate. The bubble evolution is governed by the interplay of inertial, viscous, surface tension, and buoyancy forces. Under even greater flow rates, bubble formation enters the turbulent regime. In this regime, successive bubbles coalesce with each other close to the orifice and they rise only a small distance above the orifice before they are shattered into many small bubbles of varying size by the strong turbulence (McCann and Prince, 1971).
Several phenomenological models were proposed to describe the bubble growth and detachment based on extensive experimental data sets. Earlier theoretical studies focused on the development of an analytical model to derive a scaling law for the volume of a detaching bubble assuming that the bubble retained a spherical shape (Davidson and Schüler, 1960a;Davidson and Schüler, 1960b;Kupferberg and Jameson, 1969;Gaddis and Vogelpohl, 1986). Marmur and Rubin (1976) proposed a more realistic mathematical model for non-spherical bubbles to predict the bubble volume by dividing the bubble interface into finite differential elements. In addition, a simple Young-Laplace equation could be used to describe the bubble motion at sufficiently low gas flow rates (Marmur and Rubin, 1973;Gerlach et al., 2005;Lesage and Marois, 2013). Although these approaches are in good agreement with experiments in certain regimes, the majority of them are unable to account for the viscous effects, bubble interactions and the last phase of the neck pinch-off at detachment. More advanced mathematical models are necessary to accurately describe the whole parameter space of interest.
Due to the fast developments of the computational resources and numerical methods over the last decades, it is possible to conduct detailed Direct Numerical Simulations (DNS) of multiphase flows. These simulations are performed by either interface capturing methods or interface tracking methods, since moving grid methods are often limited to moderately deforming interfaces (Baltussen et al., 2014;Weber et al., 2017). Both interface capturing and interface tracking methods assume a one fluid approximation, whose properties are determined from the interface position. The essential difference between these methods is the interface treatment. In interface capturing methods, the interface is reconstructed from an indicator function, which is advected by the fluid velocity on a fixed Eulerian grid. The most widely used interface capturing approaches are the Volume of Fluid (VOF) method (Hirt and Nichols, 1981), the Level Set (LS) method (Osher and Sethian, 1988) and a combination of these two, known as Combined Level Set and Volume of Fluid (CLSVOF) method (Sussman and Puckett, 2000). All these three methods have been extensively used to study multiphase flows including bubble formation at an orifice (Valencia et al., 2002;Gerlach et al., 2007;Chakraborty et al., 2009Chakraborty et al., , 2011Ma et al., 2012;Yujie et al., 2012;Georgoulas et al., 2015). In the front-tracking methods, the motion of the interface can be captured more accurately because the method tracks the interface with separate Lagrangian interface elements, improving also the surface tension calculation (Unverdi and Tryggvason, 1992). In spite of its attractiveness, the original front-tracking method does not allow interface merging or breakup. Therefore, it cannot be applied for the simulation of bubble growth and detachment. Several authors have proposed new techniques that enable the front-tracking method to automatically and robustly model the merging and breakup of interfaces (Shin and Juric, 2002;De Sousa et al., 2004;Quan and Schmidt, 2007;Shin et al., 2011). In the present work, we have chosen a front-tracking algorithm that allows interface merging and breakup, the Local Front Reconstruction Method (LFRM) (Shin et al., 2011). The LFRM has been validated with several benchmarking tests and the results demonstrate excellent performance in maintaining detailed interfacial shapes and local mass conservation especially at low-resolution.
In many cases of formation of bubbles from an orifice plate, the bubbles are not pinned to the orifice rim. This movement of the contact line will introduce an extra aspect in the calculation, namely the wettability. When a bubble is formed on a thinwalled nozzle, the behavior of the moving contact line has little influence on the diameter of the detached bubble (Oguz and Prosperetti, 1993). However, when forming on an orifice plate, both the apparent contact angle and the contact line diameter vary as the bubble grows in size. Unfortunately, the experimental and numerical studies on the contact line behaviors during bubble formation are very limited and as far as the authors know no universal pattern of the time-history of the contact angle and contact diameter has been determined (Chigarev and Chigareva, 1986;Kandlikar and Steinke, 2002;Gnyloskurenko et al., 2003;Chen et al., 2013). Experiments show that for a hydrophilic surface, the contact line will not recede very far from the orifice, or it may even remain pinned to the orifice rim for the whole bubble formation process Byakova et al., 2003). On the other hand, the contact line of a bubble forming on a hydrophobic surface will recede much further, thus increasing the formation time and the volume of the bubble formed .
The introduction of a moving contact line in the mathematical model results in a number of challenging issues. First of all, a complete mathematical representation of the motion of the triple contact line is still a problematic task. It is well known that the no-slip boundary condition yields stress singularity at the contact line since the fluid velocity is finite at the free surface but zero on the wall (Huh and Scriven, 1971). This singularity is usually removed by relaxing the no-slip boundary condition with a slip model. The most common approach is to use a relation between the apparent contact angle and a contact line velocity using a contact line model derived from theory (e.g. that of Cox, 1986;Blake, 2006;Kistler, 1993). The contact line position is then determined by specifying the apparent contact angle computed from the contact line model. Typically, no-slip boundary condition is still imposed on the wall for the fluid velocity since slip condition may lead to a singularity in the pressure (Sui et al., 2014). The velocity components tangential to the wall at the nearest grid node or marker point are then used as the contact line velocity (Francois and Shyy, 2003;Saha and Mitra, 2009;Yokoi et al., 2009;Muradoglu and Tasoglu, 2010). Although, the unresolved contact line region is somewhat large for DNS, Kafka and Dussan (1979) showed that for a nanometer slip length, an interfacial angle at a distance to the contact line ranging from Oð10 nmÞ to Oð10 lmÞ leads to no significant differences in the outer region.
Several attempts have been made to provide macroscopic models of the contact line dynamics based on the microscopic physics for droplet simulations (Sui et al., 2014). However, only a few numerical studies on bubble formation have taken into account the effects of moving contact line. The majority of the authors simulated the bubble formation process with the contact angle kept equal to the static contact angle (Higuera, 2005;Higuera and Medina, 2006;Gerlach et al., 2007). Chen et al. (2013) studied the contact line behaviors during bubble formation using two kinds of contact angle models in their LS method: a contact line velocity dependent model and a stick-slip model. They showed that the stick-slip model is more accurate in describing the contact line dynamics.
In the present work, we improve the 3D-LFRM by incorporating contact angle dynamics and apply the combined method to the simulations of bubble formation from submerged orifice. The paper is organized as follows: The numerical model is presented in Section 2. Section 3 presents the verification of the contact angle boundary condition implementation. In Section 4, the predictions of the improved numerical model are validated against experimental data on adiabatic bubble formation available in the literature. Next, the performance of the numerical model is assessed by simulating bubble growth and detachment under various volumetric flow rates. Finally, the summary of the main conclusions of the present work is provided in Section 5.

Numerical model
The numerical model used in this paper is based on the LFRM, established by Shin et al. (2011). In the sections below, we will describe the main characteristics of the numerical model.

Governing equations and solution methodology
In this study, the fluids are assumed to be incompressible, immiscible and Newtonian. Furthermore, a one fluid formulation is used to describe the fluid flow for all phases. The governing mass and momentum conservation equations are expressed as follows: where u is the fluid velocity, p is the pressure, and s is the stress tensor given by l ru þ ru . The local averaged density q and dynamic viscosity l depend on the local fluid phase distribution and hence are evaluated from the local phase fraction F between fluid phases. The local volumetric force accounting for the effect of surface tension, F r , is obtained by employing the hybrid Lagrangian-Eulerian formulation representation of Shin et al. (2005) to reduce unphysical parasitic currents in the vicinity of the interface using: where r is the surface tension coefficient and j H is twice the mean interface curvature field calculated on the Eulerian grid using the information from the Lagrangian interface.
The equations are spatially discretized on a regular threedimensional staggered grid using the finite volume approach. The solution of the governing equations is achieved using a two-step projection-correction method. In the first step, an explicit estimate of the velocity u nÃ is calculated, using information from the previous time step (indicated with the superscript n): The convective term is discretized using a second order fluxdelimited Barton scheme and the diffusion term with a second order central scheme. All the contributions in Eq. (1) are treated explicitly, except for the diffusion term, which is treated semiexplicit. The implicit part of the diffusion term is chosen such that the velocities in the different directions can be solved separately.
After the projected velocities are calculated, velocity field is corrected such that it satisfies the continuity equation via the following equation: r Á Dt q n r dp where dp ¼ p nþ1 À p n is the pressure correction. The equations in both the projection step and the correction step are solved with an efficient parallel Block Incomplete Cholesky Conjugate Gradient (BICCG) method (Saad, 2003).
After the fluid flow has been calculated, the Lagrangian marker points, which are used to track the interface (Section 2.1.1), are moved using locally interpolated fluid velocities. The fluid velocities are interpolated from Eulerian to Lagrangian grid points using cubic spline interpolation and the marker points are displaced with the interpolated fluid velocities using a 4th order Runge-Kutta time stepping scheme. Finally, the phase fraction in each Eulerian cell is updated using a geometrical method based on the location of the marker elements (Dijkhuizen et al., 2010). The density of each Eulerian cell is calculated by weighted averaging with the phase fraction. Following Prosperetti (2002), the local viscosity is obtained by harmonic averaging of the kinematic viscosities, based on the phase fraction.

Interface tracking in LFRM
Interface tracking is an essential part of the front-tracking method. However, due to the advection of the interface, the size of the marker elements changes, decreasing the quality of the interface mesh. To overcome this, a reconstruction procedure is required to ensure that each part of the interface is represented with an adequate resolution. In LFRM, the reconstruction procedure consists of the following steps as illustrated in Fig. 1:

(a) Localization to reconstruction cell
The interface of the discrete phase is cut by a reconstruction grid such that each part of the interface lies completely in one reconstruction cell. The reconstruction cell is independent of the Eulerian cell. In this paper, the size is usually set to two times smaller than the Eulerian cell size to ensure adequate interface resolution. This localization process makes the interface elements exclusive to each reconstruction cell.

(b) Edge line reconstruction on cell faces
An area conservative reconstruction is applied on each reconstruction cell face that contains an original interface edge line. The reconstructed edge lines consist of two original edge crossing points and an area fitting point. For cases where there is more than one edge crossing point at one of the twelve edges of the reconstruction cell, a marker reduction algorithm is applied to reduce the number of edge points to less than or equal to one (Shin et al., 2011). Even though this procedure is done on each reconstruction cell independently, all elements are still connected implicitly because the same edge crossing points will be assigned to neighboring reconstruction cells.

(c) Centroid calculation
After performing edge line reconstruction on the six faces of the reconstruction cell, intermediate interface elements with closed edge lines will be generated. Then, edge crossing points and area fitting points are connected to a centroid, which is determined by simply averaging those points.
The final step is moving the centroid in the normal direction of the intermediate interface such that the original volume of the dispersed phase in the given cell is conserved.
The reconstruction procedure above deals with an ordinary case without interface merging or breakup. For the case of merging or breakup of interface elements, four edge crossing points are located at any of the six faces of the reconstruction cell even after performing marker reduction algorithm. These cells are identified as merging/breakup cells and reconstructed using the so-called tetra-marching procedure that will ensure continuous interface elements after reconstruction (Yoon and Shin, 2010).
In the original LFRM, the marker points are not unique, i.e. marker points are stored separately for each marker elements. To reduce the memory usage and increase robustness, marker points are renumbered using prespecified point numbers distributed in the reconstruction grid creating a unique set of marker points, as illustrated in Fig. 2. Another advantage of unique marker points is that it enables the implementation of volume conserving smoothing procedure of Kuprat et al. (2001). The implementation of this procedure will prevent volume errors due to advection of marker points and interface reconstruction, which can accumulate significantly during a simulation. A detailed explanation of its implementation in the front-tracking method can be found in the work of Roghair et al. (2016).

Contact line dynamics
In the formation of bubbles at an orifice, the dynamics of the contact line plays a major role in wetting-dewetting phenomena. As mentioned before, a stress singularity arises due to the combination of the moving contact line and the no-slip boundary. This singularity can be relieved by relaxing the no slip boundary condition in the neighborhood of the contact line. The most common approach is to use a contact line model derived from theory or experiment that relates the contact line velocity to an apparent contact angle.
In the LFRM, it is relatively straightforward to implement a contact line model since the interface location is explicitly tracked. In the present work, we impose the apparent contact angle during the interface reconstruction by modifying the locations of contact points, as sketched in Fig. 3. Specifically, the position of the contact line, which is represented by the position of marker point located on the wall, cp, is extrapolated based on the position of the marker point crossing the threshold line, tp, contact angle value, h D , and interface normal parallel to the solid wall at the threshold point, n tp . In the simulations, the threshold length, h, is taken as the size of the reconstruction cell. Therefore, marker points crossing the threshold line will always be generated during the interface reconstruction when the distance between interface and the solid wall is closer than the prespecified threshold length. Following Francois and Shyy (2003), Yokoi et al. (2009), Saha and Mitra (2009), Muradoglu and Tasoglu (2010), we still impose the no-slip boundary condition on the wall. The velocity tangential to the wall at the marker point crossing the threshold line is then taken as the contact line velocity.
The apparent contact angle may have a single value (static) or different values (dynamic) depending on the contact line model used in the simulation. For the static contact angle model, a constant value is imposed during the entire course of the computations. For this study, we have implemented two kinds of dynamic contact line models, a stick-slip model (Model-A) and a velocity dependent model (Model-B).

Model-A (stick-slip model)
This model is an experimentally inspired contact line model proposed by Chen et al. (2013). In this model, there is no explicit relationship between the apparent contact angle and contact line velocity. The apparent contact angle h D is limited by the prescribed values of the static receding and advancing contact angles as follows: The variations of the contact angle and contact line position are shown schematically in Fig. 4. Consider a-b as the initial interface (Fig. 4, left). When the interfacial velocity at a distance h from the wall is positive, i.e., interface point b moves to c after a time step. Due to hysteresis, the contact point a will remain pinned, changing the contact angle from h b to h c . When h c < h rec , the contact line position will move from a to d to set the angle to h rec . The contact angle and contact line position for an advancing interface is determined in a similar way (Fig. 4, right).
In the current framework, the contact angle value associated with the contact point is identified by calculating the interface normal for the near wall interfaces. If the current contact angle is less than the prescribed receding angle h rec or greater than This schematic is a simplification of real interface, which is typically composed of more than two reconstruction cells). the advancing angle h adv , then the interface will be reconstructed with the given receding or advancing contact angle, respectively. If the contact angle is between the advancing and receding angles, we use the current contact angle without modification.

Model-B (contact line velocity dependent model)
In this model, Kistler's empirical correlation is used to dynamically relate the apparent contact angle h D to the contact line velocity, v cl , via the capillary number (Ca cl ¼ l l v cl =r) (Kistler, 1993). Although Kistler's correlation can treat non-wetting, partially wetting and fully wetting cases, the correlation is only valid for small capillary numbers. Therefore, the model has been slightly modified to deal with larger Ca following Muradoglu and Tasoglu (2010): where f À1 Hoff is the inverse of Hoffman's function defined as f Hoff ðxÞ ¼ arccos 1 À 2 tanh 5:16 In the equation above, h Di is the intermediate apparent contact angle computed using Kistler's empirical correlation, h e is the equilibrium (static) contact angle and Ca clm is defined as Ca clm ¼ minðCa max ; Ca cl Þ where Ca max is the cut-off capillary number. The apparent contact angle is computed depending on the sign of contact line velocity as follows (Muradoglu and Tasoglu, 2010): The contact line velocity v cl can be specified as the velocity of the marker point on the solid surface. However, since this marker point velocity depends on the position and thus on the apparent contact angle, the calculation of the contact line velocity has to be performed iteratively. This iterative procedure has convergence issues (Muradoglu and Tasoglu, 2010). Therefore, as mentioned before, we use the contact line velocity of the marker point crossing the threshold line as the contact line velocity. Once the apparent contact angle is determined, the marker point crossing the threshold line is connected to the solid surface such that the angle is equal to the apparent contact angle (Fig. 3).
The interface reconstruction is performed at each time step after the advection of marker points to adjust the contact line position according to the contact line model. The volume loss or gain due to the modification of the geometry near the contact region is then corrected by distributing the volume correction globally using the aforementioned smoothing procedure. However, marker elements connected to the solid wall are excluded from this procedure so that the apparent contact angle is preserved. It was verified that the volume of the bubble is conserved within machine precision by applying this global volume correction.

Verification and validation of the contact angle modeling
In this section the verification and validation test cases for the LFRM with contact angle boundary condition are presented. To test the static contact angle boundary condition, the problems of equilibrium shape of a 3D droplet on a horizontal surface, with and without the presence of gravity are investigated. The simulations are compared with known analytical solutions. To check the dynamic contact angle boundary condition, the impact and spreading of a 3D droplet on a horizontal surface are compared to experimental and numerical results available in literature. For all simulations presented in this paper, the time step Dt has been chosen such that it satisfies both Courant-Friedrichs-Lewy (CFL) and capillary time step restrictions as follows (Brackbill et al., 1992): Dt < Dt r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Here, D is the size of the computational grid, v max is the maximum fluid velocity in the computational domain, Dt CFL and Dt r are the time step restrictions according to the CFL and capillary criteria, respectively. The threshold distance, h, is taken as the size of reconstruction cell, which is half the size of Eulerian cell.

Equilibrium shape of a droplet on a horizontal surface
In this section, we consider the shape at equilibrium of a 3D droplet released on a horizontal wall to check the implementation of the static contact angle. The droplet is initially a semicircle with initial radius R 0 ¼ 1 Â 10 À2 m and initial contact angle h ¼ 90 (Fig. 5). If the static contact angle, h s , is different from the initial angle, the contact line moves and the droplet deforms to respect h s . Considering the density ratio between the liquid and the surrounding gas, the final shape of the droplet depends on two parameters, namely h s and the following Eötvös number: Fig. 4. Schematic illustration of stick-slip model (Chen et al., 2013).
The physical properties used in the simulations are summarized in Table 1. The simulations are carried out using grid size D ¼ 5 Â 10 À4 m and the typical time step used is Dt ¼ 1 Â 10 À4 s. No-slip boundary conditions are applied at all computational domain boundaries.

Absence of gravity
First, simulations are performed for the case where gravity is absent (Eo ¼ 0). In this case the equilibrium shape is a spherical cap. By conservation of the volume of the droplet, it is possible to geometrically derive at equilibrium the radius of curvature R f , the contact radius r f , the height of the droplet h f and the pressure difference DP between the drop and the surrounding gas as follows: Fig. 6 shows an excellent agreement between the numerically obtained contact radius and droplet height with Eqs. (14) and (15), respectively. The numerical and analytical equilibrium droplet shapes compared in Fig. 7, also are in good agreement. The expected value of pressure difference DP analytical is also compared with the numerical pressure difference DP numerical (Fig. 8). Here, DP numerical is Table 1 Fluid physical properties for simulations of equilibrium shape of a droplet on a horizontal surface.

Spreading effect due to gravity
When there is gravity in the simulation, the gravity will counteract the surface tension. At low Eötvös numbers (Eo ( 1), the droplet maintains its original spherical cap shape as the surface tension forces are dominant and the droplet height in this case is given by Eq. (15). At higher Eötvös numbers (Eo ) 1), the droplet spreads under the gravity effect and the thickness is directly proportional to the capillary length (Dupont and Legendre, 2010): Fig . 9 shows the equilibrium droplet shapes and dimensionless droplet height for a wide range of Eötvös numbers and a static contact angle of 90°. The figure clearly shows that the obtained results agree well with the two asymptotic solutions given by Eqs. (15) and (17). As expected, for the intermediate values of Eötvös number, the transition between a spherical cap and a puddle shape can be observed.

Impact and spreading of droplet on a horizontal surface
To validate the implementation of the dynamic contact angle model, we simulate the dynamics of a 3D glycerin droplet impacting on different flat surfaces, for which experimental measurements are available (Šikalo et al., 2005). Physical properties used for the simulations are summarized in Table 2. Schematic diagrams of the initial and equilibrium droplet states are presented in Fig. 10.
The simulations are carried out for a 2.45 mm diameter glycerin droplet impacting at 1.41 m/s on a horizontal wax and glass surfaces with equilibrium contact angles equal to 93.5°and 15°, respectively. Only Kistler's empirical correlation (Model-B) is used in the simulations because no distinct advancing and receding angles are considered. We first test the convergence of the simulation as we increase the grid resolution. The test is performed using Table 2 Fluid physical properties and wettability of the surfaces for simulations of impact and spreading of a droplet on a flat surface (Šikalo et al., 2005 ig. 10. Schematic representation of the initial (-) and equilibrium (--) shapes of a droplet impact on a flat surface.  five different grid resolutions, D = 2.0 Â 10 À4 m, 1.5 Â 10 À4 m, 1.0 Â 10 À4 m, 7.0 Â 10 À5 m, and 5.0 Â 10 À5 m. The effect of grid resolution on the time evolution of the dimensionless contact radius is shown in Fig. 11. Here, the dimensionless contact radius is defined as the radius of the wetted spot r f normalized by the equivalent droplet radius R 0 . It can be seen that the solution converges with increasing grid resolution. In this section, grid size D = 7 Â 10 À5 m and time step Dt = 5 Â 10 À6 s, for which the results are grid independent, are used for the comparison of the present numerical results with reference results.
The experimental and both literature and LFRM numerical results for the dimensionless contact radius are shown in Fig. 12 (a) and (b) for droplet spreading on wax surface and glass surface, respectively. It can be observed that the obtained result for the case of droplet impact on wax surface shows a very satisfactory match with the results reported in literature. In addition, the dynamic contact angle model performs better than the static contact angle model in the receding phase but there is no significant difference in the advancing phase. However, Fig. 12 (b) shows that both the present simulation and VOF simulation, which also uses Kistler's empirical correlation, over-predict the spreading rate when the droplet approaches the equilibrium in case of droplet impact on a glass surface. This may indicate that the empirical correlation used in the simulations is less accurate in describing the receding phase of droplet on hydrophilic surface.

Bubble formation from a submerged orifice
In this section we present 3D numerical simulations of bubble formation from a submerged orifice. The schematic of the computational domain is given in Fig. 13. The air bubble is injected through an orifice of radius R o submerged in initially quiescent water. The circular orifice is located at the bottom center of the numerical domain and represented with a staircase approximation. The gas injection is assumed to be under constant flow rate Q. The flow in the gas inlet is assumed fully developed laminar and a parabolic inflow velocity is imposed. At the top of the simulation domain, the pressure-prescribed outlet boundary condition is imposed and the outlet pressure is set to atmospheric pressure, while at the side and lower walls, the no-slip boundary condition is used. At the lower wall, a dynamic contact angle is used for the simulation of bubble formation with a moving contact line. When the contact line coincides with the orifice rim, the apparent contact angle is allowed to have any value between h rec and 180   (h rec 6 h D 6 180 ), ensuring that the contact angle will start to spread outwards when h D < h rec and preventing the contact line to recede into the orifice, respectively. For the simulations of bubble formation without a moving contact line, h rec is set to 10°to ensure that the bubble will remain at the orifice rim. The gravitational acceleration is imposed in the vertical direction. The numerical domain has a width and height of five equivalent bubble diameters. The domain width is similar to that of Gerlach et al. (2007) and  who have reported that the bubble growth is not influenced by any liquid circulations close to the wall. At the initial time step, a semi-circular bubble (h D = 90°) is positioned at the inlet with an initial radius equal to R o .

Bubble formation with a pinned contact line
To validate the numerical model, the results of the model with pinned contact line are compared to the experiments on quasistatic bubble growth of . In the experiments, the gas phase is injected through the orifice using a small and constant volumetric flow rate Q, so that the bubble is subjected only to the static forces, which includes buoyancy, surface tension, viscous drag, and inertial forces, during its growth. The physical properties that are used for the simulations, imitating the corresponding experimental conditions, are summarized in Table 3.
We first examine the effect of the grid resolution on simulation results by simulating the growth of bubble under volumetric flow rate Q ¼ 5:56 Â 10 À8 m 3 =s. The grid resolution tested are D = 2.5 Â 10 À4 m, 2.0 Â 10 À4 m, and 1.5 Â 10 À4 m. Fig. 14 shows the effect of grid resolution on bubble shape at the final instant prior to bubble pinch-off. It can be seen that the bubble shapes at the detachment for all grid resolutions are almost identical. In addition, the differences in detachment volume V det and detachment time t det are also computed. The differences in V det and t det for D = 2.5 Â 10 À4 m are less than 0.4% and for D = 2.0 Â 10 À4 m are less than 0.2% compared Fig. 14. Effect of grid resolution on bubble shape at t=t det ¼ 1 for Q ¼ 5:56 Â 10 À8 m 3 =s and Ro ¼ 0:8 mm. Fig. 15. Bubble shape predictions at six time frames t=t det $ 0, 0.2, 0.4, 0.6, 0.8, 1 for experimental  and present numerical simulation: (a) Q = 4:17 Â 10 À8 m 3 =s and (b) Q = 5:56 Â 10 À8 m 3 =s. The orifice radius is 0.8 mm.

Table 4
Comparison of time detachment t det and bubble volume at detachment V det for four different flow rates .  Table 5 Fluid physical properties and initial conditions for simulation of bubble growth of Thoroddsen et al. (2007). to the results with a grid resolution of D = 1.5 Â 10 À4 m. Based on the considerations of accuracy and computational time, the grid resolution of D = 2.0 Â 10 À4 m is employed in remainder of this section. The time step is Dt ¼ 5 Â 10 À5 s. Next, we qualitatively assess the present numerical method by comparing the temporal evolution of the bubble shape against the experimental observations (Fig. 15). The six interface shapes shown here correspond to t=t det $ 0, 0.2, 0.4, 0.6, 0.8, 1, where t=t det = 1 is the detachment time predicted by the present method or measured experimentally and t=t det = 0 is the initial time of the simulation where the bubble is assumed to have a hemispherical shape. The spatial evolution of the numerically predicted interface is very close to the corresponding experimental results. Table 4 shows a good quantitative agreement of the detachment time and volume, with maximum errors of 2.6% and 3.3% for t det and V det , respectively.  (Thoroddsen et al., 2007) and numerical (present study) interface evolution of the generated bubble, before the time of detachment (t = 0 s). The orifice radius is 1.35 mm. Fig. 17. Comparison of experimental (Thoroddsen et al., 2007) and numerical (present study) interface evolution of the generated bubble, after the time of detachment (t = 0 s). The orifice radius is 1.35 mm.

Flow rate Experimental Present study Error
Since the work of  only shows the bubble dynamics until the bubble detachment, we further validate the present numerical method against experimental results of Thoroddsen et al. (2007). In the experiments of Thoroddsen et al., 2007, the gas bubble is slowly grown from a vertical thin-walled nozzle. Although a difference in the bubble volume between the detachment from a nozzle compared to an orifice can be expected, one can compare their result with the case where the bubble is pinned at the orifice since the liquid circulation near the nozzle is expected to be negligible due to very low gas flow rate used in the experiment (Q ¼ 1:67 Â 10 À8 m 3 =s). The same simulation characteristics as well as the same computational domain and mesh are used, as in the previously presented validation case. The only difference is the size of the orifice and the fluids physical properties, which are summarized in Table 5.
In Figs. 16 and 17, a quantitative comparison is made by extracting the point coordinates of the interface position with respect to time, at several time instances before (Fig. 16) and after (Fig. 17) the detachment of the bubble from the orifice. The numerical results are in excellent agreement with the experimental results, considering the interface shape evolution of the generated bubble, both before and after the bubble detachment from the orifice. In the numerical results, an air thread forms as the bubble neck contracts, leaving a central micro bubble after the collapsing process. Although this phenomenon was also found in the experiment, the size of the micro bubble is more pronounced in the simulation due to the low interface resolution.

Bubble formation with a moving contact line
In this section, we present simulations of bubble formation with a moving contact line and compare the results to the experimental results of Chen et al. (2013). Chen et al. (2013) carried out experiments on bubble formation from a paraffin coated stainless-steel (hydrophobic) orifice and studied the macroscopic contact line behaviors during the growth. In addition, Chen et al. (2013) also performed numerical simulations using a 2D-LS method with two kinds of contact angle models, viz. contact line velocity dependent model and stick-slip model. Because Chen et al. (2013) did not observe an explicit relationship between the contact angle and the contact line velocity in their experiments, their simulation results obtained using stick-slip model show better agreement with the experiments. The fluids physical properties used in the present simulations are summarized in Table 3. The orifice radius R o is 1.5 mm and the orifice flow rate Q is 6.79 cm 3 /s. The prescribed receding and advancing contact angles for the simulations are h rec ¼ 90 and h adv ¼ 110 , respectively. These values are taken from the contact angle hysteresis observed in the experiment.
We again test the convergence of the simulation as we increase the grid resolution. The test is performed for Model-A (stick-slip) and Model-B (velocity dependent model) using five different grid resolutions, D 1 = 3.0 Â 10 À4 m, D 2 = 2.5 Â 10 À4 m, D 3 = 2.0 Â 10 À4 m, D 4 = 1.5 Â 10 À4 m, and D 5 = 1.25 Â 10 À4 m. We assume that the contact line remains pinned when the contact angle value changes from h rec to h adv during the advancing phase. The effect of grid resolution on the evolution of contact line radius is shown in Fig. 18. It can be seen that the solution converges with increasing grid resolution for both contact angle models. The L 2 relative error norm (defined with respect to the finest grid) is used to quantify the difference as follows: where r stands for the radius of the contact line, n is the index to represent nth sampling point and N is total number of sampling points. For Model-A, calculated L 2 for grids D 1 ; D 2 ; D 3 , and D 4 are 1:25%, 0:67%, 0:35%, and 0:15%, respectively. For Model-B, calculated L 2 for grids D 1 ; D 2 ; D 3 , and D 4 are 0:93%, 0:68%, 0:37%, and 0:15%, respectively. These indicate that the grid D 2 is good enough to produce accurate results. First, we compare the temporal evolution of the contact angle obtained using both contact angle models to the reference results. As can be seen from Fig. 19, the variations of the contact angle agree reasonably well with the experiments. However, larger contact angles are seen within the initial 5 ms experimentally. Chen et al. (2013) reported that this was attributed to the measurement errors in the contact angle, because the contact line might initially be located inside the orifice. In all simulations, the initial contact angle is 90°since the bubble shape is assumed to be a semicircle. The figure also shows that the variation of the contact angle obtained using Model-B is relatively similar to Model-A. This is due to the fact that the motion of the contact line in the present simulation is relatively slow, with maximum Capillary number in the order of 10 À3 .
Next, we compare the temporal evolution of the contact line radius observed in the experiment with the simulations in Fig. 20. It can be seen that the predicted contact line radius of Model-A is similar to Model-B since both models predicted similar contact angle variations. Furthermore, it is also evident that both LFRM and LS method over-predict the contact line radius. However, the present simulations obtained using dynamic contact angle models are slightly better in describing the maximum contact line radius compared to the LS method. Despite this fact, the bubble departure time obtained by these simulations agree well with the experiments. When a static contact angle is employed, the spreading length immediately shrinks towards the orifice rim once it reaches the maximum value, while in the experiment the bubble base sticks on the surface for a relatively long period. As a consequence, the predicted departure time is faster compared Fig. 19. Comparison of the temporal evolution of contact angle for bubble formation on the paraffin-coated stainless-steel surface at an orifice flow rate of 6.79 cm 3 /s and orifice radius of 1.5 mm (Chen et al., 2013). Fig. 20. Comparison of the temporal evolution of contact line radius for bubble formation on the paraffin-coated stainless-steel surface at an orifice flow rate of 6.79 cm 3 /s and orifice radius of 1.5 mm (Chen et al., 2013).
to the experiment. This indicates that dynamic contact angle model is essential for obtaining good agreement with the experiments.
We also make a comparison of experimental and numerical bubble shape evolution in Fig. 21. This figure clearly shows that bubble shapes obtained using both dynamic contact angle models are almost identical and they agree relatively well with the reference numerical results. However, the experimental bubble shapes are more elongated in the vertical direction compared to all numerical results, which may be attributed to the slower contact line movement in the experiment.

Bubble formation under various gas flow rates
In this section, we present simulations of bubble formation under various gas injection rates, ranging from 1 ml/min to 400 ml/min, in order to test the accuracy and robustness of the present numerical method. The orifice Reynolds number, defined as Re O ¼ 2Q q g =ðpR o l g Þ, is usually employed in the literature to characterize the orifice flow rate. In the present work, the values are in the range of 0:73 < Re O < 290. The orifice radius R o of 1 mm is used and the triple contact line of the bubble is pinned at the orifice rim. The simulations are performed with the fluid properties presented in Table 3. Two different grid sizes D = 1.5 Â 10 À4 m and D = 2.0 Â 10 À4 m were tested for the case of Re O ¼ 290, where the difference in the detached bubble volume is around 3%. We expect that the error will be less at lower gas flow rates. Therefore, grid resolution of D = 2.0 Â 10 À4 m is employed in this section and the typical time step is Dt ¼ 1 Â 10 À5 s.
The influence of gas flow rate on the formation time and the detached bubble volume is shown in Fig. 22. An important feature of the data in Fig. 22 is that for very low flow rates or static conditions, the bubble volume is almost independent of gas flow rate, i.e. the bubble formation time increases with decreasing flow rate. The reason is that at low flow rates, capillary and buoyancy forces govern the bubble volume at detachment. Excellent agreement of the present simulation results with the simulation results of Gerlach et al. (2007) in the case of low gas flow rates can be seen in the figure. However, there is a significant difference in the case of Q = 150 ml/min. Gerlach et al. (2007) predicted critical bubble volume and time period, at which the transition from single bubbling regime (static regime) to double pairing regime (dynamic regime) takes place for several contact angles. This prediction is also shown in Fig. 22 (b) (red dashed line). They reported that in the case where the bubble base is pinned to the orifice, the transition occurs at Q = 150 ml/min. For Q P150 ml/min, the wake of the departing bubble affects the formation time and bubble volume of the following bubble. Therefore, the discrepancies in the bubble volume and formation time may be attributed to the fact that the present work consider bubble detachment characteristics of the first bubble, whereas the work of Gerlach et al. (2007) considered the growth of a continuous chain of bubbles. At higher gas flow rates, the formation time is almost independent of the flow rate, but the bubble volume increases significantly with Q. For liquids of low viscosities Fig. 22. Influence of the gas flow rate on the formation time (a) and bubble volume (b) for orifice radius 1 mm (Gerlach et al., 2007;Davidson and Schüler, 1960a). Fig. 21. Comparison of the evolution of bubble shapes for bubble formation on the paraffin-coated stainless-steel surface at an orifice flow rate of 6.79 cm 3 /s and orifice radius of 1.5 mm (Chen et al., 2013). and for high flow rates, the bubble volume can be predicted by the one stage model of Davidson and Schüler (1960a) for an inviscid liquid given by the following equation: Fig. 22 shows that the obtained results are converging towards this equation.
The temporal evolution of the interface shape under various orifice flow rates Q is shown in Fig. 23. At low gas flow rates, the bubble grows quasi-statically and takes on roughly the static pendant shape, which results from a local balance of capillary and hydrostatic pressures (Q = 50 ml/min). When the bubble volume grows beyond a maximum value, the buoyancy force exceeds the anchoring surface tension force, and the bubble pinches off. At high gas flow rates, viscous and inertial forces distort an emerging bubble from the static shape. In addition, the wake of the departing bubble influences the growth of the emerging bubble (Q P 150 ml/min). Finally, a liquid jet phenomenon is found in the simulations for high gas flow rate (Q = 350 ml/min). Under higher flow rates, a jet of liquid is formed following the breaking of the bubble neck. This liquid jet directs upwards into the bubble and it can penetrate the bubble top surface under certain conditions. The formation of liquid jets could be mainly attributed to large surface tension forces right after the separation of the bubble neck. This interesting phenomenon has been reported in the experimental work of Mitrovic (1997).

Conclusions
In the present work, we extend the Local Front Reconstruction Method (LFRM), a front-tracking method that enables interface merging and breakup, by incorporating the contact angle dynamics. The implementation of the contact angle dynamics are verified and validated by performing simulations of droplet impact and spreading. The obtained results using the improved method show a satisfactory match with analytical solutions, as well as experimental and numerical data available in the literatures. The numerical method is then applied for the simulations of bubble formation from a submerged orifice. The bubble detachment characteristics predicted by the method, are extensively compared with available experimental and numerical results, showing good agreement. In the case of bubble formation with a moving contact line, the LFRM with dynamic contact angle models are better in predicting the bubble departure time and the temporal evolution of spreading length in comparison with a static contact angle model.
The problem of bubble formation under various volumetric flow rates is also considered to assess the accuracy and robustness of the numerical method. It is demonstrated that the present method is able to accurately simulate bubble formation both under low gas flow rates, where the accuracy of surface tension calculation is important, and under high gas flow rates, where the interface deformation becomes complex due to the interplay among inertial, viscous, surface tension, and buoyancy forces. The present work provides a framework to further investigate the effects of various operating conditions (e.g., gas flow rate, cross-flow), fluid physical properties, and also system properties (e.g., orifice dimensions, wettability of the orifice) on the dynamics of bubble formation.