Efficient simulation of bubble dispersion and resulting interaction

In this work, an efficient model for simulating bubble dispersion and coalescence due to turbulence is developed in the Euler-Lagrange framework. The primary liquid phase is solved on the Euler grid with the RANS turbulence model. Bubble motion is computed with the force balance equations. One-way coupling between two phases is assumed and the framework is designed for the computation of disperse bubbly flows at low Eötvös number. The turbulent dispersion of the dispersed phase is reconstructed with the continuous random walk (CRW) model. Bubble-bubble collisions and coalescence are accounted for deterministically. To accelerate the time-consuming search for potential collision partners in dense bubbly flows, the sweep and prune algorithm is employed, which can be utilized in arbitrary mesh types and sizes. Validation against experiments of turbulent pipe flows demonstrates that the one-way coupled EL-CRW dispersion model can well reproduce the bubble distribution in a typical dense bubbly pipe flow. Good agreement of the bubble size distribution at the pipe outlet between the simulation and the experiment is obtained.


Introduction
Dispersed bubbly flows appear in various engineering systems including chemical, mechanical, and biological applications. The modeling of bubble motion, turbulent dispersion, and interaction with other bubbles plays a significant role in the design process of industrial productions.
Much work has been carried out to develop numerical methods for computing unsteady bubbly flows. The point mass approach in an Euler-Lagrange framework is promising for the large-scale industrial applications due to its flexibility and computational efficiency. Compared with the Euler-Euler representation, the Euler-Lagrange approach provides the advantage to track each numerical bubble and model the interaction between bubbles deterministically (Mattson and Mahesh, 2012).
According to the review of Liao and Lucas (2010), the interaction between bubbles may be caused by several mechanisms: random motion due to turbulent fluctuations, shear flow because of gradients of the nearby fluid velocity, different rise velocities resulting from different bubble sizes, and the wake entrainment of bubbles moving in line. Compared to the turbulence induced interaction, other mechanisms are normally of minor significance. In order to predict bubble motion due to turbulence, it is ideal to describe all the temporal and spatial scales of the fluid phase with direct numerical simulation (DNS). The Large Eddy Simulation (LES) approach is conceptually similar to DNS but requires reduced computational effort by modeling instead of resolving the smallest length scales. However, both methods still need sufficiently fine grids and the computational resources increase rapidly with Reynolds number. The simulation of bubble dispersion with DNS or LES for industrial applications of high Reynolds numbers is valid but remains very expensive (Santarelli and Fröhlich, 2015;Fang et al., 2018). Furthermore, not resolving the smallest length scales in an LES poses problems for the prediction of bubble trajectories similar to the ones addressed in this paper for RANS.
Alternatively, the Reynolds-averaged Navier-Stokes (RANS) approach solves for the time-averaged velocity, which is more affordable for a wide range of complex industrial applications. In order to model bubble dispersion, the RANS statistics need to be enhanced by instantaneous fluid velocity fluctuations seen by the Lagrange phase. This is usually done in a stochastic manner. Iliopoulos and Vol. 3, No. 3, 2021, 152-170 Experimental and Computational Multiphase Flow https://doi.org/10.1007/s42757-020-0082-2

Nomenclature
Greek Letters α b Bubble volume fraction α l Liquid volume fraction μ l Liquid molecular viscosity (kg/(m·s)) μ t Turbulent viscosity (kg/(m·s)) ρ b Bubble density (kg/m 3 ) ρ l Liquid density ( Bubble volume (m/s) x b Bubble position vector (m) Δt L Lagrange time step size (s) Hanratty (1999) studied the dispersion of particles with zero inertia in a fully developed turbulent channel flow. The fluid velocity fluctuations along the particle trajectory obtained from a Lagrange dispersion model were compared to the DNS diffusion results, and good agreement was found. Bocksell and Loth (2001) investigated different dispersion models in a nearly homogeneous turbulent wake and in an inhomogeneous axisymmetric jet. The predicted statistical concentration of particles agreed well with the experimental data. Indeed, stochastic dispersion models are able to give realistic predictions in various applications (Dehbi, 2008). Therefore, one focus of this paper is to apply a RANS dispersion model for modeling turbulenceinduced interaction between the elements of the Lagrange phase, a prerequisite to address turbulence-induced bubble collisions. Another focus of this paper lies on modeling coalescence, occuring in a certain percentage of collision events. Darmana et al. (2006) adopted the Euler-Lagrange method to investigate bubble coalescence with the model of Chesters (1991) and the model of Prince and Blanch (1990). Recently, Hoppe and Breuer (2018) and Kamp et al. (2001) developed two coalescence models based on the studies of Jeelani and Hartland (1991) and Chesters (1991), respectively. Mattson and Mahesh (2012) applied the model of Kamp et al. (2001) in a turbulent pipe flow to predict the variation of the bubble size with the streamwise position. Good agreement between the predicted size distribution and the experimental results was obtained. The model of Hoppe and Breuer (2018) was validated against experiments of single bubble coalescence with a free surface for different bubble sizes and fluid conditions. In the present paper, both models are applied to turbulent bubbly pipe flows. To the best of the authors' knowledge, it is the first time that the model of Hoppe and Breuer (2018) is validated in bubble laden flows.
In all cases, contact and collision is the premise of coalescence. The time required for computing collisions between Lagrange bubbles can become extremely large for dense two-phase flows if no special algorithm for the computation of bubble-bubble distances is used. Owing to the small time step size used in the dispersion models to capture the necessary eddy properties, the computational burden is increased further. To overcome this problem, several algorithms have been proposed to accelerate the collision modeling (Shams et al., 2011;Breuer and Alletto, 2012). However, such algorithms have certain limitations. For example, the bubble size should not be larger than the cell size (Shams et al., 2011). In this work, the sweep and prune method is employed to determine the distance between Lagrange bubbles.
The paper is structured as follows. In Section 2 the computational method employed for the simulations is briefly introduced. The bubble dispersion model and the approach for modeling bubble-bubble interactions are presented in Sections 3 and 4, respectively. In Section 5, the solver is presented, together with a verification based on the well-mixed criterion. Section 6 reports the validation results for several turbulent pipe flows.

Governing equations of the continuous phase
The liquid phase is computed by solving the steady incompressible Reynolds-averaged conservation equations for mass and momentum in a 3D domain: where u and p are the Reynolds-averaged fluid velocity and average pressure, respectively. Note that the mean fluid velocity u is related to the total liquid velocity u by u u u ¢ =with u ¢ the instantanuous liquid velocity fluctuation. Statistical modeling of u ¢ will be discussed in Section 3. The density in Eq. (2) is the liquid density l ρ ρ = and assumed constant throughout. The effective viscosity eff μ is composed of the liquid viscosity and the turbulent viscosity t μ with eff t l μ μ μ = + . The turbulent viscosity t μ is obtained from two different turbulence models, the realizable k ε model or the Reynolds stress model (RSM). The governing equations of these models are given in the Appendix. Eqs. (1) and (2), combined with the transport equations for the selected turbulence model, are solved in a three-dimensional domain, covered with a stationary grid, the Euler grid.

Governing equations for bubbles
The bubble velocity is computed by solving an ordinary differential equation for the linear momentum considering different forces on the bubble. Generally, these forces may be divided into fluid dynamic forces (e.g., drag, pressure, virtual mass, Saffman forces, etc.) and field forces (gravity, electrostatic forces, etc.). Hence, the equations of motion for each individual bubble read: where x b is the bubble position vector, and u b the bubble velocity vector. The forces on the right-hand side of Eq. (4) are drag force, lift force, virtual mass force, pressure gradient force, gravity, and the collision force, respectively. For solid particles it may be required to account for the Magnus force resulting from rotation which, in turn, requires additionally solving an equation for the angular momentum. Due to the extremely small mass of bubbles and their different behavior during collisions this is not modeled in the present study.
The drag force D F results from the resistance experienced by the bubble moving through the liquid, which is a function of the velocity difference between the bubble and the liquid nearby. It takes the form (Crowe et al., 1977): where D C is the empirical drag coefficient, and b D is the bubble diameter. The correlation of Tomiyama et al. (1998) for clean bubbles in pure liquid: is used here. The bubble Reynolds number b Re and the Eötvös number Eo appearing in Eq. (6) are defined as u u with σ the surface tension coefficient. The correlation (6) was proposed by Tomiyama et al. (1998) The lift force F L is caused by a shear flow around the bubble. It can be expressed as (Saffman, 1965): According to Mei and Klausner (1994), the shear lift coefficient L C is positive for a spherical bubble so that the lift force acts perpendicularly to the flow direction and points in the side with larger relative velocity between bubble and liquid. The lift coefficient changes its sign if a substantial deformation of the bubble occurs. From the observation of the trajectories of a single air bubble rising in simple shear flow of a glycerol-water solution, Tomiyama et al. (2002)  -.
The bubbles studied in the present work are mainly spherical. So, for simplicity, the constant positive lift coefficient L 0 5 C = . is used here.
The virtual mass force F vm in Eq. (4) accounts for the additional work that the bubble performs in order to accelerate or decelerate the surrounding fluid: with vm C the virtual mass coefficient. The assumption of vm 0 5 C = . is valid for isolated spherical bubbles (Auton et al., 1988). Although vm C can vary if an additional bubble is located in close vicinity of the considered bubble (Kamp et al., 2001), the bubbly flows studied in this work are not contact dominated and the average inter-bubble distance is still large. The same holds for bubble-wall distances. Therefore, a constant value vm 0 5 C = . is selected here. The sum of the pressure force and the gravity force can be expressed as (Laín et al., 2002): The last term in Eq. (4) F coll , describes the collision force which is discussed in Section 4.1.2. Eqs. (5), (7), and (10) contain the instantaneous liquid velocity u which is the sum of the mean velocity u and the instantaneous fluctuation u ¢ . With a RANS approach according to Eqs. (1) and (2), only u is available, however. Setting u u = yields wrong results, as this disregards turbulent dispersion, for example. Hence, an adequate modeling approach is needed to provide a surrogate for the fluctuations, such that u u u mod ( ) ¢ = + can be used in the above equations. In the following discussion, the index mod is dropped and all quantities with a prime, like u ¢ , are understood to be models for the unavailable true fluctuations. Observe that in Eqs. (10), (13), and (14), the mean velocity is employed. This will be addressed in Section 3.3.

Different approaches to dispersion modeling
According to Bocksell and Loth (2001), models for generating turbulent fluctuations can be divided into discontinuous random walk (DRW) models, continuous random walk (CRW) models, and stochastic differential equation (SDE) methods in ascending order of computational effort. The first discontinuous random walk (DRW) model was proposed by Gosman and Loannides (1983) and later extended to other variants like the one of Ormancey and Martinon (1984), for example. Turbulent dispersion in DRW is a succession of interactions between a bubble and eddies with finite lengths and lifetime. It is assumed that a bubble can be captured by an eddy when it moves through the fluid. This phenomenon induces a random instantanuous velocity component, which is piecewise constant in time. After the eddy lifetime is reached or when the bubble escapes from the eddy, a new random velocity fluctuation will be generated.
Compared with DRW models, progress in predicting turbulent dispersion especially in non-homogeneous turbulence can be obtained with the continuous random walk (CRW) models, which produce velocity fluctuations that are continuous in time (Sommerfeld, 2001). To account for the non-homogeneous drift of the diperse phase, a drift correction term associated with gradients of the Reynolds stresses or the turbulent kinetic energy is included in the approach (Bocksell and Loth, 2006).
The SDE approaches have many properties in common with the CRW methods, like the continuous velocity fluctuations in time. But a major difference is that the SDE technique does not demand additional drift corrections to account for the effects of inhomogenities (Reynolds and Lo Iacono, 2004). Instead, SDE models require models for the triple-velocity moments which makes them more complex compared with other two. This technique is still in the development stage, because some problems need to be further investigated, such as, e.g., the robust specification of particle-wall interactions (Bocksell and Loth, 2001). For these reasons the approach is not widely used in engineering applications.
Despite their deficiencies, CRW models have been shown to provide realistic predictions of turbulent bubble dispersion (Laín et al., 2002) and to reach a good compromise between computational effort and accuracy. Hence, a CRW model will be adopted in this work.

Normalized continuous random walk model
With the CRW approach, a separate Langvin equation is employed for each bubble to model the liquid velocity fluctuation u ¢ seen by this bubble. The normalized CRW model of Bocksell and Loth (2006): for the y-component, and so on. The first term on the right-hand side of Eq. (16) represents the correlated part and the second term the random contribution to the velocity fluctuation. The correlation function has an exponential form: . In Pope (2000), τ 0 952 C » .
(due to the different formulation of L τ in that reference, the value given there was converted) was found to provide a good estimation in a turbulent thermal wake experiment. But as the Reynolds number of the flow increases, smaller τ C down to 0 33 . gives better results (Dehbi, 2008). In the present work τ 0 4 C = . is employed. The last term D in Eq. (16) is the drift term, which accounts for the inhomogeneity of turbulence and avoids non-physical dispersion. The theoretical derivation of the full expression of the drift term goes back to Bocksell and Loth (2006). Under the assumption of isotropic turbulence, the expression provided by these authors can be simplified to is employed. The numerical implementation of the CRW model (16) leads to where the superscripts n and 1 ndenote the current and the previous time step, respectively.
Differing from the normalized CRW (22), the conventional CRW model assumes n n σ σ -/ = on the righthand side. According to Bocksell and Loth (2006), in the case a tracer particle is heading to a wall, the conventional form cannot correctly decrease the "history" effect, i.e., reduce the first term on the right-hand side. This may cause a non-physical drive of particles into the wall and result in a large amount of wall collisions (Bocksell and Loth, 2006). Therefore, a normalized CRW model with f b ( ) is employed here.

Employing the mean velocity field to determine forces
The models for the forces acting on the bubbles, i.e., the expressions in Eqs. (5), (7), (10), (13), and (14), contain the fluid velocity. In some of the terms the velocity u u u ¢ = + is employed, in other terms only the mean velocity u . This is discussed now. First, the term u u u ¢ ´= ´+ ´ cannot be evaluated, because u ¢ is defined in a Lagrange manner along the trajectory of each bubble. Transferring in each time step the information from the set of bubbles to the Euler grid, where the curl could be evaluated, would involve overwhelming computational cost. Most of all, the Lagrange model for u ¢ is exempt of information on the spatial correlation of realistic turbulent fluctuations, so that spatial derivatives cannot be captured and would be entirely unrealistic if modeled in this way.
Second, Eqs. (13) and (14) contain the material derivative along the bubble trajectory so that is preferred, because the derivative in time would generate unrealistically large values due to the random number ξ employed in Eq. (22), without improving the realism of the model.

Previous studies and requirement in this work
Beyond transport and turbulent dispersion turbulenceinduced coalescence is the most relevant phenomenon for the bubbly flows studied in this work. The random bubble motion described by the CRW model leads to collisions, which are the prerequisite for subsequent coalescence.
One major challenge related to the collision modeling is the expensive search for potential collision partners in the Euler-Lagrange representation. If there are b N bubbles in the domain, computation of binary collision is an 2 b ( ) O N problem, which strongly restricts the simulation to small bubble numbers. Sommerfeld (2001) proposed a stochastic collision model for particulate flows. Instead of searching collision partners, the momentum change due to collision is computed by generating a fictitious partner and assuming that the collision of a real particle takes place with the fictitious partner. However, from a physical point of view, contact and collision of real bubbles is the premise of coalescence. Modeling of bubble coalescence based on a fictitious collision partner is non-trivial. For these reasons, deterministic collision modeling proposed by Breuer and Alletto (2012) and Shams et al. (2011), for example, is preferred here.
To improve upon the time-consuming search, Shams et al. (2011) employed a node-based linked-list of bubbles for a cell-centered mesh. First, every cell corner node, e.g., some red points in Fig. 1, on the Euler mesh is assigned with an empty list. Then a loop going through all bubbles will find the nearest node to the center of that bubble. The corresponding bubble index will be appended into the list of that node. Collision detection is then carried out for bubbles inside the same list. This procedure reduces the number of bubbles that are searched. However, if the bubble size is of the same order as the grid size or if cells of the Euler grid are strongly anisotropic, this approach may fail to identify some collision events as illustrated in Fig. 1.   Fig. 1 Node-based linked-list may fail to detect collisions in certain cases. The arrows indicate the nearest nodes (in red) for corresponding bubbles. Breuer and Alletto (2012) developed an efficient collision model based on their tracking algorithm for the Lagrange points (Breuer et al., 2006). The computational domain is divided into secondary virtual cells. The size of these cells is chosen based on some criterion so that an appropriate number of particles are located inside. Substantial computational savings are achieved by restricting the collision detection to particles in the same virtual cell or in the neighbouring cells. Similar virtual cell methods are employed in smooth particle hydrodynamics or granular flows (Hopkins and Louge, 1991) as well. The computational efficiency of such methods largely depends on the selected virtual cell size (Rousset et al., 2018). It could be delicate to choose the cell size in such a way that the number of bubbles per cell is sufficiently low, in particular if the bubble sizes vary substantially in time and space due to coalescence, which is the situation considered in the present study.
In the present work, an efficient algorithm which can be applied on arbitrary cell types, cell sizes, and bubble sizes is demanded. To do so the sweep and prune (SaP) approach proposed by Baraff (1992) is employed to model bubblebubble interactions. Rousset et al. (2018) have evaluated and compared several different collision detection algorithms for the Discrete Element Method (DEM). It was found that the choice of the best algorithm is a trade-off between many criteria including the size of the search space, the number of particles and the memory usage. For the situation considered there, the SaP method was found to yield competitive execution time up to about 10 5 particles and the lowest memory consumption of all investigated algorithms. Baraff (1992) originally developed the SaP method to check the overlapping of multiple axis-aligned bounding boxes whose edges are parallel to the Cartesian coordinate axes. The basic idea is illustrated here with a two-dimensional example. Figure 2 shows three axis-aligned bounding boxes labeled A, B, C. At first, the coordinates of each box corner are sorted in ascending order, yielding the two vectors:

Sweep and prune with bounding boxes
in this 2D case. A loop sweeps through these sorted arrays to identify overlapping boxes according to the following criterion: If two boxes overlap with each other, like box B with box C in Fig. 2, 0 b x or 1 b x must be placed between 0 c x and 1 c x , or vice versa. This criterion must be fulfilled also for the y-coordinates of the boxes. On the contrary, boxes A and B do not overlap, because no intersection can be found with respect to the y-axis. The same is done in the three-dimensional case, just with a third vector of corner coordinates.
Since it only works on one-dimensional vectors of length proportional to the number of elements to be checked, and since fast sorting schemes can be employed, the SaP method requires only operations. For three-dimensional problems the original SaP was improved by Capannini and Larsson (2016) who proposed a succinct data structure to store the information of overlapping boxes during the sweep process. This algorithm is selected here for accelerating the computation of bubble interactions.
The link between box overlapping and collision detection is established by a bounding box around each bubble with edges in the direction of the Cartesian coordinates. During the simulation, such a bounding box is determined for each bubble based on its position and diameter. Then, with the help of the SaP algorithm, the collision of boxes is identified. If box overlapping is confirmed, the distance between the two bubbles involved is computed based on their true, spherical shape.

Soft-sphere collision model
After the determination of the distance between the bubble surfaces a collision model has to be applied. For the present work this is a soft-sphere model. Such a model is based on a repulsive spring-like force. As a result, it usually allows a very small overlap of the two bodies with the force, depending on the amount of overlap, such that the bubbles are driven apart (Fig. 3). The same is done if more than two bubbles overlap with each other. Here, the repulsive force of the model of  is used ① , with F n 2 coll eq coll eq eq In the equation, D is the overlap distance shown in Fig. 3, and eq b1 b2 b1 b2 is the equivalent bubble radius computed with the radii of the colliding pair. Furthermore, σ is the surface tension coefficient and n coll a normal vector pointing in the direction from the center of the collision partner to the current bubble center. The complete model of  also contains viscous contributions which are neglected here, because clean bubbles are considered.

Condition of bubble interface
During the bubble collision, two bubbles may merge with each other if the thin film of liquid trapped between them is drained out under the influence of external driving forces. According to the review of Chesters (1991), various regimes of drainage may be distinguished depending on the conditions of the bubble interfaces. Important regimes of film drainage are those of immobile, partially mobile, and fully mobile interfaces. The surface of clean bubbles is often assumed to be fully mobile, while the surface of contaminated bubbles is considered immobile (Chesters, 1991). From a physical point of view, this means that at the corresponding surfaces, slip or no-slip boundary conditions hold for clean or contaminated bubbles, respectively. Fully mobile interfaces are known to accelerate coalescence by orders of magnitude (Chesters, 1991).
The surface mobility can be decreased by the presence of different kinds of surfactants in the fluid. Besides, the material properties of the interfaces like Gibbs elasticity can play other important roles in the bubble coalescence (Danov et al., 1999). In this work, clean bubbles with mobile surfaces are assumed. Liao and Lucas (2010) listed three kinds of models for computing coalescence: (i) the energy model, which is related to the kinetic collision energy and the interfacial energy of the bubble, (ii) the critical velocity model based on the relative approach velocity before collision, and (iii) the film drainage model. The last one determines the possibility of coalescence from two characteristic time scales, the contact time and the drainage time. The drainage time is the time required for the thin film between the two bubbles to reduce its thickness down to a critical ① Note the recently published corrigendum to this paper (Heitkam Et Al., 2020). thickness. The contact time is obtained from the velocity of the bubbles relative to each other. Because the film drainage model is derived from the physical mechanisms of bubble coalescence and can be extended to different interface conditions, many varieties of this model have been proposed in the literature (Chesters, 1991).
In the following the two film drainage models used in this work are briefly presented. Based on the bounding box method introduced before, the respective coalescence criterion is evaluated when overlapping between two bubbles is found at the current time step but not found at the previous time step. This is to ensure that during each collision the coalescence is determined only once and only at the first contact of the two bubbles. The reason is as follows. The whole film drainage process for clean bubbles takes place over a very short lap of time and involves complex variations of the bubble shape which cannot be resolved by the point mass Euler-Lagrange approach. For this reason, the two coalescence models determine the probability of coalescence based on the relative collision velocity at the first contact without attempting to resolve the film drainage process. Thus, when bubbles keep overlapping with each other in several subsequent time steps in the present approach, the coalescence criterion should not be evaluated repeatedly.
If no coalescence takes place, bubbles stay in touch and the elastic collision force will be added to the right-hand side of the force equation (4). Kamp et al. (2001)

Model of
The interaction time i t is defined as the interval between the onset of film formation and the moment at which the bubbles begin to separate in a collision process due to the elasticity of the bubbles. Without viscous dissipation and body forces during the collision, a balance between the increasing surface free energy and the corresponding reduction in the kinetic energy of the systems can be established. was found in the case of two equally sized bubbles just touching each other (Kamp et al., 2001). For simplicity, vm 0 8 C ¢ = . is applied here throughout. The drainage time d t in the model of Kamp et al. (2001) is obtained from the numerical solutions of Chesters and Hofman (1982) for the case of two colliding bubbles with constant approach velocity yielding: where 0 u is the normal component of the relative velocity at the start of the collision.
In many instances the assumption of a constant approach velocity in deriving d t is invalid. Thus, the final formulation of the coalescence time scale ratio d i t t / is multiplied with an additional correction factor 1 k to take into account the effect of various approximations made in the derivation (Kamp et al., 2001), thus leading to The coefficient 1 k was found to be equal to 2 5 . when comparing to experimental results for bubble coalescence with a free surface in highly purified water (Kamp et al., 2001). Mattson and Mahesh (2012) applied the model of Kamp et al. (2001) to simulate bubble coalescence using an Euler-Lagrange approach with LES. They suggested to compare c P with a random number ζ uniformly distributed in the interval [0, 1] to determine if coalescence occurs when a collision is detected. The condition proposed is that c ζ P < , with otherwise collision only and no coalescence.

Model of Hoppe and Breuer (2018)
Hoppe and Breuer (2018) have developed film drainage models in the context of the Euler-Lagrange approach based on the work of Jeelani and Hartland (1991). Their coalescence model for clean bubbles is distinguished from the one of contaminated bubbles. In the case of clean bubbles, the film drainage is an inertia-controlled process, where the viscous force during the film drainage is neglected.
One advantage of this model is the consideration of a temporally evolving contact surface to avoid the simplification of a constant contact area. The varying contact surface f A is accounted for in the force balance of the bubble motion during the collision by setting: The drainage equation describing the reduction of the film thickness h for clean bubbles employed in the model reads: where h is the film thickness. By inserting the initial film thickness i h at 0 t = and the final film thickness f h at c 2 t t = / and integrating Eq. (32) in time, an expression The model of Hoppe and Breuer (2018)  .´ corresponding to water molecules according to Israelachvili (2011) and the surface tension was set to 0 072 N/m σ = .
which is representative of air-water interfaces (Hoppe and Breuer, 2018).

Comparison of coalescence models
The probabilities of coalescence delivered by both models in the case of air bubbles in pure water at room temperature are illustrated in Fig. 4. The stochastic contribution in the model of Kamp et al. (2001) leads to a continuous variation of the probability coal P . In contrast, the model of Hoppe and Breuer (2018) is based on a critical state, thus leading to a binary result of the model. On the other hand, it is exempt from any random number generation. For both models, as the approach velocity 0 u or the diameter eq D increases, the coalescence of clean bubbles may be inhibited. Indeed, in the experiment of Orvalho et al. (2015), it was found that the bubble contact time monotonously decreases with the bubble approach velocity obeying a power law. Further increase in 0 u can lead to a regime without coalescence. Ribeiro and Mewes (2006) recorded many random bubble collisions in polydisperse dilute bubble swarms suspended in a 2D flat downflow cell. Two classes of collisions were observed, those leading to coalescence for low 0 u and those with pure collision and rebouncing for high 0 u . Small values of eq D are found to facilitate bubble coalescence. According to Orvalho et al. (2015), the reason is that small bubbles with large curvature and less deformable surface can quickly squeeze out the liquid film. On the contrary, between two large bubbles the liquid film occupies more area resulting in a longer and more complex drainage process.

Three bubbles
The probability that three bubbles collide is very low in the flows investigated here. For this reason a simple strategy is selected for efficiency, which is sufficient here. If two bubbles stay in contact without coalescence and a third bubble collides with one of them, coalescence is not initiated.

Flow solver and Lagrange solver
The governing equations of the liquid phase were solved using the code Ansys Fluent v19.2 based on a cell-centered Finite Volume method for unstructured grids. The spatial derivatives were discretized using a third-order MUSCL and the pressure interpolation scheme was PRESTO!. The convergence criterion for all equations was to require that the scaled residuals have decreased to 5 10 -. The disperse bubble phase was computed by solving Eqs. (3), (4), and (22) using the in-house solver bTrack (standing for "bubble tracking tool") developed in the present study. It was implemented in C and parallelized with OpenMPI. The bubble velocity was obtained with the 2nd order implicit Broyden method. Only one-way coupling is considered in this work. The steady solution of the liquid phase is exported from Ansys Fluent to bTrack, where the transient bubble motion and interactions are calculated, as illustrated in Fig. 5.

Model verification with well-mixed criterion
To verify the implementation of the CRW dispersion model in bTrack, the well-mixed criterion is used. It states the requirement that very low inertia particles which are initially uniformly mixed in the domain should remain well mixed as time evolves and is commonly employed in the literature to verify dispersion models (Thomson, 1987;Dehbi, 2008). The random motion of tracer particles from the CRW model should fulfill this criterion in turbulence.
The numerical setup for the verification case is given in Fig. 6. The domain is periodic in all three directions for the particles and gravity is absent. The flow field with its statistics is given explicitly with turbulence assumed to be isotropic and inhomogeneous, so that k and ε suffice for its definition. Both are non-uniform according to the relations given in Fig. 6.
A number of 60,000 tracer particles with very low mass are uniformly and simultaneously released on a line along the x-axis in the middle of the domain at the start of the simulation at 0 t = . Hence, the particles need to travel in y and z to reach the statistically steady state. The time step size for computing CRW was chosen to be   Figure 7 shows the particle distribution in the domain at 15 s t = resulting from these simulations. Beyond this time, the average distribution of the particles remains unchanged. Figure 8 reports the normalized particle concentration obtained from averaging in y and z. The label "no normalization" stands for a simulation with in the first term on the right-hand side of Eq. (22). The drift correction refers to the last term on the right-hand side of Eq. (22). The effect of different terms in the CRW equation is clearly seen in these graphs. It is apparent that only when employing both the normalization and the drift correction the well-mixed criterion can be satisfied. If only the drift correction is used, the tracer particles still tend to migrate to the area with low turbulent kinetic energy.

Model validation
6.1 Overview of the cases studied Table 1 shows a summary of the cases employed for validation of the present bubble dispersion and coalescence modeling. These cases are based on the experiments of   Colin (1990) and Colin et al. (1991Colin et al. ( , 2012. In the experiment of the case 0, coalescence was inhibited by adding 3 0 0016 mol/m . of sodium dodecyl sulphate to the tap water. Therefore, coalescence modeling is switched off in the case 0 and the focus is placed on the validation of the bubble laden dispersion in the one-way coupled simulation. In the absence of surfactant (cases A to E), it was observed that the bubble size drastically increases downstream (Colin et al., 1991). Because no break-up was observed and the changes in bubble size due to mass transfer or pressure variation were found negligible in the experiment (Kamp et al., 2001), the effect of coalescence dominates here. Thus, bubble coalescence due to the turbulent fluctuations is studied in the cases A to E and the assumption of clean bubbles for the coalescence models discussed in Section 4.2 is not violated.
It might be argued that the transport of the surfactant at the bubble surface could promote the Marangoni effect so that the drag model (6) for clean bubbles might not be suitable in the case 0. According to Colin et al. (2012), however, this concentration is still smaller than the minimal concentration at which the Marangoni effect is observed for a 2 mm diameter air bubble rising in quiescent water (Duineveld, 1995). This justifies to employ the coefficient according to (6) for all cases. Table 1 Overview of the physical cases studied in this work. The superficial velocity of liquid, l j , and of the gas, b j , are defined by the respective volume flow rate divided by the cross-sectional area. The average void fraction for the case B is not provided in the reference  Colin (1990) All the experiments were carried out in a pipe of the same diameter, 40 mm D = . Gas was injected and mixed with the liquid phase through a Venturi mixer at the inlet. For the case 0, the local measurement of gas volume fraction and axial bubble velocities were performed near the pipe outlet. The bubble size spectra in the cases A to E were obtained by an optical measuring procedure with a camera at 500 frames/s. The experiments were performed in microgravity conditions on an aircraft that flew a parabolic path to achieve a gravity smaller than 2% of the earth's gravity for 20-30 s. The most important difference between these conditions and an up-or downward bubbly flow on earth is the fact that the bubble-induced turbulence (BIT) is negligible, so that the flow is dominated by the shear-induced turbulence alone (Colin et al., 2012). As a result, the flow is unaffected by BIT, giving a clearer picture and removing the need for any BIT models in the present work.
The same Euler mesh was used for all cases. Information on turbulence modeling, boundary conditions, mesh and solver parameters is assembled in Fig. 9. The symbol x D and w r D stand for the grid spacing in x and the height of the first cell adjacent to the wall, respectively. For the liquid phase, the properties of water at room temperature are used. The pipe Reynolds numbers of the single phase flows D l l l Re j ρ D μ = / are listed in Table 1.

Single phase simulation
The single phase simulation with the k ε realizable model was carried out for the case 0. Because the inlet and outlet were defined as periodic, the single phase simulation after convergence generates the same velocity profile in each cross section. A mesh independence study showed that this grid of 2 46 . M grid points provides the same results as a grid with 8 3 . M cells. Validation of the simulation against the experimental measurement is shown in Fig. 10. The RMS of the liquid velocity fluctuation was computed as 2 3 u k ¢ = / . Good agreement between the experimental data and the single phase simulation is observed. The liquid data in the single phase experiment and the same data in the multiphase experiment are both plotted in the figures. Comparing these, the difference is  found small which means that the presence of the bubbles has not significantly affected the liquid phase velocity, as stated in Colin et al. (2012). In the absence of buoyancy, the bubble slip velocity was smaller than the measurement uncertainty in the experiment. Due to a near zero slip velocity, the momentum transfer from the bubbles to the liquid phase is negligible in this flow. Based on this observation, the assumption of one-way coupling is valid in the present work.

Bubble dispersion without modeling of collisions
Based on the single phase steady RANS results, the random bubble motion was computed as a next step using the CRW model. Bubbles with a constant diameter b 1 2 mm D = . were released uniformly at the pipe inlet at 0 x = every 0.01 s . The number of released bubbles was determined according to the experimental gas volume flow rate. During the simulation a maximum of about 120,000 bubbles were simultaneously present in the domain. The Lagrange time step size L t D was set to 0.00025 s and the total physical simulation time was 10 s .
The bubble volume fraction b α was sampled along the radius, with averaging in circumferential direction, at 2 8 m x = .
in each time step. The bubble velocity in the axial direction was also collected every 0.01 s at this position to compare the mean velocity fluctuations of bubbles with the experiment. The bubble dispersion near the outlet is illustrated qualitatively in Fig. 11, where the bubbles are scaled to their real size. Recall that bubble-bubble collisions were not accounted for in this simulation. Figure 12 illustrates the time averaged bubble volume fraction b α along the radius at 2.8 m x = with averaging in circumferential direction. The standard model employs  the forces described in Section 2.2 and the dispersion model of Section 3. The simulation results show that with the modeling employed in this study the radial distribution of the bubbles can be predicted very well for this case. Indeed, if no lift force model is employed, there is an erroneous peak of volume fraction near the wall as seen in Fig. 12. An explanation for this may be that the turbophoresis, as the force balance between the drag force and the effect of turbulence, leads to a buildup of concentration near the wall. Comparing the results from the two lift models it is seen that the particular value of the lift coefficient has no large influence on the results. If no dispersion model is used, peaks of b α appear in the figure, because bubbles near the wall are constantly pushed towards the pipe axis by the lift force and they may overlap with each other as collision modeling is neglected. The lift force decreases as the bubbles move from the wall to the center, so that the shift of the peak in b α is also becoming slower with increasing streamwise position. Figure 13 shows the results for b α and for the mean velocity fluctuations of bubbles in the axial direction b,x u ¢ with different values of the model parameter τ C in Eq. little. But overall the fluctuation is somewhat underestimated. The bubble dispersion based on the liquid phase solution with the Reynolds stress model was also computed. As mentioned in Section 3.2, for a RSM, the mean velocity fluctuation f σ can become anisotropic in Eq. (22). The liquid phase velocity fluctuations in the three directions obtained from the RSM simulation are shown in Fig. 14. Nearly no difference between them is observed, except in the near-wall region. The anisotropic effect, hence, plays a minor role in this case.
Overall, Figs. 12 and 13 demonstrate that in this validation case the bubble distribution can be well estimated by the CRW model. The underestimation of b u ¢ is caused by the fact that bubble collisions are neglected in the simulations discussed so far here.

Bubble dispersion with modeling of collisions
Figure 13 contains one result with collision modeling activated. It is apparent that the bubble distribution in the domain is not affected by this step. However, the statistics of bubble velocity fluctuations are improved markedly. Finally, it is worth taking a glance at the performance enhancement by the SaP method presented in Section 4.1.2 reported in Fig. 15. The SaP implementation used is only partly parallelized. It consists of parallelized sorting of box coordinates, sequential sweep and prune on the root CPU according to Capannini and Larsson (2016), and parallelized force computing. More improvement can be achieved by parallelizing the sweep and prune step using, e.g., the algorithm from Capannini and Larsson (2018). But already  First, the single phase simulation for the case A is validated. Unfortunately, Colin et al. (1991) did not provide measurement data of the single phase flow, so that the experimental data of a pipe flow at a similar Reynolds number from Lawn (1971) are adopted here for comparison. In this reference, the anisotropic fluctuating velocities of the single phase flows, when normalized by the wall friction velocity τ u , were found invariant with respect to Reynolds number Re from around 30,000 to 250,000.
The continous phase of the case A was computed with the Reynolds stress model and the results are plotted in Fig.  16. The computed liquid velocity fluctuations are normalized by τ u obtained from the simulation. The x-axis represents the streamwise direction. Basically, good agreement between the single phase simulation and experiment can be regarded.

Results of simulating bubble coalescence
In the simulation, a constant Lagrange time step size L t D = 0 00001 s . for the disperse phase was employed over the total simulation time of 20 s. Bubbles were released randomly at the inlet every 0.006 s with their number chosen to match the gas flow rate in the experiment. The bubble size was randomly chosen from the experimental probability distribution function (PDF). Coalescence modeling was switched off until the bubbles moved to 10 x D > . This was done for two reasons. First, this inflow region allows the bubble PDF to equilibrate. Second, the experimental bubble size PDF at the inlet was actually measured at 10 x D = of the simulation domain. The simulated bubble size near the outlet in the region 70 80 D x D < < was determined to compare with the corresponding experimental data at the outlet.
An overview of the bubble dispersion and the size distributions is given in Fig. 17. It is apparent that quite large bubbles can be obtained, with b D up to 15 mm in the case A, when bubble injection is performed with an initial mean diameter of 2.48 mm which is a factor of 6 . Furthermore, this amounts to b D D / as large as 0 375 . . An influence of the Reynolds number and the average void fraction is seen when comparing the case A and C. The snapshots also reveal that, as expected, the bubbles are not distributed uniformly over the cross section.
A more quantitative analysis of the results is now provided with appropriate statistical data. Results of the case A are shown in Fig. 18. The bubble size PDF in the simulations near the outlet was computed by the 1D kernel density estimation function in Python (SciPy, 2019) and it was evaluated for bins of width 1 mm , centered at b D = 1 5 2 5 . , . ,... to 20.5 mm , and normalized to 100% . The results show that both coalescence models can reproduce the size distribution for the case A fairly well, with the model of Kamp et al. (2001) yielding slightly smaller bubbles. A simulation result for the same case obtained with LES by Mattson and Mahesh (2012) employing the model of Kamp et al. (2001) is also plotted as a reference. These data are very close to the present RANS result. The computational effort, however, is much larger with the LES approach, so that the present method is more cost effective for this case. Figure 19 illustrates the bubble size PDF at the outlet, shown in Fig. 18(b), by displaying a snapshot of bubble positions and velocity (not size, as in Fig. 17) near inlet and outlet. Very large diameters up to 2 D / are observed, with bubbles concentrating near the center. A diameter of this magnitude yields substantial velocity differences of the fluid between the wall-oriented side of the bubble and the center-oriented side. Hence, one might expect that bubbles of this size undergo deformations of their shape which would have to be accounted for in the simulation. However,   in the microgravity experiments of Kamp et al. (2001) and Colin et al. (1991), bubbles stay nearly spherical up to bubble diameters of 20 mm . This can be seen also in the small Eötvös number for bubbles. If the experiment was performed at a gravity equal to 1% of the Earth's gravity, the Eötvös number for bubbles of 20 mm size would only be 0 54 . . So the approximation of a spherical shape in this study is acceptable. The bubble size distribution near inlet and outlet in the cases B and C are presented in Figs. 20 and 21, respectively. In these cases coalescence is not so pronouced, resulting in a lower average probability of large bubbles than in the case A. Cases B and C have the same liquid flow conditions and actually the same gas volume flow rate. The mean bubble size at the inlet, however, differs by more than 25% resulting in a different size at the outlet. In both cases the model of Kamp et al. (2001) yields smaller bubbles, with the model of Hoppe and Breuer (2018) generating the same maximum size as observed in the experiment. The experimental data  are narrow, so that the maximum value is more pronounced. Figures 22 and 23 show the predictions for the cases D and E. Similar to the case A, the simulation results at the outlet display a smooth variation from smaller bubbles to larger ones. In contrast, the experimental data demonstrate irregular shapes in the PDF. As discussed in Kamp et al. (2001), a log-normal law should be expected to represent the PDF of bubble diameters and the irregularity might be caused by accuracy issues of the measurement principle. It was reported that for both cases only around 100 bubbles were recorded for analysis during the measurement. The differences observed with both coalescence models is comparatively small and exhibits the same tendency as noticed before with the Hoppe and Breuer (2018) model yielding slightly larger bubble sizes.
Several parameters were calculated to characterize the bubble size distribution near the outlet in Fig. 24. The characterizes the standard deviation of the size distribution. Figure 24 displays the correlation of the parameter values obtained from the simulations with the corresponding experimental data. Good agreement between predicted and measured statistical values is obtained for both models. With respect to the difference between both coalescence models, the model of Hoppe and Breuer (2018) yields slightly larger bubble sizes than those from Kamp et al. (2001). Nevertheless, their outcomes are quite similar in the studied cases.

Conclusions
A simulation framework for the computation of disperse bubbly flows at low Eötvös number was presented which can be employed with Eulerian grids for the liquid of arbitrary cell size and shape. Specially, a CRW Lagrange dispersion model was employed to model turbulence with time-averaged fluid parameters and found to be suitable and efficient. By applying the sweep and prune method to detect bubble collision, the overall performance was enhanced furthermore.
Simulations of bubble collision and coalescence in turbulent pipe flows were undertaken for validation purpose. Computed bubble size distributions from the two coalescence models of Hoppe and Breuer (2018) and Kamp et al. (2001) were in good agreement with the experimental measurement. All in all, the proposed combination of physical and numerical models was successfully validated. Compared with the model in Mattson and Mahesh (2012), the current framework was able to predict bubble dispersion and collision with a stationary primary phase solution, which is proved to be an efficient simulation strategy. Further investigations will be concerned with applying and assessing the proposed modeling framework to more complex configurations. With regard to the numerical development, further focus will be placed on the two-way coupling. Large bubbles appear due to coalescence and their disturbance of the liquid phase may not be negligible. However, in the two-way coupled Euler-Lagrange simulation with point mass assumption, the accuracy and stability strongly depend on the ratio of mesh to bubble size. A point source assumption may not be valid for finite size bubble whose size is larger than the local grid size. This requires a proper interpolation scheme for finite size bubbles and appropriate distribution of forces, both being topics of ongoing work.

Appendix
In the present study turbulence is described by two different turbulence models. One is the two-equation realizable k ε model based on the transport equations: where the terms P k and P ε represent the production of k and ε , Y k and Y ε the dissipation terms of k and ε , and k σ and ε σ the turbulent Prandtl numbers, respectively (Fluent Theory Guide, 2019).
To account for anisotropic effects of turbulence the Reynolds stress model (RSM) with linear pressure strain ( ) ρ φ ¢ ¢  ⋅ = + -+ uu u D P Y (A.3) was adopted. In this equation D , P , and Y are diffusion, production, and dissipation term, respectively. The pressurestrain term φ is modeled by a linear combination of different pressure-strain effects (Fluent Theory Guide, 2019). Bubble-induced turbulence is not considered in this study, since the flows investigated are micro-gravity flows, where this effect is known to be small.

Funding note
Open access funding provided by Projekt DEAL.