University of Birmingham An Eulerian- Lagrangian numerical method to predict bubbly flows

. A large-eddy simulation based Eulerian-Lagrangian model is employed to study bubble plumes in an open channel with crossflow. The numerical results are validated with PIV experimental data. Good agreement between simulated and observed velocities is found. The impact of the crossflow on the structure of the plume and the resulting turbulent structures are described.


Introduction
Gas bubble plumes are widely used in environmental applications such as destratification of lakes and reservoirs (e.g. [1], [2]), ice prevention, accidental deep sea blowouts of gas, or gas leakage from natural vents (e.g. [3], [4], [5]). Crossflow, i.e. ambient current, is usually present in many of the aforementioned applications. Lima Neto [6] reported a crossflow of 0.3m/sec for an aeration project in the iced-covered Athabasca River, Canada. Although the effect of crossflows on single-phase plumes and jets has been studied extensively, little is known about their effect on multi-phase plumes.
The understanding of the underlying physics of the interaction between the gas and liquid phases has been a challenging task for both Computational Fluid Dynamics modellers and experimentalists. In order to capture the three-dimensional structures of multi-phase plumes, advanced 2D and 3D numerical simulations have been performed extensively on stagnant and/or vertically stratified ambient fluid. Three main approaches have been used in order to treat the continuous phase: Reynolds-Averaged Navier Stokes simulation (RANS) with different closures, mostly κ-ε model, ( [7]); Large-eddy simulation (LES) models that resolve directly large-scale turbulent motions and modelling only the unresolved subgrid-scale (SGS) effects ( [8]); and direct numerical simulation (DNS) approach models that directly resolve the entire turbulent spectra. Dispersed phase may be described in the same static reference frame as the continuous, leading to the Eulerian-Eulerian (EE) approach ( [3]) or in a dynamic (Lagrangian) frame, leading to the Eulerian-Lagrangian (EL) approach ( [9]).
Limited numerical studies are published on multi-phase plumes in crossflow. Yapa and Zheng [10], Johansen [11] and Zheng et al. [12] have developed integral models to simulate oil and gas blowouts in crossflow and stratified environments. Le Moullec et al. [13] applied a three-dimensional EE two-phase numerical approach for the simulation of a crossflow gasliquid wastewater treatment reactor using the CFD code FLUENT, testing two different turbulence models: κ-ε and Reynolds Stress Model (RSM). Even though time-averaged horizontal velocity along the length of the reactor and vertical velocities at specific heights were presented, the study is focused on the mixing and flow within the reactor and not on the gas-liquid physics. Wang et al. [14] presented instantaneous snapshots of the interface of deforming liquid bubbles and drops in a gas cross-flow, using a second-order hybrid level set-volume constraint method for numerically simulating a co-flowing liquid jet in gas deforming bubbles. Chen et al. [15] applied a LES based model in order to simulate doubleplume formation by oceanic CO2, including seawater entrainment, gravity waves and peeling out of the plume. Momentum, mass and heat transfer phenomena are described using empirical formulae. The simulated double-plume structure after the CO2 injection and temperature field contours are shown. Decrop et al. [16] extensively studied negativelybuoyant, vertically downward injected turbidity plumes, i.e. mixtures of water and sediment particles, in crossflow using two-phase Large-eddy simulations.
Eulerian-Lagrangian based Large-Eddy Simulation (EL-LES) employs Lagrangian Particle Tracking (LPT) to simulate the dispersed phase. Each bubble is represented by a Lagrangian marker which describes a trajectory across the Eulerian mesh according to its own motion equations. EL gives detailed information about every bubble's position, force and velocity. It is more expensive than EE because each particle requires the calculation of a set of equations and a mapping procedure between the Lagrangian and Eulerian coordinates.
In this paper the refinements to an in-house large-eddy simulation (LES) based CFD code to allow its application to study bubbly flows are reported. This is a very challenging case for CFD modellers due to the interaction between two strong sources of advection: the bulk velocity of the channel flow and the buoyancy of the plume. The effect of the lateral current on the plume is shown by the turbulent structures as a result of the interaction of the plume with the surrounding liquid.

Numerical Framework
In this study the in-house finite-difference-based Large-Eddy Simulation code Hydro3D is employed. Hydro3D solves the filtered Navier-Stokes equations on a staggered grid for the continuous (liquid) phase and has been validated thoroughly for many different flows (e.g. [17], [18], [19], [20]). The code is equipped with a Lagrangian Particle Tracking algorithm to allow for accurate predictions of the dispersed (bubbles) phase. This approach has been successfully applied to bubble plumes and validated in [9] and [21].
Hydro3D solves the space-filtered mass and momentum conservation equations for an incompressible fluid: (1) where ui refers to the filtered velocity component in the spatial direction i, t is the time, p the filtered pressure, ν the dynamic viscosity and Sij the strain rate tensor. The term τij accounts for the unresolved, subgrid-scale (SGS) turbulence, which is calculated through the turbulent viscosity νt using the Wall-Adapting Local Eddy-viscosity (WALE) sub-grid scale (SGS) model, with WALE constant Cw=0.45 for all simulations. ξi is a source term and accounts for the contribution of the dispersed phase to the flow. The derivatives in the governing equations are discretized with a two-step Runge-Kutta algorithm for the time derivative and secondorder central differencing schemes for both convective and diffusive terms. The code is based on a predictor-corrector fractional step method with the solution of the Poisson pressure equation using a multi-grid method as the corrector. The code is parallelized via domain decomposition and sub-domain communicate with the standard Message Passing Interface (MPI). The bubbles are represented by volumeless Lagrangian points/markers. The physical effect of the liquid-gas mixture is accounted for through forcing terms. The assumptions made are that the bubbles are rigid and spheric and that there is no direct interaction between them (due to the relatively dilute gas mixture). Also, only linear interaction between interfacial forces is considered. The motion of individual bubbles (from here onwards called particles) is computed by Newton's second law: where mp is the particle's mass, up,i is the particle's velocity in spatial direction i and Fp,i is the sum of the interfacial liquid forces acting on the particle in direction i. The integral forces acting on each particle are approximated by semi-empirical formulae. According to the most commonly accepted formulation (see [8]), five forces are considered: buoyancy, fluid stress, added mass, drag and lift. The expressions for each force are detailed in Fraga et al. [9] and should not be repeated for brevity. These forces are calculated based on the size of the bubble, its velocity on the previous time step and the velocity field of the surrounding liquid.
In a two-way coupling approach, as it is proposed here, exchange of information is required twice. First, the interfacial particle forces are calculated and the particles' velocities obtained. This is called forward coupling. Second, the contribution of the dispersed phase to the continuous phase is computed and added as a source term, ξi, to the liquid's momentum equations (Eq. 2). This is called reverse or backward coupling. Forward and backward coupling are achieved by connecting randomly placed Lagrangian particles with fixed locations of the Eulerian framework through the PSI-cell (Particle-Source-In Cell) mapping technique. By employing the PSI-cell technique, only the fluid nodes of the cells in which the Lagrangian particle's center is located receive the momentum source. For the forward coupling, the interpolation and/or transfer of quantities between phases is done through the smoothed delta functions developed by Yang et al. [22]. The momentum source term, ξi is computed as: (4) where M the number of bubbles in the influence volume of the fluid node under consideration, δ is the delta function at the node's location xj for a stencil length h and Fp,i* is the sum of the forces acting on that cell. The calculation of the backward forces is analogous to the forward ones, except that in this case the buoyancy force is not included.

Computational Setup
The bubble plume simulations reported here are an extension of a sets of experiments conducted at the Fluid Dynamics Laboratory in the Zachry Department of Civil Engineering at Texas A&M University (details can be found in [23]). The setup of the experiments is sketched in Fig. 1. In the experiments compressed air was injected at a constant gas flow rate at standard conditions through an aquarium airstone located 0.12m above the bottom of a glass-walled flume 35m long, 0.9m wide and 1.2m deep. The gas flow rate was set to Qg=0.5l/min and the bubble size to Dp=2.5mm with a standard deviation of 0.28mm. The approach flow velocity was 0.04m/s. Standard Particle Image Velocimeter (PIV) technique was utilized to measure the velocity field of the continuous phase in the plume center plane. The numerical simulations are carried out under analogous geometrical conditions to both experiments. Bubbles are initially randomly distributed over the diffuser and once they make their way through the tank and reach the water surface they are removed from the computational domain. Boundary conditions for the fluid phase include the use of the no-slip boundary condition at all solid walls and a rigid lid at the water surface with a free slip condition. Grid size sensitivity was studied previously ( [9]) and a constant grid spacing of Δx=7.5mm, Δy=4.5mm and Δz=4.1mm was chosen (22.4 Mio gridpoints). In order to provide a realistic, fully-turbulent inflow while avoiding extreme computational costs, the velocity inlet boundary conditions are generated using a Synthetic Eddy Method (SEM). The outflow is treated with a convective outflow boundary condition.      Figure 5 illustrates their origin and fate. Near the bottom of the plume a horseshoe vortex is formed but this vortex is rather weak compared to the vortex structure near the water surface. The strength and size of the vortex systems depend on a number of parameters including the gas flow rate and the strength of the crossflow. More details and some animations of this flow will be presented at the conference.

Conclusions
Large eddy simulations of a bubble plume were conducted by using a LES Eulerian-Lagrangian two-way coupling algorithm in an open channel with turbulent crossflow.
The velocities have been compared with PIV data at two different vertical profiles located at the plume's wake. The numerical results are in good agreement with experimental PIV measurements.
The turbulence structures inducted by the plume in the surrounding liquid are represented using Q-criterion isosurfaces. In the presence of a crossflow, the bubbles constitute an obstacle for the horizontal motion of the liquid, which is confirmed by the negative values upstream the bubble column in the horizontal velocity contour. A complex secondary flow in the wake of the plume is triggered by the presence of the plume. The resulting turbulent structure is characterised by the formation of two strong counter-rotating vortices at the upper half of the channel and weaker horse-shoe vortices close to the bottom.