Experimental measurement and CFD simulation on the hydrodynamics of an internal-loop airlift reactor

This paper concerns with the experimental measurement and computational fluid dynamics simulation on local hydrodynamics of a gas-liquid internal-loop airlift reactor. The aim of this work is to study the sensitivity of the drag models and the significance of considering the lift force on the predictive accuracy of the simulation. The experimental analysis was carried out using laser Doppler anemometry at three different heights (i.e. Y = 0.20 m, 0.30 m and 0.38 m) across the riser and downcomer at volumetric flow rate of 0.30 m3/h to provide validation for the simulation results. A transient three-dimensional gasliquid internal-loop airlift reactor was carried out using FLUENT 16.2 by implementing the two-fluid model approach. The Eulerian-Eulerian multiphase and standard kdispersed turbulence model were employed in this study. Results suggest that the spherical drag model performed poorly and that the drag model governed by Rayleigh-Taylor shows promising accuracy in the prediction of overall mean axial liquid velocity. On the other hand, the consideration of lift model shows slightly improvement in accuracy. These findings may serve as a guidance for future scale-up and design of airlift reactor studies


Introduction
The airlift reactor is a class of pneumatic reactor that is designed to address the shortcomings found in the bubble column [1].This device is attractive for its lower likelihood for mechanical failures, relatively lower shear rate in the absence of mechanical agitation and minimal loss of sterility [2].Scaling up an airlift reactor is still currently challenging as the complexity of the hydrodynamics is influenced by the geometry of the reactor.
Many empirical correlations have been proposed in the past relating to riser-to-downcomer ratio, gas holdup and superficial gas velocity which were used to aid in current industrial scale-up procedures [3][4][5].These correlations were however derived based on global hydrodynamics measurements which may not reflect the local intrinsic information such as Kolmogorov effect and flow patterns of the fluid flow.Recent breakthrough in optical imaging techniques has made local hydrodynamics studies easier in turbulent gas-liquid dispersion systems.The laser Doppler anemometry (LDA) is one of the non-intrusive methods available to measure local velocities of the flow field.The LDA is capable of studying transient flow relatively close to the wall, providing single point time-averaged measurements [6][7].This technique requires seeding particles which are assumed to track the flow dynamics following the principles of Stokes number [8].LDA has been successfully employed in stirred tanks and bubble columns to measure mean axial liquid velocity profiles despite the challenges of capturing accurate two-phase flow field [6,[9][10][11][12][13].Though experimental using LDA as a measurement method in airlift reactors are still limited [14][15].In addition, both global and local hydrodynamics measurements used to derive empirical correlations for scale-up purposes is accurate only for similar geometrical designs, operating conditions and liquid properties.
Thus numerical approach becomes an attractive alternative method that allows flexible design solutions in a virtual experimentation environment where it screens potential design and configurations.This would reduce fabrication cost and time for re-optimization due to scale-up failures.Computational fluid dynamics (CFD) has garnered interested in the field of numerical studies for its capabilities to provide an in-depth understanding into the hydrodynamics and mass transfer of airlift reactors with reasonable accuracy.The modelling approach of an airlift reactor depends on its flow regime.In bubbly flow, the two-fluid model is widely employed across literature whereby the bubble sizes are assume to be mono-dispersed [16][17][18].The twofluid model implements the Eulerian framework whereby the conservation equations for both the liquid and gas phases are ensemble-averaged reducing the computational cost to predict the gas-liquid flow fields in an airlift reactor.However, this procedure inevitably eliminates microscopic information describing the interacting forces between the gas-liquid phases.Detailed information that describes the contact between the bubbles with liquid phase for effective mass transfer are important for a successful scale-up and design of an airlift reactor.Hence, additional closure models are introduced to account the two-phase momentum exchanges.
As of current, no universal set of closure models have been concluded for bubbly flow in an airlift reactor although a diverse set of closure models have been proposed across literature [19][20][21][22].These closure models which composed of interfacial momentum forces can be divided down further as drag and non-drag forces.
The drag force is known to be the pre-dominant contributing force whereby it determines the velocity of the gas phase and overall flow pattern in the gas-liquid system [23].Past airlift reactor studies have favoured Schiller and Naumann (1935) model over other proposed drag models to predict the airlift reactor hydrodynamics [16,24].This model is however meant for spherical bubbles which may not be suitable for more turbulent flows where bubble may deform due to turbulent forces.Hence, comparison studies have been conducted over the decade to evaluate the sensitivity of these drag models on local hydrodynamics.A study carried out in an internal-loop airlift reactor have shown superior performance through the drag model Karamanev and Nikolov (1992) as a function of Reynolds number when compared with other models governed by Eötvös number [19].Another study showed that the combination of both Reynolds and Eötvös number (i.e.Dijkhuizen et al., 2010) gave improvements in radial profile of the axial liquid velocity prediction [20].However, both the these studies were conducted in a two-dimensional spatial simulation whereby its accuracy can be limited to lower superficial gas velocity and sensitive to gas inlet geometry [22,25].Comparison on drag model performance in a three-dimensional airlift reactor showed better performance in predicting gas velocity when the drag coefficient as a function of superficial phase velocity was used [21].Presently, there is absent in literature on the evaluation of the performance of the Universal drag model which is governed by different flow regimes (i.e.function of Reynolds number, Rayleigh-Taylor instability and gas holdup) in a threedimensional airlift reactor.
On the other hand, the lift force is a non-drag force that is not obvious in comparison of the drag force but has a role in radial gas holdup.A lack of clarity on the effects of considering the lift model has result some studies to neglect it [21,26].A past study found that the lift force assist in correcting the radial profiles of the gas holdup, axial liquid velocity and turbulent kinetic energy [27].In addition, another study pointed that lift should be considered in order to capture the bubble plume oscillations [28].Plume oscillations are responsible in predicting the instantaneous bubble velocity.However, another study conducted in a pipe showed that the effect of lift and turbulent dispersion force cancels out each other near walls in the gas holdup and liquid velocity [29].Due to the discrepancy, much clarity I still needed on the significance of considering the lift model in airlift reactors, hence is studied in our work as well.
The outline of this paper consists of the experimental measurement on local hydrodynamics of an internal-loop airlift reactor as the first part.The measurements are obtained using LDA taken across the riser and downcomer at three different heights and will be used to validate with the simulation sensitivity study of the interfacial forces.The second part concerns with a transient simulation of the internal-loop airlift reactor to evaluate the performance of different drag models and the effect of considering the lift force for a closure model of a two-phase internal-loop airlift reactor.The drag models considered in this study consists of a drag coefficient governed by Reynolds number, Eötvös number and bubble swarm which accommodates the bubble swarm correction model.In addition, drift velocity was considered in this simulation to take into account the effects of fluctuating liquids on the gas phase [30].Air is introduced into the airlift reactor through a 4" Classica air stone with the dimensions of 0.10 m in length and 0.20 m in height located at the base of the reactor.The reactor was filled with water at a height of 0.46 m and aerated with a volumetric flow rate of Qg = 0.30 m 3 /h using a vacuum pump.In order to obtain a non-coalescence system, 4 kg/m 3 of salt was added into the reactor [31].

Experimental procedure
Measurements of the axial liquid velocity component is taken using two component FlowExplorer.LDA system from Dantec Dynamics whereby the system is operating in back-scattered mode.The gas-liquid system was seeded with polyamide seeding particles of diameter size 20 ȝm and density of 1030 kg/m 3 from Dantec Dynamics.The measurements were taken at heights 0.20 m, 0.30 m and 0.38 m of the column as shown in Figure 1(b).These measurements were taken radially at an interval of 10 mm between each points with a 2.5 mm gap away from the wall.Both 15 mW He-Ne and 5 W Ar-ion lasers manufactured by Spectra-Physics were used as light sources emitting wavelengths of 532 nm and 561 nm.The laser intensity is set at 40 mW for both lasers.A 500 mm focal length lens with a probe volume of 2.1 mm 3 was used as a receiving lens in this study.The laser and focal lens were mounted on a bidirectional traverse controlled using a workstation.
Data validation and signal processing were carried out using BSA3 Processor v6.10 software collected from the site using burst mode data collection method.It is worth noting that the data rate of acquisition drops as the laser approaches the bubble region.This is especially observed when the laser moves from the wall to the bubble plume region in the riser.When light hits the bubble, the intensity of the light reflected will be higher [32].However, these lights are reflected and scattered, lowering the probability of the lights hitting the detector.This results in the reduction in sample count obtained.
Several measurements were performed at these locations repetitively by adjusting the gain, sensitivity, and anode to improve the accuracy of the time-averaged values.The duration of time specified to collect the samples is set at 60 s which is sufficient enough to obtain time-averaged value with an accuracy of ±0.01 m/s.Within that timeframe and setup as describe in section 2, 10 2 to 10 3 amount of samples were obtained.Although longer sampling would be encouraged, it was found large data quantities would not result to significant improvements in accuracy [26].

Geometry and computational grid
A transient simulation of a three-dimensional internalloop airlift reactor was carried out using ANSYS FLUENT 16.2.The reactor was prepared in GAMBIT 2.4.6 according to the dimensions as described in section 2 as illustrated in Figure 2 using hexahedral cells.The headspace was extended with an additional 0.30 m to prevent recirculation and complex flow features from affecting flow field results.The governing laws and interfacial forces adopted in this simulation were further elaborated in section 3.2.No-slip boundary conditions were applied on the walls and the top of the headspace region was set as pressure outlet.A gas with a volumetric flow rate, Qg = 0.30 m 3 /h was injected and treated as a continuous source of gas at the sparger.The pressure and temperature were set to ambient.The CFD simulation was performed using mean bubble size correlation [33].All residuals were set below 1x10 -8 for each time step to achieve good convergence.The axial liquid velocity at height 0.30 m was monitored and the iterations were halted once a constant value was observed.The data was averaged once a pseudo-steady condition was achieved.The time step was set at 0.0001 s.The CFD simulation was performed using six units of HP Compaq Pro 6300 MT with a quad core processor (3.2 GHz i7-3770, 4 GB of RAM).

Multiphase modelling
The Eulerian-Eulerian framework was adopted in this simulation work, whereby the mass and momentum equations for the continuous and dispersed phases were weighted by their corresponding volume fractions.Each of the phases are ensemble-averaged and that the total sum of the volume fractions should satisfy to a unity as described in the following governing equations: where α is the volume fraction, ρ is the density, t is the time and u is the velocity vector, respectively.The interphase mass transfer is not considered in this study, hence the mass source term in the right hand side of (1) is neglected.Meanwhile in the first part of the right hand-side of expression (2) describes the pressure gradient, followed by the tensor stress between the phases, gravity and lastly the interfacial momentum transfer between the two-phases, respectively. , The interfacial momentum forces as shown in (3) compose of drag, D F and non-drag forces whereby, the non-drag forces that are considered in this study consists of lift, L F and turbulent dispersion force, TD F .Virtual mass force is neglected in this simulation study as it is computationally costly to simulate with no significant improvements [22].

Turbulence modelling
This work employs a standard k-İ dispersed turbulence model based on Tchen-theory for systems where the concentration of the dispersed phase is dilute.The transport equation for k and İ in the dispersed model are given as: , k l G is the production rate of the turbulent kinetic energy and its form is similar to the one applied for single flow.The values of the constant in (4) and ( 5) are given by 1 C ε = and 2, 1.92 C ε = .

Drag force
Drag force is the resistance force acting on the body (i.e.bubble) moving in a surrounding fluid (i.e.liquid phase).The drag force is governed by inertia whereby the force becomes increasing significant at larger Reynolds number.A layer of viscous stress and pressure distribution forms around the moving body, affecting the bubble size and shape resulting to the drag force embodying the bubble.The drag force can be describe as a function of local velocity as described in (6): where Re b is the bubble Reynolds number defined as

Eo
is the Eötvös number defined as where σ is the surface tension.The Eötvös number measures the distortion of the bubbles based on the gravitational force to surface tension ratio.When the Eötvös number is Eo < 1, the surface tension governs which retain the sphere of the bubble which is likely to be true in disperse homogeneous regime.As the Eötvös number increases, the gravitational force is more prominent resulting to bubble deformation.On the other hand, the Universal drag considers a wide range of flow regimes.The model divides the flow regimes into viscous, distorted and capped bubble regimes.In addition, the model also incorporates the bubble swarm correction as described below [36]: ( ) 0.75 , 24

Lift force
Lift force is the interfacial force acting on the lateral direction, perpendicular to the direction of the flow.This results to the migration of the bubbles in a shear field through the traverse lift force.The migration of the bubbles can be translated through the shear-induced lift force model, (12) as given by [37][38][39]: ( ) where, L C represents the lift coefficient.The lift coefficient can be integrated either as a non-dimensional constant or through a lift model as a function of the bubble Eötvös and bubble Reynolds number.As lift coefficient is affected by bubble size, this study opt to incorporate lift coefficient using lift model instead of assigning a fixed constant value as it tends to be influenced by bubble size and the proposed constant values varies across literature [40][41].Tomiyama et al. (2002) lift model is implemented to calculate the lift coefficient in this study as shown in (13).The principle behind this correlation was derived from a single bubble moving in a glycerol-water solution.This model has also shown success in air-water system and thus is being implemented in this study [42].
The modified Eötvös number is given by

Turbulent dispersion force
The turbulent dispersion force describes the effect of the turbulent fluctuations of liquid velocity on the bubbles.In this study, the Simonin and Viollet (1990) model on drift velocity is being implemented in both the sensitivity study in the drag models and consideration of the lift force.The correlation assumes the bubbles are to be caught up in turbulent eddies of the liquid phase and carried from high concentration regions to lower concentration regions as described in the following: , , where d v is the drift velocity, D is elaborated in [19] and will not be repeated here.

Measurement of mean axial liquid velocity from airlift reactor
The mean axial liquid velocity was measured using LDA in an internal-loop airlift reactor and the measurements were illustrated in Figure 3. 0.00 m < XR < 0.25 m represents the horizontal position of the riser meanwhile -0.15 m < XD < 0.00 m represents the horizontal position of the downcomer.

FluidsChE 2017
At height 0.20 m and 0.30 m the axial liquid velocity shows negative in the downcomer capturing the downward flow direction.Meanwhile, the riser at all heights showed an upward flow as shown through the positive values of the mean axial liquid velocity.This indicates there is a liquid circulation occurring in the airlift reactor which is vital for mixing.The axial velocity in the riser region shows greater fluctuation in the riser than those in the downcomer which is similar in trend with [14] study.Nevertheless, results obtained from the LDA is satisfactory within a ±0.01 m/s range of discrepancy between the repetitive results as recommended [15].

Comparison between different drag models
Prior to the CFD simulation, a grid sensitivity study was conducted by preparing a coarse grid (236k cells), intermediate grid (373k cells) and fine grid (717k cells).Figure 4 shows the results from the grid independence study, whereby the Universal drag model, Tomiyama lift model and drift velocity were used for the evaluation.At height Y1 = 0.20 m, all grids were in close agreement with the experimental data especially in the downcomer.At higher heights where the liquid velocity of more driven by the buoyancy of the bubble, both the intermediate and fine grid achieved good agreement with experimental data obtained from LDA.Meanwhile, large discrepancy occurs with experimental data using the coarse grid at higher heights.The computational time spent by each grids respectively are; coarse grid (0.12 s/iteration), intermediate grid (0.26 s/iteration) and fine grid (0.35 s/iteration), respectively.From the overall result, the intermediate grid (373k cells) yielded a grid independent solution as considering the fine grid (717k cells) will be computationally costly with no significant improvements to the result.Thus, the intermediate grid was considered sufficient enough to be used for the remainder of this work.
The effect of different drag models (i.e.

The effect of lift model
The effect of lift model on the predictive accuracy of the mean axial liquid velocity is studied with Universal drag and turbulent dispersion force enabled as illustrated in Figure 6.When the lift effect is neglected, it is observed that the discrepancy between the simulation results with experimental result is minimal.This indicates that drag force plays a dominant role amongst the interfacial momentum forces.From the figure, the radial gas profile peaks are much higher, when lift force is neglected.
When lift force is enabled, the radial gas profile peak subsides.This is also similarly reported in a two dimensional airlift reactor study where they observed flatter radial profiles on the mean axial liquid velocity in a bubble column at different superficial gas velocity [45].This is attributed to the consideration of lift force, drifting the bubbles tangentially from its flow direction, heading towards the wall.From the results it can be seen the role of the lift force more prominent with increasing height across the airlift reactor.Although the difference in peaks here are not that substantial, it is worth noting that the effect of the lift model is deem more significant at higher superficial gas velocity as the bubble sizes are polydisperse in nature [46].It is observable for all case studies that the riser at height Y2 = 0.30 m and both the riser and downcomer height at Y3 = 0.38 m, the prediction accuracy of the mean axial liquid velocity in the riser is far off although the trend is present.The combination of drag force, lift force and drift velocity leaves room for further improvement in the accuracy at those heights.Those are the position at which a swarm of bubbles are dispersed at.Hence, it is desired that in the future studies to include the bubble-induced model to improve the predictive accuracy.It is also interesting to observe the performance between the drag models, Tomiyama et al. (1999) and Universal drag model at different gas volume flow rate especially at higher gas flow rates.A study conducted in the bubble column have deem the Tomiyama et al. (1999) unsuitable at heterogeneous flow [47].The study here indicates Tomiyama et al. (1999) to be on par with Universal drag model.Nevertheless, the results of this study will be useful to understand the roles of each interfacial forces especially the role of lift force at different heights across the reactor.

Conclusion
A simulation case study for a closure model was carried out in a gas-liquid internal-loop airlift reactor using twofluid model approach.The study of the closure model consists of the evaluation of the effects between different drag models and the consideration of the lift model on the predictive accuracy.The simulation results were validated with experimental local measurements taken using LDA, obtained across three different heights (Y = 0.20 m, 0.30 m and 0.38 m) of the internal-loop airlift reactor.
Results suggest that Schiller and Naumann (1935) drag model is poor at predicting the mean axial liquid velocity, giving major discrepancy with increasing height.The Universal drag performed the best in comparison to the Tomiyama et al. (1999), with closer agreement with experimental data.The Universal drag model was then implemented to evaluate the significance of considering the lift force in the gas-liquid internalloop airlift simulation.The outcome reflects that the lift model subsides mean axial liquid velocity peaks, giving better improvements to the accuracy of the simulation results.The role of the lift force is seen more prominent with increasing height across the airlift reactor.This work reflects the importance of the closure model for better accuracy in the gas-liquid simulation results and LDA's capability at capturing the local hydrodynamics in an internal-loop airlift reactor.Future works in heterogeneous flow will allow further evaluation on the potential of this closure model.

Fig. 1 .
Fig. 1.Laboratory scale setup of an internal-loop airlift reactor, (a) Dimensions of the airlift reactor and (b) LDA measurements taken at three positions.

Figure 1 (
Figure1(a) shows a laboratory scale internal-loop airlift reactor made of transparent Perspex materials with measurements of 0.40 m (width) x 0.60 m (height) x 0.25 m (depth).An internal baffle with a height of 0.26 m is positioned 0.10 m from the bottom of the reactor.The top of the reactor was left open.Air is introduced into the airlift reactor through a 4" Classica air stone with the dimensions of 0.10 m in length and 0.20 m in height located at the base of the reactor.The reactor was filled with water at a height of 0.46 m and aerated with a volumetric flow rate of Qg = 0.30 m 3 /h using a vacuum pump.In order to obtain a non-coalescence system, 4 kg/m 3 of salt was added into the reactor[31].Measurements of the axial liquid velocity component is taken using two component FlowExplorer.LDA system from Dantec Dynamics whereby the system is operating in back-scattered mode.The gas-liquid system was seeded with polyamide seeding particles of diameter size 20 ȝm and density of 1030 kg/m 3 from Dantec Dynamics.The measurements were taken at heights 0.20 m, 0.30 m and 0.38 m of the column as shown in Figure1(b).These measurements were taken radially at an interval of 10 mm between each points with a 2.5 mm gap away from the wall.Both 15 mW He-Ne and 5 W Ar-ion lasers manufactured by Spectra-Physics were used as light sources emitting wavelengths of 532 nm and 561 nm.The laser intensity is set at 40 mW for both lasers.A 500 mm focal length lens with a probe volume of 2.1 mm 3 was used as a receiving lens in this study.The laser and focal lens were mounted on a bidirectional traverse controlled using a workstation.Data validation and signal processing were carried out using BSA3 Processor v6.10 software collected from the site using burst mode data collection method.It is worth noting that the data rate of acquisition drops as the laser approaches the bubble region.This is especially observed when the laser moves from the wall to the bubble plume region in the riser.When light hits the bubble, the intensity of the light reflected will be higher[32].However, these lights are reflected and scattered, lowering the probability of the lights hitting the detector.This results in the reduction in sample count obtained.Several measurements were performed at these locations repetitively by adjusting the gain, sensitivity, and anode to improve the accuracy of the time-averaged values.The duration of time specified to collect the samples is set at 60 s which is sufficient enough to obtain time-averaged value with an accuracy of ±0.01 m/s.Within that timeframe and setup as describe in section 2, 10 2 to 10 3 amount of samples were obtained.Although longer sampling would be encouraged, it was found large data quantities would not result to significant improvements in accuracy[26].

d
represents the Sauter mean bubble diameter, D C is the drag coefficient which can be retrieved from the drag model and the relative velocity, Simonin and Viollet (1990) model as shown in(17).N this study, the drag models Schiller and Naumann (1935),Tomiyama et al. (1999) and Universal drag were considered for comparison.Schiller and Naumann (1935) drag model is derived from a single spherical bubble moving in a stagnant liquid.It has been a favourable drag model applied widely in airlift reactor studies even in those which consist of polydisperse bubbles in their study[23][24][25]35].Hence, it is chosen as one of the drag models for evaluation: al. (1999) drag model is derived from a single bubble in an infinite stagnant liquid studied under a wide range of bubble sizes, gravity acceleration, different degrees of contamination and different fluid properties.The expression of the model is given by : is the major dimension of the bubble which can be calculated through an empirical correlation of the aspect ratio as given[43]:
Schiller and Naumann, 1935; Tomiyama et al. 1999 and Universal drag) is being studied by comparing the predicted mean axial liquid velocities with experimental measured data obtained from LDA as shown in Figure 5.The lift model and the turbulent dispersion were enabled for this case study.It can be seen Schiller and Naumann (1935), a drag coefficient as a function of Reynolds number performed poorly in the overall results.It only performs fairly at height, Y1 = 0.20 m which is the region most nearest to the sparger.Beyond that, the discrepancy became even larger with increasing height.It is worth noting that Schiller and Naumann (1935) model assumes spherical bubbles which neglects deformation of the bubble surface.On the other hand, both drag coefficients by Tomiyama et al. (1999) and Universal drag governed by the Eötvös number and different flow regimes, respectively were in closer agreement with experimental data obtained from LDA.The Tomiyama et al. (1999) included the bubble Reynolds number and an additional Eötvös number in its expression whereby it considers the terminal velocities that are affected by gravity and surface tensions due to bubble deformation.The Tomiyama et al. (1999) drag model is expected to be more accurate when paired up with Tomiyama et al. (2002) lift model as the lift correlation was developed while accounting the Tomiyama drag coefficient[44].However, the result especially in the downcomer in our study shows otherwise.

Fig. 4 .
Fig. 4. Results on grid sensitivity study on the mean axial liquid velocity at three different heights, (a) Y1 = 0.20 m, (b) Y2 = 0.30 m and (c) Y3 = 0.38 m.Amongst the drag models, the Universal drag is best fitted with the experimental results in both the riser and downcomer although it only gives a slight improvement in accuracy in comparison to the Tomiyama et al. (1999) model.This result suggests that the drag model as a function of Rayleigh-Taylor instability is more accurate.Hence, the Universal drag model is being implemented

Fig. 5 .
Fig. 5.The effect of different drag models on the mean axial liquid velocity compared with experimental data obtained from LDA at different heights, (a) Y1 = 0.20 m, (b) Y2 = 0.30 m and (c) Y3 = 0.38 m.

Fig. 6 .
Fig. 6.The effect of lift force on the mean axial liquid velocity compared with experimental data obtained from LDA at different heights, (a) Y1 = 0.20 m, (b) Y2 = 0.30 m and (c) Y3 = 0.38 m.