Bubble columns can be found in a various set of operations in chemical and biological applications. In order to reach reliable operation, an optimal Bubble Size Distribution (BSD) in any reactor type has to be achieved. This markedly depends on the local hydrodynamics [1]. Computational Fluid Dynamics (CFD) simulation is one possibility to obtain 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 [2-4].

This work focuses on the comparison of Euler-lagrange (EL) and the Euler-euler (EE) approach. The EL approach yields a high level of detailed information on the bubble scale. In contrast to 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 EE approach is less dependent on the bubble number and can outperform the EL approach when large industrial scale reactors are simulated.

Although both simulation techniques build on fundamentally different assumptions, they still try to resolve the same issue, namely a transient bubble population. Be it a faster computation or a higher level of detail, each approach has its advantages over the other. This is why numerous references mention both different approaches and their application area [5-8]. A consensus of scientific work is that the principles of EL and EE simulations lead to a point, where one of the methods should be preferred in terms of computational time. This point is mostly dependent on the problem size, or rather on the bubble number and hold-up. Fact is that the simulation cost scales stronger with the bubble number in the EL than in EE approach. This is well known in literature, many authors state that the EL approach is faster than EE simulations until a certain number of bubbles is exceeded [9]. However, one important question is not answered: What is this certain number of bubbles? Is it of order hundred, ten thousand or even a million bubbles?

This work tries to give an answer to this question by comparing state of the art EE simulations with an EL approach. Our EL approach makes use of stochastic modeling and improved collision algorithms to compete against the other simulation techniques. All results are compared with solutions coming from different Method of Moment (MoM) based EE simulations. We investigate the solution quality and computational cost for both (EL & EE) and identify critical computation steps.

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.

Note, that this detail is achieved by resolving each single bubble. Thus, information of each individual bubble can be tracked. Bubble coordinates are not bound to the Eulerian mesh; trajectories are calculated with double precision floating-point manner. However, the underlying eulerian mesh does not contribute any more details than it would do in EE simulations. The opposite is rather the case, since a lower mesh resolution has to be chosen for the EL approach.

The original Open FOAM 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 bubbles – gaseous particles – require a rework of this original solver to account for the additional physical effects.

2.1. 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, 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:

(1)

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 k – epsilon 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, F_D, 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:

(2)

(3)

The open literature reveals different constants, C_g, and turbulent time scales, τ, according to specific physical effects in the bubbly flow [10-12]. In the presented EL simulation, the model of [11] was used.

(4)

2.2. 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 X_.b. Local mass can change due to break-up and coalescence events. The new position of any specific bubble is calculated discretely with the velocity .

(5)

Forces acting on the bubble will change its velocity. Here m_b resembles the bubble’s mass and ΣF stands for the sum of all acting forces on a single bubble.

(6)

The sum of forces ΣF consists of the buoyancy and weight force F_B, the drag force F_D, the lift force F_L, the virtual mass force F_VM, the wall lubrication force F_W and the dispersion force F_TD. Here the subscripts b and C stand for the bubble and the continuous phase accordingly, the subscript identifies the relative differences between them. Furthermore, g denotes the gravitational acceleration, p stands for densities, u for velocity, d_b for the bubble’s diameter, k for the turbulent kinetic energy and α for the phase fraction. Adjustable model parameters are denoted with C_i- .

(7)

These forces are in particular:

(8)

(9)

(10)

(11)

(12)

(13)

In (11), V_b stands for the bubble’s volume and Di/Dt denotes the material derivative, meaning that the derivative is made while following the bubble. In (12) the wall force parameter C_W 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 [13] several models for the drag and lift force coefficients C_D, C_L have been tested. Based on this, the models of [14] were chosen. The virtual mass coefficient is set to C_VM = 0.5 according to [15], the coefficient for the dispersion force is set to C_TD = 0.1 [16].

Besides movement and position, each bubble is carrying information of further properties. This is in particular information that cannot be obtained in EE simulations, e.g. shape, orientation, age and thermophysical properties. Of course, these properties react to surrounding conditions, e.g. the bubble density depends on the surrounding pressure and bubble temperature. Species concentration and chemical reactions can also be simulated. However, in this simulations chemical reactions and mass transfer are neglected (no solute), which is a topic for further studies.

2.2.1. Stochastic Modeling

Many aspects of bubble behavior are simulated using stochastic Monte Carlo (MC) models [17-21]. 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.

2.2.1.1. Bubble Collision

The first bubble interaction, which is described through a stochastic model, is the collision between bubbles. Here the collision probability by [22] is used. The actual position of bubbles is not considered, but rather the relative velocity and sizes of the bubbles are taken into account:

(14)

Note, that the volume of the computational cell C_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 will be drawn. The constraint for collision will therefore be:

(15)

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 [23] was customized to work also with the stochastic collision algorithm gaining a major speed-up 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.

2.2.1.2. Bubble Coalescence

After a successful collision, the next step is to calculate the coalescence probability. Previous validation simulations [13] have shown that the model of Coulalouglou and Tavlarides [24] is an appropriate choice. Basis of the calculus is the combination of bubble contact time t_contact and film drainage time t_contact. 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:

(16)

(17)

(18)

Here σ represents the surface tension, h_f and h_i are critical and initial film thicknesses, respectively, and d_eq 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.

2.2.1.3. 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 [24] into a break-up frequency, ω_b, and probability, P_break. Again, a random number will be drawn and compared in each time step for each bubble.

(19)

(20)

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 f of the two bubbles. This is a controversial point in break-up modeling, since many researchers claim that f should follow various different distributions. The model assumption of [24] is a normal distribution of f, where the ratio f = 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 < f < 0.99.

2.2.1.4. 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.

(21)

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.

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, 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 [25].

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 [26] 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) [27] and Direct Quadrature Method of Moments (DQMoM) [28], 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) [29] and Interfacial Area Transport Equation (IATE) [30] are based on only one equation to be solved, which improves the CPU speed even further.

For the comparison, the simulations of [25] and [31] will be used.

3.1. 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 d_3x - calculated via the tracked moments - the phase hold-up a_c & a_b and the mean relative bubble velocity u_rel the resulting forces can be calculated for each computational cell. For the drag force, the following equation is used.

(22)

Same procedure is done for the buoyancy force. Other forces were neglected in the given simulations by [25] and [31]. For the drag coefficient C_D, [25] used the model of [14], [31] used the model of [32]. 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).

3.2. 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:

(23)

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 [25] is adopted from [24], the break-up model originates from [33] 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 [31] are both adopted from [33]. 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 , the mean diameter 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.

There exist experimental data [34-36] and simulations ([25], [31]) for a set-up described in [36] consisting of a rectangular bubble column depicted in Fig. (1): An optimized mesh is used (Fig. 2).

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 [25] has been used for the EL approach in a first simulation. To improve simulation time and solution quality further, 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.

A mesh resolution study showed that the B&C kernels need to be fitted to the specific mesh resolution. Lowering the cell size led to a reduction of the mean bubble diameter. Thus, B&C model parameters and recommended mesh resolution have been adopted from literature [18]. A detailed grid convergence analysis and B&C parameter fitting was therefore disregarded for our EL simulation. 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 (d_b,max ≈ 1 cm) and the grid cell length L_c. 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, mostly following the theory of [37] that tells d_b/L_c < 0.67, while [36] even detected a numerical deviation when using a finer mesh. We followed a similar regulation with d_b,max/L_c ≈ 1 and with the resulting Sauter diameter between 5.5 mm and 7 mm, the theory of [37] is also fulfilled. 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.

Fig. (1). Simulation case set-up

Fig. (2). Different meshes in comparison: left – original EE mesh; middle – simplified EE mesh; right – optimized coarse EL mesh.

In experiment and simulation, different gas flow rates have been tested. Corresponding velocities and bubble injection sizes have been calculated with equations from [38] 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.

Note, that the numerical tolerances for the calculation of the eulerian fields (e.g. velocity, pressure etc.) were chosen identical in the EL and EE simulations. All simulations have been carried out using the OpenFOAM framework with the same numerical solver properties, wherever applicable. Unfortunately, there are field values, which are only calculated in some of the methods (e.g. moments) and can therefore not be compared in terms of numerical tolerance. In addition, the lagrangian velocities and positions are calculated using Newtonian dynamics, where there is no need for iterative matrix solving. Thus, a definition of a numerical tolerance is not possible; all lagrangian variables are calculated in double precision floating-point manner. Furthermore, particles are solved with a lagrangian time step that is smaller than the eulerian time step. Thus, even with the same tolerance levels no equivalent solution convergence could be defined.

Experimental and simulation results (s. Table 2) for Sauter diameter and hold-up are taken from [25], [31] and [36]. 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. The results from the EL simulation presented here are mere characteristic values calculated through averaging of the whole bubble size distribution to enable comparability to EE simulations. In general, the explicit simulation of each bubble will always yield higher detailed results than it is possible in any EE simulation.

All simulations have been computed on the same machine (Core i7 4790, 16 GiB RAM, Ubuntu 14.04) 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 (3) 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. It 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. Resulting hold-up for this setup was 6.3%, 6.7% and 14%, respectively. Note, that high hold-up values are not realistic anymore, since the EL method presented is not designated for a hold-up larger than 10%.

Fig. (3). Simulation CPU time for different numbers of bubbles in the domain.

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 [23], the optimal collision algorithm would need n log (n) steps to calculate for collisions of n particles. Since this can never be reached, the quality of an algorithm has to be quantified in terms of exponential notation n^a, where a smaller value for a qualifies a better algorithm. Fig. (4) compares the computational time of the original (old) and improved algorithm (new) that has been used.

Lowering the complexity from n²³ to n^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 algorithm 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.

Fig. (4). CPU time for colliding bubbles.

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, 0000 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 (d₃₂ < 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.