Simulating spiraling bubble movement in the EL approach

Simulating the detailed movement of a rising bubble can be challenging, especially when it comes to bubble path instabilities. Free rising ellipsoidal bubbles not only move in straight lines but can describe sinusoidal, zigzag or spiraling paths. The common Euler Euler (EE) simulation techniques can no longer resolve the actual movement patterns and Direct Numerical Simulations (DNS) tend to be very costly when simulating a larger number of bubbles. A solution based on the Euler Lagrange (EL) approach is presented, where the bubbles show oscillating shape and/or instable paths while computational cost are at a far lower level than in DNS. The model calculates direction, shape and rotation of the bubbles to finally create typical instable path lines. This is embedded in an EL simulation, which can resolve bubble size distribution, mass transfer and chemical reactions. To ensure realistic solution, validation against experimental data of single rising bubbles and bubble swarms are presented. References with 2D and also 3D analysis are taken into account to compare simulative data in terms of typical geometrical parameters and average field values.


Introduction
Understanding bubbles path instabilities is a major challenge since the 1960s. Early measurements were performed by Aybers & Tapucu (1969), where a single bubble was rising in a fluid showing different path lines, like zigzag and spiraling. Nowadays, a 3D camera setup enables a more detailed view (Shew & Pinton, 2006b) and sophisticated simulations (Mougin & Magnaudet, 2006) could help to explain the complex rising behavior. The origin of this behavior is in the turbulent eddies induced behind the bubble during its rise. Instabilities in these eddies lead to an eccentric force on the bubble, which then leads to a deviation from a straight rise resulting into a spiraling or zigzagging path. The development of the turbulent eddy behind a bubble is connected also to its shape and size. Different types of bubble paths can be experienced, based on the bubble size or, in terms of turbulence intensity, on the bubble Reynolds number.
Simulating bubble hydrodynamics is an issue since many years. The most common approach to do this is the simplified Euler Euler (EE) method, where single bubbles are no longer resolved but a bubble number density is taken into account instead. This leads to a much simpler way to calculate for an instable bubble path, where a diffusion term is used, which gives a high loss of detail while the overall spatial distribution of bubbles can be forecast with adequate precision. The other extreme is a high detail direct numerical simulation (DNS) of a single bubble or rather a small bubble swarm. Bubbles are resolved in full detail including the hydrodynamics inside. In doing so, path instabilities can be simulated consequently on the lowest level of scale but at a high computational load.
Free rising ellipsoidal bubbles not only move in straight lines but can describe sinusoidal, zigzag or spiraling paths. The common Euler Euler (EE) simulation techniques can no longer resolve the actual movement patterns and Direct Numerical Simulations (DNS) tend to be very costly when simulating a larger number of bubbles. This work presents a solution to calculate the orientation and shape of bubbles using the Euler Lagrange (EL) approach. Advantages lie within the fast computation and the high level of detail. In comparison to DNS the insides of the bubbles are not calculated in full detail but macroscopic models are employed. Every bubble is calculated individually, having its own size, direction and shape. The surrounding fluid will influence not only the bubble's movement but also the rotation and shape. The actual calculation of the turbulent eddies behind the bubble will not be carried out, but an oscillation orientation model is used and model parameters are calibrated from experimental data.

EL modeling
The general EL model describes bubbles as a point volume acting under Newtonian dynamics. Forces created by the surrounding fluid and neighboring bubbles accelerate these Lagrangian points through the domain. The continuous phase itself is calculated using Navier-Stokes equations and can be coupled to the interaction forces by a source term. All forces are calculated for each bubble individually, which produces an individual path for each bubble. Therein lies one of the advantage of the EL approach in comparison to the EE methods. Bubbles can coalesce and break, which gives a bubble size distribution, with even more detail than common method of moments approaches. Downside of the EL approach is the higher computational load, which is strongly dependent on the number of bubbles simulated. Nevertheless, the EL approach has been used frequently in several different simulations of bubbly flow (Gruber et al., 2015).

Liquid phase hydrodynamics
The continuous phase is assumed to be incompressible, basis for calculation is a modified Navier-Stokes equation.
Given is the continuous velocity uc, the pressure p, the density the viscosity and the source term f.
The source term f depicts the forces of the bubbles and will thereby serve as a coupling of the phases. Turbulence is computed using the standard RANS k-epsilon model (Launder & Spalding, 1972) with additional bubble induce turbulence (BIT) (Rzehak & Krepper, 2013). Each bubbles' drag force induces turbulent energy and dissipation in the associated computational grid cell. For each cell this sums up and is added to the source term in the turbulence model.
Here, Sk denotes the source term for the turbulent energy k, FD stands for the drag force. Sε is the turbulent dissipation with τ, the turbulent time scale.

Bubble hydrodynamics
Bubbles are modeled as point volumes acting under Newtonian dynamics. Their movement is calculated using a number of different forces.
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 bubble dispersion force FBD. Here the subscripts b and c stand for the bubble and the continuous phase accordingly, the subscript rel identifies the relative differences between them. Furthermore, g denotes the gravitational acceleration, stands for densities, u for velocity, V for the bubble's volume, k for the turbulent kinetic energy and α for the phase fraction. Appropriate model parameters are denoted with Ci.
Di/Dt in eq. (8) denotes the material derivative, meaning that the derivative is made while following the bubble. Drag and lift force coefficients CD,CL are calculated using the models of Tomiyama (2004). The virtual mass coefficient is set to CVM = 0.5 according to Delnoij et al. (1997), the coefficient for the dispersion force is set to CTD = 0.1 (Lahey et al., 1993). This dispersion describes the ambition of bubbles to spread due to collisions with other bubbles. Additionally, a second turbulent dispersion is used to model the collision of bubbles and turbulent eddies. The Random Dispersion Model (Smith & Milelli, 1998) is used to calculate eddies according to the surrounding level of turbulence. Assuming an isotropic turbulence, eddies are traveling through the liquid in a uniformly random direction with a specific life time. In the model, the turbulent eddy lifetime tE is evaluated, after which a new eddy is calculated.
Then the movement direction of the eddy is uniformly chosen while its velocity follows a normal distribution with a variance dependent on the turbulent energy k.
This turbulent velocity is added to the underlying continuous phase velocity for each bubble individually, which is then used to calculate the different bubble forces. Especially the drag force calculation is a crucial step to produce dispersion of the bubbles. In order to achieve the same amount of bubble dispersion like in the experiment, the turbulent velocity for drag force calculation had to be slightly increased.
Change in bubble volume due to the pressure drop while rising is calculated assuming an ideal gas. Mass transfer and/or chemical reactions are not considered in this simulation.

Ellipsoidal bubble model
According to the state of the art of bubble simulations almost any model assumptions are based on bubbles having a spherical shape. Simplest example is the Sauter diameter d32, which maps the mean volume/surface area ratio of a bubble population to one spherical bubble size. Collision frequency or rather probability, e.g. by O'Rourke (1981), is based on the calculation of the overlapping volume of spheres. Other examples are mass transfer calculation (spherical transport area), coalescence and breakup (interfacial energy of spheres) or simply the distance calculation between bubbles or wall and bubble. Only exception is the drag force, which is often modeled with respect to the shape of the bubble e.g. by Tomiyama (2004) but uses no information about orientation or rotation of the bubble. One has to admit, that the assumption of bubbles being spheres is adequate for many problems. For example, the deviation of the collision probability is negligible if a spatially uniform random orientation of the bubbles is assumed. Nevertheless, for the computation of realistic bubble movement, a better model for the bubble shape has to be chosen.
In this work, the deformed bubble will be approximated with an oblate spheroid, an ellipsoid with two different axes a and c, where c > a . For simplicities sake, all following figures will contain a simplified 2D version of the spheroid in Figure 1. The ratio of the axes is chosen to describe the bubble shape via the shape factor sf.
A shape factor of sf = 1 describes a perfect spherical shape, while a lower value stands for a more deformed sphere. sf = 0 would describe an infinitely thinned spheroid. Some models are based on empirical measurements using dimensionless numbers for the calculation, others are physical models derived from the interfacial tension and pressure distribution. It turned out that a good result for the simulation problem was achieved by Moore (1965) with the underlying equation: Since the shape factor has to be computed as a function of the Weber number, an approximation has been used: Note that this equation is based on the dimensionless Weber number which can represent changes in the shape induced by fluctuating relative velocities. The critical Weber Number Wecrit = 3.745 describes the transition to irregular bubble shapes. If the current Weber number is higher than Wecrit a shape factor of sf = 0.2 is chosen. Differently than with common rotation calculation, using the moment of inertia and torque action (Shew & Pinton, 2006a), an approach originally developed by Taylor (1923) is used. The direction of the bubble is described by a vector p pointing in the direction of axis . The change in direction is given as the vector , which is deduced from the rotation vector ω (s. Figure 2). This notation is beneficial because there is no need for a change to spherical coordinates (Θ, Φ), which would lead to higher computational load. The only requirement is a sufficiently small change of the vector p, which is fulfilled within the computed time step of the simulation.

̇= ⃗ ⃗ ×
(20) Figure 2: left: orientation in spherical coordinates; right: orientation with p and change in direction ̇ This rotation notation has been used by Junk & Illner (2007) to calculate the orientation of rigid ellipsoidal bodies in a Stokes flow. Based on these equations, a modified model for ellipsoidal bubbles was developed. In general, the bubble rotation is calculated using an explicit Euler algorithm: Just as the bubble's position and velocity, the orientation p changes due to the orientation change . The bubble's rotation relaxes against the outer rotation G with the factor ̅ . Here, ̅ denotes a simplified interaction of the bubble's moment of inertia and torque. A low ̅ implies a high bubble mass/inertia and thus a slower rotation due to outer forces. The additional term ̅̅̅ resembles a random rotation due to turbulence, where is used to scale the effect of turbulent randomness. It is generated similarly to the turbulent dispersion velocity ̅̅̅̅. Since the effective moment of inertia of a gas bubble inside a fluid is unknown, ̅ has to be derived from experimental data.
The outer rotation G is a summation of three main mechanics acting on the body (eq. 23). The first line implies that the rotation of the surrounding fluid is transmitted to the bubble itself. The second line describes the rotation induced by shear stress in the surrounding fluid, [⃗⃗⃗⃗ ] is the symmetric part of the Jacobian matrix.
[⃗⃗⃗⃗ ] = 1 2 The third line is an addition to the original model (Junk & Illner, 2007) and describes the ambition of the gas bubble to orient itself along the gravitation direction g. Since the original equation was deduced for rigid bodies only, all of the above-mentioned mechanics are weighted by the parameters Ji to fit experimental data of (non-rigid) bubbles. This also compensates for the fact, that there is no Stokes flow around the bubble and its interface is mobile (slip condition).
Especially the third line of eq. 23 leads to an oscillatory system, which enables the bubble to describe sinusoidal or helical orientation characteristics. While the orientation of the bubble is changing, the vector of orientation change ̇ will eventually point to a direction not perpendicular to p, thus rising in size and finally damping and stopping the oscillation. To prevent this, the direction change vector will be moved to the plane normal to the orientation vector by subtracting the part parallel to it: In total, this will lead to a slight damping of the oscillation which can be eliminated by preserving the magnitude of the orientation change vector: To preserve robust behavior a slight damping (δ = 0.2) is executed for the orientation change: Finally, with an interaction force based on the orientation, the bubble experiences a drift perpendicular to the main movement direction. This results in a bubble trajectory describing helical and sinusoidal (zigzag) paths. This perpendicular force has the same direction as the change in direction (Mougin & Magnaudet, 2006), which gives: This force is modeled as a modified drag force. Also, to imitate real bubble behavior, we limit the force to only occur on bubbles within a certain range of Reynolds number (Aybers & Tapucu, 1969). This will stop the oscillating motion for very small (d < 0.8 mm) and large bubbles (d > 3.5 mm) as observed in experiments.
The additional side force due to bubble rotation is modeled with the magnitude of the current drag force and the bubble rotation but will point in the direction of ̇ instead. (It turned out, that a scaling linear to the magnitude of the direction change || will lead to instabilities easily, which can be stabilized by using the square root √| | instead.) The bubble path amplitude is calibrated with the parameter β, where a higher value implies a larger amplitude of the resulting oscillating path.

Parameter estimation
Parameters ̅ , R and Ji (eqs. (22), (23)) are derived from a parameter study and with analysis of the oscillatory equations. The characteristic bubble path amplitudes, frequency and wavelength shown in Shew & Pinton (2006b) were taken to set first limits to the parameter study. In a more detailed analysis, the parameters are correlated to experimental data showing bubble swarms rather than single rising bubbles. The resulting simulated bubble paths are then compared to experimental data from the Bubble Column Reactor Database of the University of Magdeburg 1 . Characteristic oscillatory dimensions can be derived from eq. (22) and the third line of eq. (23). This will lead to a simple harmonic oscillator equation with the following form: Without damping of the oscillation, the frequency f0 can be calculated with the parameters J3 and ̅ .
It turned out, that this is the case for our oscillatory system, but only if the random rotation is set to zero. A case without random rotation was set up for this reason, the results are shown in Figure 4. Over a wide range, the simulated frequency is identical to the analytical solution for the undamped oscillation. Most references are made on a basis of a stagnant liquid phase, such that the bubble rise velocity ub can be easily estimated. This makes it easy to calculate an appropriate wavelength of the generated/observed bubble path using its characteristic frequency. = In case of a non-stagnant liquid, like in a bubble column, the bubble velocity has to be identified at first. It is therefore important to mention, that the wavelength of a bubble path in a dynamic system will differ from most reference experiments, which are made using a well-defined surrounding. Also, in some areas of the column a downward flow occurs, lowering the bubble's rising velocity. To overcome this problem, the bubble velocity and/or wavelengths are averaged over a large number of bubble paths.

Experimental setup
For the parameter study, simulation results were compared to experimental measurements from the Institut für Strömungstechnik und Thermodynamik, University of Magdeburg (Hlawitschka et al., 2016;Zähringer et al., 2014). The experimental setup consists of a cylindrical air-water bubble column with an inner diameter of 14.22 cm and a filling height of 73 cm. The air inlets are positioned at the bottom of the column and consist of four nozzles in a row. The distance between nozzles is 2.2 cm. The measurements were made with a camera system facing the four nozzles side by side (front view) and were also done a second time when facing the nozzles behind each other (side view). With this setup, the bubble positions and velocities were captured so that the trajectories could be reconstructed. With an air throughput of 10 l/h a mean bubble size of d = 3 mm was estimated using the same camera setup. To ensure a correct analysis, only those trajectories were considered, that were of appropriate length (minimum of 30 measured positions). Unfortunately, a 3D reconstruction was not possible for the bubble paths, because the side view measurements were not made simultaneously with the front view.

Computational setup
The corresponding computational mesh was created using a rectangular grid of 28 x 28 x 146 cells with overlapping cells being removed and reshaped to fit the cylindrical shape. This results in a mesh with 90000 cells, shown in Figure 6. Mesh resolution was intentionally left at approximately 0.125 cm³ per cell (mean cell length of 0.5 cm), because the bubble size must be smaller than the cell size in our EL approach. Cells near the walls are not rectangular anymore, but are also solemnly populated by bubbles due to wall lubrication force. The main bubble flow therefore happens in the center part, where the cells are in perfect rectangular shape.  Table 1), the simplest possible case could be adopted. All walls except the top outlet patch are combined and have the same properties. Since bubbles are injected using the Lagrangian approach, there is no need for an inlet patch. Instead of a free surface at the top, a slip condition has been set. This approach enables to use a single phase solver but is surely not suitable for estimating hydrodynamics near the surface. Since our interest lies within the lower part of the column, this error is acceptable. On the basis of the experimental measurements, the mean inlet diameter of the bubbles was set to d = 3 mm. With a throughput of 10 l/h this leads to a total number of 196.4 bubbles per second to inject from the four nozzles. Bubble breakup, coalescence and mass transfer was turned off in this case, since only the bubble movement is of concern and measurements showed only a minimum of bubble size variations during rise. The experimental measurement area is positioned in the lower 312 mm of the column only, but the whole height of the column was simulated to achieve the correct flow pattern.
Since there are 4 parameters in total to be calibrated, an automated parameter study software has been used (Dakota 2 ) to simulate a variety of combinations. The target function was the characteristic bubble flow path parameters, which were compared by screenshots at a first sight. After choosing the most promising parameter interval, a detailed analysis of bubble path frequency and amplitude followed. Therefore, an automated analysis of the generated path lines was used to calculate for mean wavelength and amplitude (~600 trajectories per simulation). The same method has been used to characterize the experimental measurements. Since the experimental data was collected using a 2D visual acquisition, simulation results were also calculated using a 2D mapping. Additionally, a 3D analysis was made and compared to other reference measurements.

Velocity profiles
Simulated liquid velocity and its deviation are compared to the experimental measurements in Figure  Calculation of the simulated fluctuating velocity is made using the k-epsilon turbulence model. Values descripted here (s. Figure 8) are solution of the turbulent energy k with assumed isotropic turbulence velocities: The upper figure shows the fluctuating velocity deviation of experiment (left) and simulation (right) on the middle plane of the column. Again, the scale on the left side shows the line plot positions A-D. In the lower sections, the fluctuation velocity is similar to the measured values from the experiments, in higher sections the simulation underestimates the level of turbulence, leading to low fluctuating velocities in comparison to experimental data. Here again, a standard parameter set was used for the BIT model, which was not optimized on this particular case as investigations on bubble induced turbulence was also not in the focus. However, as a consequence the impact of turbulence on bubble rotation is adjusted via the parameter study.

Bubble path characteristics
According to literature (Shew & Pinton, 2006a), the bubble path oscillations of a single air bubble (d = 3 mm) rising in tap water is approximately f = 5 s -1 and with its mean free rising velocity of ub = 35 cm/s the bubble path wavelength is λ = 7 cm. According to eq. (32) this would recommend setting the parameter ̅ 3 ≈ 1000.
Anyhow, analysis of the experimental data from the University of Magdeburg yielded a mean bubble path wavelength of λ = 4.14 cm and an amplitude of A = 1.87 mm. Mean bubble rise velocity is ub = 29 cm/s in this case which gives a frequency of f = 7 s -1 and a product of parameters of ̅ 3 ≈ 2000. This was only a first estimation, since the random rotation is not be set to zero in the final simulation, which changes the frequency slightly.
Experimental data from 4 different measurements on the same properties were used to reconstruct and analyze bubble trajectories (s. Table 2). The trailing numbers in the Case name stand for the frame rate and the number of frames used in the analysis. A given frame rate of 0.1 kHz and 100,250 or 500 frames implies a measuring time of 1, 2.5 or 5 seconds accordingly. Deviation of wavelength and especially of amplitude is quite high, which is partly owed to the 2D analysis of the bubble paths. Since a flat zigzag path can only be seen correctly from one perspective, the calculated 2D amplitude will hardly correspond to the true 3D amplitude, it will generally underestimate the true value. This is why the standard deviation of the analyzed amplitude show rather high values in all experimental and simulation cases. With the appropriate set of parameters, the simulated bubble path characteristics properly match the experimental values. Final simulations show that the assumed parameter ̅ 3 ≈ 1000 to 2000 is a quite good starting value, the most promising results were done with parameters in the range of ̅ 3 = 1250 to 1350. It turned out that the values of the parameters for direct influence of liquid rotation (γJ1) and shear (γJ2) are very low in contrast to the oscillation (γJ3), while the random factor (R) is of the same magnitude. Thus, the impact of rotational/shear flow around bubbles plays a minor role on the rotation while it is dictated mostly by turbulent eddies hitting the bubbles. A visual comparison of the three most promising bubble paths in Figure 9 reveal only few evident differences. The overall shape and distribution of path lines is almost identical. Most noticeable differences can be observed directly at the bubble inlets at the very bottom of the column. In the simulation bubbles are spread earlier than it is observed in the experiments. In the experiment bubbles describe a straight line directly after being injected to the column, while the simulation shows a rotation of the bubbles at the very beginning of the path, pushing the bubbles into a turn after the injection. In the first 100 mm a pattern can be seen in the flow paths of the simulation. The flow in the middle of the column pushes the bubbles away from the middle line, creating an area where almost no bubble cross. Some bubbles tend to take the same path multiple times, while in the experiment a slightly more chaotic distribution is present. In the simulation, a parameter for random rotation is used, but the analysis does not concern the overall distribution of bubbles. The spread/diffusion of bubbles in the upper area however (height > 150 mm) yields no (qualitatively) visible difference between experiment and simulation. Another point of concern is the behavior of the simulated orientation in comparison with the velocity of the bubble. In order to describe this relations, different deduced angles are used; the movement angle lies between the vertical axis and movement direction, the orientation angle is between the short axis direction and vertical axis and the drift angle lies between orientation and movement vectors. With 2D analysis of the simulation data, the orientation angles show a normal distribution with mean value near to zero (Figure 10, left), just like the experimental results (Sommerfeld & Bröder, 2009). A 3D analysis however proves the oversimplification when using 2D analysis because the simulated mean 3D orientation is 18.3°. This gets clear when figuring a perfectly spiraling path seen from only one side, where there is a constant orientation angle in reality but an alternating angle is seen from a fixed observer. In a worst case scenario, a flat sinusoidal path oscillates parallel to the line of sight, thus making it impossible to see any oscillations at all. Figure 10: Characteristic angles between velocity and orientation vectors of simulated bubbles, left: orientation angle using 2D data, middle: orientation angle using 3D data, right: drift angle vs. orientation angle (3D data) Bubble velocity and orientation vectors roughly point to the same direction, this has been shown in several references (Sommerfeld & Bröder, 2009;Ellingsen & Risso, 2001;Mougin & Magnaudet, 2006), which implies a small drift angle. Simulation results also show mostly small drift angles ( Figure  10, right) with a mean value of 19.2°. Experimental data shows, that the largest drift angles occur at the most outer points of the oscillation, especially when a planar zigzag path is fulfilled (Mougin & Magnaudet, 2006). This means a high drift occurs when the orientation angle is near zero. Simulation data shows other characteristics; a high drift angle correlates with a high orientation angle (s. Figure 10, right). In DNS simulation, the drift and orientation angles have shown sinusoidal characteristics, while both were 90° out of phase. In the EL simulation these angles are in phase. Unfortunately, there is no experimental proof for this correlation with 3D measurements in a bubble swarm.
When taking a closer look at single trajectories of one of the bubbles rising in the swarm it will reveal characteristic parts of the bubble's movement. As shown in Ellingsen & Risso (2001), the bubble paths can be described as 'flattened helixes' which become less flattened while rising. Exactly this behavior can be found in the simulated bubble paths depicted in Figure 11. The bubble path shown here is not picked randomly, since characteristic movement is not achieved with every bubble. Note that this simulation considers bubbles in a large swarm rather than single rising bubbles (like in the experimental reference). This leads to a more chaotic liquid and bubble movement which can abruptly change bubble movement direction and orientation. Near the walls, a downward flow occurs in which bubbles performed differently, mostly not showing steady oscillations. As most bubbles rising in the middle part of the column, oscillation slowly starts within the lower 100 mm where it evolves to a sinusoidal and finally to a helical path (at approx. 400 mm height). We can observe influence of the flow in terms of a strong drift in the bubble paths (s. Figure 11, top row). In order to give a three dimensional description of the paths, a polynomial fit (n = 3) is used to smooth out the bubbles drift movement. In Figure 11 (top row) the original bubble paths are shown from three sides with the polynomial fit (dotted line). The bubble paths are normalized using the polynomial fit and the resulting path data (s. Figure 11, bottom row) describes a spiral, which is more comparable to experimental work (where the drift movement is almost zero). Assuming that our path describes a perfect flattened spiral, its coordinates perpendicular to the rising direction can be taken to calculate for a 'dynamic radius' of the path. A flattened spiral would look like an ellipsoid in the X-Y view (s. Figure  11, bottom left). The maximum of this radius will give the amplitude of the major oscillation mode, while the minimum will give the minor modes amplitude. When the minor amplitude is zero, a completely flat spiral, a planar zigzag path, is present. Also the wavelength in rising direction can be extracted from this data by measuring the distance for one spiral rotation. In Simulation case 1 the mean (3D) wavelength is λ = 4.79 cm and mean amplitude is A = 4.04 mm. This is a much larger amplitude, than it has been calculated for the 2D analysis, but comparison with Ellingsen & Risso (2001) shows a good match. They analyzed a single rising bubble with equivalent diameter d = 2.48 mm and could measure a major mode amplitude of A = 4.3 mm. Measured path frequency was ω = 39 rad/s, which is ω = 35.4 rad/s in our simulation case (mean rise velocity ub = 0.264 m/s).

Conclusions
The presented EL simulation is capable of simulating unstable sinusoidal/spiraling bubble paths using macroscopic models. Bubble orientation, rotation and shape are calculated to achieve characteristic movement. Due to the assumption of bubbles describing rotational spheroids, the additional parameters that have to be calculated reduce to a shape factor, rotation and orientation vectors. A force induced by the bubbles rotation produces the lateral force leading to an oscillation movement.
A parameter study was used to fit the model constants to experimental data for mean bubble size of d = 3 mm. Characterization of the bubble paths was made using the amplitude and wavelength of the typical spiral movement. Evaluation of the 3D bubble path in a bubble swarm is difficult, since most references only supply 2D camera setups or single bubble trajectories in 3D analysis systems. It turned out, that most 2D measurements cannot reflect characteristic path parameters entirely, especially orientation angles are problematic. Nevertheless, the simulation results were mapped to a 2D point of view and compared to the experimental data. After parameter fitting, comparison to reference bubble path data was made using 2D and 3D analysis and could prove correct reproduction of unstable bubble paths. Comparison to DNS simulation and single rising bubble path data could also show good agreement. Amplitude and wavelength of the simulated bubble path are in unison with the measurements. Detailed comparison of DNS results of the drift angle reveal slight disagreement.
For further improvement of the model, a predictive parameter approach should be used to also cover different bubble sizes. Interaction of deformed bubbles are not considered in the EL model shown here, this could include collision, break-up, mass transfer and other shape dependent processes.