Complex bubble deformation and break-up dynamics studies using interface capturing approach

The dynamics of bubble deformation has significant impacts on two-phase flow fundamentals such as bubble induced turbulence and flow regime transition. Despite the significant progress achieved by experimental studies on bubble deformation, certain limitations still exist especially for wide-range datasets. To significantly expand the flow conditions available from experiments, direct numerical simulation (DNS) is utilized to study the bubble-liquid interactions using finite-element solver with level-set interface capturing method. Different from conventional investigations of bubble rising and deforming in stagnant liquids, a proportional-integral-derivative (PID) bubble controller is leveraged to maintain the bubble location in uniform liquid flow. This paper evaluates the reliability and reproducibility of the PID bubble controller for complex bubble deformation studies through a comprehensive set of verification and validation studies. An improved bubble deformation map is developed, based on Weber number and bubble Reynolds number, showing six zones for different deformation and break-up mechanisms. This research aims at producing virtual experiment level data source using interface resolved DNS and shedding light into the physics of interface dynamics. The insights obtained can be further incorporated in improved multiphase CFD models to guide the engineering designs and industrial processes where bubble deformation and break-up play a pivotal role.


Introduction
Bubble deformation is a ubiquitous phenomenon in gasliquid two-phase flow and has significant impacts on various interface related flow physics. For example, bubble induced turbulence has distinct characteristics from particle induced turbulence. The existence of bubbles can either augment or suppress local liquid turbulence given specific bubble deformability (Gore and Crowe, 1989); in other words, the turbulent intensity change depends on the bubble distortion level (Feng and Bolotnov, 2017a). In addition, the bubble deformation also plays an important role in determining bubble interfacial forces, such as drag force, lift force, and virtual mass force (Bhaga and Weber, 1981). Severe deformation may result in bubble break-up, which leads to more complex flow regimes. One such example is the slug bubble during the slug-to-churn flow regime transition (Zimmer and Bolotnov, 2019). As the gas flow rate increases, satellite bubbles are teared off from the slug bubble which contribute to the interface instability and eventually trigger the transition. Fundamental studies on bubble deformation and break-up could help reveal the underneath mechanisms and produce improved predictive models. This in turn will lead to optimized engineering designs (e.g., nuclear, thermal, and chemical reactors) and industrial processes (e.g., purification of polluted water, and the charging of plasma bubbles) where bubble deformation and break-up play a pivotal role (Yamatake et al., 2007;Guillen et al., 2018).
A considerable amount of literature has been published on the experimental studies concerning bubble deformation (Bhaga and Weber, 1981;Sharaf et al., 2017;Dueñas, 2019). Bhaga and Weber (1981) injected the air bubble in aqueous sugar solution and the various liquid viscosities were obtained by solutions of different concentrations. The photos of bubble deformation and the streamline in the wake region have become classical references for two-phase flow model development and validation (Sussman et al., 2007;Tripathi et al., 2015). As an extension of the numerical study conducted by Tripathi et al. (2015), Sharaf et al. (2017) performed an experimental investigation with a high-speed camera and a satisfying agreement was reported between experiments and the corresponding simulations. The measurement techniques are also evolving for two-phase flow experiments, and Dueñas (2019) utilized particle image velocimetry (PIV) to track the velocity of rising ellipsoidal bubbles, and developed an improved drag coefficient model. Despite the significant progress achieved by experiments, certain limitations persist. Examples include the availability of working fluid properties, fine control of experimental conditions as well as measurement uncertainties.
Although Bhaga and Weber's experiments covered an adequate range of flow conditions (e.g., Eotvos number, Morton number, and Reynolds number) to produce a bubble deformation map, the parametric exploration is primarily achieved by altering the liquid viscosity with varying concentration of the solutions. With almost unchanged surface tension, liquid density, and gas parameters, they were not able to capture the bubble break-up phenomenon which is more sensitive to surface tension (Bhaga and Weber, 1981). The experimental conditions are hard to control at constant values in two-phase studies. To simplify the control of fluids, all of the aforementioned experiments require liquid to be stagnant when bubble rises up, which limits the research scope into purely gravity-driven. However, it is rarely the case in engineering applications. As bubbles and liquid are moving together, the inertia of the bubble could be sometimes more dominant than gravity (Ishii et al., 2004). Although experimental data remains the major knowledge source on two-phase flow physics, the uncertainty quantification (UQ) of measurements is still challenging and the UQ analysis is often lacked (especially in legacy experimental studies). Bubble columns are usually built with fixed height, and thus for fast-rising bubbles it is hard to determine if bubbles reached steady state at the observation section. Consequently, the measured bubble terminal velocity generally has higher uncertainty than other experimental parameters. Furthermore, accurate three-dimensional (3D) bubble shape is extremely difficult to obtain in experiments. Cameras could only take two-dimensional (2D) photos, and thus most studies have to assume bubbles to be axisymmetric which is only applicable for simply deformed bubbles (Dueñas, 2019). Even with tomography technique to reconstruct the 3D bubble shape, the accuracy still depends on the rotation speed of camera, the strength of received signal from the sensors, and the selection of back projection algorithms (Utomo et al., 2001).
Meanwhile, numerical simulations have grown into an attractive investigation approach to compensate the shortcomings of experiments. Any working fluid and physical condition could be readily represented by changing material properties and initial/boundary conditions. In particular, 3D interface resolved simulations can provide all the quantities of interest in bubble dynamics studies, such as bubble topology changes, velocity and pressure distributions, flow field streamlines, etc. Sussman et al. (2007) reproduced nearly identical deformable bubbles in the experiments with a high-order interface capturing method. Tripathi et al. (2015) computationally investigated bubble deformation and break-up with wide-ranging fluid properties. Their simulations revealed that bubble topology is not only affected by fluid properties, but also strongly dependent on bubble sizes. Those studies were all carried out with direct numerical simulation (DNS). DNS is a first-principle based simulation approach, which directly solves Navier-Stokes equations (NSE) with adequate spatial and temporal resolutions to capture all flow scales of interest. Coupled with interface capturing/tracking approaches, DNS can be an ideal tool to study the bubble liquid interactions (sometimes referred as "virtual experiment"). For example, in complex bubble deformation studies, DNS is able to capture the extremely thin film and even the satellite bubbles after break-up.
With the rapid development of high-performance computing, DNS is becoming increasingly affordable to study the complex two-phase phenomenon with high fidelity. As the DNS flow solver adopted in the current study, PHASTA is developed for Parallel, Hierarchic, higher-order accurate, Adaptive, Stabilized, Transient Analysis (Whiting and Jansen, 2001). It is a finite element code with level-set method implemented to track the topological evolution of gas-liquid interface (Sussman et al., 1998). Two important features have made PHASTA particularly suitable to investigate the bubble deformation and break-up dynamics: a novel proportionalintegral-derivative (PID) bubble controller and the support of unstructured meshes. Instead of having bubbles rise in a stagnant liquid tank, the PID bubble controller would control the bubble's position in moving liquid with or without local shear rate (Thomas et al., 2015;Feng and Bolotnov, 2017b). The bubble will be anchored at a fixed location under steady state to facilitate the bubble deformation study and interfacial force quantification (based on a force balance assumption). It is reported by previous DNS and experimental studies that some bubbles with severe deformation also have complicated travel trajectories (Tripathi et al., 2015;Sharaf et al., 2017). The trajectories vary with fluid conditions thus determining the size of computation domain and the mesh resolution is a process of trial-and-error. A more common practice is to refine all the regions where bubbles could possibly travel, which may lead to exceedingly large computational costs. However, with the PID bubble controller, bubbles would only vibrate initially in the vicinity of and eventually settle down at the prescribed control location, which significantly reduces the size of mesh refinement zone and thus the computational costs. Moreover, due to the bubble movement, a conventional investigation would have to realign the bubble center before comparing the relative interface topology change at different time stages (Tripathi et al., 2015;Sharaf et al., 2017). With the PID bubble controller, this burdensome task can be totally avoided because the bubble has a fixed steady-state location. With the support of unstructured mesh, an optimized computational mesh can be easily created with additional mesh refinement only around the bubble control location.
With the advantages of PID bubble controller and unstructured mesh in PHASTA, only moderate bubble deformation under PID control and complex bubble deformation without using PID bubble controller were studied (Feng and Bolotnov, 2017b;Guillen et al., 2018). It is necessary to prove the capability and reproductivity of PHASTA on PID-controlled complex bubble deformation, which will shed light into the physics of interface dynamics including the flow regime transition, bubble induced turbulence, and interfacial force closure development. The application of PID bubble controller and unstructured mesh for complex bubble deformation study will be validated in this paper.
The rest of the paper is structured as follows: the numerical methods are to be introduced in Section 2 including the basics of PHASTA code, the level-set method, and the PID bubble controller; simulation setups are presented in Section 3; and the result discussions (Section 4) consist of the grid convergence study of unstructured mesh, validation and verification of the current interface capturing approach in studying the bubble deformation and break-up dynamics, and an improved bubble deformation map based on Weber number and bubble Reynolds number. As a fundamental numerical study, this research aims at producing virtual experiment level data source and shedding light onto the underneath physics of bubble deformation and break-up dynamics.

Flow solver and interface capturing approach
As a multiphase flow solver, PHASTA can solve both incompressible and compressible flows governed by Navier-Stokes equations (Whiting and Jansen, 2001). In the bubble deformation studies, incompressible Navier-Stokes equations (INS) are solved in PHASTA, and the "one-fluid" approach is used to describe the two-phase flow fields (Prosperetti and Tryggvason, 2007).
where , , ρ p u , and μ are the velocity, density, pressure, and viscosity of either liquid or gas. ρg represents the most common body force, gravity force. For problems where electric field is resolved, the relevant forces can be added as another source term (Wang et al., 2018). A finite interface thickness is assumed and the surface tension force x ( ) γκδ n is added in the interface region as an additional body force. γ is the surface tension coefficient, κ is the curvature, x is the interface location, and n is the unit normal vector at the interface.
Level-set method is a well-known interface capturing approach and was implemented in PHASTA for two-phase computations (Sussman et al., 1998;Nagrath et al., 2005;Nagrath et al., 2006). In the level-set distance field, the interface is modeled by the zero level-set 0 =  . The gas phase and liquid phase have negative and positive  correspondingly and the magnitude of  is the distance to the nearest interface. To improve the numerical robustness, the transitions of fluid density and viscosity are described by a smoothed Heaviside function where g l 1 , , , ρ ρ μ and g μ are the densities and viscosities for liquid and gas phases respectively and  is the level-set value. The solving of level-set field in each time iteration consists of the advection of  and then a redistancing process.
PHASTA is a proven tool to simulate the two-phase flows with highly distorted interface. Behafarid et al. (2015) simulated the cap bubble rising in a narrow near-vertical channel and excellent agreement was noticed when comparing with the experimental data. Code-to-code verification was also successfully conducted between PHASTA and STAR-CCM+ on deformed bubbles in waste glass melter (Guillen et al., 2018). The adaptive mesh refinement capabilities used in PHASTA were also verified and it was shown to successfully capture wavy interface structures of annular flow as well as reduce the computational costs (Rodriguez et al., 2013). More recently, a slug-to-churn flow regime transition study reported favorable agreement between PHASTA results and the existing experimental and analytical results for complex two-phase flow regimes (Zimmer and Bolotnov, 2019). The two-phase flow simulation capability can also be extended to even larger engineering problems, such as nuclear reactor subchannels (Fang et al., 2017(Fang et al., , 2018(Fang et al., , 2020.

Bubble control approach
The proportional-integral-derivative (PID) controller was implemented in PHASTA to switch the reference system from a container-based frame to a bubble-based frame of reference (Thomas et al., 2015). The container-based frame is stationary, while the bubble-based frame moves at the same speed as the bubble, and thus the bubble in the latter always has a zero velocity, equivalent to being controlled at a fixed position. Figure 1 illustrates the bubble and liquid velocity conversion as the reference system is switched. In the container-based frame, the rising bubble has a terminal velocity g t ( ) u u = at steady state in the stagnant liquid l ( 0); u = in the bubble-based frame, bubble velocity is 0 g ( 0) u = and the surrounding liquid will be assigned a velocity of l t u u = -. Simple math shows that the relative velocities in the two frames are both t u as long as steady state is achieved: Container-based reference frame: Bubble-based reference frame: Under Galilean transformation theory, as long as the relative velocity is maintained identical, all the physics involved (including mass, momentum, and energy) are invariant in the two inertial frames (McComb, 1999). Therefore, bubble motions in the container-based and bubble-based frames are invariant with equivalent bubble physics as the bubble rises in the liquid. In addition, the interfacial forces that the bubble experiences are also invariant in two systems, so the PID controller can further make use of the steady state force balance to infer the interfacial forces.
Similar to particles, at steady state a bubble in a laminar flow only experiences drag force and lift force if local shear flow exists (Ejeh et al., 2020). The existence of the walls will result in wall forces on the bubble, but if the bubble is at the centerline of the domain, wall forces pointing from the two sides to the bubble are opposite thus will negate each other. Therefore, the PID bubble controller will only need to balance drag force D ( ) F and lift force L ( ) F on bubble by streamwise and lateral control forces, cf x F and cf y F (Fig. 2). The PID control model has three components: the proportional, integral, and derivative terms respectively. In bubble control scenario, the proportional component takes  into consideration the bubble location change with respect to the initial location. Next, the integral component averages the historical control forces. Finally, the derivative component will conduct derivatives on the bubble location and bubble velocity, which are velocity and acceleration respectively. The PID bubble controller was previously implemented in PHASTA with the following control force equations (Thomas et al., 2015): where ( ) 1 cf n i F + is the control force at the new time step (all terms with ( ) n are older time steps) and cf i F is the control force along the i direction. In the current study, control forces are enabled in all 3 directions ( , , ) and derivative ( i v and d i v ) component. When the simulation reaches a steady state, the control force can be extracted as the force counteracting against interfacial forces. For bubble rising in liquid, the counterforce of drag force D F is buoyancy force B F , and thus the streamwise control force in bubble-based reference system is equivalent to the buoyancy force in container-based reference system.
Container-based reference frame: Bubble-based reference frame: where g is the earth gravity, g V is the volume of the bubble, D C is the drag coefficient, and b A is the cross sectional area of the bubble. The drag coefficient can be then computed accordingly.
PID bubble controller coupled with unstructured mesh provides significant advantages in bubble deformation and break-up studies. Since the bubble's location is fixed, the locally refined unstructured mesh around the interface is sufficient to resolve the deformed topology, simplifying the overall geometry and eliminating the need for adaptive meshes. Also, a relatively small computation domain can be utilized for all deformation and break-up regimes for each bubble size. In addition, postprocessing the bubble topology evolution does not require additional procedures in experiments or other DNS to align bubble centers from different trajectories. With the mass center of bubble maintained at one location using PID controller, the bubble deformation relative to the initial shape can be observed and compared directly.
In addition, the bubble terminal velocity generally has the largest uncertainty in experiments compared to other parameters. With this reverse force determination approach leveraging the PID controller, the uncertainty of bubble terminal velocity measurement can be eliminated since it becomes an input. The liquid will maintain this relative velocity once a steady state is achieved, improving the overall accuracy on bubble topology prediction. The accurate input of the terminal velocity will also enhance the reproductivity of deformation studies, with more precise classification of regimes for different bubble deformations. As the mechanisms of two-phase interaction are more thoroughly resolved, a detailed bubble deformation map can be built and utilized in scientific and engineering applications.
PID bubble controller in PHASTA has been verified and validated on the interfacial force evaluation in laminar and turbulent flows and achieved reliable agreement with experimental data (Thomas et al., 2015;Feng and Bolotnov, 2017b). This paper will attempt to justify the reliability of the PID bubble controller for complicated interface distortions.

Spatial discretization
The computational fluid domain discretized by unstructured mesh with the deformed bubble under PID control is illustrated by Fig. 3. x axis is the streamwise direction, y axis is the lateral direction, and z axis is the spanwise direction.
The fluid has a uniform inflow velocity distribution towards the bubble in the controlled environment, equivalent to the bubble rising in stagnant liquid. The computation domain size and the distance from bubble center to the inlet plane is proportional to the initial bubble radius R . The edge of the cube-shaped computation domain has a length of 10R and the bubble is located at in all the simulations. From the surrounding liquid to the bubble, the mesh is gradually refined. Since the bubble shape will most likely expand against the inflow, refinement regions are wider on lateral and spanwise directions ( y z plane) than the streamwise direction ( x axis). Four different refinement zones are designed to better capture the interface distortion. Cylindrical refinement approach is utilized over box refinement as it better approximates the bubble shape while also reducing the computational cost. The near bubble wake region has the finest resolution due to the potential of forming thin film and satellite bubbles after complex deformation and break-up.

Simulation parameters
Simulation parameters were selected based on two principles: the need to match dimensionless numbers provided in published studies for topology reproducibility (Bhaga and Weber, 1981;Tripathi et al., 2015;Sharaf et al., 2017) and the need to expand the range of each parameter for classification of deformation regimes through the variation of bubble topology.
To reproduce the bubble deformation in the classical experiment, the bubble size was selected to be 13.0 mm , and the viscosities were designed to yield identical Eotvos number ( ) Eo , Morton number ( ) Mo , and bubble Reynolds number ( b Re ) (Bhaga and Weber, 1981). As for the reproduction of the bubble deformation and break-up in the more recent numerical and experimental studies, the bubble radii were selected to be 2.0 mm and 20.0 mm , and the ratio of liquid and gas densities is fixed at l g : 1000 ρ ρ = and l g : 100 μ μ = (Tripathi et al., 2015;Sharaf et al., 2017). Some key parameters of bubble and liquid properties are summarized in Table 1.  10, 0.16, 0.17, 0.198, 0.218, 0.238, 0.263, 0.27, 0.80, 0.83, 0.92

Results
The results of PID controlled bubbles in uniform liquid flow simulations are equivalent to those of rising bubbles in stagnant liquids driven by buoyancy. In the beginning of this section, a grid convergence study is performed to quantify the uncertainty from mesh discretization (Roache, 1998). The film size of the deformed bubble in each simulation is measured and compared to support the choice of mesh configuration. Followed by a verification and validation (V&V) of the presented approach solutions, with published results, the capability and precision of the PID bubble controller are assessed in revealing complex bubble deformation and break-up behaviors. The bubble deformation and breakup mechanisms are classified into six different zones on a bubble deformation map and the corresponding dynamics of each zone is discussed. The impact of different parameters on bubble dynamics is elaborated as well to shed light onto the underneath physics of deformation and break-up phenomena.

Grid convergence study
Five mesh configurations depicted in Figs. 4(a)-4(e) are selected to conduct the grid convergence study. The element numbers across initial bubble diameter are 34, 39, 45, 52, and 60 in each configuration, with a refinement ratio of 1.2 r = . The bubble has a radius of 20 mm R = and the relative velocity between the bubble and liquid is r u = 0.8 m / s . Severe deformation will happen at and l g : μ μ = 100, and the wake side of the bubble will form a thin film ( Fig. 3(a) and Fig. 4(f)). It has traditionally been a challenge for simulations to resolve thin films but by utilizing levelset method in PHASTA it is possible to capture such severe interface distortions.
A grid convergence study can be used to estimate the uncertainty of spatial discretization via the well-developed theory of grid convergence index (GCI) (Roache, 1998). As GCI reduces to 0, the numerical result theoretically approaches to the exact solution. For a group of meshes with a uniform refinement ratio r , GCI can be computed using a factor of safety s ( ) F , the results from each mesh ( i f ), and the order of convergence ( ) p .
where s 1.25 F = when more than three meshes were involved. i f and 1 i f + are results of finer and coarser meshes. The peripheral size of the thin film is a crucial parameter to describe deformability. The film has a ring shape on y z plane (Fig. 4(f)). The larger the ring, the more severe the bubble deformation is with respect to the initial state. Therefore, the outer radius of the film at the same x location is selected as the GCI metric in the grid convergence study and is computed with different elements across the initial bubble diameter (#/D). From Fig. 5, GCI curve is asymptotic to 0, which indicates that the numerical solution gradually approaches to the exact solution. The bubble deformation will be the most accurate if using 60 #/D. However, Table 2 shows that the difference in GCI of 52 #/D and 60 #/D is only 0.4% , while (a) film outer radius varied with mesh resolutions (b) GCI reduced as the mesh is refined  60 #/D requires an increase in mesh size of 43.7%. The minor gain in solution accuracy comes at the price of much larger computational cost, and therefore 52 #/D is selected to conduct all the deformation studies as it is a more affordable mesh without significantly deteriorating the result accuracy.

Verification
The simulation results of PHASTA are compared with published bubble deformation DNS studies (Tripathi et al., 2015). The differences between the presented simulation and Tripathi's research are summarized in Table 3.
The reference study provides rough ranges of Eotvos number ( ) Eo and Galilei number ( ), Ga and representative bubble shapes (Table 5) for different deformation dynamics (Tripathi et al., 2015). The simulation parameters in Table 4 were designed to accommodate the given ranges of Eo and Ga for each representative bubble, and the corresponding simulations reproduce the bubble topology (Table 5). Table 3 Difference between presented simulation and the reference research Tripathi et al. (2015) Simulations in this paper

Interface capturing method
Volume-of-fluid Level-set Reference system Container-based ( Fig. 1(a)) Bubble-based (under PID control, Fig. 1   The first four bubble deformation patterns (columns 2-5 in Table 5) can achieve a quasi-steady state, but the fifth bubble (last column) corresponds to an early stage deformation which will develop into an unsteady break-up, classified as a central break-up (Tripathi et al., 2015). The temporal topology change is also considered in the verification. Table 6 compares the predicted bubble deformations under identical normalized computation time ( ) τ . Regardless of the bubble topology change or time varying shape change, from dimpled ellipsoid to uneven toroid, the generated results from the PHASTA PID controller exhibit high consistency with those of Tripathi et al., further testifying its accuracy and versatility.

Validation
Two experiments were selected as validation references in this paper, Bhaga and Weber's classical experiments on bubble deformation in viscous liquids, and Sharaf 's experiment which was motivated by previously published DNS studies (Bhaga and Weber, 1981;Sharaf et al., 2017). This section reveals the advancement of experimental techniques benefiting the validations of numerical studies.   Table 7 lists the unchanged parameters in each simulation case, and Table 8 compares bubble shapes in the experiments and simulations with viscosity monotonically decreased from bubble (a) to bubble (f). Table 9 summarized the corresponding dimensionless    Table 8, the aspect ratios of numerical bubbles are found to slightly deviate from the experiments. Since the experiment only provides a range for densities and surface tensions as clues, the material properties tentatively assigned in simulations may not be consistent with the selections of Bhaga and Weber. Comparing the dimensionless numbers in Table 9, numerical Eo and Mo are close to the experimental data, with relative differences between 1%  and 8%  . However, the discrepancy of b Re can be observed to grow with bubble deformation level. As experimental measurements on bubble terminal velocity generally have higher uncertainty than other parameters, it is expected that b Re (proportional to relative velocity r u ) will consequently have larger discrepancies. In conclusion, the trends of bubbles flattening can be easily observed in both results, and the experimental and numerical Eo and Mo match well, therefore the PID controlled bubble deformation is validated by the classical experiments.
It is worthwhile to mention that due to visualization from the simulation, the existence of a cavity inside the bubble was confirmed. The experimental figures are either too dark and opaque to demonstrate the interface or are poorly clarified by streamlines to depict the cavity. Although Bhaga and Weber used "oblate ellipsoidal cap" to describe the bubble shape, the "cap" was only confirmed with halftransparent numerical results.

Validation by simulation-motivated experiment
This section will validate complex bubble deformation and break-up simulations carried out by PHASTA and PID bubble controller using parallel experiments and simulations in more recent studies (Tripathi et al., 2015;Sharaf et al., 2017). Table 10 compares the bubble deformation with a thin film forming at the wake side of the bubble, which satisfies (30,500) Eo Î and (0,50) GaÎ . Table 11 demon- Table 10 Comparison of complex bubble deformation with a film formed at the wake side of the bubble Experiment (Sharaf et al., 2017) DNS simulation (Sharaf et al., 2017) PHASTA with PID bubble controller  (Sharaf et al., 2017) DNS simulation (Tripathi et al., 2015) PHASTA with PID bubble controller strates the bubble break-up process with satellite bubbles shedding away from the leading bubble, which satisfies (1,500) Eo Î and (10,500) GaÎ . The high-speed camera both provides valuable reference on temporal interface change and increases photo resolutions dramatically. Table 10 describes the deformation process for a helmetshape bubble. Compared with the cavity shape and film thickness in the reference simulation, PHASTA results match with the experiment better. The numerical satellite bubbles in Table 11 are not identical with the experimental findings, and the satellite bubbles produced by PHASTA are with different distributions compared with the existing simulations. It is worth to mention that the thin film and satellite bubble formation in Table 11 is a recurring dynamic process so it is almost not possible to reproduce the identical satellite bubbles as the ones in the experiment or in the published simulations. However, the thin film formation and interface topology change in the break-up process agree with the existing findings. Both Table 10 and Table 11 confirm the accuracy of PHASTA code and reinforce the precision of PID bubble controller to simulate complex interface distortions. The validation also demonstrates that two-phase DNS can be as reliable as a virtual experiment. Specifically, DNS has a better visualization at interface compared with both classical and advanced experimental photos (Table 8, Table 10, and Table 11).

Bubble deformation map
This is the first bubble deformation map established utilizing PID bubble controller with local grid refinement.
The key difference is that, while all other published maps use bubble terminal velocity (also r u ) to be the output of a bubble rising and deforming in stagnant flow (Bhaga and Weber, 1981;Tripathi et al., 2015;Guillen et al., 2018), the map shown here uses r u as an input to produce equivalent bubble deformations. Most of the cases involved in this deformation map are already verified and validated adding to the creditability of the proposed deformation map. Replicated tests also showed excellent reproducibility for bubble topology under the same r u . To construct the map, the dimensionless numbers involving r u were selected as the two axes of the deformation map. Weber number (We ) is chosen to represent the bubble deformability and bubble Reynolds number ( b Re ) is used to describe the effects of viscosity and inertia.
where D is the initial bubble diameter (or the volume equivalent diameter). We and b Re take into consideration of all potential variables in bubble deformation studies, which is applicable for comprehensive deformation map establishment. Another advantage of We and b Re axes is that their common variables l r ( , , ) u D ρ are all in the numerators with similar math forms. Therefore, any parameter change in r l , , ρ u D will simultaneously increase or decrease We and b , Re which yields a deformation map with different deformation behaviors distributed elegantly (Fig. 6).
The bubble deformation map depicted in Fig. 6 includes some of the cases described in the V&V studies as well as representative cases in additional explorations. ( b Re , We ) values for each bubble are roughly denoted by the bubble center, to compare the shape change with b Re or We . Each deformation zones in Fig. 6 is lettered and further described in Table 12. Since the shapes of the zones are regular squares, it will be more convenient to clarify the involved physical phenomenon, and will provide more distinguishable criteria for application purposes.

Deformation and break-up mechanisms
This section will thoroughly discuss and analyze the mechanisms of bubble deformation and break-up. Each zone in Fig. 6 will be further explained with the help of streamlines listed in Fig. 7. As the next step, the effects of physical parameters will be clarified including surface tension , σ liquid density l , ρ initial bubble diameter , D and relative velocity r u to better understand the underneath physics of bubble deformation and break-up. Figure 7 shows the streamlines for the deformed or breakup bubble of each zone in the deformation map (Fig. 6).

Individual deformation zone analysis
Zone A represents the most common bubble shape, where the surface tension is large enough to restrain the deformation. Regardless if the bubble is rising slowly or quickly, a spherical or ellipsoidal shape can be maintained. The streamline (Fig. 7(a)) is evenly distributed around the bubble and no vortex is formed in the wake. This deformation is very easy to simulate and to achieve a steady state.
Zone B includes disk shape and mildly deformed ellipsoidal bubbles. The surface tension in Zone B is just adequate to hold a convex shape for bubbles. In the wake of the bubble, streamline plot ( Fig. 7(b)) shows that vortices are formed which reflects the bubble motions. It is found that although the bubble is primarily symmetrical, the liquid vortices have different sizes, where the left and right vortices will alternate in length over time. This is an indication that if a PID bubble controller is not utilized, the bubble will both rise and swirl in the liquid (Sharaf et al., 2017). The corresponding post-processing will require aligning the bubbles on the swirling trajectories onto the same bubble center before analyzing the deformation. With the leverage of a PID bubble controller, post-processing is negligible as the bubble center is already fixed until reaching a quasi-steady state.
Zone C corresponds to all cap bubbles seen under Table  8 and all other bubbles with cavity formed on the wake side of the bubble. The viscosity in the liquid is low so multiple vortices form denoted by the streamline (Fig. 7(c)). As the cavity in the bubble grows, the vortices formed inside the cavity will enlarge along the peripheral. The cavity will continue to expand against the liquid in streamwise direction and as a result, the bubble will become elongated until steady state.
Zone D is the peripheral break-up conducted in the V&V study of Section 4.3.2 (see Table 11). If the cap bubbles in Zone C have reduced surface tensions or are traveling faster, either satellite bubbles or films will form at the periphery of the cap and will shed away from the leading bubble. The streamline (Fig. 7(d)) shows the perturbation and asymmetry resulted from the satellite bubbles. The mass loss in this process is mild and the vortex motion is moderate, and thus the leading bubble could maintain a quasi-steady state in each shedding process. The film thickness or satellite bubble size is less than 7% of the leading bubble diameter, and therefore the peripheral break-up can sustain.
Zone E depicts helmet bubbles as seen under V&V of Section 4.3.2 (see Table 10). Different from Zones B, C, and D, the vortices are formed completely inside the bubble (Fig.  7(e)). The liquid vortices will mildly and gradually shed the gas away from the center to the rim of the bubble, and thus a film can form and grow. The surface tension is comparatively low in Zone E, and thus the bubble deformation is severe, gradually developing a helmet shape. The cavity inside the bubble grows more spacious while vortices become more unstable. However, as the viscosity is comparatively higher, the helmet bubble can maintain integrity even under great deformation. Furthermore, the relative velocity between the gas and the liquid is small such that the overall bubble shape is primarily spherical and a steady state will be achieved.
Zone F is the only unstable pattern among the six zones. The relative velocity is so high that the liquid breaks through the center of the bubble, which forms a toroid shape. The opening allows the liquid on the front and wake sides of the bubble to interact with each other, which creates vortices encircling the toroid like magnetic field lines ( Fig. 7(f)). The toroid bubble will continue to expand which widens the opening. Mass flux of the liquid going through the front side of the opening gradually increases, which distorts the vortices inside the bubble and forms a very thin and wide film. When the surface tension force can no longer balance the bubble expansion, the entire film will be torn off from the bubble and the vortices will encircle the film. This process results in destructive topology change.

Parameter analysis
Surface tension σ plays a very important role in all deformation and break-up schemes. From the bottom to the top in Fig. 6, it is clear that deformability increases with 1 σor We . For slightly or mildly deformed zones (A, B, and C), as 1 σincreases, the bubble gradually loses the convex shape and becomes half-convex and half-concave. Under high deformability (Zones E, D, and F), the bubble will form a film on the wake side. As b Re increases, the film will be more unsustainable. Zones E, D, and F demonstrate the sequence of the film to be stable, quasi-steady, and torn off.
Liquid density l ρ affects both the inertia of the liquid and the deformation of bubbles. The most affected zones are E and F. In Zone E, as l ρ grows, the liquid vortex circulation inside the bubble cavity is enhanced. As a result, the cavity will grow deeper, and the film will grow thinner and longer. In Zone F, the main reason for central break-up to occur is the inertia of the liquid. The heavier the liquid, the more chaotic the break-up will be.
Initial bubble diameter D is more important in complex deformation and break up schemes (Zones D and F). Only large bubbles (i.e., 20 mm R = ) were found to generate the deformation in Zone D as the vortex needs space to grow inside the bubble. Large D also corresponds to a small local surface tension, which is essential for helmet-shape deformation to be achieved. The majority of PHASTA simulations yielding Zone F break-up are large bubbles. However, there are indeed instances where central break up occurs with 1.76 mm R = using very heavy liquids. This may be a new discovery since in the published studies all presented bubble sizes with the same break-up dynamics are of a magnitude of 20 mm R = (Tripathi et al., 2015). Relative velocity r u is also analyzed as a consequence of buoyancy force. For less deformable bubbles, the bubble shape changing from Zone B to C is dependent on either the reduction of surface tension, or the increase of r u . For break-up regimes, low r u is essential in order to maintain stability in Zone D when satellite bubbles/film may form and shed away. At high r u the majority of energy is channeled into the vortices forming inside the bubble wake which eventually will break through the bubble center.

Conclusions
Understanding of bubble deformation and break-up dynamics is important for modeling two-phase flow behavior in a variety of conditions. Compared with published experiments and DNS efforts where deformation is studied with bubble rising in stagnant flows, this paper changed the reference system and performed DNS of liquid flowing around a fixed bubble. This approach allowed performing extreme mesh refinements around the regions of interest without the need for adaptive mesh refinement which would follow a moving bubble. With the support of unstructured mesh, grid convergence study demonstrates that the complex bubble deformation studies only have an uncertainty of 3.8% with the selected mesh configuration. V&V of simulation results agree well with both previously published experimental and numerical results, which demonstrate the accuracy of interface capturing using PID bubble controller and local grid refinement. Multiple simulations were carried out covering a wide range of bubble sizes and fluid properties, eventually establishing a bubble deformation map. There are six deformation zones in the map, denoting bubble shapes as sphere, disk, cap, skirt, helmet, and toroid. Each zone distinctly represents a deformation or break-up mechanism and can be a reference for future research and engineering applications.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.