Optimal Design of Double Stage Internal Loop Air-Lift Bioreactor

: Bioreﬁnery systems play a critical role in the transition towards a sustainable bioeconomy, and bioreactors are a key component in these systems. While mechanically stirred reactors have been extensively studied, there is a lack of research on pneumatically driven systems like air-lift reactors (ALRs). This study aims to address this gap by examining the hydrodynamic behavior of a double draft tube airlift bioreactor using Computational ﬂuid dynamics simulations. Ten different geometric conﬁgurations were investigated, with variations in draft tube placement, liquid height, distance between draft tubes and draft tube diameters. Results showed that the placement of the draft tubes heavily inﬂuenced hydrodynamic behavior, with smaller distances between draft tubes and a funnel conﬁguration leading to higher velocities. Stable downcomer velocities were achieved by maintaining a consistent distance between the bottom clearance and the sum of the distance between draft tubes and the bottom clearance on the top clearance. The model was validated against literature experimental data. This study provides insight into the optimal design of ALRs, which can contribute to the development of more efﬁcient and effective bioreactor systems. The ﬁndings can be used to forecast the most optimal conﬁgurations of airlift bioreactors and have signiﬁcant value for the development of more efﬁcient bioreﬁning concepts in light of the increasing importance of studying bioreﬁneries and their components in the shift towards a biomass-based economy.


Introduction
In order to transition to the bioeconomy, state-of-the-art biorefinery systems are needed.Biorefineries offer potential pathways for producing food, feed, organic chemicals, fuels, polymers, and electricity from biomass with complex processing technologies.Among these complex processing technologies, the design and analysis of bioreactors are crucial to the successful operation of biological processes.Effective mixing is crucial in the design of bioreactors as it enhances mass transfer within and between phases, which is essential for their optimal performance.There are many studies focused on mechanically stirred systems and barely any on pneumatically driven ones.Among the pneumatical stirred reactors the Air Lift Reactor (ALR) has withdrawn attention due to the low shear forces while mixing, low energy consumption, and reasonable heat and mass transfer.
The ALR was first introduced by Lefrancois, Mariller, and Mejane [1].This type of reactor is widely utilized in the wastewater treatment, biochemical, chemical, and pharmaceutical industries.There are two types of ALRs.The first type is the internal loop ALR, where an insert (baffle or draft tube) is utilized to separate the riser from the downcomer.The second type is the external loop ALR, where the riser and downcomer are two independent vertical compartments that are connected by horizontal sections on the top and bottom.In this study, the focus is centered on a modified Internal Loop Air Lift Reactor (ILALR).
Energies 2023, 16, 3267 3 of 21 predict the local flow structure, the riser was modeled correctly, but an underprediction was found on the downcomer.
Computational fluid dynamics (CFD) is an interdisciplinary approach that combines physics, numerical mathematics, and computer science to model fluid flow phenomena [19].The first CFD simulations on ALRs were performed by Becker, Sokolichin and Eigenberger [20] as a split reactor, utilizing the commercial software Fluent.Cockx et al. [21] utilized a 1D CFD model based on mass and momentum two-phase flow balances which predicted correctly the liquid circulation and gas hold-up.
Padial et al. [22] performed the first 3D CFD simulation of an ALR with conical bottom and reported that their model gave good qualitative results and the liquid circulation time was similar to the experiments performed by Pironti et al. [15] when using only the downward velocity.On other geometries, a 2D CFD simulation with an Eulerian-Eulerian approach was performed by Oey et al. [23] to a rectangular ALR.Mudde and van den Akker [24] performed two and three-dimensional simulations of an internal loop ALR on the basis of a two-fluid model and compared the simulations to Laser Doppler Anemometry (LDA) experiments.At lower superficial velocities, LDA findings match the outcomes of 2D simulations, while at higher gas flow rates, the LDA results shift towards the 3D findings [24].
Different ALR geometries have been proposed.Starting with dividing the single draft tube geometries into two or more parts, known as multi-stage ALR.The first appearance in the literature of the multi-stage internal loop airlift reactor was in the works of Blenke [2] and Chistit and Young [3].In the review of Petersen and Margaritis [25] it appears as a staged internal loop gas lift bioreactor.The division of the draft tube to communicate the riser and downcomer can be categorized under multi-stage in-series ALR.A patent was filed in 2003 [26] but then abandoned.On the patent by Ding et al. [26], the multi-stage ALR also appears in parallel.Tunthikul, Wongsuchoto and Pavasant [27] compared a single draft tube with triple and quadruple parallel draft tube multi-stage ALRs and found that a more efficient recirculation of fluid was achieved when increasing the number of draft tubes.The first study on parallel two-stage ALR reported on literature compared the ALR without, with one, and with two draft tubes for investigating the mixing time and oxygen transfer characteristics [28].Behin [29] performed experiments on a double draft tube ILALR in parallel and found that with this configuration mixing was reduced by up to 48% and circulation time by 35.5%.This type of parallel draft tube configuration was patented by Oldani [30], with the geometrical variation of a diameter increase at the top of the draft tube.Another interesting configuration was studied by Bakker et al. [31] which consisted of a three-concentric compartment multi-layer air-lift reactor.Finally, Bun et al. [32] proposed a baffled modified ALR, which resulted in improved oxygen transfer, achieved by maintaining bubble size and prolonging the gas-liquid transfer period.
There are different applications for ALRS.Karim, Thoma, and Al-Dahhan [33] and Coughtrie, Borman, and Sleigh [34] utilized an internal loop ALR to simulate with CFD the hydrodynamics of an anaerobic digester.Aslanbay Guler et al. [35] utilized CFD to simulate an ALR as a photobioreactor for astaxanthin production.Ruitenberg, Schultz and Buisman [36] designed a novel ALR for biological processes for the mining industry to improve conditions for biological processes.
Luo and Al-Dahhan [47] studied the local characteristics of hydrodynamics in an ALR with a draft tube and found that the bottom and top clearances have significant effects on the reactor's flow.Luo and Al-Dahhan [48] studied the pressure fluctuations in an internal loop airlift bioreactor and found that their CFD model provided a reasonable prediction of microorganism recirculation around the draft tube, but it overestimated the residence time of cells in the wall regions.
Recent investigations have not only been limited to CFD studies.Fadavi and Chisti [8], studied mixing and gas hold-up in a forced circulation loop reactor for various combinations of liquid (water) and gas (air) flow rates.The authors [8] found that when operating the reactor under forced flow conditions, a bubble flow regimen was observed and, when working as an airlift reactor, it produced a heterogeneous churn-turbulent flow regimen.
Karim, Thoma and Al-Dahhan [33] simulated a gas lift anaerobic digester with CFD considering a single phase in a 2D grid.The poorly mixed zones in the bottom of the digester were reduced from 33.6% to 31.9% and 29.6% respectively by varying the bottom from flat to 25 • and 45 • and adding a hanging baffle [33].A similar study performed by Zhang et al. [49] also utilized an angle on the bottom and a funnel instead of a baffle.Zhang et al. [49] used CFD to forecast the hydrodynamics of an improved ALR.In their study a funnel is added on top of the draft tube before the surface, with this the gas hold-up was increased by up to 15% and the turbulent kinetic energy reduce by up to 7.5%, resulting ideal for shear-sensitive biological applications [49].
Pawar [50] demonstrated with CFD that the Lagrangian approach can represent the global flow characteristics of an ALR.The author [50] also identified that the high shear zones were on top of the draft tube and near the sparger.Coughtrie, Borman and Sleigh [34] performed a study with a 2D CFD approach to compare experimental values with the commercial software Fluent and the open-source software OpenFOAM ® .Their results showed that single-phase modeling could be utilized to give reliable predictions of the reactor's flow [34].
Gemello et al. [51] performed a CFD-based scale-up of hydrodynamics and mixing in bubble columns and mentioned that carrying out CFD simulations utilizing LES turbulence modeling could result in more comprehensive conclusions.McClure et al. [52], utilized CFD to analyze the fluid flows within an airlift system and found that a single bubble-size model provided fair predictions of the experimental results, except for underestimating the amount of air dragged in the downcomer section.Aslanbay Guler et al. [35] performed a CFD study on ALRs for microalgae cultivation and found that the most important parameter for the scale-up of airlift photobioreactors is height to diameter ratio of the reactor and draft tube structure.In the case of airlift photobioreactors, a more uniform distribution can be achieved by reducing the height or increasing the diameter of the draft tube [35].
In recent years, multi-stage internal loop ALRs have gained research interest due to the growing demand for more efficient reactors.Li and Qi [53] conducted research on the flow patterns and hydrodynamics of a multi-stage internal loop ALR, while Li et al. [54] investigated the behavior of bubbles and the hydrodynamics in a three-phase two-stage internal loop ALR.Tao et al. [55] experimentally investigated the hydrodynamics and mass transfer in a three-stage internal loop ALR, highlighting the differences between gas-liquid bubbly flow and gas-liquid-solid slurry flow.
A new CFD model for anaerobic digestion has been published recently by Dabiri et al. [56] coupling mixing with biochemical reactions based on the benchmark Anaerobic Digestion Model No 1 (ADM1) by Batstone et al. [57].The new OpenFOAM ® solver passiveScalarADM-1Foam developed by Dabiri et al. [56], was inspired by their previous work [58].Their previous work was validated by the first CFD Smooth Particle Hydrodynamics (SPH) of Rezavand et al. [59] that couples mixing and biochemical reactions in AD using a fully Lagrangian method.
Shi et al. [60] used CFD to investigate a two-stage internal loop ALR equipped with a contraction-expansion vane.Their double draft tube geometry without an expansion and contraction vane was taken to compare the performance of a double and single draft tube in our previous publication [61].In our previous work [61], three reactor geometries with single and double draft tube configurations were evaluated with CFD to investigate the flow behavior.This study utilizes the Eulerian-Eulerian two-fluid numerical model to resolve the hydrodynamics of different double-stage internal loop airlift reactor configurations.The study uses the OpenFOAM ® multiphaseEulerFoam solver, which considers two Eulerian phases to analyze the effects of different internal geometries on bioreactor performance.

Materials and Methods
Computational fluid dynamics (CFD) was utilized to resolve the hydrodynamics of 10 different ALR configurations.The particulars of the CFD model will be discussed in the following section, followed by the mathematical models, the geometries and the meshes.

CFD Model
The open-source toolbox OpenFOAM ® version 9 was utilized to resolve the CFD hydrodynamic calculations with the two-fluid model solver multiphaseEulerFoam.The multiphaseEulerFoam solver is based on the Euler-Euler model.The Euler-Euler model is a type of two-fluid model that utilizes the Lagrangian framework to track the motion of individual parcels.To resolve the flow of the studied internal loop airlift reactors, the Large Eddy Simulation (LES) method was coupled with a two-fluid model (Euler-Euler), using the Schiller-Naumman drag model and the NicenoKEqn and ContinuousGasKEqn one equation sub-grid turbulence models.The NicenoKEqn Large-Eddy Simulation (LES) sub-grid model has been utilized to determine the turbulence of the fluid phase.This model consists of a one-equation sub-grid stress model that resolves the continuous phase in a two-phase system, making it ideal for modeling turbulence generated by bubbles [62].The ContinuosgasKEqn LES sub-grid model has been utilized to determine the turbulence of the gas phase.This model consists of a one-equation sub-grid stress model to resolve the dispersed phase (gas) in a two-phase system accounting for phase-inversion [63].For the phase properties, the drag of each phase is solved by the Schiller-Neuman model [64].A single bubble-size is used since a single bubble-size model provides fair predictions of experimental results according to a study performed by McClure [52].Regarding the chosen bubble size, according to van Baten, Ellenberger and Krishna [65], the rise velocity is practically independent of the bubble size in the 3 to 8 mm range.The diameter of the bubbles was set to 4.5 mm since literature values are in the range of 3 mm to 5 mm.While the diameter of the water particles was left as the solver's default (0.1 mm).

Mathematical Models
Large Eddy Simulations (LES) applied to a Two-Phase Fluid Model (Euler-Euler) for an Air Lift Reactor (ALR) is a simulation approach used to study the behavior of turbulent multiphase flows in ALRs.In an ALR, a gas is utilized to lift and mix a liquid phase, in the process creating a highly turbulent flow.The Two-Fluid Model separates the liquid and gas phases and considers their interactions through a set of governing equations and additional source terms [66].In the Euler-Euler model, the fluid equations for each phase are solved independently, taking into account interphase transfer through the introduction of interphase transfer terms.LES is utilized to resolve the large-scale motions of the fluid while filtering out the small-scale turbulence.This provides a more accurate representation of the turbulence, making LES well-suited for simulating complex, multi-scale phenomena in ALRs.The Schiller-Naumman Drag, the turbulence sub-grid models NicenoKEqn and ContinuousGasKEqn one equation are all commonly used models for describing the drag force between the gas and liquid phases in pneumatical stirred reactors [66,67].These models provide the means to account for the transfer of momentum between the phases and are crucial for obtaining accurate predictions of the flow and transport processes in the reactor.
The continuity equation describes the conservation of mass in each phase (Equation ( 1)) and the momentum equation describes the conservation of momentum (Equation ( 2)) for a two-phase flow in the multiphaseEulerFoam solver: where, α q , ρ q , and u q are the volume fraction, density and phase average velocity of the liquid phase (q = l) and the gas phase (q = g); τ q is the interphase momentum transfer, u p is the velocity of the p-phase, and the interfacial forces between the phases are represented in Equations ( 3)-( 5).The Schiller-Naumann drag model [64] is based on the empirical relationship between the interphase transfer rate and the relative velocity between the two fluids.It provides an efficient estimate of the interphase transfer rate.The Schiller-Naumann drag force, the Reynolds number (bubble based) and the lift can be calculated with the following equations: where, C D is the drag coefficient, Re the Reynolds number, d B is the bubble diameter, U R is the relative velocity between the phases, and v L is the kinematic viscosity of the liquid phase, C L is the lift coefficient, ∝ is the volume fraction of the dispersed phase, ρ q is the density of the dispersed phase, and U L is the local liquid velocity.The Niceno one-equation turbulent kinetic energy LES sub-grid scale model (NicenoK-Eqn) [62] is a turbulence-based model that describes the interphase transfer rate in terms of the turbulence kinetic energy and turbulence dissipation rate of the gas phase.It is well-suited for high-resolution simulations, such as LES, as it provides a more accurate representation of interphase transfer in turbulent flows.The transport equation for k SGS for multiphase flows is the following: where, P k = µ q,SGS S : S (11) Energies 2023, 16, 3267 where, α inv is the alpha inversion and C k , C p , C mub , C ε are model constants, ∝ q is the volume fraction of the q-phase, ρ q is the density of the q-phase, k SGS is the sub-grid scale turbulent kinetic energy, u q is the q-phase velocity, D k is the eddy viscosity, µ q is the molecular viscosity of the q-phase, µ q,SGS is the SGS viscosity of the q-phase, P k is the SGS production term, which is related to the SGS stresses, and S is the strain rate tensor and SGS is a dissipation term.The continuous gas one-equation turbulent kinetic energy LES sub-grid scale model (ContinuousGasKEqn) [63] is based on a set of differential equations that describe the interphase transfer rate in terms of the turbulence kinetic energy and turbulence dissipation rate of the gas phase.This model provides a detailed description of interphase transfer, making it ideal for simulations that require a higher level of accuracy.The ContinuousGasKEqn One-Equation Model in the multiphaseEulerFoam solver is as follows: where, k s,g is the sub-grid scale turbulent kinetic energy, k l is the liquid-phase turbulent kinetic energy, k g is the gas-phase turbulent kinetic energy, C ptr is a function of the gasphase density, µ t,q is the turbulent viscosity for the continuous gas kinetic energy equation, ∆ is the filter width, ∂ t is the time step, C k , C , ∝ inv.g , and C mub are constants.By combining LES with a Two-Phase Fluid Model (Euler-Euler), it is possible to obtain a more detailed and accurate representation of the flow and transport processes in an ALR.This can be useful for optimizing the design and performance of these reactors, as well as for understanding and predicting the behavior of multiphase flows under different operating conditions.

Geometries
In our previous publication [61], three different vessel geometries were evaluated with a single and double draft tube configuration.Since the double draft tube geometries resulted with better mixing performance than the ones with a single draft tube, this study focuses on varying the distance between draft tubes and the diameters of these inserts.
The ten air lift reactor geometries on this study are based on the two geometries shown in Figure 1a,b.Figure 1a is based on Shi et al. [60] and Figure 1b consists of the same cylindrical shape with a coned bottom to transition from the inlet to the main body of the reactor.
All internal geometrical variations, the liquid height, the bottom and top clearance distances are presented on Table 1.
Figure 2 shows all ten geometries utilized in this study.Geometries one and three were taken from Ramonet et al. [61], variations were made to them to investigate how the variation of liquid height, distance between draft tubes and acylindricity by varying the diameters of the draft tubes affects the internal hydrodynamics of an ALR.Geometry three is based on [60], while geometry one is an adaptation of the same adding a coned bottom inspired on the reactors utilized in waste water treatment.Geometries four and five are the only ones where the top draft tube moves.In all the other geometries, when the distance between draft tubes is reduced, the bottom draft tube changes origin.All internal geometrical variations, the liquid height, the bottom and top clearance distances are presented on Table 1. Figure 2 shows all ten geometries utilized in this study.Geometries one and three were taken from Ramonet et al. [61], variations were made to them to investigate how the variation of liquid height, distance between draft tubes and acylindricity by varying the diameters of the draft tubes affects the internal hydrodynamics of an ALR.Geometry three is based on [60], while geometry one is an adaptation of the same adding a coned bottom inspired on the reactors utilized in waste water treatment.Geometries four and five are the only ones where the top draft tube moves.In all the other geometries, when the distance between draft tubes is reduced, the bottom draft tube changes origin.A mesh independency study has been performed on a similar geometry without any draft tubes.A refinement level (r) of 1.2 was utilized between meshes.Table 2 presents the four meshes with their respective cell count, maximum skewness detected and the extracted air velocity magnitude from the first peak marked with red square on Figure 3.A mesh independency study has been performed on a similar geometry without any draft tubes.A refinement level (r) of 1.2 was utilized between meshes.Table 2

presents
Energies 2023, 16, 3267 9 of 21 the four meshes with their respective cell count, maximum skewness detected and the extracted air velocity magnitude from the first peak marked with red square on Figure 3.The obtained value for M0 was 0.363951 when compared to mesh 4 (0.365633), a relative error of discretization (Equation ( 20)) of 0.46% is obtained.The Relative Discretization Error (RDE) is the difference between the numerical solution and the exact solution [70].On our previous publication [61], an error of 0.38% was obtained for a similar geometry with a refinement level of 1.5.Geometries 1, 3 and 6 from this study were taken from that previous work, making all the grids utilized for this study with a RDE of less than 0.5%.
Once validated that meshes over 300 thousand cells provide accurate results, the snappyHexMesh mesh generation utility from OpenFOAM ® was utilized to mesh the 10 geometries.To find a good combination of the mesh generating parameters for the Open-FOAM ® utilities blockMesh and snappyHexMesh, the trial-and-error methodology is used [71].Due to the internal geometry the cell refinement varies for the coned bottom cases from 600 thousand to 1 million cells to keep a maximum skewness of under 4.
Table 3, shows the resulting number of cells, the maximum skewness and maximum non-orthogonality of the meshes for the 10 geometries on this study.According to [72], the recommended mesh quality in OpenFOAM ® involves a maximum skewness of approximately 4 and a maximum non-orthogonality of 70°.The geometry utilized for cases 3 and 6 has a skew cell on the bottom draft tube which brings the maximum skewness of the geometry to 4.48.Since the location of this skew cell is not on a critical path, this maximum skewness of 4.48 can be tolerated.To compare the order of convergence the air velocity magnitude was extracted with a diagonal line from the inlet to the outlet (plotted in Figure 3).The air velocity magnitude values were extracted from the first peak (marked with the red rectangle in Figure 3)and are presented in Table 2.
The next step in a mesh study is to calculate first the order of convergence, also known as order of accuracy, which is the uncertainty of the numerical solution caused by the spatial discretization errors [68].To calculate the order of convergence the first three meshes are utilized, then to corroborate the Richardson extrapolation is used.Equation ( 18) is utilized to calculate the order of convergence [68].Where M 1 is the finest mesh (3), M 2 is the medium mesh (2) and M 1 is the coarser mesh (1).
The order of accuracy (p) equals −10.42, then with Equation ( 19) the Richardson extrapolation [69] was utilized: The obtained value for M 0 was 0.363951 when compared to mesh 4 (0.365633), a relative error of discretization (Equation ( 20)) of 0.46% is obtained.The Relative Discretization Error (RDE) is the difference between the numerical solution and the exact solution [70].On our previous publication [61], an error of 0.38% was obtained for a similar geometry with a refinement level of 1.5.Geometries 1, 3 and 6 from this study were taken from that previous work, making all the grids utilized for this study with a RDE of less than 0.5%.
Once validated that meshes over 300 thousand cells provide accurate results, the snappyHexMesh mesh generation utility from OpenFOAM ® was utilized to mesh the 10 geometries.To find a good combination of the mesh generating parameters for the OpenFOAM ® utilities blockMesh and snappyHexMesh, the trial-and-error methodology is used [71].Due to the internal geometry the cell refinement varies for the coned bottom cases from 600 thousand to 1 million cells to keep a maximum skewness of under 4.
Table 3, shows the resulting number of cells, the maximum skewness and maximum non-orthogonality of the meshes for the 10 geometries on this study.According to [72], the recommended mesh quality in OpenFOAM ® involves a maximum skewness of approximately 4 and a maximum non-orthogonality of 70 • .The geometry utilized for cases 3 and 6 has a skew cell on the bottom draft tube which brings the maximum skewness of the geometry to 4.48.Since the location of this skew cell is not on a critical path, this maximum skewness of 4.48 can be tolerated.The time step of the simulations was set to adjustable, which means that the time step is self-calculated by the algorithm depending on the timestep's Courant number.The Courant number [73], also known as the Courant-Friedrichs-Lewy (CFL) condition, consists in imposing a stability condition when solving a numerical equation that limits the size of the time step based on the chosen spatial discretization.
The equations are resolved in OpenFOAM ® using the finite volume method, and the coupling between pressure and velocity is resolved using the PIMPLE algorithm [74].To assure numerical stability and convergency the following criteria has been applied.To guarantee stability the time step is set to adjustable, which is regulated by the algorithm to have a CFL under 0.8.The convergence criteria for the matrix solver were set to 1 × 10 −8 for each variable.
The ten meshes are shown in Figure 4.All meshes have wall refinement except for geometries 4 and 5. Geometries 1, 2, 3, 6, 7 and 9 have 600 thousand cells; geometries 4 and 5 have 300 thousand cells and geometries 8 and 10 have 1 million cells.Geometries 4 and 5 were the only two without wall refinement.For the other geometries a wall refinement of 1 level and 2 cells in between levels.To illustrate the mesh, geometry 8 has been chosen for a close up to view the wall refinement at the draft tube section, inlet, outlet and walls, as shown in Figure 5. Geometries 4 and 5 were the only two without wall refinement.For the other geometries a wall refinement of 1 level and 2 cells in between levels.To illustrate the mesh, geometry 8 has been chosen for a close up to view the wall refinement at the draft tube section, inlet, outlet and walls, as shown in Figure 5.To calculate the loop circulation time, values were extracted from the riser and downcomer.Nine points were sampled on the riser (marked with a U on Figure 6b) and other nine on the downcomer (marked with D on Figure 6b).Based on extracted velocity values, the velocity was interpolated from one point to another until forming a loop.A similar method was reported in literature, where the total circulation time was calculated by adding the ascending and descending time for a single draft tube geometry [22].In our case, the connection between the upcomer and downcomer is also calculated by interpolating the velocity of the vector moving in this direction.To calculate the loop circulation time, values were extracted from the riser and downcomer.Nine points were sampled on the riser (marked with a U on Figure 6b) and other nine on the downcomer (marked with D on Figure 6b).Based on extracted velocity values, the velocity was interpolated from one point to another until forming a loop.A similar method was reported in literature, where the total circulation time was calculated by adding the ascending and descending time for a single draft tube geometry [22].In our case, the connection between the upcomer and downcomer is also calculated by interpolating the velocity of the vector moving in this direction.
comer.Nine points were sampled on the riser (marked with a U on Figure 6b) and other nine on the downcomer (marked with D on Figure 6b).Based on extracted velocity values, the velocity was interpolated from one point to another until forming a loop.A similar method was reported in literature, where the total circulation time was calculated by adding the ascending and descending time for a single draft tube geometry [22].In our case, the connection between the upcomer and downcomer is also calculated by interpolating the velocity of the vector moving in this direction.

Results
Two points were utilized on the top and bottom of the downcomer side to extract data for each case.These points were selected at approximately the same height; the

Results
Two points were utilized on the top and bottom of the downcomer side to extract data for each case.These points were selected at approximately the same height; the pressure varies from one case to another depending on the liquid height.The axial velocity is presented on m/s as an absolute value.The turbulent kinetic energy is also presented in Table 4.The draft tube circulation time was calculated for each case.All simulations were performed on 4 cores, the computation time varied from 81 to 438 h, and the adjustable time step was used in all cases.In Table 4, the column marked as first draft tube indicates the time it takes for a particle to go from the inlet all the way up, travel down until the top draft tube finishes, and incorporates itself back into the flow forming a first loop as shown in Figure 6a with the blue arrows.The column marked as second draft tube indicates the time it takes for a particle to travel from the inlet all the way up to the interphase section and back down to the bottom of the lower draft tube forming a loop as shown in Figure 6a with the black arrows.
From geometry 1 to 2, the draft tube was moved up by 8 cm, which was also done with geometry 7 to 9. In literature, this phenomenon was studied experimentally with a single draft tube geometry with conical bottom by Pironti et al. [15], who found that the draft tube placement heavily affects the hydrodynamic behavior.When analyzing the data in Table 4, it can be seen that the top velocity on geometry 9 is more than double from its reference geometry (7), which can relate to the double velocity gradient achieved by the smaller distance between the draft tubes and the double funnel configuration.On the other hand, the downcomer velocity is one third of its counterpart, making the bottom of the reactor a less active zone.Out of the six configurations with coned bottom, geometry 9 is the one with the fastest recirculation loop on the upper loop of the reactor (top draft tube).The fastest circulation for the second loop on the conical bottom geometries was seen on geometry 7.For the cylindrical draft tube (geometries 1 and 2), the only remark is the decrease of velocity on the bottom of the reactor on the measuring point caused by the movement of the draft tube.
When comparing by reactor type, the conical bottom reactor shows the following.Cases 1, 2, 7, 8, 9, and 10 all have similar values for pressure and draft tube loop circulation time, with pressure ranging from 100,374 to 100,402 Pa and draft tube loop circulation time on the first loop ranging from 21.4 to 36.97 s.However, there are noticeable differences in axial velocity and turbulent kinetic energy.Cases 1 and 2 have the highest values for turbulent kinetic energy, while cases 9 and 10 have the highest values for axial velocity.Case 8 has the highest value for turbulent kinetic energy and the second-highest value for axial velocity.Case 7 has the lowest values for axial velocity and turbulent kinetic energy.Overall, it appears that the design parameters used in cases 1, 2, 8, 9, and 10 were more effective at achieving high axial velocities and turbulent kinetic energies than the design parameters used in case 7.
From cases 3, 4, 5, and 6, cases 3 and 4 have the same liquid height, and cases 5 and 6 have the same liquid height.Cases 3 and 6 have identical geometry, and cases 4 and 5 have identical geometry.Looking at the results, case 3 has the highest axial velocity, turbulent kinetic energy, and circulation time for the first draft tube loop, indicating that it has the most efficient fluid circulation.Case 4 has the lowest axial velocity and turbulent kinetic energy for both the top and bottom sections, indicating less efficient mixing.Case 5 has the highest turbulent kinetic energy for the top section, but the lowest for the bottom section, indicating a non-uniform flow field.Case 6 has the highest pressure and circulation time for the second draft tube loop, indicating that it has the most efficient flow in this loop.In summary, case 3 has the most efficient fluid circulation and mixing, followed by case 6, case 5, and finally case 4.
The following section shows with Figures 7-9 the averaged velocities on all geometries at the extracted points.On the downcomer section, the first six seconds were not considered since there was no flow, the values were taken since the moment the downcomer flow showed a stabilized movement.
Figure 7 compares the axial velocity on the riser and downcomer on geometries 1 and 2. From this figure it can be deduced that reducing the distance between draft tubes by half by moving the bottom draft tube does not have a significant effect on the overall hydrodynamics of the reactor since the curves present very similar behavior.
On Figure 8a, geometries 5 and 6 show a higher rising velocity on the bottom of the draft tube.The axial velocity on the riser, becomes the approximately the same on the four cases on the draft tube section.Further, ones the flow leaves this section, geometries 5 and 6 show higher velocities again.Geometries 3 and 4 share the same liquid height and show a similar tendency.
Geometry 3 and 6 have the same geometrical characteristics, what changes is the liquid height.Geometry 4 and 5 have the same geometrical characteristics, what changes is the liquid height.In Figure 8a, it can be observed that the riser is affected mainly by the liquid height, since the geometries with the same liquid height have the same tendency.
On the other hand, the downcomer velocities shown in Figure 8b have the same behavior at the bottom of the reactor, but this changes on the top part.Geometry three shows a linear decrease in velocity as the flow goes downwards.Geometry 5, shows a more stable velocity on the upper part of the reactor, which could be beneficial for microorganism cultures.
Considering that in these four cases the variations are the liquid height and the distance between draft tubes, the most stable downcomer velocity case was when having approximate the same distance between the bottom clearance and distance between the draft tubes and approximate the sum of these two values on the top clearance.
From Figure 9a it can be clearly seen that geometry 7 performs worse than the other 3 cases on the upper half of the reactor, since it has a lower riser velocity.Geometries 8 and 10 show a similar behavior on the riser, although geometry 8 shows a higher velocity at the top of the reactor which can be due to the contraction-expansion effect.
On the downcomer velocities (Figure 9b), all geometries show a similar behavior at the bottom of the reactor.At the top of the reactor, geometry 9 doubles the velocity of geometries 8 and 10, and almost triples the velocity of geometry 7. Geometry 7 shows the more stable on the top half of the reactor's liquid phase for the downcomer velocity, but the lowest on this same area on the riser.The following section shows with Figures 7-9 the averaged velocities on all geometries at the extracted points.On the downcomer section, the first six seconds were not considered since there was no flow, the values were taken since the moment the downcomer flow showed a stabilized movement.

Top
Figure 7 compares the axial velocity on the riser and downcomer on geometries 1 and 2. From this figure it can be deduced that reducing the distance between draft tubes by half by moving the bottom draft tube does not have a significant effect on the overall hydrodynamics of the reactor since the curves present very similar behavior.On Figure 8a, geometries 5 and 6 show a higher rising velocity on the bottom of the draft tube.The axial velocity on the riser, becomes the approximately the same on the four cases on the draft tube section.Further, ones the flow leaves this section, geometries 5 and 6 show higher velocities again.Geometries 3 and 4 share the same liquid height and show a similar tendency.
Geometry 3 and 6 have the same geometrical characteristics, what changes is the liquid height.Geometry 4 and 5 have the same geometrical characteristics, what changes is the liquid height.In Figure 8a, it can be observed that the riser is affected mainly by the liquid height, since the geometries with the same liquid height have the same tendency.On the other hand, the downcomer velocities shown in Figure 8b have the same behavior at the bottom of the reactor, but this changes on the top part.Geometry three shows a linear decrease in velocity as the flow goes downwards.Geometry 5, shows a more stable velocity on the upper part of the reactor, which could be beneficial for microorganism cultures.
Considering that in these four cases the variations are the liquid height and the distance between draft tubes, the most stable downcomer velocity case was when having approximate the same distance between the bottom clearance and distance between the draft tubes and approximate the sum of these two values on the top clearance.From Figure 9a it can be clearly seen that geometry 7 performs worse than the other 3 cases on the upper half of the reactor, since it has a lower riser velocity.Geometries 8 and 10 show a similar behavior on the riser, although geometry 8 shows a higher velocity at the top of the reactor which can be due to the contraction-expansion effect.
On the downcomer velocities (Figure 9b), all geometries show a similar behavior at the bottom of the reactor.At the top of the reactor, geometry 9 doubles the velocity of  Figure 10 colors the axial velocity of all ten cases.The comparison is divided as in the previous section, a first comparison of geometries 1 and 2 (Figure 10a), a second of geometries 3-6 (Figure 10b) and, a third of geometries 7-10 (Figure 10c).The riser velocity is represented with the positive values from orange to red, while the downcomer velocity is represented with negative values from green to blue.
On geometries 1 and 2 (Figure 10a), it can be seen that the downcomer enters the bottom of the draft tube.This can also be seen on most of the geometries on Figure 10b.This counterflow has been eliminated by transforming the bottom draft tube into a funnel (reducing the diameter of the top part of the draft tube) as shown in Figure 10c.To represent the flow movement the axial water velocity of geometry 9 is captured every 5 s in Figure 11, colored with a maximum riser velocity of 0.05 m/s and a maximum downcomer velocity of 0.03 m/s. Figure 10 colors the axial velocity of all ten cases.The comparison is divided as in the previous section, a first comparison of geometries 1 and 2 (Figure 10a), a second of geometries 3-6 (Figure 10b) and, a third of geometries 7-10 (Figure 10c).The riser velocity is represented with the positive values from orange to red, while the downcomer velocity is represented with negative values from green to blue. Figure 10 colors the axial velocity of all ten cases.The comparison is divided as in the previous section, a first comparison of geometries 1 and 2 (Figure 10a), a second of geometries 3-6 (Figure 10b) and, a third of geometries 7-10 (Figure 10c).The riser velocity is represented with the positive values from orange to red, while the downcomer velocity is represented with negative values from green to blue.
On geometries 1 and 2 (Figure 10a), it can be seen that the downcomer enters the bottom of the draft tube.This can also be seen on most of the geometries on Figure 10b.This counterflow has been eliminated by transforming the bottom draft tube into a funnel (reducing the diameter of the top part of the draft tube) as shown in Figure 10c.To represent the flow movement the axial water velocity of geometry 9 is captured every 5 s in Figure 11, colored with a maximum riser velocity of 0.05 m/s and a maximum downcomer velocity of 0.03 m/s.On geometries 1 and 2 (Figure 10a), it can be seen that the downcomer enters the bottom of the draft tube.This can also be seen on most of the geometries on Figure 10b.This counterflow has been eliminated by transforming the bottom draft tube into a funnel (reducing the diameter of the top part of the draft tube) as shown in Figure 10c.
To represent the flow movement the axial water velocity of geometry 9 is captured every 5 s in Figure 11, colored with a maximum riser velocity of 0.05 m/s and a maximum downcomer velocity of 0.03 m/s. Figure 12, colors the downcomer axial velocity of geometries 1, 2, 7, 8, 9 and 10.To ignore the riser velocity, the maximum velocity was set to zero and the minimum velocity was set to −0.01 m/s.Selecting a darker color for the minimum velocity gives an adequate graphical representation of the maximum downcomer velocity.On the left of Figure 12, it can be seen that the downcomer invades the riser section of the bottom draft tube geometries.Reducing diameter of the top section of the bottom draft tube eliminates this phenomenon, as shown on geometries 7-10 of Figure 12.On Figure 12, it can be seen that when there is the presence of a funnel draft tube on the bottom of the reactor a dead zone is created on the corners before the conical bottom transition.Figure 12, colors the downcomer axial velocity of geometries 1, 2, 7, 8, 9 and 10.To ignore the riser velocity, the maximum velocity was set to zero and the minimum velocity was set to −0.01 m/s.Selecting a darker color for the minimum velocity gives an adequate graphical representation of the maximum downcomer velocity.On the left of Figure 12, it can be seen that the downcomer invades the riser section of the bottom draft tube geometries.Reducing diameter of the top section of the bottom draft tube eliminates this phenomenon, as shown on geometries 7-10 of Figure 12.On Figure 12, it can be seen that when there is the presence of a funnel draft tube on the bottom of the reactor a dead zone is created on the corners before the conical bottom transition.Figure 12, colors the downcomer axial velocity of geometries 1, 2, 7, 8, 9 and 10.T ignore the riser velocity, the maximum velocity was set to zero and the minimum veloci was set to −0.01 m/s.Selecting a darker color for the minimum velocity gives an adequa graphical representation of the maximum downcomer velocity.On the left of Figure 12, can be seen that the downcomer invades the riser section of the bottom draft tube geom etries.Reducing diameter of the top section of the bottom draft tube eliminates this ph nomenon, as shown on geometries 7-10 of Figure 12.On Figure 12, it can be seen th when there is the presence of a funnel draft tube on the bottom of the reactor a dead zon is created on the corners before the conical bottom transition.To represent the flow inside air lift bioreactors 25 particles were taken as a sample to color streamlines with the velocity magnitude in Figure 13.To represent the flow inside air lift bioreactors 25 particles were taken as a sample to color streamlines with the velocity magnitude in Figure 13.To validate our CFD model the squared single draft tube geometry of Ramonet et al. [61] has been compared to the literature it was extracted from Jia et al. [74].As reported on Ramonet et al. [61] for the case with 0.01 m/s inlet velocity the loop circulation time is 16.56 s and for the case with 0.02 m/s inlet velocity the loop circulation time is 12.28 s.Axial velocity values were extracted and averaged from the last 20 s of the simulation.Figure 14 compares these extracted averaged axial velocity values to the ones reported by Jia et al. [74].To validate our CFD model the squared single draft tube geometry of Ramonet et al. [61] has been compared to the literature it was extracted from Jia et al. [74].As reported on Ramonet et al. [61] for the case with 0.01 m/s inlet velocity the loop circulation time is 16.56 s and for the case with 0.02 m/s inlet velocity the loop circulation time is 12.28 s.Axial velocity values were extracted and averaged from the last 20 s of the simulation.Figure 14 compares these extracted averaged axial velocity values to the ones reported by Jia et al. [74].

Discussion
The investigation of the hydrodynamic behavior of a draft tube airlift bioreactor using Computational fluid dynamics simulations was addressed in this paper.Ten different geometric configurations with variations in the placement of the draft tube, liquid height, and distance between draft tubes have been studied.The simulations were performed on

Figure 1 .
Figure 1.Multi stage air-lift reactors: (a) cylindrical double draft tube geometry based on Shi et al. [60] and, (b) coned bottom cylindrical double draft tube geometry.

Figure 1 .
Figure 1.Multi stage air-lift reactors: (a) cylindrical double draft tube geometry based on Shi et al. [60] and, (b) coned bottom cylindrical double draft tube geometry.

Figure 2 .
Figure 2. Geometries utilized in this study.

Figure 2 .
Figure 2. Geometries utilized in this study.

Figure 3 .
Figure 3. Air velocity magnitude from inlet to outlet.Values extracted from first peak (marked in red) for mesh study.

Figure 3 .
Figure 3. Air velocity magnitude from inlet to outlet.Values extracted from first peak (marked in red) for mesh study.

22 Figure 5 .
Figure 5. Visualization of the mesh of geometry 8. Inlet and coned transition marked with 1, bottom draft tube marked with 2 and reactor's top and outlet marked with 3.

Figure 5 .
Figure 5. Visualization of the mesh of geometry 8. Inlet and coned transition marked with 1, bottom draft tube marked with 2 and reactor's top and outlet marked with 3.

Figure 7 .
Figure 7. (a) Axial velocity in riser of geometries 1 and 2 based on 9 sampling points; (b) Axial velocity (absolute value) in downcomer of geometries 1 and 2 based on 9 sampling points.

Figure 7 .
Figure 7. (a) Axial velocity in riser of geometries 1 and 2 based on 9 sampling points; (b) Axial velocity (absolute value) in downcomer of geometries 1 and 2 based on 9 sampling points.

Figure 11 .
Figure 11.Axial water velocity of geometry 9 with 5 s interval for 60 s.

Figure 11 .
Figure 11.Axial water velocity of geometry 9 with 5 s interval for 60 s.

Figure 11 .
Figure 11.Axial water velocity of geometry 9 with 5 s interval for 60 s.

Figure 13 .
Figure 13.Streamlines coloring the velocity magnitude of a sample of 25 particles representing the flow field inside of geometries 1, 2, 7, 8, 9 and 10.
Jia et al. 2007 inlet 0.02 m/s Simulation inlet 0.02 m/s

Figure 13 .
Figure 13.Streamlines coloring the velocity magnitude of a sample of 25 particles representing the flow field inside of geometries 1, 2, 7, 8, 9 and 10.

Figure 14 .
Figure 14.Comparison of CFD model to literature values of Jia et al.Ref. [74]: (a) inlet velocity of 0.01 m/s and (b) inlet velocity of 0.02 m/s.

Figure 14 .
Figure 14.Comparison of CFD model to literature values of Jia et al.Ref. [74]: (a) inlet velocity of 0.01 m/s and (b) inlet velocity of 0.02 m/s.

Table 4 .
Results from the simulations.