Interfacial Area Transport Equation for Bubble Coalescence and Breakup: Developments and Comparisons

Bubble coalescence and breakup play important roles in physical-chemical processes and bubbles are treated in two groups in the interfacial area transport equation (IATE). This paper presents a review of IATE for bubble coalescence and breakup to model five bubble interaction mechanisms: bubble coalescence due to random collision, bubble coalescence due to wake entrainment, bubble breakup due to turbulent impact, bubble breakup due to shearing-off, and bubble breakup due to surface instability. In bubble coalescence, bubble size, velocity and collision frequency are dominant. In bubble breakup, the influence of viscous shear, shearing-off, and surface instability are neglected, and their corresponding theory and modelling are rare in the literature. Furthermore, combining turbulent kinetic energy and inertial force together is the best choice for the bubble breakup criterion. The reviewed one-group constitutive models include the one developed by Wu et al., Ishii and Kim, Hibiki and Ishii, Yao and Morel, and Nguyen et al. To extend the IATE prediction capability beyond bubbly flow, two-group IATE is needed and its performance is strongly dependent on the channel size and geometry. Therefore, constitutive models for two-group IATE in a three-type channel (i.e., narrow confined channel, round pipe and relatively larger pipe) are summarized. Although great progress in extending the IATE beyond churn-turbulent flow to churn-annual flow was made, there are still some issues in their modelling and experiments due to the highly distorted interface measurement. Regarded as the challenges to be addressed in the further study, some limitations of IATE general applicability and the directions for future development are highlighted.


Introduction
Gas-liquid two-phase flow plays a critical role in natural and industrial systems, which exists in air bubbles in the ocean [1], steam in the water of an electric power plant [2], carbon dioxide in methylamine of the chemical engineering [3], and oxygen gas in the hot metal of steelmaking [4]. The gas-liquid flow generates the deformation in two phases and leads to a constant change of the interphase and the medium distribution, which greatly affects the flow, heat, and mass transfer characteristics. Different flow rates, pressure, channel layout, and channel geometry can result in a different distribution of phase interface in gas-liquid two-phase flow, that is, different flow structure modes are formed. Therefore, depending on the gas, flow properties, and container geometry, the gas-liquid flow patterns are classified: stratified (horizontal channels), slug, churn, annular, and bubbly flow.
Under certain conditions, the fluctuation of a gas-liquid interface is unstable in stratified flow, leading to a hydrodynamic slug. Slug flow is characterized by alternating slugs of liquid and large gas bubbles. In annular flow, the liquid flows around the inner circumference of the pipe as a film, while gas flows in the center of the pipe. The film becomes thinner and more liquid is carried by the gas as a mist when the gas velocity increases in annular flow. Dispersed bubbly flow is the opposite of annular flow. Bubbly flow is distributed in the liquid phase as variable size, deformable bubbles moving upward, and undergo a transition to slug flow when the gas superficial velocity increases due to the corresponding rise in the volume fraction and the coalescence of bubbles. It is readily acknowledged that the formation and motion of bubbles and interaction mechanisms of bubbles are very complicated and the gas-liquid interfacial area is the deciding factor in controlling mass and heat transfer in gas-liquid two-phase flow. As the basis of bubble motion, bubble coalescence and breakup are often regarded as a game, just like playing soap bubbles at early ages, but have a wide range of applications in natural and industrial systems. In the mass flux between the sea surface and the atmosphere, bubble breakup is considered a source of aerosol droplets in the air [5]. In metallurgical processes, gas bubbles are injected at the bottom of bulk liquid metal to mix and homogenize the metal, where the importance of bubble coalescence and breakup cannot be emphasized enough [6]. In transdermal drug delivery, the bubble with controlled size can elevate the efficiency of drug delivery and reduce medical waste [7]. As for the biotechnology industry, bubbles that rise to the surface of cell suspension can damage cells, when they break up [8]. In the food industry, a large of bubbles and foams are needed to form cavities and produce low-density fillers in the fields like adhesives. In the coating industry, the defoamer which prevents bubble formation can improve the quality of paint film. Therefore, the principle of bubble coalescence and breakup is an important problem to be solved urgently in industry, and bubble interactions are studied experimentally and calculated by scientists and engineers in various industries in the past and in the future.
As the determination of flow patterns in a pipe is the precondition in the computation, many approaches were developed for flow pattern maps and flow transition criterion according to own term, which can be used to predict the specific bubble behaviors with a given set of local variables. In the beginning, 1D system codes are frequently used for the description of two-phase flow, which is based on the empirical work and can predict the gas-liquid flow in the case of the given volume flow rates of liquid and gas. However, these flow patterns can neither predict the flow pattern changes along the flow path nor the flow structures with time and space. Then, a 3D computational fluid dynamics (CFD) program has been tried to solve the above issue by combining it with additional transport equations for fluid particle number density or its corresponding parameters in the advanced two-phase flow model. On account of generous theoretical and simulation analysis, the definition of bubble interaction mechanisms is very important to obtain accurate and reliable two-phase flow phenomena by CFD simulation. A better method than algebraic approximation is to directly use transport equation to simulate the interfacial concentration, especially for the 3D of two-phase flow. Therefore, the interfacial area transport equation (IATE) developed from the population balance equation is one of the prospective ways to dynamically describe the temporal and spatial evolutions. A major advantage of IATE is that the evolutions of the interfacial area concentration (IAC) can be modeled without flow patterns.
The paper presents a review of constitutive models for gas-liquid two-phase flow in the publications, especially for fundamental issues. The transition of flow patterns of bubble coalescence and breakup are introduced in Section 2, transport equations of bubbles are provided in Section 3. Then, bubble interaction mechanisms and comparisons of constitutive models for bubble coalescence and breakup are generalized in Sections 4 and 5, respectively. In Section 6, conclusions and future developments are highlighted.

Transition of Bubbly to Churn-Turbulent Flows
IAC is an essential parameter of the phase interaction terms in the two-phase model, and its dynamic prediction has been realized by IATE [9]. It is consistent with the change mechanism of the interfacial area and able to predict the patterns without static flow. So, discontinuity caused by the transition between flow patterns can be avoided. In order to predict IAC, the physical mechanisms of bubbly to churn-turbulent flows should be studied first.

Bubbly to Slug Transition
In bubbly flow, fluid turbulence is an important cause of bubble interactions. Transfer motion and interaction of bubbles affect the flowing path, and the flow evolution is determined by the interactions of bubble migration, coalescence, and breakup as shown in Figure 1 [10]. In Figure 1 (left), a stable bubbly flow is observed far away from the transition region to slug flow and bubbles are nearly spherical. Upward bubbly flow has the wallpeaked void profile, while downward bubbly flow concentrates some small bubbles at the pipe center. Coalescence and breakup occur among bubbles with approximately the same diameter. In Figure 1 (right), the bubbly flow changes to slug flow obviously, where there are both nearly spherical bubbles and elliptical or cap bubbles. In such a case, larger bubbles in the form of ellipses or cap bubbles concentrate in the central region. Finally, two or more cap bubbles merge into a bigger bubble. If the bubble is sufficiently large, it will fill most of the cross-section of the pipe, and slug flow may result. For bubbles with different sizes have significantly different decomposition volume fraction and turbulent energy dissipation rates, bubble coalescence, and breakup should be considered locally. As the above discussion, the transition from dispersed bubbly to slug flow observed at low gas flow rates requires a coalescence process, where discrete bubbles with a diameter are closed to the pipe form a larger vapor space. As the gas rate increases, so does the amount of gas in the pipe. In such a case, tighter bubbles lead to an increase in the coalescence rate. Meanwhile, turbulence fluctuation can lead to the breakup of large bubbles in the coalescence process, as the liquid rate increases. If the breakup rate is high enough to avoid recoalescence, dispersed bubbles can be maintained. Therefore, a critical size of 15 mm for air-water is observed, bubbles deform and move randomly along the zigzag path [11]. They collide and coalesce, forming larger-cap bubbles similar to Taylor bubbles but with a diameter smaller than the pipe. In a high gas flow rate and a low liquid flow rate, the bubble packing intensifies. The value of volume fraction at which the transition to slug flow discussed above is about 0.25 to 0.3 [11].

Slug to Churn Transition
If small bubbles coalesce to form Taylor bubbles, a slug flow can form from the bubbly flow. If two serial Taylor bubbles are too small to be stable, a churn flow occurs. As Taitel et al. mentioned [11], churn flow is a developing length region for slug flow in essence. If a position along the pipe length is shorter than the developing length, churn or slug flow can be observed. On the contrary, only the slug flow exists.

Slug/Churn to Anular Flow Transition
For a high gas flow rate, the flow goes annular and the liquid film flows upward along the wall. The gas flows in the center, carrying droplets. The film flows uphill against gravity because the rapid-moving gas core exerts a force on it. This film has a wavy interface, which can be broken easily and enter the gas core in the form of droplets. In fact, when the gas flow rate is low, droplets will fall back and accumulate, resulting in churn or slug flow.

Interfacial Area Transport Equation (IATE) for Bubble Coalescence and Breakup
As mentioned in Section 2, bubble coalescence and breakup should be considered locally for various bubbles in bubbly, slug, churn, and annular flow patterns. IATE can handle these bubbles in two groups: the spherical/distorted bubble group and the cap/slug bubble group, leading to four-type interactions as shown in Figure 2. Two kinds of IATE for all bubbles are possible and needed to cover the flow patterns including bubbly, slug, churn, and annular flow as seen in Figure 3. Group I bubbles (spherical and distorted bubbles) in bubbly flow are treated in one-group IATE, while others (cap-shaped bubbles) are managed into two-group IATE in slug, churn, and annular flow. If there are not Group II bubbles, two-group IATE can be simplified to one-group IATE [12]. In addition, the main challenge of the two-group IATE comes from the complex modeling of inner-and inter-group interactions (as seen in Figure 2), and two-group bubbles provide slight contributions at low volume fraction which may be ignored for practical modeling [13]. Constitutive models of one-group IATE which can discuss the dynamic characteristics of interfacial area enough is our main focus. Based on Boltzmann transport equation, an IATE was formulated by Kocamustafaogullari and Ishii [9] as: where f (V, x , t) is the bubble number density distribution function. The right-hand side (RHS) of this equation is the source and sink terms. The first term S j represents the net rate of change in the number density distribution function due to the bubble breakup and coalescence processes, while S ph is the bubble number sink or source rate due to phase change. Even though the Boltzmann transport equation is much more macroscopic, a better version is the bubble number density transport equation. As the IAC of bubbles is related to the bubble number and size [14], the IATE can be derived from Equation (1). Applying the Leibnitz integration rule, the bubble number density transport equation can be written as: where v pm , the average local particle velocity, n is the number of particles of all sizes per unit mixture volume and R j are the source and sink rates.

One-group IATE for Bubble Coalescence and Breakup
Similar to the above, Ishii and Hibiki [14] developed their IATE as follows.
where α i is the average interfacial area density of all bubbles of volume V min and V max and v i is the interfacial velocity. For the expression of volume source one obtained: The third term on the left-hand side (LHS) is the change of interfacial area density due to bubble volume change, the source term expresses V max V min ∑ j S j dV = ∑ j R j (in the bubble number density transport equation) and the sink term is V max V min ∑ j S j A i dV = ∑ j Φ j (in the interfacial area density transport equation), where Φ ph = πD 2 bc R ph , and D bc is the bubble critical size due to phase change.
Bubble-bubble and bubbles-turbulence interactions lead to the temporal and spatial evolution of gas-liquid two-phase flow, which has aroused extensive attention. Models of interaction area transfer proposed by Wu et al. [13,15], Ishii and Kim [16], Hibiki and Ishii [17], Yao and Morel [18], and Nguyen et al. [19] were used widely. Of these numerous studies pertaining to bubble coalescence and breakup for bubbly flow, the main phenomenon and mechanism have been summarized: (a) bubble breakup due to turbulent impact (TI); (b) bubble coalescence due to random collision (RC); (c) bubble coalescence due to wake entrainment (WE). Major bubble interaction mechanisms in bubbly flow are shown in Figure 4. With bubble Reynolds number increase, the shape of bubbles, mechanism, and constitutive equations changes. Taking three bubble interactions above into consideration, Wu et al. [15], Ishii et al. [20], and Kim [21], formulated the one-group IATE for bubbly flow. In their studies, a one-group IATE is developed: where R TI , R RC and R WE are the bubble number source/sink rate due to turbulent impact, random collision, wake entrainment, respectively. The LHS in Equation (5) represents the time-dependent and the convective rate of change of the interfacial area density. In the RHS, the first term represents bubble expansion or compression due to pressure effects, the second term is the rate of change due to bubble interaction mechanisms and, the third, due to phase change.
The mechanism of the IATE is likely to depend on the bubble mixing length, turbulence intensity, and volume fraction. Therefore, the evolution of the number density distribution is described by the destruction and creation of the population due to bubble coalescence and breakup. In order to understand the mechanisms of bubble transport and the transient change of the flow patterns, the mechanisms and dynamics of coalescence and breakup must be studied systemically as summarized in Section 4.

Two-Group IATE for Bubble Coalescence and Breakup
A wide variety of bubbles sizes and shapes result in various bubble-turbulence and bubble-bubble interaction mechanisms. Considering bubbles with various sizes and shapes under a given set of gas-liquid two-phase flow conditions, an IATE that can manage the transport characteristics of various bubbles is necessary. As shown in Figure 3, Group I bubbles only exist in the spherical and distorted bubble regime, while Group II bubbles exist in the cap bubble regime. The boundary of two regimes is given by Ishii and Zuber [22] as the following Equations (6) and (7), respectively: where D distorted,max is the maximum critical size of distorted bubbles in Group I, D cap,max is the maximum stable bubble size limit, σ is the surface tension, and ∆ρ is the density difference.
Averaging the volume range of each group in Group II, a two-group IATE was obtained by Ishii and Kim [12].
where χ is the inter-group transfer coefficient, V cr is the critical bubble volume, D sc and D sm1 are the surface equivalent diameter and Sauter mean diameter, and subscripts 1 and 2 stand for Group I andII, respectively. The LHS of Equations (8) and (9) are the gross changing rate for IAC of each group and the terms in the RHS are the changing rate for IAC due to particle volume change, phase change, and intergroup transfer due to particle volume change and various particle interactions.
In the latest study proposed by Worosz [23], two-group IATE was statistically obtained as follows.
where Φ nuc is interfacial area concentration source/sink rate due to nucleation. Therefore, two sets of equations are used to describe the generation and destruction rates of bubble number density, volume fraction, and IAC for two-group bubbles due to expansion and compression, coalescence and disintegration, and phase change.
The one-group IATE proposed by Worosz is consistent with the one by Ishii and Kim [12], but a difference in the view on two-group IATE (Equations (10) and (11)). In the two-group IATE of Ishii and Kim (Equations (8) and (9)), the effects of the intergroup transfer due to particle interactions and particles expansion across the group boundary are considered into particle volume changing rate and they are double-calculated.
For the sink term and source term concerning the bubble coalescence and breakup mechanism are the main components of IATE, the constitutive models of bubble coalescence, and breakup which are usually established from experimental data based on various flow conditions is largely dependent on the performance of IATE. Thus, there are a large number of bubble coalescence and breakup models which are usually supported by the experimental data. However, recent findings [23,24] indicated that the performance of IATE in transition flow is not satisfactory. There are still some issues and difficulties: the establishment of bubble behavior turbulence model, bubble coupling of mechanisms, the lack of reliable experimental data, and the simultaneous occurrence of coalescence and breakup. Therefore, an evaluation of existing bubble coalescence and breakup constitutive models in the literature is needed as a basis for further developments, which is reviewed in Section 5.

Identification of Bubble Coalescence and Breakup Interaction Mechanisms
As previously mentioned in Section 3, most coalescence and breakup mechanisms are based on the assumption of the interactions among spherical bubbles. Bubble interaction and transfer mechanisms depend on bubble shapes. Due to the wide range of bubble shapes and sizes, a set of transport equations describing bubble transport in a two-phase flow is needed. Fully understanding the dynamics of a single bubble is the initial form of multi bubble system modeling and is a fundamental issue of bubbly flow. Therefore, most publication papers study the behavior of a single isolated bubble and the transfer phenomenon of dispersed bubbles in bubbly flow. The process is shown schematically in Figure 5.

Bubble Coalescence Mechanisms
Bubbles collide with each other constantly in involving gas-liquid dispersion. The phenomenon of coalescence entails three stages: initial contact between bubbles, controlled essentially by the hydrodynamics of the liquid, which leads to a film of the thickness of a few microns separating the two bubbles. The second step is the thinning to a few Angstroms, of which the rate is decided by the hydrodynamics of thin films. The final step is the rupturing of this film resulting in the coalescence. Obviously, the rate of bubble drainage and thinning in the second stage determines whether the bubbles coalesce or not. If the drainage time to rupture is longer than the contact time, the two bubbles will separate rather than coalesce. The final step, or the rupturing step, is usually faster than the other two stages.
The large spherical cap bubble with a strong wake region may lead to the collision and coalescence of other bubbles when others enter the wake region and accelerate. Bubble coalescence depends on bubble size, velocity, and collision frequency. Two main coalescence mechanisms have been identified: bubble coalescence due to random collision (RC) and bubble coalescence due to wake entrainment (WE) as shown in Figure 6.

Bubble Breakup Mechanisms
A bubble disintegrates into two or more new bubbles in the breakup process, and the interfacial area increases and so does the interfacial transfer between phases. Group I bubbles (spherical/distorted bubbles) break up due to the turbulent impact (TI) of the eddies against the bubble as shown in the left of Figure 7. As for Group II bubbles (caps or slugs), smaller bubbles shear off from them with a low viscosity fluid as shown in the intermedia of Figure 7, called shearing-off (SO). When the properties of the bubble surface are unstable, smaller bubbles are formed from the unstable bubbles as shown in the right of Figure 7. This breakup process is due to surface instability (SI).

Breakup due to Turbulent Impact
When the amplitude of the oscillation is close to the amplitude needed to destabilize the surface of the bubble, it begins to deform and stretch in a direction. This leads to the further contraction of the neck and finally breaks into two or more bubbles. Furthermore, bubble breakup is attributed mainly to pressure fluctuations along its surface and by bubble-eddy collisions in the moving fluid or turbulent liquid. In the view of the force required to cause a breakup, the disintegration process can be expressed as the balance between dynamic pressure and its surface stress. Hence, the criterions in the literature are summarized due to different conditions and various bubble properties (see Table 1). There are three categories for the breakup: turbulent kinetic energy, fluid velocity, or inertial force hitting on the bubbles. As turbulent kinetic energy is the function of fluid velocity, (a) and (b) in Table 1 have a similar effect on breakup behaviors. If there is no turbulent flow in fluid and the bubble with a low velocity is broken up by an object, behaviors of the bubble can be explained by the two criterions above. Therefore, d) is a solid option. In brief, combining turbulent kinetic energy (velocity) and inertial force together was proposed by Wang et al. [25] and Zhao and Ge [26] is the best choice for the bubble breakup criterion. Table 1. Breakup criterion of several literature [27].

Criterion Authors
(a) Turbulent kinetic energy of the particle greater than a critical value Coulaloglou and Tavlaride [28] Chatzi et al. [29] (b) Velocity fluctuation around the particle surface greater than a critical value Narsimhan and Ghupta [30] Alopaeus et al. [31] (c) Turbulent kinetic energy of the hitting eddy greater than a critical value Lee et al. [32,33] Prince and Blanch [34] Tsouris and Tavlarides [35] Luo and Svendsen [36] Martinez Bazan et al. [37] (d) Inertial force of the hitting eddy greater than the interfacial force of the smallest daughter particle Lehr and Mewes [38] Lehr et al. [39] (e) Combination of (c) and (d) Wang et al. [25] Zhao and Ge [26]

Breakup due to Shearing-Off
With the increasing bubble size, breakup mechanisms become much more complicated owing to additional phenomena such as shearing-off and interfacial instability, both of which are induced by a velocity difference across the interface. During the shearing-off process, that many small bubbles are sheared off from a large bubble is called erosive breakage. When the relative velocity is high enough to make the bubble skirt unstable, the bubble skirt can fall off from the large one with the generation of small bubbles at the edge. However, experiments and theoretical studies have revealed that the interfacial viscous shear force can be neglected due to the low viscosity of water in the air-water flow, in which the shearing-off process is caused by the gas distribution inside the cap bubble [40]. In a word, the shearing-off is determined by the balance between the viscous shear force and the surface tension at skirts of the cap/slug bubbles in highly viscous flows

Breakup due to Viscous Shear Forces
When the fluid has a high viscosity, viscous shear forces in the continuous phase lead to a velocity gradient near the interface, resulting in bubble deformation and final breakup. If a large part of a trailing bubble locates outside the wake region, the shear stress across the wake boundary may split the bubble due to extension, and the bubble surface is depressed and necked. In the first stage, the bubble is elongated into two lumps by a thread and then divided into two equal-sized daughter bubbles corresponding to the lumps and a series of smaller bubbles called satellites. The bubble may also break into a cylindrical thread which may break into many smaller bubbles, where it is called thorough the breakage [9,41,42]. In terms of the viscous case, the breakup mechanism is expressed as the force balance between surface tension forces and external viscous stresses, which is usually formulated as a dimensionless Capillary number.
Expressed as the ratio of viscous stress over the surface tension, it can be formulated as If the Capillary number or bubbles size surpasses the critical value, that is, Ca ≥ Ca cr , interfacial tension forces cannot keep the form of the bubble, leading to the breakup of two or more daughter bubbles. The critical diameter is determined by As for the case of turbulent flows, the viscous force includes the laminar shear in the bulk flow and the turbulent one ε ν . The latter one is responsible for the collision of bubbles with the size smaller than the Kolmogorov length scale. Briefly, the breakup mechanism due to viscous shear force is expressed by the force balance between surface tension forces and external viscous stresses.

Breakup due to Surface Instability
The mechanisms above depend on the dynamic characteristics of continuous flow. However, it is proved that even if there is no net flow in the continuous phase, the interface instability will lead to breakup. The breakup process due to SI includes two instabilities: Rayleigh-Taylor instability (RTI) and Kelvin-Helmholtz instability (KHI). RTI is common in natural and industrial processes. It describes the formation of irregularities along the gasliquid interface, as it is accelerated in a direction perpendicular to the gas-liquid interface of bubbles. That is, when there is a density difference in two fluids, RTI occurs. KHI arises due to shear along an interface between the gas and liquid. Thus, when the density ratio is approximately 1, the breakup process is dominated by RTI.
In general, five bubble interaction mechanisms in the general two-phase flow are summarized. They are (1) RC: bubble coalescence motivated by turbulent eddies; (2) WE: bubble coalescence due to acceleration of the following bubbles in the wake of a leading bubble; (3) TI: bubble breakup due to the impact of turbulent eddies; (4) SO: small bubbles shearing-off from the rim of a large bubble; (5) SI: breakup due to surface instability at the interface. Table 2 summarizes the inner-and inter-group interactions based on the discussion above. Superscript (1) stands for Group I bubbles and (2) for Group II bubbles. Table 2. Major bubble interaction mechanisms [43].

Mechanism
Interaction Source in 1 (multiple number); sink in 2 (no number change) Source in 2

Evaluation of Frequency Models and Constitutive Models for Bubble Coalescence and Breakup
The bubble interaction mechanisms of coalescence and breakup discussed in Section 4 should be pursued as constitutive relations to solve the IATEs. In this section, several coalescence and breakup frequency expressions are reviewed, classified, and compared. Furthermore, a literature review regarding the bubble coalescence and breakup constitutive models is presented, followed by discussions and conclusions.

Breakup Frequency Models
In the past few years, great progress has been made in the analysis and modeling of the breakup process, and a large number of breakup frequency and daughter size distribution models have been proposed. In Figure 8, there are a large number of breakup frequency models, which are classified into four parts according to the discussion in Section 4.2, namely, turbulent impact, shearing-off, viscous shear force, and surface instability.
Turbulence is the most common situation, on which breakup models are focused. Compared to the turbulent impact, the influence of viscous shear, shearing-off, and surface instability in a turbulent flow is usually neglected and the corresponding theory and modeling are rare in the literature. In the classification of turbulent impact, Prince and Blanch [34] set the upper limit of the integration as 10π/D arbitrarily, rather than the lower limit of integration. They believed that there is not enough energy to break up bubbles in the case of eddies with characteristic dimensions less than 20% of the bubble size. Finally, it is appropriate to point out here that the models by Prince and Blanch [34] and by Luo and Svendsen [36] are sensitive to the upper and lower limit of the integration, so they cannot be chosen arbitrarily. Different values of the upper integration limit and of the lower integration limit have been tested and the results are shown in Figure 9.  Most of these breakup frequency models are based on a theory fundamentally similar to the kinetic theory of gases, and the models above assume that turbulence is considered as an array of eddies with well-defined sizes and densities. After defining a collision cross-section and obtaining an eddy arrival frequency, closure parameters that change the conditions, such as integration limits, have to be set. However, it is impossible to prove the models based on the eddy concept, and it is too difficult to obtain correct and reliable data concerning the density or size of eddies that take part in interactions with bubbles.

Coalescence Frequency Models
Coalescence is considered to be more complex than breakup, for the interactions of bubbles and between bubbles and surrounding liquid are involved, once they are gathered together by outflow or external force. The coalescence frequency is calculated by the colli-sion frequency and the coalescence efficiency which are based on the physical quantities, and most of them considered only the inertial collision caused by turbulent fluctuations (see Figure 6). Collision frequency models are classified into RC and WE, but both frequencies derive the mathematical functions with bubbles size, fluid properties (viscosity and turbulent energy dissipation rate) by experimental data or theoretical derivation in Figure 10. Collision efficiency models are divided into three models, namely, the energy model, film drainage model, as well as the critical velocity model. Figure 10. Classification of bubble breakup models.
As previously stated, there are various mechanisms that promote the collision. Notwithstanding the fact that five mechanisms are listed [27], one main thing of the one-group interfacial area density transport equation is the turbulent random motion-induced collisions. The bubble collision is induced by the turbulent fluctuation around them. The concepts of collision cross-sectional area and bubble relative velocity u rel are of leading importance. Kennard explained that collision efficiency is the effective volume swept by particles in unit time [45]. In order to determine the bubble relative velocity, Coulaloglou and Tavlarides [28], Lee et al. [32,33], Prince and Blanch [34], and Luo and Svendsen [36] regarded the colliding bubble velocity as the velocity of an equally sized eddy, in which way bubble relative velocity and collision frequency are obtained. It is necessary to apply three modification factors (see Table 3) into collision frequency. Although these modification factors are of help to obtain collision frequency, the necessity of Π and the formulation of Υ, Π still need for a further study, because the main difference in the models presents the results from the calculation of u rel and modification factors Υ, Π. For example, Chesters [46]

Λ
The effect of the size ratio between bubbles and eddies.

Y
The effect of the bubble that reduces the free space for bubble movement and causes an increase in the collision frequency.

Π
The ratio of the mean distance between bubbles to their average relative turbulent path length.
Turbulent collision frequencies of the dependent on the bubble sizes from the literature are shown in Figure 11. With the increase in bubble size (D 1 and D 2 ), the collision frequencies from all authors increase. Even though all authors take into account the buoyancy, shear rate, and wake interaction, the collision frequency of Wang et al.'s (the purple curve) [47,48] is the highest and twice the value of Prince and Blanch's [34] (the blue one) due to an identical efficiency for all collision mechanisms by Prince

Constitutive Models for the One-Group IATE
One-group IATE formulations in Section 3.1 and its bubble coalescence and breakup frequency models in Section 5.1 have been proposed. Though a lot of efforts have been made to develop them in the past two decades, some obvious discrepancies exist in their constitutive models. Therefore, it is necessary to compare the typical expressions for bubble coalescence and breakup frequency and efficiency terms that have been commonly used in the literature. Analyzed on the constitutive models, frequency and efficiency of bubble coalescence and breakup obtained by these models are compared in this section.
The bubble coalescence and breakup constitutive models for the one-group IATE are usually developed from the experimental data taken in various flow conditions and flow channels. The salient findings of these studies are summarized in Table 4.
K TI = 6.0, C TI = 0.021 Yao  For the consideration in the modeling of bubble interaction mechanisms, Wu et al. set up bubble interaction mechanism terms, considering the contribution of RC, WE, and TI of the eddies against the bubbles [13]. In order to determine the adjustable parameters in the source and sink terms, experimental data of air-water upward flow in a 50.8 mm diameter pipe were used. Rather than WE, only RC and TI of the eddies against the bubbles were considered in the models of Hibiki and Ishii [49], and local flow measurements of vertical upward air-water flow in a round tube with an inner diameter of 50.8 mm were performed. Yao and Morel [18] took the contribution of RC and the TI of the eddies against the bubbles into consideration. Compared to previous studies, a significant modification made by Yao and Morel arises in two points: (a) bubble coalescence time = free traveling time + interaction time; (b) the breakup interaction time is modified by the bubble-eddy resonance mechanism. Furthermore, the contact time can be obtained from the characteristic time of eddies with the same radius of bubbles. Considering turbulent suppression phenomena in the turbulent bubbly flow with high velocity, Nguyen et al. [19] improved the RC model and the TI model proposed by Yao and Morel [18] and made a dynamical prediction of IAC evolution by the following RC model and TI model. The experimental conditions of Nguyen et al.'s work is the adiabatic air-water bubbly flow taken in vertical upward round pipes with inner diameters of 80 mm [50].
For the work on empirical coefficients, Ishii and Kim [16] developed a micro foursensor conductivity probe and obtained the main parameters in radial distribution, providing the basis for studying the evolution of the two-phase flow. They enlarged the database for assessment of the source and sink terms proposed by Wu et al., simplified the evaluation procedure, and assumed that all parameters along the radial direction have uniform profiles, in which way the covariance of the terms under study can be neglected. Hence, a new set of empirical coefficients were presented by them. Applying the model by Wu et al. [13] into a 3D simulation of vertical upward two-phase flow in a pipeline using Fluent of CFD code, Wang [51] tested the primary model and proposed a set of new 3D simulation coefficients. The method for determining these coefficients is based on the experimental observation that different bubble interaction mechanisms dominate under different flow conditions. According to this thought, each coefficient can be estimated independently.
Although a lot of studies have been conducted on the IATE in the past two decades, further efforts should be made to improve the prediction accuracy and general applicability of the IATE. Some points can be summarized as follows:

1.
Flow orientation and channel size: As summarized in Table 4, nearly all the interfacial area transport data is derived from the vertical two-phase flows, hence the constitutive relations in vertical flow are fully-fledged compared to other flow directions. The latest applicability for horizontal flow was proposed by Talley [52]. Only one mechanism (TI) was observed in horizontal flow [52,53], while RC, TI, and WE were managed in vertical bubbly flow. It is necessary to model additional covariance parameters to consider the effects of non-uniform and asymmetric gas distribution in horizontal flow. Covariance parameters are coincident in vertical flow, but different for horizontal flow, owing to a skewed bubble distribution [54]. In order to calculate the covariance parameters, Talley proposed a void fraction profile reconstruction method to estimate the local void fraction [52], in which a 1D one-group IATE for horizontal flow was proved by the experiment with a small diameter horizontal pipe (38 mm) [52]. Yadav [55], Yadav et al. [56,57], Qiao [58], Qiao and Kim [59,60], and Qiao et al. [61] tried to study the interfacial structure in other orientation flows, but it was not enough. Although most of the work has focused on the straight pipe, two-phase flow undergoes various limitations such as elbows and tees in practice. Therefore, bubble coalescence and breakup IATE should be benchmarked by consistent data at various flow directions in a future study.
Kong [62,63] and Talley [52,53] thought that their work could be applied to different diameter pipes (38 mm, 51 mm, and 102 mm). As discussed above, the current interface area transport models are applied to a relatively small and limited channel with a diameter of 19.2~102 mm. With the decreasing of the channel size, RC may be restrained and WE may also become much more important. In short, channel size should be studied in a wide period.

2.
Pressure and liquid type: Water interfacial transport data derived from standard atmospheric pressure (0.1 MPa) is used for benchmarking existing sink and source terms, except for the boiling refrigerants and SF6-ethanol, where they can be used for the simulation of high pressure [64,65]. It is probably due to bubbles in water are common and used everywhere. As is known, bubbly flow with high pressure in industries, such as steam in the water of an electric power plant, oxygen gas in the hot metal of steelmaking, chlorine in sodium chloride molten salt system of titanium metallurgy, and carbon dioxide in methylamine of the chemical engineering, are widely used. Therefore, bubble studies under high pressure, metal liquid, and molten salt liquid should be carried out in the future for validation [66].

3.
Simulation of developing and transient flow [65]: Most of the interface area transport data is obtained at a certain distance from the inlet of the experiments section in steady flow as seen in Table 4. It means that the existing interface transport models have been verified under stable quasi-developing conditions, but not under transient flow conditions. The applicability of the interface transport model in developing and transient flow needs to be further studied.

4.
Bubble coalescence and breakup: Compared with the dimensions in Table 4, bubble coalescence models were developed from 1D to 3D, from a simple constant (C RC ) to a function of the Weber number. Furthermore, bubble breakup models were evaluated from the function of We We cr and kept at a similar form for almost all the models. The steady-state shape of bubbles depends on the balance between drag and surface tension. For this reason, the Weber number has been used to map possible bubble shapes as shown in Figures 12 and 13 [55]. Considering the coupling of the structure of air-water bubbly flow in a vertical pipe (volume fraction, energy dissipation rate, and bubble size distribution) [67], the evaluation of one-group constitutive models is presented in Sections 5.2.1 and 5.2.2.

Bubble Coalescence due to RC
As shown in Figure 14a, four random collision frequency models are compared. As summarized in Table 4  The trend of the coalescence efficiency as a function of the turbulent energy dissipation can be seen in Figure 14b, which varies significantly amongst different authors. In general, bubble coalescence occurs in the following process: (a) bubbles collide, (b) a little liquid between the bubbles is rejected and gradually drain, (c) when the film between bubbles reaches the critical thickness, the film breaks and bubble coalescence occurs. The results of the three models above are constant or decrease with the increases of turbulent energy dissipation. Significantly, Wu et al. did not apply a coalescence model, because the coalescence rate decreases exponentially with respect to the turbulent fluctuation velocity and different liquid flow conditions will lead to great differences. In order to describe the randomness of coalescence after each collision, a constant coefficient is used as the collision frequency η c , which is not given explicitly but contained in C RC . Regarding the coalescence efficiency as a constant, the models proposed by Wu et al. and Ishii and Kim are impractical. Nguyen et al. took turbulent suppression phenomena into consideration. They valued the adjustable coefficient K c3 of 0.913, resulting in a little bit higher coalescence efficiency than the one proposed by Yao and Morel, which is regarded as a solid consideration.
The turbulent energy dissipation is a function of the random collision frequency multiplied by a constant coalescence efficiency. As can be seen in Figure 14c

Bubble Breakup due to TI
Bubble breakup is considered to be due to the collision of turbulent eddies with bubbles. As for the effect of turbulent eddies on bubbles, all available bubble-eddy collision frequency and efficiency models can be compared, because the calculation formulation and almost all the coefficients have been given clearly. In Figure 15a Figure 15b). In short, Wu et al.'s model cannot be applied to the flow with low turbulent dissipation (less than 9).
As a function of turbulent energy dissipation, bubble breakup frequency models are compared in Figure 16a and the highest value is again obtained by Nguyen et al.'s model. It should be noted that there are the same trends in terms of bubble coalescence and breakup. The model delivering the second-highest rates is the one proposed by Yao and Morel. The curves of Ishii and Kim and Hibiki and Ishii's models in Figure 16a are almost zero when turbulent energy dissipation is from zero to ten. It indicates turbulent energy dissipation has less effect on bubble breakup frequency (see η b in Table 4) and is opposite with the facts.
As seen in Figure 16b However, the assumption-based values and the evaluation of model capability are not correlated. As a matter of fact, the source and sink terms explained in the previous sections are applied together, for each author, to modify the transported value in the IATE. Thus, the resultant value from the balance between creation and destruction of bubbles is what matters-the values of the single interaction terms are of secondary importance. Based on the discussion above, there are no perfect models for bubble coalescence and breakup. The best model for bubble coalescence due to RC is from Yao and Morel and Nguyen et al. Table 4) and the best one for bubble breakup due to TI is from Yao and Morel and Nguyen et al. ( Table 4), considering the sensibility of turbulent energy dissipation and efficiency in bubbly flow.

Constitutive Models for the Two-Group IATE
With the increase in bubble size, transport mechanisms become much more complicated because of the extra inner-and inter-group interaction mechanisms. Hence, two-group IATE for bubble coalescence and breakup is needed. In the gas-liquid flow beyond bubbly flow, the length scale of the Group II bubble can be comparable to the length scale of the flow channel. Therefore, two-group flow can be influenced by some geometry structures such as elbows and tees and the channel structure should be considered in the modeling and simulation. As the bubble size expands and reaches a limit point (D cap,max in Equation (6)), the gas-liquid interface will become unstable and no longer be sustained. Before the bubble size goes to D cap,max , bubbles will break up due to SI. D cap,max proved by Miller et al. [68], is approximately 10 cm in an ambient air-water system. Table 5 shows the existing constitutive models for two-group IATE in various conditions. All of them presented in the vertical upward air-water system study the modeling of five bubble interaction mechanisms (i.e., RC, WE, TI, SO, and SI) in two-group IATE.         In Table 5, in a narrow confined channel (10 mm × 200 mm), Sun [69] and Sun et al. [43] proposed a model of bubble interaction mechanisms in the two-group IATE for gas-liquid flow, including five major bubble interaction mechanisms (i.e., RC, WE, TI, SI, and SO). Bubbles are confined in the depth direction rather than the width direction due to the relatively high aspect ratio of the cross-section. There cannot be stable slug bubbles in such a channel, because the channel width is larger than D cap,max [43,69]. Later, they presented the evaluation approach and results of two-group IATE on experimental data obtained in confined upward flow [77]. The data derived from their results [77] and Ozar et al. [70] in an annular channel at elevated pressures (0.58 MPa) were evaluated for 1D two-group IATE. Yang et al. [71] developed Sun et al.'s model [43] and focused on the benchmark of IATE performance in a rod bundle geometry. The newly developed two-group model showed great prediction in bubbly flow to churn-turbulent flows. Results indicated that the average IAC and void fraction prediction error is less than 15% in most of the flow conditions [71].
In a round channel, pipe size affects the bubble interactions significantly. Fu [72] and Fu and Ishii [40] proposed a set of two-group IATE for moderate diameter round pipes (25 mm~100 mm). Solid databases for a two-inch air-water loop have been built [72] and experiments indicated that WE and SO are dominant in the bubble interaction mechanisms. Doup [78] and Worosz [23] improved Fu's model and their predictions were in good agreement with experimental data. Wang et al. [24] established a database in a round pipe with a diameter of 25.4 mm (D H ≈ 2D cap,max ). Their model could not predict the drastic inter-group transfer and overestimated the IAC, hence extra experiments and modeling are necessary for the bubbly to slug transition flow in order to explain the sharp inter-group transfer in the two-phase flow. Furthermore, a more accurate model has been performed after the improvements including re-deriving the coalescence terms and optimizing the experimental coefficients.
In a relatively larger channel (D H > 2D cap,max ), Smith [73] developed a two-group IATE, and the model prediction was in good agreement with the experiments using a four-sensor conductivity probe in large diameter test sections, 101.6 mm and 152.4 mm diameter. Large bubbles will break up due to SI before they fill the pipe cross-section, hence, stable slug flow cannot exist in such a large pipe [79]. A modified model applicable to large pipes was developed by Smith [75] and Smith et al. [74], and two vertical test sections, one with a diameter of 102 mm and one with a diameter of 152 mm, were used in their work. In order to manage the poor predicting performance in interfacial area transport beyond the flow condition range in the original benchmarking effort, Schlegel et al. [76] modified the constitutive relations using an expanded range of experimental conditions extending to pipe sizes of 304 mm and the combined RMS error is below 15%.
In the past 20 years, substantial progress has been made in the simulation of bubbly, slug, and churn flow, but troubles in the experiment and simulation of churn-turbulent and annular flow. In the terms of modeling, the gas phase in churn-turbulent and annular flow is continuous, rather than discrete bubbles. Hence it is unreliable for investigators to set up models based on the previous work of bubble interaction mechanisms, and their experiments are difficult due to the highly distorted interface measurement. As discussed above, in order to expand the IATE to annular flow, the interface between the gas core and liquid film should be considered into IAC. Furthermore, the prediction ability of the 3D two-fluid model is a comprehensive coverage from bubbly to annular flow [54].

Conclusions and Outlook
Bubble interaction mechanisms play critical roles in both natural and industrial processes, of which constitutive models help to reveal, design, improve and optimize the bubble performance. To extend the IATE prediction capability beyond Group I bubbles in bubbly flow, two-group IATE which can deal with five bubble interaction mechanisms (i.e., RC, WE, TI, SO, and SI) is needed. However, to our knowledge, there is not yet a comprehensive review on both one-group and two-group IATE for bubble coalescence and breakup in the two-phase flow until now. Therefore, in the present work, bubble interaction mechanisms and their constitutive models of gas-liquid two-phase for both the one-group and two-group IATE are summarized and analyzed extensively from the literature. This includes transition phenomena from bubbly to turbulent-churn flows for bubble coalescence and breakup, IATE formulations, identification of bubble interaction mechanisms, bubble coalescence and breakup frequency models, and constitutive model development and evaluation. The following conclusions are highlighted.

1.
Models for bubble coalescence and breakup processes of interfacial area density transport equation are based on the kinetic theory of gases. (1) The breakup frequency models are inconsistent and different models show completely different behaviors, where closure parameters change the conditions, such as integration limits, can have to be set. (2) In coalescence frequency models, the main difference in the models present results from the calculation of u rel and modification factors Y, Π which need further study to obtain collision frequency. Wang et al. gave a small collision frequency for small bubbles due to the larger mean distance between small bubbles than big ones in the case of the same bubble number [47,48].

2.
Some important limitations of constitutive models for one-group IATE general applicability, such as flow orientation, channel size, pressure, liquid type, simulation of developing, and transient flow and bubble coalescence and breakup are summarized. Based on the summary and evaluation of the constitutive models for the one-group IATE in the literature, the best model for bubble coalescence due to RC is

3.
Constitutive models for two-group IATE in a three-type channel (i.e., narrow confined channel, round pipe, and relatively larger pipe) are reviewed, including five bubble interaction mechanisms, and their corresponding experimental conditions are summarized. Although great progress in extending the IATE beyond churn-turbulent flow to churn-annual flow was made, there is still some trouble in their modeling and experiments due to the highly distorted interface measurement.
For the future perspective of the IATE, its evolutions could be addressed in two different directions. On the one hand, it would be possible to complement the one-group IATE by means of the phase change and nucleation terms which makes the simulation of phase transition of subcooled boiling and saturated boiling regions possible, and the transition zone from bubbly to slug is avoided. On the other hand, as the adiabatic approximation is maintained, the two-group IATE should be considered and the modeling of other bubble interaction mechanisms is needed.
Significant efforts have been made to develop IATE over the past two decades, but this still requires further research to be conducted to improve the performance and make them applicable to a wider range of flow conditions. Therefore, inter-and inner-group terms should be proposed and the modification of the two-phase models is necessary. The inter-group transfers of mass momentum and eventually energy between group I and II bubbles need to be taken into account. Furthermore, the modeling of bubble interactions for two-group IATE is strongly dependent on the channel size and geometry. Hence different flow restrictions (elbows, tees, etc.) and flow orientations (horizontal and inclined) need to be studied in the future for the application of real life, where extensive measurements and experimental data are essential for the validation under different flow conditions. The smooth transition from one-group to two-group IATE, that is, from bubbly to turbulentchurn flow, needs to be further improved. shearing-off mechanism TI turbulent impact mechanism WE wake entrainment mechanism SI surface instability mechanism Superscripts * non-dimensional value (1) interactions within Group I bubbles (11,2) coalescence of a Group I bubble with another Group I bubble to generate a Group II bubble (12,2) coalescence of a Group I bubble with a Group II bubble to generate a new Group II bubble (2) interactions within Group II bubbles (2,1) Group I bubbles generated from breakup of a Group II bubble (2,2) Group II bubbles generated from breakup of a Group II bubble (2,11) breakup of a Group II bubble to generate two Group I bubbles (2,12) breakup of a Group II bubble to generate a (or multiple) Group I bubble(s) and a Group II bubble