Flow simulation in a 2D bubble column with the Euler-Lagrange and Euler-Euler method

Bubbly flows, as present in bubble column reactors, can be simulated using a variety of simulation techniques. In order to gain high resolution CFD methods are used to simulate a pseudo 2D bubble column using EL and EE techniques. The forces on bubble dynamics are solved within open access software OpenFOAM with bubble interactions computed via Monte Carlo methods. The estimated bubble size distribution and the predicted hold-up are compared to experimental data and other simulative work using EE approach and show reasonable consensus for both. Benchmarks with state of the art EE simulations shows that the EL approach is advantageous if the bubble number stays at a certain level, as the EL approach scales linearly with the number of bubbles simulated. Therefore, different computational meshes have been used to also account for influence of the resolution quality. The EL approach indicated faster solution for all realistic cases, only deliberate decrease of coalescence rates could push CPU time to the limits. Critical bubble number - when EE becomes advantageous over the EL approach - was estimated to be 40.000 in this particular case.


Introduction
Bubble columns can be found in a various set of operations in chemical and biological applications. In order to reach a high efficiency, an optimal bubble size distribution (BSD) in any reactor type has to be achieved, which markedly depends on the local hydrodynamics (Akita and Yoshida, 1974). Computational fluid dynamics (CFD) simulation is one possibility of obtaining time and space resolved information in respect to complex bubbly flows. In general, they can be simulated in various ways, from detailed resolved single bubble behavior up to large-scale averaging methods with reactors with billions of bubbles inside. Special classes are based on particle population balances (PPB) relying on different break-up and coalescence kernels (B&C), that have been developed in the last few years to characterize the particle interactions Lucas, 2009, 2010;Martıńez-Bazán, 1998).
This work focuses on the Euler-Lagrange (EL) approach, giving a high level of detailed information on the bubble scale. In contrast to Euler-Euler (EE) simulations where the dispersed phase is considered as pseudo-continuous, every single bubble is resolved and can be tracked on its way through the domain. The EL approach makes use of stochastic modeling and improved collision algorithms to compete against other simulation techniques. All results are compared to solutions coming from state of the art EE simulations. Further, we investigate the computational cost of both (EL & EE) and identify critical parameters.

EL modeling
The EL approach for simulation of bubbly flows can be ranked between the direct numerical simulation (DNS) and the EE approach. While the EE model resembles the bubbles as a density distribution, in the EL model every single bubble is calculated as a point volume acting under Newtonian dynamics, while the movement of the surrounding fluid is solved via Navier-Stokes equations on an Eulerian grid. The internal coordinates of each bubble (size, concentration, etc.) are calculated using macroscopic models, which is the main difference to DNS simulations. However, the combination of this Eulerian grid and Lagrangian points yields a higher level of detail with a moderate increase in computational cost.
The original OpenFOAM solver sprayFoam provided the basis for the bubbly flow simulation shown in this work. sprayFoam uses an EL approach to simulate liquid or solid particles inside a continuous phase, e.g. a spray injected through a nozzle. The simulation of bubblesgaseous particlesrequire a rework of this original solver to account for the additional physical effects.

Continuous phase hydrodynamics
Continuous phase and dispersed bubbles can be calculated predominantly independent from each other, which allows a solution in two steps. To achieve coupling of phases in the second, the exchange of momentum and energy from the bubbles is calculated and condensed into the Navier-Stokes source term . With the continuous phase velocity , the pressure , the density and the viscosity , the modified Navier-Stokes equation reads as: Turbulence is computed with Reynolds-averaged Navier-Stokes (RANS) modeling, which requires significantly less computational effort than solving for a Large Eddy Simulation (LES) approach. In the RANS model, only the averaged turbulence properties, like turbulent energy and dissipation, are calculated, while the LES approach will give a more detailed shape and movement of the actual turbulent eddies in the domain. A single phase kepsilon turbulence model was chosen to calculate for the turbulence in the continuous phase and in addition to the original functionality, bubble induced turbulence (BIT) models have been added. The drag force, FD, of a rising bubble with velocity induces turbulent energy and dissipation into the surrounding phase. This can be calculated using a source term for the turbulent and dissipative energy equations: The open literature reveals different constants, , and turbulent time scales, , according to specific physical effects in the bubbly flow (Pfleger and Becker, 2001;Rzehak and Krepper, 2013;Yao and Morel, 2004). In the presented EL simulation, the model of Rzehak and Krepper (2013) was used.

Bubble Dynamics
Since the continuous phase solution has been vastly modified, the dispersed phase hydrodynamics yield a much larger rework of the original OpenFOAM solver. Nevertheless, the basic model of Lagrangian particle tracking persists. The bubbles are assumed to be points of mass at the position . Local mass can change due to break-up and coalescence events. The new position of any specific bubble is calculated discretely with the velocity .
Forces acting on the bubble will change its velocity. Here resembles the bubble's mass and Σ stands for the sum of all acting forces on a single bubble.
The sum of forces Σ consists of the buoyancy and weight force , the drag force , the lift force , the virtual mass force , the wall lubrication force and the dispersion force . Here the subscripts and stand for the bubble and the continuous phase accordingly, the subscript identifies the relative differences between them. Furthermore, g denotes the gravitational acceleration, stands for densities, u for velocity, for the bubble's diameter, for the turbulent kinetic energy and for the phase fraction. Adjustable model parameters are denoted with .
These forces are in particular: In (11), stands for the bubble's volume and / denotes the material derivative, meaning that the derivative is made while following the bubble. In (12) the wall force parameter is dependent on the distance and relative velocity of the bubble to the next wall, stands for the normal vector on the wall. In prior validation simulations (Weber and Bart, 2016) several models for the drag and lift force coefficients , have been tested. Based on this, the models of Tomiyama (2004) were chosen. The virtual mass coefficient is set to = 0.5 according to Delnoij et al. (1997), the coefficient for the dispersion force is set to = 0.1 (Lahey et al., 1993).
Besides the movement and position information, each bubble is carrying an extra set of inner coordinates, in our case age, diameter, orientation, temperature and concentration of a solute. Additionally, physical properties like density, heat capacity or surface tension are calculated, which change with the surrounding pressure, temperature and solute concentration. However, in this initial simulations chemical reactions and mass transfer are neglected (no solute), which is a topic for further studies.

Stochastic modeling
Many aspects of bubble behavior are simulated using stochastic Monte Carlo (MC) models. On the one hand, it resembles the actual random character of bubble movement or other unpredictable interactions, and on the other hand it yields an immense benefit in computational effort. The procedure is similar in all situations; a calculated probability is compared to a random number which will trigger the event or neglect it. Through the law of large numbers, this algorithm will converge to the expectation value step by step. This seems a bit oversimplified but it represents a reliable method to calculate very complex systems, without the need of high dimensional differential equations.

Bubble collision
The first bubble interaction, which is described through a stochastic model, is the collision between bubbles. Here the collision probability by O'Rourke (1981) is used. The actual position of bubbles is not considered, but rather the relative velocity and sizes of the bubbles are taken into account: Note, that the volume of the computational cell is also a parameter in this equation. This implies, that only bubbles belonging to the same cell can actually collide. Also, this has to be evaluated for every possible pair of bubbles in the interesting volume. For each pair, a random number ∼ (0,1) will be drawn. The constraint for collision will therefore be: Compared to deterministic models, where the precise trajectory of every bubble has to be predicted, this stochastic approach is much faster. Because the collision step is one of the largest computational efforts in the code, it is crucial to use an efficient algorithm. One promising approach shown by Sigurgeirsson et al. (2001) was customized to work also with the stochastic collision algorithm gaining a major speedup in comparison to the original code. The algorithm first creates a list of bubbles per each computational cell. In the second step, all combinations of bubble pairs for each cell are evaluated. Thereby, only bubbles in relative proximity will be observed. Moreover, the optimal cell size for this algorithm has been identified as the maximum bubble diameter. Smaller cells would presume a bubble, that has more volume than the enclosing cell, leading to numerical instabilities. However, a larger cell leads to more possible bubble pair combinations necessary to be calculated and slows down the process.

Bubble coalescence
After a successful collision, the next step is to calculate the coalescence probability. Previous validation simulations (Weber and Bart, 2016) have shown that the model of Coulaloglou and Tavlarides (1977) is an appropriate choice. Basis of the calculus is the combination of bubble contact time and film drainage time . The drainage time can be seen as the time it takes for a thin liquid film between two bubbles to flow out. Only when bubbles are in long enough contact their surfaces can merge and coalescence happens. The time for contact is evaluated based on turbulence energy, . A normal distributed turbulent energy is assumed, which gives the average time that the bubble gets hit by a turbulent eddy, thus pushing them away from each other again. Analog to the collision, a new random number will be used for the constraint of coalescence: Here represents the surface tension, hf and hi are critical and initial film thicknesses, respectively, and deq is the equivalent diameter defined as: If no coalescence happens, the bubbles will bounce of each other and continue their movement. Successful coalescence will lead to the deletion of one bubble, while its mass is added to the other one.

Bubble break-up
When a turbulent eddy with enough energy hits a bubble, it is deformed and eventually splits into smaller bubbles. The smaller a bubble is the more energy is necessary to deform it. Thus, the break-up frequency of small bubbles is low. This has also been modeled by Coulaloglou and Tavlarides (1977) into a breakup frequency, , and probability, . Again, a random number will be drawn and compared in each time step for each bubble.
When a break-up event has been evaluated as positive, a new bubble has to be created. This new daughter bubble has the same inner dimensions like its mother bubble. Only volume, diameter and mass will be changed on a basis of the volume ratio of the two bubbles. This is a controversial point in break-up modeling, since many researchers claim that should follow various different distributions. The model assumption of Coulaloglou and Tavlarides (1977) is a normal distribution of , where the ratio = 0.5 has the greatest probability. To sustain numerical stability, the creation of bubbles with an extremely small diameter is excluded by restricting 0.01 < < 0.99.

Turbulent bubble movement
In order to achieve chaotic bubble movement trough turbulent eddies, the bubble velocity is overlaid with a Brownian motion. The magnitude of this turbulent velocity is capped by the current turbulent energy. The direction is chosen uniformly random.
This will lead to a trembling movement of bubbles positioned in a high turbulence area representing a diffusional movement. The higher the turbulence, the faster the bubbles will spread. Note that this additional velocity is not coupled to the forces and will not affect the creation of turbulent energy or the probability of collision.

EE modeling
The common approach to simulate bubbly flows is the usage of multiphase system solvers. Here, the dispersed phase is also modeled as a continuum and no explicit bubble positions will be evaluated anymore. Thereby, it is not possible to track single pathways of bubbles or to even assign precise diameters to bubbles in a certain volume. Anyhow, when incorporating a suitable model for coalescence and break-up, a EE simulation can still give valuable information about the mean diameter, surface area and phase hold-up of the gas phase inside a bubble column. The EE simulation is also capable of simulating a much higher dispersed phase hold-up like it would be possible in the EL approach (Hlawitschka et al., 2016).
There are many subcategories of EE solvers, while this work will concentrate on the comparison of EL simulations to EE simulations, which are combined with different moment based PPB models. The method of moments (MoM) (Hulburt and Katz, 1964) describes the dispersed phase based on the characteristic moments of the BSD: number, diameter, surface area and volume. In principle, the moments movement is solved like an incompressible isothermal fluid, also respecting the conservation of mass and impulse. Further improvements of this basic method have led to the quadrature method of moments (QMoM) (McGraw, 1997) and direct quadrature method of moments (DQMoM) (Marchisio and Fox, 2005), where closure problems have been overcome. Those two models are solved with a multiple equation system. In contrast to that, the one primary one secondary particle method (OPOSPM) (Drumm et al., 2010) and interfacial area transport equation (IATE) (Ishii et al., 2005) are based on only one equation to be solved, which improves the CPU speed even further.
For the comparison, the simulations of Hlawitschka et al. (2016) and Buffo et al. (2013) will be used.

Bubble Dynamics
In contrast to EL simulations, no explicit bubble position or velocity is calculated in the EE approach. Instead of that, the dispersed phase moments are used and their velocity is supposed to be equal in a computational cell. Transport of the moments is similar to a scalar transport, no interface tracking/sharpening steps are used. Nevertheless, the bubble forces are included by using volumetric equivalents of the explicit equations. With the mean bubble diameter d3x -calculated via the tracked moments -the phase hold-up & and the mean relative bubble velocity the resulting forces can be calculated for each computational cell. For the drag force, the following equation is used.  (1956). Bubble induced turbulence and turbulent dispersion were neglected in the EE simulations, the standard mixture k-epsilon RANS model was used. The coupling of phases was realized with a source term similar to eq. (1).

Bubble interaction
In general, the interaction of bubbles can be calculated using the same model equations like in the EL simulation. The difference is a conversion from explicit Monte Carlo simulation of a single or a pair of bubbles to a discrete application on the moments. The calculated event probability is no longer compared to random numbers but will be converted to an event frequency instead. This implies that it is possible to get fractional values (e.g. 0.5 bubbles per computational cell), which is impossible in the EL approach. While bubble volume stays constant the combined loss and gain in bubble number, size and surface area is recalculated during the bubble interaction step. The source term for the population balance can thereby be written as: Where B stands for birth and D for death of bubbles, exponents C and B stand for coalescence and break-up respectively. While the coalescence model in Hlawitschka et al. (2016) is adopted from Coulaloglou and Tavlarides (1977), the break-up model originates from Laakkonen et al. (2006) which uses a similar approach based on the turbulent energy dissipation. Break-up of bubbles is supposed to be always binary and resulting daughter bubbles are supposed to be equal in size. Models for break-up and coalescence used by Buffo et al. (2013) are both adopted from Laakkonen et al. (2006). Here, the daughter bubbles created by break-up are supposed to follow a beta distributed volume ratio .
While in the EL approach every single bubble diameter is inserted into the MC simulation, the EE approaches use representative diameters for all bubbles in the considered computational cell. This can either be the Sauter diameter 32 , the mean diameter 30 or nodes of the quadrature points when using multi equation MoM models. This simplification enables EE simulation to be less dependent on the bubble number, making it a feasible calculation even with high gas hold-up.

Simulation case
There exist experimental data (Cachaza et al., 2009;Díaz et al., 2008b;Díaz et al., 2008a) and simulations (Hlawitschka et al., (2016), Buffo et al.,(2013) for a set-up described in Díaz et al. (2008b) consisting of a rectangular bubble column depicted in Fig. 1: In the original experiment, the column is filled with tap water at room temperature while air is injected trough a sparger plate at the bottom. The sparger consists of eight holes with 1 mm diameter. They are arranged in two rows with a distance of 6 mm between holes.
In order to compare the two types (EE vs. EL) of simulation, the original EE mesh (Hlawitschka et al., 2016) has been used for the EL approach in a first simulation. To further improve simulation time and solution quality, an optimized mesh is used. The EE mesh consists of 24640 hexahedral cells with a length between 4 and 6.7 mm. Cells around the gas inlet zone are finer to account for a better solution.
Since the EL simulation does not compute the free surface, the mesh is shortened to the filling height. The mesh refinement in the inlet area was turned off to achieve a more isotropic mesh. The resulting fine EL mesh consists of 11745 cells with a length between 6 and 7.5 mm.
The optimized coarse EL simulation case resembles the experimental volume within 20x4x45 computational cells, which gives uniformly sized cells of 1x1x1 cm. Note, that this coarse grid is inevitable in order to maintain a suitable ratio of maximum bubble size (≈ 1 cm) and the grid cell length. This ensures that bubbles are always smaller in volume than the surrounding computational cell and that the collision of bubbles within a cell performs well. Other simulative work was built on a comparable coarse mesh, while Díaz et al. (2008b) even detected a numerical deviation when using a finer mesh. A grid convergence analysis was therefore disregarded for our EL simulation. Unlike EE Simulations, there is no need to define an inlet patch. The bubbles are injected on predefined points, depicting the eight sparger holes at the bottom. In experiment and simulation different gas flow rates have been tested. Corresponding velocities and bubble injection sizes have been calculated with equations from Geary and Rice (1991) and can be found in Table 1. Additionally, the max+ case was simulated for the most promising models, although no experimental data exists for this case.

Results and discussion
Experimental and simulation results (s. Table 3) for Sauter diameter and hold-up are taken from Hlawitschka et al. (2016), Buffo et al. (2013) and (Díaz et al., 2008b). Note that in all cases, the hold-up is slightly under predicted in the EL in comparison to the EE simulations. The predicted Sauter diameter fits to the experimental values, except for the 'max' case, where it is markedly under predicted. The overall results of the EL simulation can be considered accurate enough in comparison with the experimental data and the other simulations. Values for the max+ case are not listed, because no experimental results exist. to ensure equal conditions. Also, the simulations were done in single core mode, due to issues with some of the solvers in parallel mode. For all models, a total of 120 s experimental time has been simulated. Table 2 shows the resulting CPU time (in hours) in direct comparison. As expected, the more complex two equation models QMOM and DQMOM need significantly more time, while the one equation models IATE and OPOSPM yield a faster solution. The EL approach could easily beat the two equation EE models, while the one equation models were faster at lower gas throughput. When switching to the coarse mesh, the EL simulation could achieve solution in an even shorter duration, beating the simulation time of all EE simulations. Note the strong increase in computational time for the EE simulations from the min to the med case while the increase from med to max is not that significant. Simulations for the max+ case showed no further significant increase in the computational time for higher gas velocities. Still, the EL simulation could achieve a faster solution for all listed cases and it does not seem useful to compare the models based on superficial velocity vs. CPU time.
A better approach (at least for the EL simulation) is the comparison of average bubble number and computational cost. I turns out that the computational time scales almost linearly to the average number of bubbles in the domain for the original coarse EL mesh (s. Fig. 3). Additional simulations were performed to prove this behavior even for a larger mesh. Therefore, the mesh dimensions have been doubled in every direction while resolution stayed the same, which gives a mesh with eight times more cells (mesh x 8 in fig. 3). The superficial gas velocity has been further increased from 21.3 mm/s to 42.6 mm/s and 85.2 mm/s for this larger mesh. To find the critical number of bubbles, where the EE simulation finally becomes faster than the EL simulation, the coalescence efficiency was successively lowered. This leads to smaller bubbles and thus a higher number of bubbles. With an average bubble number of 36.000, the achieved computational time for the EL simulation reached a value of about 14 h (s. Fig. 4, grey dashed line), which is the time that the OPOSPM solver needed. A further increase to a bubble number of about 50.000 led to a simulation time of 23 h, which is even slower than the IATE solver. Thus, at least for this case geometry, the critical bubble number can be pinned to 40.000.
As mentioned above, a major speed up for the EL approach was achieved by using an efficient collision algorithm. According to literature (Sigurgeirsson et al., 2001), the optimal collision algorithm would need ( ) steps to calculate for collisions of particles. Since this can never be reached, the quality of an algorithm has to be quantified in terms of exponential notation , where a smaller value for qualifies a better algorithm. Fig.3 compares the computational time of the original and improved algorithm that has been used. Lowering the complexity from 2.3 to 1.5 is a large step: doubling the bubble number would have meant a five times higher CPU time but will now take only about 2.8 times higher effort with the new algorithm. The new one shows a slower solution for bubble numbers under 100 due to memory preallocation but its order for higher numbers shows a significantly faster calculation. Taking into account that in a bubble column are usually more than 100 bubbles, it is clear that the new approach is advantageous.

Conclusions
The EE and EL simulations are capable of solving the flow inside a bubble column, showing similar results in comparison with the experimental data. The computational time needed differs a lot from approximately one hour up to 40 hours, where a higher gas throughput corresponds to higher computational time. This dependency is existent even for the EE simulations but much more observable in the EL simulation. The EL algorithm is almost linear dependent on the number of bubbles and tends to be very slow, when the number increases to higher ranges ( > 40.000 bubbles). Here it is clearly advantageous to use an EE approach for an efficient simulation in cases, when the gas hold-up is high ( > 5%) or when the bubble size is small (d32 < 4 mm). Still, there are many applications with a medium bubble size and/or low gas hold-up. In this case, the EL simulation technique can easily outperform the EE simulations in CPU time and level of detail. This is due to the scaling of the EL simulation to the number of bubbles, which is almost linear, while the EE technique has already high computational cost even at a relatively low hold-up. Computational cells with no bubbles present have almost no impact on CPU time in the EL approach. In the EE approach, even those cells have to be computed concerning transport and interaction of moments.
Further advantage of the EL approach is the easy visualization of trajectories and positioning of single bubbles on their way through the reactor. The EE simulation results are treated in post processing to reproduce such a meaningful visual representation although no "real" path of single bubbles were calculated. When concerning chemical reactions, detailed information about its temperature and concentration progress can also be combined with particle residence time, path length and other explicit single bubble properties within the EL approach.