Department of Chemical Engineering,College of Chemistry and Chemical Engineering,State Key Laboratory of Metal Matrix Composites,Shanghai Jiao Tong University,Shanghai 200240,China

Keywords:Bubble column reactor Computational fluid dynamics Volume of fluid method Surface tension models Gas distributor

ABSTRACT This work aims at comparing surface tension models in VOF(Volume of Fluid)modeling and investigating the effects of gas distributor and gas velocity.Hydrodynamics of a continuous chain of bubbles inside a bubble column reactor was simulated.The grid independence study was first conducted and a grid size of 1.0 mm was adopted in order to minimize the computing time without compromising the accuracy of the results.The predictions were validated by comparing the experimental studies reported in the literature.It was found that all surface tension models can describe the bubble rise and bubble plume in a column with slight deviations.

1.Introduction

Bubble column reactors,with its extensive applications in the biochemical and petrochemical industries,have been studied for decades with theoretical and experimental endeavors as an interesting and significant subject[1-5].Besides these approaches,computational fluid dynamics(CFD)has received much attention as a powerful tool in modeling fluid hydrodynamics in bubble columns since the rapid improvement in computer science in recent years[6-12].Among different modeling techniques for bubble column reactor,the Euler-Euler[13,14]and the Euler-Lagrange [15-19]approaches are mostly used.The superiority of the Euler-Euler approach is that macro-scale average hydrodynamics properties could be obtained,which is of practical values.Nevertheless,constitutive equations with closure models are required in this approach and dynamics of the dispersed phase are ensemble averaged.The Euler-Lagrange model,with external forces accounted for individual particle,can be used to track the individual bubble with the motion equations.However,the bubble breakage and coalescence are difficult to be modeled by the Euler-Lagrange approach due to its limitation in bubble interface tracking[20,21].

Recently,the interface tracking methods (e.g.the front-tracking approach and the volume-tracking approach)are applied to solve this problem and some satisfactory results have been obtained,among which the volume of fluid (VOF)method is most commonly used due to its reliability and stability [18,19,22-28].Compared with the two fluid model (TFM)and the Euler-Lagrange method,VOF method has not yet received much attention and most of the available investigations usually deal with the single bubble behaviors [5,24,28-32].Notwithstanding the wide applications of the TFM and the Euler-Lagrange and other approaches [13,33-37],it is interesting to simulate bubbly flow with VOF approach,especially in modeling the behaviors of bubble swarms.During the previous studies where VOF approach had been used [18,19,25-28],the continuum surface force (CSF)model [38]was mostly applied to deal with the surface tension in the interface.To our best knowledge,there is hardly any literature focused on the comparison of surface tension models.In this study,the VOF method was used to simulate the hydrodynamics of a continuous chain of bubbles inside a typical column with a square cross-section.Different surface tension (e.g.the continuum surface force (CSF)model and the continuum surface stress (CSS)model)and adhesion models (e.g.wall adhesion model and jump adhesion model)were thoroughly compared.Subsequently,the influences of gas distributor and inlet gas velocity on the bubble rise velocity,bubble size distribution and bubble trajectory were studied.The predictions were validated with the experimental studies reported in the literature.

2.Model Developments

The VOF method was adopted to simulate the continuous chain of bubbles hydrodynamics inside a typical column with a square crosssection[18,19,22,25-28],in which the continuity equation of volume fraction can be expressed as follows:

a single momentum equation was solved throughout the domain,and the resulting velocity field was shared among the phases.The momentum equation,shown below,is dependent on the volume fractions of both phases through the properties:

where ρ is the density,μ is the viscosity,in the two-dimensional Cartesian coordinate system u is the velocity vector which is comprised of uxand uy,g is the gravitational acceleration vector which is comprised of gxequaling to zero and gyequaling to 9.81 m·s?2.Fbis bubble-liquid interface coupling interaction.In this study,the continuum surface force(CSF)model[38]can be written as following:

where σ is the surface tension coefficient.κ is the free surface curvature which has the following expression:

in the equation,the normal unit vector^n has the definition as following:

in modeling surface tension in a conservative manner,the continuum surface stress (CSS)model is an alternative method,which can be expressed as follows:

where I represents the unit tensor and σ represents the surface tension coefficient.In the VOF model,an alternative option to ascertain a jump and wall adhesion angle together with the surface tension model is also accessible.The surface normal next to the wall at the live cell or the neighboring cell to the porous jump can be expressed as follows:

For the wall adhesion model,is the unit vector to the wall.is the unit vector tangential to the wall.While for the jump adhesion model,is the unit vector normal to the porous jumpis the unit vector tangential to the porous jump.In the rest parts of this paper,V,CF,CFJ,CFJW,CS and CSW were used to represent different surface tension and adhesion models,whose meanings are listed in nomenclature.In this study,the mixture density and viscosity should obey the mixing rule as following:

3.Simulation Strategy

A two-dimensional rectangle bubble column reactor was adopted as the simulation domain with 0.15 m in length and 0.45 m in height that is sketched in Fig.1.Firstly,the rectangular reactor was meshed into structured grids that are of different amounts.In this step,GAMBIT 2.4.6(Ansys Inc.,US)was applied,which is a commercial software package.Subsequently,the equations formulated in Section 2 were calculated with the computational fluid dynamics code which is named FLUENT 15.0(Ansys Inc.,US).All the cases were simulated in an unsteady flow manner,that is to say,a transient technology was applied in this study.Liquid was the primary phase and bubble phase was set as the second phase,which is modeled with the multiphase VOF approach.In the multiphase VOF approach,implicit body force was applied.Outlet condition was set with outflow.For all stationary reactor walls,wall conditions were set as no-slip conditions.For spatial discretization of volume fraction,the geo-reconstruct scheme was adopted.In order to ensure accuracy,convergence criterion was set as 1×10?6for all the simulated cases.Under-relaxation values were the default for density,pressure as well as momentum in order to achieve good astringency.The detailed information of simulation settings is listed in Table 1.

Fig.1.The 2D rectangular simulation domain.

Table 1 Simulation parameters in this work

In modeling the bubble column,a total actual flow time of 5.0 s was adopted for each simulation case.The simulation tests were conducted on a platform of 2.6 GHz Intel CPU(8 cores)with 16 GB of RAM.

Fig.2.Checking work of grid size independence analysis.

Fig.3.Comparisons simulation results with experimental data.

4.Results and Discussion

4.1.Grid independence and numerical validation

The checking work of mesh size independence here was focused on a bubble column reactor.Three types of grid resolution including 0.50 mm,1.00 mm,and 2.00 mm were designed.The instantaneous velocity profiles were curved in Fig.2.As can be seen in Fig.2,the velocity curves are almost identical for the three types of mesh size in most region of bubble column.However,deviations can be observed in the central region.The predicted velocity magnitude in 0.50 mm grid size case is slightly higher than the cases of 1.00 mm and 2.00 mm.For grid resolution of 1.00 mm and 2.00 mm,the predicted velocity results are identical in the whole region of bubble column reactor,which means that the grid resolution of 1.00 mm and 2.00 mm are fine enough to obtain reasonable calculation accuracy.Therefore,the grid resolution of 1.00 mm was adopted for rest cases in order to minimize the computing time without compromising the accuracy of the results.The total number of grid elements is 67500.

The validation of the method was carried out by comparing simulation results with experimental data reported in Deen's work[33].As can be seen in Fig.3 where the instantaneous axial mixture velocities are curved,the simulation results agree well with the experimental data.In the center of the bubble column,the axial mixture velocity is slightly overestimated by the VOF model.While in the region between the center and column wall,the axial velocity is slightly underestimated.In the near wall region the velocity is overestimated again which is contributing to the wall effects.The experimental data curve is more smooth than the simulation curve.It is noteworthy that in the region with liquid as the main phase,the simulation and experimental data can be compared.However,in the region with both gas and liquid phases,the data comparison is unreasonable except when there's no velocity slip between the two phases.The data deviation in the central region may be due to the difference in the slip condition in the simulation and experimental cases.Generally speaking,the tendencies of axial velocity are the same for simulation as well as experimental results,which indicates that the VOF model can predict bubble behavior in a liquid system with adequate precision in most cases.

Fig.4.Comparisons of different surface tension models in terms of axial velocity magnitude.(A)H=0.25 m;(B)H=0.352 m.

Fig.5.Comparisons of different surface tension models in terms of bubble trajectory.

Fig.6.Velocity vector distributions under different hole size distributors.

4.2.Comparison of surface tension models

The comparisons of different surface tension models in terms of axial velocity magnitude are curved in Fig.4.The adopted grid resolution is 1.0 mm.Velocity profiles in two different altitudes (e.g.H=0.25 m and H=0.352 m)are presented.As can be seen in Fig.4,all surface tension models show the same tendency when simulating axial velocity.There exists the maximum axial velocity magnitude in the center of the bubble column.While in the near-wall region,the axial velocity reduces to the minimum.What's more,there also exists backmixing phenomenon in the near-wall region since the magnitude of axial velocity becomes negative.When compared the velocity magnitude of the two different altitudes,it can be seen that the maximum velocity in the magnitude of 0.25 m is slightly higher than the magnitude of 0.352 m,which indicate that the axial velocity is reduced as bubble column height rises.When compared with different surface tension models,it is found that the axial velocity curves simulated with combined continuum surface force model or continuum surface stress model are slightly lean to the right.That is to say,the bubble plume centers simulated with combined continuum surface force model or continuum surface stress model are slightly lean to the right.

Fig.7.Simulated velocity profiles of a different distributor and column height.(A),(B)VOF model;(C),(D)VOF model combined with the CSF model;(E),(F)VOF model combined with the CSS model.

A further confirmed comparisons of different surface tension models are contoured in Fig.5,where instantaneous bubble ascending trajectories are presented.As can be seen,all surface tension models can describe the bubble rise and bubble plume in a column with slight deviations.Bubbles rise straightly in the center of the column.At the same time,bubble breaks up and coalescent with each other since bubbles with different size have different rise velocity.The shape of the leading bubble is mostly like a cap while the trailing bubbles are like circles.The edge of the leading bubble would be sheared by liquid and break up with the bulk.When compared different surface tension models in terms of gas volume fraction contours,it is found that a bigger leading bubble would be obtained with the combined continuum surface force model and continuum surface stress model,which is attributed to the coalescence phenomenon of the leading bubble and following trailing bubbles.However,the daughter bubble plume is longer and the leading bubble is correspondingly smaller when the VOF model is adopted.

4.3.Effect of the gas distributor

Effect of the distributor on bubble plumes and hydrodynamics were studied at constant superficial gas velocity.Two distributors with different hole sizes were conducted in this study.The bigger distributor with the hole size of 0.03 m is recorded as distributor-1,while the smaller distributor with the hole size of 0.015 m is recorded as distributor-2.The simulated velocity vector distributions are contoured in Fig.6,where three different surface tension models are also compared.As can be seen in Fig.6 where instantaneous bubble ascending trajectories are presented,bubble plumes in the column with distributor-1 is zigzag.In the column with distributor-2,bubble plumne is more straight and hardly any vortex exists.It is attributed to the violent shearing effect in the gas inlet region of distributor-1,where bubbles with large size are easily broken.Consequently,the bubble trajectory is influenced and become wandering.In the column with a smaller distributor,the simulated result by VOF model is obviously different from the results simulated by the combined continuum surface force model and continuum surface stress model.As can be seen,the plume is more zigzag and the bubbles are more dispersed when compared with the combined CF and CS model.The simulated results are the same in the column with distributor-2,which confirms that more straight bubble plume could be obtained with the combined CF and CS model.

Fig.8.Velocity magnitude distributions under different gas inlet velocities.(A)gas inlet velocity is 5.0×10?3 m·s?1;(B)gas inlet velocity is 7.5×10?3 m·s?1.

The simulated instantaneous velocity profiles of a different distributor and column height are curved in Fig.7.When comparing the velocity profiles of the two different altitudes,it can be seen that the velocity distribution in the magnitude of 0.25 m is higher than the magnitude of 0.352 m,which indicates that the axial velocity is reduced as bubble column height rises.Generally,the velocity profiles in the column with distributor-1 are higher than in the column with distributor-2.The results are the same with the velocity vector distributions which contoured in Fig.6.When comparing the velocity profiles interms of different simulating methods,it can be found that the velocity profiles are overall small obtained by VOF model,especially in the column with distributor-2.It can be concluded that zigzag bubble plume and overall small velocity profiles are obtained with the VOF model than the combined CF model and CS model.

4.4.Effect of the gas velocity

Effect of the gas velocity on bubble plumes and hydrodynamics were studied in the column with a smaller distributor.Two different inlet gas velocities(e.g.5.0×10?3m·s?1and 7.5×10?3m·s?1)were tested in this study.The velocity magnitude distribution is contoured in Fig.8,where three different surface tension models are also compared.As can be seen in Fig.8,the effects of inlet gas velocity are not obvious if the combined CF model and CS model are adopted.In the case where the VOF model was applied,the influences of inlet gas velocity on velocity distribution are obvious.It can be found that the bubble plume in low inlet gas velocity is more straight than the bubble plume in high inlet gas velocity.It is because that the shearing effect of liquid is violent in the high inlet gas velocity case.When comparing the simulated results obtained by the VOF model and the combined CF and CS models,it can be concluded that the simulated velocity distribution obtained by VOF model is different from combined CF and CS models.The velocity are overall smaller if the VOF model is adopted.This conclusion is the same with the velocity vector distributions and velocity profiles which contoured in Figs.6 and 7,respectively.

The simulated instantaneous velocity profiles of a different distributor and column height are curved in Fig.9.Generally speaking,the effects of inlet gas velocity are not obvious in terms of velocity profiles.When comparing the velocity profiles of the two different altitudes,it can be seen that the velocity distribution in the magnitude of 0.25 m is higher than the magnitude of 0.352 m,which indicates that the axial velocity is reduced as bubble column height rises.When comparing the simulated results obtained by the VOF model and the combined CF and CS models,it can be concluded that the velocity is overall smaller if the VOF model is adopted.This conclusion is the same with the velocity vector distributions and velocity profiles which contoured in previous figures.

5.Concluding Remarks

Fig.9.Simulated velocity profiles of different gas inlet velocity and column height.(A),(B)VOF model;(C),(D)VOF model combined with the CSF model;(E),(F)VOF model combined with the CSS model.

In this study,the VOF method was used to simulate the hydrodynamics of a continuous chain of bubbles inside a typical column with a square cross-section.The checking work of grid independence was first conducted.The grid resolution of 1.0 mm was adopted in order to minimize the computing time without compromising the accuracy of the results.Then the predictions were validated with the experimental studies reported in the literature.After that,different surface tension models (e.g.the continuum surface force model and the continuum surface stress model)and adhesion models (e.g.the wall adhesion model and the jump adhesion model)were thoroughly compared.All surface tension models show the same tendency when simulating axial velocity.Subsequently,the influences of the gas distributor and the inlet gas velocity on the bubble rise velocity,bubble size distribution and bubble trajectory were studied.It can be concluded that bubble plumes in the column with distributor-1 is zigzag.While in the column with distributor-2,bubble plume is straighter and hardly any vortex exists.What's more,it can be found that zigzag bubble plume and overall small velocity profiles are obtained with the VOF model than the combined CF model and CS model.Finally the effects of inlet gas velocity were tested.It can be found that the effects of inlet gas velocity are not obvious if the combined CF model and CS model are adopted.In the case where the VOF model was applied,the influences of the inlet gas velocity on velocity distribution are obvious.As for the applicable scope of the model,it should be based on pressure solver.Besides,the calculation of surface tension effects on triangular and tetrahedral meshes is not as accurate as on a quadrilateral and hexahedral meshes.The region where surface tension effects are most important should therefore be meshed with quadrilaterals or hexahedra.

Acknowledgements

The authors thank the Center for High Performance Computing,Shanghai Jiao Tong University for supporting this work.