Simulation of the rising of gas bubbles in a pilot-scale external loop airlift photobioreactor

Air-lift reactors offer an interesting option as a microalgae cultivation system, especially for biorefineries. To optimize this application, a precise description of the moving interfaces formed by the liquid and gas phase is critical. In this paper, a coupled Level Set Method (LSM) and finite difference method is used to simulate gas bubbles dynamics in a pilot-scale external loop air-lift photobioreactor in which microalgae are used to capture CO2 from flue gas and to treat wastewater. Numerical simulations are carried out on a rectangular domain representing the section on the vertical axis of the riser. The data employed were either acquired from previous experimental campaigns carried out in the ALR or found in the literature. The rise, shape dynamics and coalescence process of the bubbles of flue gas are studied. The issue of volume loss characteristic of standard LSM is dealt with the conservative level set method. Computation results show good correspondence with the experimental ones.


Introduction
The attention to renewable energy sources, from topic of interest solicited by the necessary awareness towards sustainable development, is also starting to be seen as an economic opportunity, thus embracing the triple paradigm of the properly-defined sustainability (economic, environmental, social). Closely related to this discussion is the concept of circular economy, which aligns well with that of biorefinery. Biorefineries are fully integrated systems that exploit biomass to generate energy, biofuels, and/or products with high added value. By resorting to a holistic view, these systems aim at minimizing the waste flows outcoming from the closed-loop. In recent years, interest in microalgae has been spiked since they represent a multipurpose environmental tool, as they can be used in a variety of ways, including mainly wastewater treatment, CO2 biofixation and biofuels production. Due to such characteristics, microalgae growth constitutes a widely studied process to be incorporated in biorefineries. An example is offered by the integrated system presented in [1] that makes use of waste frying oils, solid organic and algal biomass to produce biodiesel and energy while treating wastewater and flue gases as well as possibly bio-fixating CO2. The most critical component of the whole process is probably the algal biomass production. In fact, it is carried out in particular photobioreactors, known as Air-Lift Reactors (ALRs). A sort of evolution of bubble columns, ALRs represent a class of pneumatic reactors capable of achieving good mixing conditions without high energy demand thanks to the prompting of a swirling motion of natural circulation due to density gradients. Due to these qualities, as well as the possibility of obtaining very fast photoperiods, ALRs display superior performances in terms of microalgae cultivation with respect to most other microalgae culture systems [2]. According to their structure, ALRs are further classified in internal and External Loop (EL-ALRs) reactors. The latter, employed in the discussed biorefinery, are constituted by a transparent column, the riser, equipped with an air injector at the bottom and connected by two horizontal collectors to another vertical column, the downcomer. The reactor, filled with water, contains as the dispersed phase the microalgae. For the case at hand, it was designed as to contain bubbles formed due to flue gas insufflation from the sparger, only in the riser.
In the following, a front-capturing method known as Conservative Level Set Method (CLSM) is employed to model bubble dynamics in the riser of one of the pilot-scale EL-ALRs used in the considered biorefinery. The motivations behind the simulation of such phenomenon in detail are multiple. The shape, trajectory and size of the ascending bubbles have clear repercussions on mass transfer since they are related to permanence time and interfacial surficial area. As LSMs are based on a Eulerian framework, they allow to take into account topology variations, i.e. the possibility of bubble coalescing and splitting, rather straightforwardly.
According to what has been just said, bubble shape and trajectory optimization positively influences microalgae growth as well as their recirculation. Excessive turbulence in the flow could spell damage to the cells. A precise understanding of the bubble flow regime is thus vital to improving the efficiency of the processes happening inside ALRs and the ability to opportunely manipulate the bubble flow would represent a major achievement for microalgae cultivation.

Mathematical modelling
To simulate the hydrodynamic environment inside the riser, a physically-based mathematical model was built upon the CLSM. CLSM [3] is a front-capturing method hinging on the one-fluid formulation of physical variables. It was conceived as a revisiting of the standard Level Set Method (LSM) [4] and it is aimed at overcoming the main drawback of the original method: the intrinsic volume loss.
The idea behind level set methods is straightforward: think of the interface, e.g. in 2D, as the zero-level set of a three-dimensional surface. This latter, in LSM, is given by the level set function ϕ, a continuous map defined as a space and time-dependent signed distance field that assumes positive values on one side of the front, negative on the other and is zero at the interface. As ϕ is advected, the front is propagated implicitly, with topology changes being handled effectively and with simplicity. As said, the main downside of this method is that, for numerical reasons, the front is smeared out instead of sharp and volume conservation is therefore not achieved [5]. CLSM, to overcome this issue, employs a different phase function ψ (Fig. 1), defined as ψ(x, t) = Hε(ϕ(x, t)), where Hε is a smoothed Heaviside function [6]. It is initialized as a hyperbolic tangent (Eq. 1) of thickness ε as, for ε→0, it tends to the exact Heaviside step function H(ϕ).
CLSM is articulated in two steps. The first one, similarly to standard LSM, is the advection of the phase function, which, in the presence of a solenoidal flow field, can be expressed as: Eq. 2 is a conservation equation and as such lends this property to the algorithm.The second step is the reinitialization, needed to regularize the shape of ψ and enhance numerical robustness. It consists in the resolution of the hyperbolic Partial Differential Equation (PDE) Eq. 3, comprised by a compressive limiter intended to sharpen the profile and a diffusive one in the normal direction (given by the unit normal vector n), used to balance it and maintain the adequate interface thickness. The evolution equation Eq. 3, reported in conservative form, has to be solved to steady state in an artificial time τ framework. Courant-Friedrichs-Lewy (CFL) condition ensues for both advection and re-initialization steps. In particular, when applying the artificial compression technique, Olsson and Kreiss [3] suggested that with Courant number C = 0.25, ε interface thickness and d either 0 or very small (e.g. 0.1) for complicated flows.
To work with Eq. 2, another ingredient is needed: the knowledge of the speed vector field u. This can be achieved handling the incompressible Navier-Stokes (NS) equations, which were employed here in their non-dimensional formulation for multiphase flow [5].
where t is time, p the pressure of the fluid, g ⋆ the gravity acceleration unit vector, κ the curvature, δ the Dirac delta distribution, Fr is Froude number, We is Weber number and Re is Reynolds number. Moreover, superscript ⋆ denotes the dimensionless character of the parameters.
Eq. 5 allows to describe a multiphase flow thanks to the so-called one-fluid formulation with which material properties like density and viscosity are assimilated to discontinuous fields and are thus valid for the whole computational domain. The notation is therefore kept compact, combining the various multiphase equations in a single one including the jump conditions at the interface. Practically, the modelling of the interfacial terms concentrated on the front is carried out with Dirac delta distribution while the discontinuity in material properties is conveyed through the Heaviside step function [5].

Model implementation
The simulation framework was articulated in MATLAB environment, with a cascade model hinging on a flow solver and on a code devoted to the resolution of the CLSM equations so as to model the moving fronts related to the process described by the first PDEs system.
In particular, bidimensional NSE in non-dimensional form were solved by employing Chorin's projection on a Marker-and-Cell (MAC) staggered grid. The relative code is based on the NS solver programmed by Seibold [7] and utilizes a three-steps semi-implicit scheme for time discretization (explicit treatment of the nonlinear convective term, implicit handling of the diffusive term and pressure correction). Spatial discretization for the portion of the model concerning CSLM method is carried out with an upwind second-order Essentially Non-Oscillatory (ENO) scheme whilst time discretization is performed with a three-steps, second-order Total Variation Diminishing (TVD) Runge-Kutta scheme. The choice of upwind approximations is motivated by the hyperbolic nature of the PDE to avoid numerical instabilities. All PDEs were treated by relying on finite difference approximations.
A crucial point regards the choice of the mesh size, as a tradeoff between volume conservation and accuracy emerges: if on the one hand smaller ε implies better volume preservation and compliance with the CFL condition, on the other it entails a decrease in the order of accuracy. Exceeding in setting ε to a small value can determine the formation of spurious oscillations, which in turn may damage the normal field by switching its direction. Different ways to tackle this problem were proposed. We opted for a remapping of ψ by first restoring the steep profile of the Heaviside step function, and then smearing it out with a Gaussian filter (Fig. 2). Since this causes volume loss, in order not to waste the benefits of using CLS, we coupled it with a correction by translation of the level set function (performed through bisection method to identify the best translation value).
Simulations were carried out on a regular square grid made up by 300 × 300 points, as Ω = dr × dr, with Dirichlet boundary conditions for each edge of the computational domain. dr = 0.110 m is the riser diameter. In particular, no-slip conditions were applied in both directions to the southern and western edges, to the upper half on the eastern boundary. The definition is instead trickier on the northern boundary. We opted to impose vertical component velocity of the same value of the mean liquid velocity measured experimentally in the riser vL,r added with Gaussian White Noise (GWN, zero mean, 0.05 vL,r variance), whilst pure GWN was used to introduce small aleatory deviations in the x-direction. The same stratagem was employed to model the components of the imposed velocity profile on the eastern edge in the lower half. In particular, uy is described as GWN and ux the sum of mean liquid velocity in the collector vL,c and GWN.
Bubbles were initially imposed as perfect circles. For what regards the data used in the simulation, dynamic parameters were retrieved from measurements in the experimental campaign [1] as well as the ALR dimensions ( Table 1). The remaining information necessary for using the model was gathered from literature and is reported in Table 2.

Results and discussion
In order to maximize mass transfer, it was mentioned that the maximum surficial area of exchange is required and that the longer the transport phenomenon lasts, the better it is. Thus, hydrodynamics of the multiphase flow holds a controlling influence. Keeping in mind the aim of this study, two types of flow were considered: bubble and churn flows. The first, obtained as long as the gas inlet velocity is maintained below a threshold value, dependent on the tube geometry, is characterized a mostly rectilinear ascent path of the bubbles, which rise almost individually without significant interactions between them and with narrow bubble size distribution. In this flow condition, values of the diameter of the bubble generally fall within the range 1 ÷ 7 mm. Whenever the gas phase velocity exceeds the aforementioned threshold, the density of the gaseous fraction in the liquid gradually increases, resulting in greater interaction between the bubbles, with collisions, clusters formation and the occurrence of coalescence phenomena. The consequential appearance of larger bubbles significantly alters the hydrodynamic scenario, with the concomitant presence of large (more than 20 mm) and small bubbles. This latter rise rather fast (1 ÷ 2 m s -1 ) stirring the liquid while the larger ones tend to churn up the liquid [8]. In this state, as the corresponding Reynolds numbers prove to be higher, spiraling and zigzagging motions can be observed. As researches in this lead to an organic classification of bubble shapes (spherical, ellipsoidal, dimpled ellipsoidal-cap, skirted, spherical-cap [9]), it is possible to surmise that the preferable one is that of a spherical cap, given the high ratio between the exchange surface and the occupied volume. Proceeding with the reasoning, a zigzagging trajectory seems to be preferable, as it would extend the permanence time of the bubbles in the riser.
The simulations regarded the instantaneous puff behavior for a multiple bubbles flow inside the EL-ALR riser, taking into account the geometry, and in particular the position of the sparger and of the horizontal collectors. In the simulations, we were able to identify bubble shapes coherent with those belonging to the classes found experimentally. Coalescence and breaking phenomena were successfully modeled as well, as evidenced by the formation of spherical caps as a result of bubble merging. In addition, the correct definition of the characteristics of the flow field highlighted the zig-zagging ascending trajectories desired and obtained in the experimental tests. Finally, the employed model based on CLSM allowed overcoming volume loss limitations.