CFD-PBM Coupled Simulation of an Airlift Reactor with Non-Newtonian Fluid

Hydrodynamics of an AirLift Reactor (ALR) with tap water and non-Newtonian fluid was studied experimentally and by numerical simulations. The Population Balance Model (PBM) with multiple breakup and coalescence mechanisms was used to describe bubble size characteristics in the ALR. The interphase forces for closing the two-fluid model were formulated by considering the effect of Bubble Size Distribution (BSD). The BSD in the ALR obtained from the coupled Computational Fluid Dynamics (CFD)-PBM model was validated against results from digital imaging measurements. The simulated velocity fields of both the gas and liquid phases were compared to measured fields obtained with Particle Image Velocimetry (PIV). The simulated results show different velocity field profile features at the top of the ALR between tap water and non-Newtonian fluid, which are in agreement with experiments. In addition, good agreement between simulations and experiments was obtained in terms of overall gas holdup and bubble Sauter mean diameter.

Turbulent dissipation rate due to bubbleinduced turbulence, m 2 s À3 Subscripts m Gas-liquid mixture l,L Liquid phase g Gas phase i,j Classes of bubble size

INTRODUCTION
Airlift Reactors (ALR) attract more and more interests in many process applications, such as syngas fermentation [1], direct coal liquefaction [2] and waste water/gas treatment [3].The advantages of ALR over bubble columns include enhanced mixing efficiency and good suspension of particles due to improved liquid circulation.In comparison to mechanically stirred tanks, simple structure and low energy consumption make ALR even more attractive, particularly considering the current stringent energy and safety demands.Airlift loop reactors are an important type of pneumatically agitated multiphase reactors where the bubble phase functions as the agitation source and reactant.Bubble size characteristics, determined by the initial bubble size, bubble breakup and coalescence, is an important hydrodynamic parameter and closely related to mass transfer performance and even production yield in the reactor.Bubble size influences gas holdup, gas residence time, specific interfacial area, and mass and momentum transfer between the gas and liquid phase.Meanwhile, bubble size and its behavior are affected by several factors, such as gas holdup, liquid velocity field and liquid physical property.Extensive studies [4][5][6][7][8] have shown that bubble size characteristics have significant effects on the hydrodynamics and mass transfer performance of ALR, especially for liquids with high viscosity.For example, Hwang and Cheng [5] found that gas holdup decreased with increasing Carboxyl Methyl Cellulose (CMC) concentration due to formation of large bubbles in highly viscous liquid.Deng et al. [6] showed that Bubble Size Distribution (BSD) for highly viscous liquids became broader in heterogeneous flow regime and increase of gas holdup slowed down with aeration rate due to formation of large bubbles.Gumery et al. [8] reported that liquid circulation velocity decreased in the ALR by increasing liquid apparent viscosity due to enhanced bubble coalescence.
Computational Fluid Dynamics (CFD) simulation is a powerful tool to investigate and predict hydrodynamics in multiphase reactors.However, it is required that the CFD models describe the phenomena behind the physics accurately in order to predict results with good confidence.This is still difficult and challenging since fluid mechanics in multiphase flow is quite complex and far from fully understood.Therefore, it is necessary to validate the CFD simulations of ALR experimentally, particularly with different liquid physical properties.The two-fluid model with the assumption of a single bubble size can give reasonable predictions for multiphase reactor only in cases where the bubble size is narrowly distributed.Coupling of Population Balance Model (PBM) into CFD models could make it possible to extend the simulation and prediction of flow regimes from homogeneous to heterogeneous where the bubble size is widely distributed.PBM is an effective method to calculate BSD for gas-liquid flows on the condition that bubble coalescence and breakup models are provided as model closures.Many researchers have used the coupled CFD-PBM model to predict the hydrodynamics and mass transfer in multiphase reactors [9][10][11][12][13][14]. Their studies showed that a better agreement between simulations and experiments was obtained when BSD was considered compared to using a constant mean bubble diameter.Chen et al. [10] found that the coupled CFD-PBM simulation was able to predict the effects of aeration rate and liquid surface tension on the gas holdup in bubble columns.They pointed out that simulations made by implementing PBM produced much better comparison to experiments for gas holdup and liquid velocity than by using mean bubble diameter.Wang and Wang [12] investigated mass transfer in bubble column with the coupled CFD-PBM model and showed that the CFD-PBM coupled model could predict effectively the hydrodynamics and mass transfer.Liu et al. [13] performed CFD-PBM coupled simulations of airlift bioreactor for hairy root culture and presented good agreement between simulations and experiments in terms of bubble number density.Xing et al. [14] investigated the ability of the CFD-PBM coupled model to account for the effect of viscosity on the gas holdup in a bubble column.The feature that the increased viscosity resulted to a decrease in the volume fraction of small bubbles and an increase in the volume fraction of large bubbles was captured well using the developed CFD-PBM model [12,15].The coupled CFD-PBM simulations of ALR are relatively rare, especially with the viscous or non-Newtonian fluid.In industrial applications, many processes [16][17][18] are operated in the viscous mediums or with catalysts showing non-Newtonian behavior.
The aim of this study was to develop a CFD model for ALR with non-Newtonian fluid in order to contribute to the design and optimization of ALR in practical applications.A PBM with the coalescence mechanisms by turbulent eddies, wake entrainment and different rise velocities of bubbles is coupled into the two-fluid model.Bubble breakups due to eddy collision and eddy instability are taken into account.Gas holdup, bubble size and liquid and gas velocity fields obtained by Particle Image Velocimetry (PIV) were used to validate the CFD model.

EXPERIMENT 1.Experimental Setup and Material
The schematic diagram of the experimental setup used in this work is shown in Figure 1.
The AR_ALR is composed of a riser (outer column), downcomer (draft tube) and gas distributor.The riser is 0.15 m in inner diameter and 1.2 m in height.The downcomer is 0.60 m high with the inner diameter of 0.10 m.The downcomer was co-axially mounted with the riser locating at 0.05 m above the column bottom.A plastic O-ring gas distributor with a diameter of 0.12 m was installed in the annulus between the riser and the downcomer and it was at 0.07 m above the bottom.Gas is fed to the column through 36 pores, with hole size of 0.5 mm distributed equidistantly at the top of the O-ring.The AR_ALR is located inside a transparent rectangular box in order to diminish optical distortions during visualization measurements.
Tap water and 0.14 wt% CMC (Finnfix 50000) aqueous solution were employed to study the Newtonian and non-Newtonian fluids.The dynamic viscosity of the CMC solution was measured with the Modular Compact Rheometer 302 (Anton Paar) and the curve is shown in Figure 2. As the reference, the viscosity of tap water was also tested and is shown in Figure 2. The CMC solution exhibited shear   thinning behavior and the apparent viscosity, l app , could be described by the following power law model: where j and n are the consistency index and the flow index of the fluid, respectively.j is 0.118 Pa s n and n is 0.68 for the CMC solution obtained by fitting Equation ( 1) with the leastsquares method.Tap water and CMC solutions were operated in experiments in batch mode.The static level for both liquid phases was 0.8 m.Oil-free compressed air as the gas phase was run in once-through mode.Gas flow rate was controlled by a mass flow controller (Bronkhorst E-7000).Gas superficial velocity, U g , was calculated based on the cross-sectional area of the riser (the annulus).All experiments were carried out at atmospheric pressure and room temperature.

Measurement Techniques
The overall gas holdup in the reactor, ɛ g,o , was determined by the height expansion method: where H l is the static liquid level and H m is the gas-liquid mixture height after aeration.
Turbulent dispersion force [24] F TD ¼ ÀC TD a g q l k l oa or Population balance model for bubbles in the ALR [11] PIV method was adopted to measure flow velocity fields for both gas and liquid phases simultaneously.The PIV system (LaVision GmbH) mainly included a double pulsed Nd Yag laser, two high resolution CCD cameras (1600 pixel 9 1200 pixel) and DaVis 8.3 software.The liquid velocity field was traced by PMMA-RhB fluorescent particles with the diameter range of 20-50 lm.The measuring range of the flow field was from the wall to the centerline of the ALR in the radial direction.Each measurement window is 75 9 50 mm in physical size.Eight hundred (800) pairs of images were acquired by each camera within 10 min for each test when the flow field was in steady state.DaVis 8.3 software was utilized to analyze the acquired images for each phase.The spatial resolution of the measured velocity field was 3 9 3 mm (grid size).Ensemble-and time-averaged velocity fields were calculated and used to compare to the simulations.
Digital imaging method was employed to evaluate the bubble size and size distribution in the ALR.One thousand (1000) digital images were acquired and used to analyze the bubble size for each condition.An overlapping object recognition algorithm (Pixact Object Recognition and Analysis software, Pixact Ltd) developed by Eloranta et al. [19] was used to recognize and analyze the bubbles from the acquired digital images.More than 10 4 bubbles were detected for each condition and used to calculate the bubble Sauter mean diameter, d 32 : where d bi and n i are the bubble equivalent diameter and the number of bubbles corresponding to the equivalent diameter, and i indicates the size class that refers to the bubble equivalent diameter.

Coupled CFD-PBM Model
The governing equations of the time-averaged two-fluid model and the discrete PBM model from Wang et al. [11] are listed in Table 1.BSD was resolved from the PBM using the gas holdup and the kinetic energy dissipation rate calculated in CFD.Then BSD obtained from the PBM was used to formulate the interphase forces and the turbulence modification closure terms for the two-fluid model.A modified k-ɛ model based on the two-time constant model from Lahey et al. [20] was used to describe the liquid phase Coalescence rate due to turbulent eddies Coalescence rate due to bubble wake turbulence, where both the shear-induced turbulence and the bubble-induced turbulence were considered.The turbulent viscosity of the gas phase was formulated using the correlation proposed by Jakobsen [21].The drag, lift, wall lubrication and turbulent dispersion forces were taken into account for gas and liquid momentum transfer.The drag coefficient model from Tomiyama et al. [22] was used to calculate drag force for each bubble size group.The Tomiyama drag coefficient model, being developed based on the balance of forces acting on a bubble and available theoretical and empirical correlations of terminal rising velocity, is well suited for gas-liquid flows where the bubbles can be of different size to a certain extent.The lift force was calculated with the lift force coefficient model of Tomiyama et al. [23], which is applicable to larger-scale deformable bubbles in the ellipsoidal and spherical cap regimes.The Tomiyama wall lubrication force model [24], which modified the formulation of Antal et al. [25], was employed for wall lubrication force calculation, as listed in Table 1.
The coupling of PBM into CFD model was implemented in ANSYS Fluent 14.5 using user defined scalars that were defined for volume fractions of bubble group i in the gas holdup.Bubble breakup and coalescence rate models were loaded by the user defined function.All bubble groups were assumed to move at identical velocity that equals to the local ensemble averaged gas velocity.
To describe BSD characterization in viscous fluid, the multiple bubble breakup and coalescence models from Xing et al.'s [14] were used to close the PBM.The total bubble breakup rate was calculated as the sum of the rates from eddy collision and from the instability of large bubbles.For the total coalescence rate three mechanisms were included, namely coalescence due to turbulent eddies, bubble wake entrainment, and different bubble rise velocities.The details of the breakup and coalescence models from Xing et al. [14] used are listed in Table 2.

Numerical Details
Steady state simulations were performed for the ALR.Two liquid phases, which are the CMC solution as the non-Newtonian fluid and tap water as the reference, were simulated.The apparent viscosity of the CMC solution was described with the power law equation, as measured in the experiments.The numerical simulations were performed at two gas superficial velocities of 1.18 cm s À1 and 2.04 cm s À1 , which corresponded to the homogeneous and heterogeneous flow regimes, respectively.The bubble size was divided into 30 classes with the geometric method (v i+1 = v i q), where the smallest bubble volume v 1 was 1.0 9 10 À10 m 3 with the factor q of 1.7.
The two-dimensional axis-symmetry assumption was used in this work due to the geometric symmetry of the ALR and low computational costs.The computational domain is shown in Figure 3a.The computational grid was created with ANSYS ICEM and the structural hexahedra gridding at the bottom and top is shown in Figure 3b and 3c, respectively.The grid density and the corresponding cell size in the  three modeled zones along the ALR height were 23 9 3 mm (below the gas distributor), 21 9 1 mm (around the gas distributor), and 284 9 2.5 mm (above the gas distributor).
In the radial direction, the grid was fined near the wall with the minimum edge size of 0.5 mm and the maximum edge size of 2 mm.Grid independent solution was verified by comparing to the experiments in terms of gas holdup and liquid velocity.The wall function was used for the liquid phase and the non-slip condition was used for the gas phase.The velocity inlet boundary was applied to the upper face of the gas distributor, where the inlet velocity, gas holdup and volume fractions of bubble groups were specified according to the experiments.At the top surface of disengagement zone, pressure outlet boundary condition was used.The second order upwind scheme was used for the momentum, turbulence kinetic energy and dissipation rate equations and the first-order scheme for the volume fraction equations.

Overall Gas Holdup
The overall gas holdups for tap water and CMC solution at two different gas superficial velocities U g , are shown in Figure 4.The overall gas holdups increase with increasing gas superficial velocity U g for both in tap water and in CMC solution.The value of the overall gas holdup is lower in CMC solution than in tap water at the same U g .The simulated results for the overall gas holdup agrees well with the experimental data for both the tap water and CMC solution with the relative error being less than ±10%.

Bubble Size
The simulation results of the BSD in the ALR with different liquid phases were obtained from the discrete PBM.The comparison between the simulation results and experimentally determined bubble number densities in the ALR with different liquid phases is shown in Figures 5a-5d.It can be seen that the bubble sizes are narrowly distributed both in water and in CMC solution at low U g .With increasing U g , BSD becomes broad in water and even bimodal in CMC solution.All these features are captured by the PBM with the bubble breakup and coalescence models including several mechanisms.The relative errors between the experiments and the simulations is less than ±10% for most of the bubble groups, as shown in Figure 5.
Figure 6 shows a good agreement between the simulations and experiments for bubble Sauter mean diameter in the ALR with different fluids at two investigated U g .It can be found that the bubble Sauter mean diameter is slightly under-predicted in water and over-predicted in CMC solution.These results are obtained most probably because the breakup models in the PBM are much more sensitive to the viscosity compared to the coalescence models, as shown by Xing et al. [14].In addition, the breakup model due to turbulent eddies by Wang et al. [15] is able to give much better prediction for unequal daughter BSD in breakup.

Velocity Profiles
The velocity fields of both liquid and gas phase in the ALR were measured with PIV method to compare the velocities against simulations.Figures 7 and 8 show the comparison of the velocity fields for water and CMC solution at U g 2.04 cm s À1 .It can be seen that the simulated velocity profile generally matched well to that measured with PIV.In addition, the magnitudes of the velocities are in agreement between the simulations and experiments.
By comparing Figures 7a and 8a, it can be found that the measured liquid flow structures at the top of the ALR for water and CMC solution are different from each other.These structures are also captured by the coupled CFD-PBM model, as shown in Figures 7c and 8c.This illustrates that the CFD-PBM model is able to describe the hydrodynamics in the ALR with reasonable accuracy for non-Newtonian fluid at the investigated U g .
Figures 9 and 10 show the radial profile of the axial liquid velocity at three axial heights in the ALR located at 0.15, 0.45, and 0.7 m from the bottom of the ALR.It can be seen that the radial profile of liquid velocity is more uniform in water compared to CMC solution, particularly in the downcomer.At the top of the ALR, the liquid velocity profile in CMC solution is different from that of water at high U g , most probably due to the effect caused by movement of large bubbles and bubble wakes.

CONCLUSION
The hydrodynamic characteristics of fluid flow was studied by experiments and numerical simulations in an ALR for tap water and non-Newtonian fluid.Coupled CFD-PBM simulations were carried out taking into account multiple breakup and coalescence models.Simulated overall gas holdup, BSD, Sauter mean bubble diameter, and velocity fields for both phases were compared and validated against experimental data.PIV and digital imaging techniques were employed in experiments in order to validate the simulations.
It can be concluded that the coupled CFD-PBM model can reasonably well predict the effect of apparent viscosity on the overall gas holdup, bubble Sauter mean diameter and bubble number density in the range of the investigated U g .The overall gas holdup and bubble Sauter mean diameter are smaller in CMC solution than in water.The BSD becomes broad in water and even bimodal in CMC solution at U g of 2.04 cm s À1 .All these features were captured well by the coupled CFD-PBM model.A difference in the liquid velocity profile was observed between the CMC solution and tap water at the top of ALR with the PIV method.The same difference was also predicted well by the coupled CFD-PBM simulations.

Figure 1
Figure 1 Schematic diagram of the experimental setup.

Figure 2 Dynamic
Figure 2 Dynamic viscosity of water and CMC solution measured by rheometer.

Figure 3
Figure 3 Computational domain and grid profile: a) computational domain, b) grid profile at the top, c) grid profile at the bottom.

Figure 4 Overall
Figure 4Overall gas holdup from experiments and simulations.

Figure 5
Figure 5 Bubble volume fraction of experiment and simulation.

Figure 6 Bubble
Figure 6Bubble Sauter mean diameter from experiments and simulations.

Figure 9 Figure 10 Radial
Figure 9 Radial profile of liquid axial velocity of experiments and simulations at U g = 1.18 cm s À1 .
À1 -T (d i , d j ) Collision rate due to turbulent eddies, s À1 P T (d i , d j ) Coalescence efficiency due to turbulent eddies -W (d i , d j ) Collision rate due to bubble wake, s À1 P W (d i , d j ) Coalescence efficiency due to bubble wake -U (d i , d j ) Collision rate due to different rise velocities, s À1 P U (d i , d j ) Coalescence efficiency due to different rise velocities

TABLE 1
Governing equations and discrete PBM for the ALR.

TABLE 2
Models of bubbles breakup and coalescence rate.