Flow Simulation in a 2D Bubble Column with the Euler-lagrange and Euler-euler Method

Andreas Weber, Hans-Jörg Bart*
Chair of Separation Science and Technology, TU Kaiserslautern, PO box 3049, Germany

Article Metrics

CrossRef Citations:
Total Statistics:

Full-Text HTML Views: 1194
Abstract HTML Views: 1068
PDF Downloads: 1050
ePub Downloads: 853
Total Views/Downloads: 4165
Unique Statistics:

Full-Text HTML Views: 711
Abstract HTML Views: 378
PDF Downloads: 391
ePub Downloads: 216
Total Views/Downloads: 1696

© 2018 Andreas Weber and Hans-Jörg Bart.

open-access license: This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International Public License (CC-BY 4.0), a copy of which is available at: https://creativecommons.org/licenses/by/4.0/legalcode. This license permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

* Address correspondence to this author at the Chair of Separation Science and Technology, TU Kaiserslautern, PO box 3049, Germany, Tel: +496312052414, Fax: +496312052119; E-mail: bart@mv.uni-kl.de



Bubbly flows, as present in bubble column reactors, can be simulated using a variety of simulation techniques. It is presented, how Computational Fluid Dynamics (CFD) methods are used to simulate a pseudo 2D bubble column using Euler-Lagrange (EL) and Euler-Euler (EE) techniques.


The presented EL method uses the open access software OpenFOAM to solve bubble dynamics with bubble interactions computed via Monte Carlo methods. The estimated bubble size distribution and the predicted hold-up are compared with experimental data and other simulative EE work with a reasonable consensus for both. Benchmarks with state of the art EE simulations shows that the EL approach shows good performance 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 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 superior to the EL approach - was estimated to be 40.000 in this particular case.

Keywords: Bubble size distribution, Euler-Lagrange, Euler-Euler, Bubble column, Method of moments.


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:


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, 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, Cg, 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.


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 .


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


The sum of forces ΣF consists of the buoyancy and weight force FB, the drag force FD, the lift force FL, the virtual mass force FVM, the wall lubrication force FW and the dispersion force FTD. 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, db for the bubble’s diameter, k for the turbulent kinetic energy and α for the phase fraction. Adjustable model parameters are denoted with Ci- .


These forces are in particular:


In (11), Vb 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 CW 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 CD, CL have been tested. Based on this, the models of [14] were chosen. The virtual mass coefficient is set to CVM = 0.5 according to [15], the coefficient for the dispersion force is set to CTD = 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. 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:


Note, that the volume of the computational cell Ccell 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:


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. 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 tcontact and film drainage time tcontact. 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 [24] into a break-up frequency, ωb, and probability, Pbreak. 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 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. 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.


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


Same procedure is done for the buoyancy force. Other forces were neglected in the given simulations by [25] and [31]. For the drag coefficient CD, [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:


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 (db,max ≈ 1 cm) and the grid cell length Lc. 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 db/Lc < 0.67, while [36] even detected a numerical deviation when using a finer mesh. We followed a similar regulation with db,max/Lc ≈ 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.

Table 1. Bubble sizes and number of bubbles to be injected for the different gas velocities.
case Name Min Med Max Max+
superficial gas velocity 2.4 mm/s 11.9 mm/s 21.3 mm/s 35.5 mm/s
bubble inlet diameter 6.3 mm 8.3 mm 10 mm 10 mm
bubbles per second 147 318 325 540


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.

Table 2. Resulting Sauter diameter, hold-up and plume oscillation time in comparison. 1d30
case Model / Ref. Experiment / Diaz DQMOM / Marchisio OPOSPM / Hlawitschka EL / own sim.
(coarse mesh)
min d32 6.83 mm 6.22 mm 5.1 mm1 6.17 mm
hold-up 0.62% 0.62% 1.28% 0.52%
med d32 6.5 mm 6.69 mm 5.9 mm1 6.93 mm
hold-up 2.63% 2.27% 2.86% 1.8%
max d32 7.73 mm 8.21 mm 6.8 mm1 5.52 mm
hold-up 4.1% 4.06% 3.92% 4.2%

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.

Table 3. Direct comparison of CPU time [h]
Simulation type EE EL
model QMOM DQMOM IATE OPOSPM Fine coarse
min 10.99 4.95 3.03 3.27 4.22 0.47
case med 33.10 33.82 14.72 15.52 7.82 2.78
max 43.91 33.57 19.79 14.56 11.26 3.33
max+ - - 20.3 14.09 - 5.18

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 na, 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 n23 to n1.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 (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.


d (Bubble) diameter
d23 Sauter diameter
f Volume ratio
F Force
g Gravity
h Layer thickness
k Turbulent energy
m mass
n Number; normal vector
P Probability
u Velocity
rel Relative
t Time
p Pressure
V Volume
x Position Greek letters
α phase fraction
x Random number
ε Turbulent energy dissipation
σ Interface tension
ω Frequency
p Density
µ Viscosity
t Turbulent time scale Subscripts, superscripts, symbols
b bubble
B Buoyancy; break-up
C Continuous phase; coalescence
D Drag
eq Equivalent
L Lift
VM Virtual mass
W Wall
TD Turbulent dispersion
turb Turbulent
k Turbulent energy
rel Relative
ε Turbulent energy dissipation
Material derivative
Normal distribution


Not applicable.


The authors declare no conflict of interest, financial or otherwise.


Funding by the Deutsche Forschungsgemeinschaft (DFG) within the RTG GrK 1932 ”Stochastic Models for Innovations in the Engineering Sciences”, project area P1, is gratefully acknowledged. We want to thank Mark W. Hlawitschka, TU Kaiserslautern, for supply of the EE solvers used in this work.


[1] Akita K, Yoshida F. Bubble Size, Interfacial Area, and Liquid-Phase Mass Transfer Coefficient in Bubble Columns. Ind Eng Chem Process Des Dev 1974; 13(1): 84-91.
[2] Liao Y, Lucas D. A literature review of theoretical models for drop and bubble breakup in turbulent dispersions. Chem Eng Sci 2009; 64(15): 3389-406.
[3] Liao Y, Lucas D. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem Eng Sci 2010; 65(10): 2851-64.
[4] Martínez-Bazán C. Splitting and dispersion of bubbles by turbulence 1998.
[5] Joshi JB, Vitankar VS, Kulkarni AA, Dhotre MT, Ekambara K. Coherent flow structures in bubble column reactors. Chem Eng Sci 2002; 57(16): 3157-83.
[6] Laín S, Bröder D, Sommerfeld M, Göz MF. Modelling hydrodynamics and turbulence in a bubble column using the Euler–lagrange procedure. Int J Multiph Flow 2002; 28(8): 1381-407.
[7] Sokolichin A, Eigenberger G, Lapin A, Lübert A. Dynamic numerical simulation of gas-liquid two-phase flows Euler/Euler versus Euler/Lagrange. Chem Eng Sci 1997; 52(4): 611-26.
[8] Zhang D, Deen NG, Kuipers J. Numerical simulation of the dynamic flow behavior in a bubble column: A study of closures for turbulence and interface forces. Chem Eng Sci 2006; 61(23): 7593-608.
[9] Sokolichin A, Eigenberger G. Gas-liquid flow in bubble columns and loop reactors: Part I. Detailed modelling and numerical simulation. Chem Eng Sci 1994; 49(24): 5735-46.
[10] Pfleger D, Becker S. Modelling and simulation of the dynamic flow behaviour in a bubble column. Chem Eng Sci 2001; 56(4): 1737-47.
[11] Rzehak R, Krepper E. Bubble-induced turbulence: Comparison of CFD models. Nucl Eng Des 2013; 258: 57-65.
[12] Yao W, Morel C. Volumetric interfacial area prediction in upward bubbly two-phase flow. Int J Heat Mass Transfer 2004; 47(2): 307-28.
[13] Weber A, Bart H-J. Euler-Lagrange simulation of bubble hydrodynamics using a stochastic approach Young Researchers Symposium 2016 (YRS 2016): Proceedings: April, 14th-15th 2016, Fraunhofer-Zentrum Kaiserslautern 2016; 154-9.
[14] Tomiyama A. Drag, Lift and Virtual Mass Force Acting on a Single Bubble 3rd International Symposium on Two-Phase Flow Modelling and Experimentation Pisa, Italy. 2004.2004.
[15] Delnoij E, Lammers F, Kuipers JA, van Swaaij W. Dynamic simulation of dispersed gas-liquid two-phase flow using a discrete bubble model. Chem Eng Sci 1997; 52(9): 1429-58.
[16] Lahey RT, Lopez de Bertodano M, Jones OC. Phase distribution in complex geometry conduits. Nucl Eng Des 1993; 141(1-2): 177-201.
[17] Garcia AL, van den Broeck C, Aertsens M, Serneels R. A Monte Carlo simulation of coagulation. Physica A 1987; 143(3): 535-46.
[18] Grünewald M, et al. Coalescence and Break-Up in Bubble Columns: Euler-Lagrange Simulations Using a Stochastic Approach. Chemieingenieurtechnik (Weinh) 2013; 85(7): 1118-30.
[19] Sommerfeld M. Validation of a stochastic Lagrangian modelling approach for inter-particle collisions in homogeneous isotropic turbulence. Int J Multiph Flow 2001; 27(10): 1829-58.
[20] Vikhansky A, Kraft M. Modelling of a RDC using a combined CFD-population balance approach. Chem Eng Sci 2004; 59(13): 2597-606.
[21] Vikhansky A, Kraft M. Single-particle method for stochastic simulation of coagulation processes. Chem Eng Sci 2005; 60(4): 963-7.
[22] O'Rourke P J. Collective drop effects on vaporizing liquid sprays L 1981. os Alamos National Lab, NM (USA), no No LA-9069-T 1981.
[23] Sigurgeirsson H, Stuart A, Wan W-L. Algorithms for Particle-field Simulations with Collisions. J Comput Phys 2001; 172(2): 766-807.
[24] Coulaloglou C, Tavlarides L. Description of interaction processes in agitated liquid-liquid dispersions. Chem Eng Sci 1977; 32(11): 1289-97.
[25] Hlawitschka MW, Schäfer J, Hummel M, Garth C, Bart H-J. Populationsbilanzmodellierung mit einem Mehrphasen-CFD-Code und vergleichende Visualisierung. Chemieingenieurtechnik (Weinh) 2016; 88(10): 1480-91.
[26] Hulburt HM, Katz S. Some problems in particle technology: A statistical mechanical formulation. Chem Eng Sci 1964; 19(8): 555-74.
[27] McGraw R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci Technol 1997; 27(2): 255-65.
[28] Marchisio DL, Fox RO. Solution of population balance equations using the direct quadrature method of moments. J Aerosol Sci 2005; 36(1): 43-73.
[29] Drumm C, Attarakih MM, Hlawitschka MW, Bart H-J. One-group reduced population balance model for CFD simulation of a pilot-plant extraction column. Ind Eng Chem Res 2010; 49(7): 3442-51.
[30] Ishii M, Kim S, Kelly J. Development of Interfacial Area Transport Equation. Nucl Eng Technol 2005; 2005(37): 525-36.
[31] Buffo A, Marchisio DL, Vanni M, Renze P. Simulation of polydisperse multiphase systems using population balances and example application to bubbly flows. Chem Eng Res Des 2013; 91(10): 1859-75.
[32] Haberman WL, Morton RK. An experimental study of bubbles moving in liquids. Trans Am Soc Civ Eng 1956; 121(1): 227-50.
[33] Laakkonen M, Alopaeus V, Aittamaa J. Validation of bubble breakage, coalescence and mass transfer models for gas?: Liquid dispersion in agitated vessel. Chem Eng Sci 2006; 61(1): 218-28.
[34] Cachaza EM, Díaz ME, Montes FJ, Galán MA. Simultaneous Computational Fluid Dynamics (CFD) Simulation of the Hydrodynamics and Mass Transfer in a Partially Aerated Bubble Column. Ind Eng Chem Res 2009; 48(18): 8685-96.
[35] Díaz ME, et al. Numerical simulation of the gas–liquid flow in a laboratory scale bubble column. Chem Eng J 2008; 139(2): 363-79.
[36] Díaz ME, Montes FJ, Galán MA. Experimental study of the transition between unsteady flow regimes in a partially aerated two-dimensional bubble column. Chemical Engineering and Processing: Process Intensification 2008; 47(9-10): 1867-76.
[37] Milelli M, Smith BL, Lakehal D. Large-Eddy Simulation of Turbulent Shear Flows Laden with Bubbles. In: Geurts BJ, Friedrich R, Métais O, Eds. Direct and Large-Eddy Simulation IV 2001; 461-70.
[38] Geary NW, Rice RG. Bubble size prediction for rigid and flexible spargers. AIChE J 1991; 37(2): 161-8.