Jingcai Cheng,Qian Li,2,Chao Yang,2,*,Yongqiang Zhang,Zaisha Mao

1CAS Key Laboratory of Green Process and Engineering,Institute of Process Engineering,Chinese Academy of Sciences,Beijing 100190,China

2University of Chinese Academy of Sciences,Beijing 100049,China

3SINOPEC Research Institute of Petroleum Processing,Beijing 100083,China

Keywords:Population balance equation(PBE)Multiphase Bubble column OpenFOAM Computational fluid dynamics(CFD)

A B S T R A C T A general CFD-PBE(computational fluid dynamics-population balance equation)solver for gas–liquid poly-dispersed flows of both low and high gas volume fractions is developed in OpenFOAM(open-source field operation and manipulation)in this work.Implementation of this solver in OpenFOAM is illustrated in detail.The PBE is solved with the cell average technique.The coupling between pressure and velocity is dealt with the transient PIMPLE algorithm,which is a merged PISO-SIMPLE(pressure implicit split operator-semi-implicit method for pressure-linked equations)algorithm.Results show generally good agreement with the published experimental data,whereas the modeling precision could be improved further with more sophisticated closure models for interfacial forces,the models for the bubble-induced turbulence and those for bubble coalescence and breakage.The results also indicate that the PBE could be solved out the PIMPLE loop to save much computation time while still preserving the time information on variables.This is important for CFD-PBE modeling of many actual gas–liquid problems,which are commonly high-turbulent flows with intrinsic transient and 3D characteristics.?2017 The Chemical Industry and Engineering Society of China,and Chemical Industry Press.All rights reserved.

1.Introduction

Gas–liquid flows are often encountered in the cultivation of bacteria and mold fungi,production of cell proteins,treatment of sewage,liquid phase methanol synthesis,Fischer-Tropsch synthesis,etc.[1].Bubbles may break and coalesce due to bubble-bubble and bubble-liquid interactions,which greatly complicate these processes.Thus,in the last decades,population balance equation(PBE)is being increasingly utilized with computational fluid dynamics(CFD)to account for the effects of bubble-liquid and bubble-bubble interactions on the hydrodynamic and transport behavior in various gas–liquid configurations[2,3].

For low-turbulent bubbly flows,works on simulations of these processes are well reviewed by some researchers,for example,Joshi[2],Sokolichin et al.[4],Wang et al.[5],and Wang and Yao[6].Appropriate closure models are required for the interfacial/interphase forces,in order to close Reynolds-averaged,Favre-averaged or ensemble averaged two-phase momentum equations.Besides,closures for the interphase forces,models for the bubble-induced turbulence should also be taken into account.However,simulation results are considerably dependent on adopted models and combinations of these models.Therefore,modeling interfacial forces and bubble-induced turbulence remain an open question and a subject of considerable research activity in bubbly flow areas[7,8].For high-turbulent gas–liquid flows such as the slug or churn flow,theoretical,experimental and simulative works on these processes are well reviewed recently by Montoya et al.[1].It was concluded that there was a necessity for important further development in the closure laws in the highly-turbulent flow area.The authors also concluded that the models for interfacial forces,which have been developed and validated for predicting bubbly flows,were almost useless in most cases for slug or churn flows.Thus,besides theoretical and experimental works,further considerable modeling and simulative research on gas–liquid flows both in the bubbly and high-turbulent flow regimes are still needed.

OpenFOAM(open-source field operation and manipulation)is a free source CFD package written in C++,which uses classes and templates to manipulate and operate scalar,vectorial and tensorial fields[9–11].Combined with implementations of adequate numerical methods to the discretization of partial differential equations and to the solution of the resulting linear systems,OpenFOAM is deemed as a good choice to handle CFD problems[12].Its open-source characteristics facilitate the implementation of any addition or modification in the source code,which is very suitable for the purposes of research.And for industries,OpenFOAM could be a good alternative due to its cost-saving for parallel computation and its potential capability of dealing with specific problems.As regards the CFD packages for simulating the gas–liquid flows,most are multi-purpose commercial packages.However,the validation examples in OpenFOAM are far less than those in commercial packages.Rzehak and Kriebitzsch[8]made a comparison between OpenFOAM and the commercial code ANSYS-CFX in a mono-dispersed bubbly pipe flow.The same physical models were used in the computations with the different codes.The results show that the OpenFOAM is comparable to the commercial code,whereas different coefficients in some models should be used.The difference in the results from the two codes is somewhat considerable near the wall,which may be due to different wall treatment methods.

Very recently,some researchers have conducted the coupled CFDPBE simulations of bubble columns in OpenFOAM[7,13–19].Different interphase transfer closure models(especially the drag model)and different combinations of these closure models were adopted in these works to obtain reasonable results for the systems under consideration.For obtaining results in better agreement with the experimental data,different drag coefficients for different superficial gas velocities in the same system was also used[15].Wall lubrication and turbulent dispersion forces were ignored except that in Pena-Monferrer et al.[16].The PBE of coalescence and breakage is solved using the quadrature method of moments(QMOM)[13,14,16,17](which was first proposed by McGraw[20]and developed by Marchisio et al.[21]),the fixed pivot technique(a type of the discretization method proposed by Kumar and Ramkrishna[22])[7,13,15]or the average bubble number density(ABND)method[18,19].The coupled CFD-PBE models in these works are limited to low gas volume fractions.As the volume fraction for high-turbulent flow could be very high.A general CFD-PBE model,which is applicable to both low and high volume fractions,is desirable for future further research on gas–liquid flows.

The main contribution of this work is that a general CFD-PBE solver for gas–liquid poly-dispersed flows is developed in OpenFOAM.By using the modified pressure and the mixture turbulence model as well as taking the forms for momentum transfer terms proportional to the product of both volume fractions,this solver could be used to systems of very high gas volume fractions or even stratified flows.Other new features of this work lies in the following two aspects.(1)The transient PIMPLE algorithm[23],which is a merged PISO-SIMPLE(pressure implicit split operator[24]-semi-implicit method for pressure-linked equations[25])algorithm,is preliminarily investigated for the first.(2)The PBE is discretized and solved by the cell average technique[26],which is an improved version of the widely used fixed pivot technique[22].

2.Eulerian–Eulerian Multiphase Model

2.1.Continuity and momentum equations

The conditionally averaged continuity and momentum equations for the phase φ in a multiphase system are given by[27]

where αφis the volume fraction,Γφstands for the mass transfer among phases resulting from reaction or phase change,Mφis the momentum transfer between phase φ and other phases and satisfying∑Mφ=0 andrepresents the interface-averaged stress contribution expressed as

where index i refers to the interface.Tφis the total stress tensor usually decomposed into an isotropic pressure term and an effective deformation stress,i.e.,

where I is the identity tensor,pφis the static pressure of phase φ,andis usually modeled as

It is assumed that the average phase stress and its interfacial average equals,τφeff= τφi,and the same static pressure acts in all phases and equals to the interfacial average,pφi=pφ=p[28].Since phases may have very different densities,problems could arise from the fact that there is only a single pressure field.For example,? p has got to be different for each phase to satisfy the momentum equations at a non-vertical wall and it will suffer change across a stratified flow[11].This problem could be solved by adopting a modified pressure defined as

where x is the position vector and the mixture density is ρm= α1ρ1+α2ρ2+….

The momentum equation becomes singular for that phase when written as in the form of Eq.(2),when αφ→ 0.The phase-intensive version of the momentum equation,which is obtained dividing Eq.(2)by αφ,is introduced by Oliveira and Issa[28]to circumvent this problem.Thus,for incompressible flows without mass transfer among phases and after some mathematical manipulation,the phase-intensive momentum equation is expressed as

where Ur=Ua-Ub,CD,CLand CVMrepresent drag force coefficient,lift force coefficient,virtual mass force coefficient,respectively,KTDand KWare the functions in the corresponding turbulent dispersion and wall lubrication models adopted,and nWis the unit vector on the wall face.

In OpenFOAM,the Gauss theorem is employed to transform the volume integrals in a control volume into surface integrals of normal fluxes.Rusche[11]discussed the discretization procedure of phase-intensive momentum equations.The Reynolds stress terms are decomposed into a diffusive component and a correction:

where ?φ=S ?(Uφ)fis the volumetric phase flux,and notation ?[φ]」means the implicit discretization for the variable, and the terms without notation means they are handled explicitly.

Careful discretization of the pressure gradient prevents pressure and velocity decoupling and the occurrence of oscillations(checkerboarding)in the pressure field.And other terms such as the buoyancy term and the interphase momentum transfer term also need special treatment.The discretized phase-intensive momentum equations are expressed as

where Aφ(φ =a,b)represents the sparse matrix with coefficients resulting from implicitly handled terms and it can be decomposed into two matrices,(Aa)Dand(Aa)N,representing the diagonal and off-diagonal coefficients,respectively.(Aa)H=(Aa)S-(Aa)NUa′with(Aa)Sdenotes the matrix coefficients mainly from the terms handled explicitly and Ua′denotes the approximated solution of the linear system such as the velocity field of the former iteration or time step.Symbols Uφand Uφ′here do not denote true vectors but a list of values defined at the centers of the control volumes.

2.2.Interphase transfer models

Interfacial forces include the drag force,lift force,virtual mass force,wall lubrication force and turbulence dispersion force.There are many different models for each of these forces or force coefficients as well reviewed by Joshi[2],Sokolichin et al.[4]and Montoya et al.[1].The applicability of the interphase force models was dependent on the bubble regimes and flow conditions in gas–liquid flows.

A commonly used expression for the drag coefficient is due to Schiller&Nauman[29],which is simple in form but neglects the bubble deformity and is limited to small spherical bubbles.Some models consider the bubble deformability based on the E?tv?s number,for example the drag model by Tomiyama et al.[30].Ishii and Zuber[31]even formulated the coefficient for different bubble shape regimes.These models can get reasonable results for higher Reynolds number flows,where bubbles have evolved to the ellipsoidal regime[6].As the void fraction increases,there is a high possibility that bubbles exist in swarms,which can influence the drag forces significantly. It is highly challenging to measure the drag coefficient of a bubble swarm experimentally,which increases the difficulty to differentiate among the quality of the large number of swarm corrections in the literature[1].In this work,the popular expression for the drag coefficient proposed by Tomiyama et al.[30]is adopted:

where the bubble Reynolds number and the E?tv?s number are defined by.

and σ is the surface tension.The Sauter mean bubble diameter d32is computed by

Although the influence of the virtual mass force is generally accepted because of its physical soundness, it is rather difficult to estimate the coefficient correctly[4].The overwhelming majority of authors took the coefficient to be a theoretical 0.5.Other authors[16]believed that the influence of the virtual mass force in steady simulations can be neglected.The value of 0.5 for CVMis taken in this work.

The lift force significantly affects the bubble distribution pattern across the pipe section in simulation,whereas its existence and sign are still not validated for realistic bubble sizes [4]. Both positive and negative values for the lift coefficient have been used,for example,0.01 in Lahey[32],0.5 in Auton[33],Jakobsen et al.[34]and Gupta and Roy[35],-0.5 in Boisson and Malin[36]and Grienberger and Hofmann[37],-2.0 to-1.5 in Jakobsen et al.[38],and even-3.0 in Svendsen et al.[39].Tomiyama et al.[40]thought the coefficient is dependent on the bubble size and the flow conditions,where CLtakes positive values for small bubbles and negative values for large bubbles.Thus,large bubbles are driven towards the center of the pipe and small bubbles are towards the wall.Frank et al.[41]modified this expression to be consistent at the transition points:

The wall lubrication force pushes bubbles near a wall away from the wall.It is well known that the wall and lift forces are responsible for the determination of the gas fraction in the region close to the wall[1].However,different values of wall force coefficient were used by many authors,and there is no widely applicable choice.The expression of Tomiyama[42]for the wall lubrication function KWis used here:

where CWis the lubrication coefficient,and its value is about 0.1 for air water system.

The turbulence dispersion force accounts for the effect of liquid turbulent fluctuation on the bubbles [6]. It is derived by the ensemble average of the fluctuating component of the drag forces between bubbles and the liquid phase.In the Reynolds-like averaging procedure,the turbulent dispersion force was neglected because these mechanisms occur as a phasic diffusion term in the continuity equations[43].The model of Lopez de Bertodano[44]was widely used by many authors,but the coefficient was taken differently between 0.09 and 1.0,for example,0.1 in Pe?a-Monferrer et al.[16],0.5 in Ekambara et al.[45],0.7 in Wang et al.[5]and Wang and Wang[46],1.0 in Wang and Yao[6].Burns et al.[47]derived a formulation based on the Favre averaging of the drag term,and obtained acceptance results:

2.3.Pressure equation

The solution of the pressure equation produces the updated pressure field,which in turn is used to correct the velocities and their fluxes so that the continuity is guaranteed.Rusche[11]combined the continuity equations of both phases to obtain a volumetric continuity equation,which is used to form the pressure equation.The volumetric continuity equation for a two-phase system is

where the mixture velocity U= αaUa+ αbUb.Substituting the estimated velocities from Eq.(25)into Eq.(26),the pressure equation is expressed as

Like done for the phase-intensive momentum equation,the discretized form of the pressure equation can be obtained by formulatingEq.(27)at cell faces.Followingabove discretization procedures,it is not difficult to be extended to flows of more than two phases.

2.4.Phase volume fraction equation

Rusche [11] and Oliveira and Issa [28] discussed some techniques for solving the phase fraction equation as well as advantages and disadvantages of these methods.Here,the equation for the dispersed phase volume fraction is used:

This approach implicitly couples the two phases through the relative velocity Ur.

2.5.Turbulence model

Experimental studies have shown that when the gas holdup exceeds a certain limit(could be as small as6%),both phases tend to fluctuate as one entity[48,49],which suggests that the use of one set of turbulent equations for the mixture phase is appropriate.The form of the mixture turbulence model is identical to that in Behzadi et al.[50],though some different models for the shear-induced and bubble-induced turbulence are used in this work.

The following relations are used:

where Ctis the turbulent response function representing the ratio of the dispersed phase velocity fluctuations to those of the continuous phase[51].The turbulent properties of the mixture are expressed by

The mixture turbulent kinetic energy generation term is related to those of both phases by

where AD= αaαbKD.

For modeling of the bubble-induced turbulence,two different approaches are commonly used:the method of additional source terms and the method of additional viscosity.Many authors[4,6,52–54]presented their own models using the former method.The later method is usually referred as the Sato model[55,56].It is the simplest model for consideration of the bubble influence on the liquid,but it underestimates the influence as pointed out by Sokolichin et al.[4].The model of Malin and Spalding[57]and Davidson[58],which belongs to the method of additional source terms,is employed here:

The turbulent response function Ctis related to the volume fraction[11,50]:

where Ct0is the value of Ctat αa≈ 0 given by[51]

3.Population Balance Equation

3.1.Cell average technique

Based on the work of Kumar and Ramkrishna[22],Kumar et al.[26]presented a cell average technique,which has all the advantages of the fixed pivot technique but alleviates the problem of over-prediction to a great extent.As investigated by Kumar et al.[59],the computing time of the cell average technique is comparable to and in some cases is even less than that of the fixed pivot technique.The cell average technique has just been validated in homogeneous systems,and thus we will couple it with CFD to solve the PBE in a poly-dispersed gas–liquid flow in this work.The PBE for a spatially inhomogeneous system can be described as

where αi=xiNiis the volume fraction for the ith bin(the ith bin is bounded by its lower and upper boundaries,xi-1/2and xi+1/2),xiis the representing volume for the bin,and Nidenotes the number of bubbles per unit volume for that bin.The velocity Uais the average velocity of the gas phase,which is to be calculated from the momentum equation of the gas phase.The assumption that all bubbles travel at the same mean velocity reduces greatly the computational time and resource,otherwise N gas phase momentum equations need to be solved,in addition to the greater difficulty in resolving the gas mass conservation.The termsrepresent the cell average modified birth terms in the ith bin due to aggregation and breakage,respectively. However,the death termsdo not need to be modified and remain the same as those in the fixed pivot method.They are given as follows:

where βi,k= β(xi,xk)= β(di,dk)is the aggregation kernel of two bubbles with volumes xiand xk(or diameters diand dk),Sk=S(xk)represents the breakage rate of a bubble of size xkand b(x|xk)is the fragmentation distribution function that contains the information on fragments produced by a breakage event.For the binary breakage,Mkequals to 2.The Heaviside step function is defined by

and the function used for the redistribution of particles is given as

The average volume of all new born particles in the ith bin,is computed as

where the total volume fluxes into the ith bin as a result of aggregation and breakage are expressed by

3.2.Aggregation and breakage models

The aggregation model proposed by Prince and Blanch[60]and the breakage model of Lehr et al.[61]are used in this work.These two models are widely used in bubble columns,and there is no tuning parameter in them.According to Prince and Blanch[60],collisions between two bubbles in a turbulent environment occur due to three main mechanisms:turbulence,buoyancy and laminar shear,as shown in Fig.1.It is assumed that collisions from these various mechanisms are cumulative. Only collisions due to turbulence and buoyancy are considered in this work because of the low superficial gas velocity in the bubble column simulated in this work.The aggregation rate is expressed as

where ωt(di,dj),ωb(di,dj)andP(di,dj)are the turbulent collision rate,the buoyancy-driven collision rate and the collision efficiency,respectively.The turbulent collision rate is

Fig.1.Three bubble collision mechanisms according to Prince and Blanch[60]:(a)turbulence;(b)buoyancy;(c)laminar shear.

The buoyancy-driven collision rate is given as

Here uriis the rise velocity of bubbles,and can be expressed as a function of size:

The collision efficiency is a function of the contact time τijbetween bubbles and the time required for bubbles to aggregate,tij,as

Here τijand tijare given by.

where rijrepresents the equivalent radius,rij=2（ri-1+rj-1）-1,and h0and hfare the initial thickness and the critical thickness of the film,h0=1×10-4m and hf=1×10-7m.

The breakage rate of Lehr et al.[61]is calculated by

The daughter bubble size distribution is described by

where x' and d' represent the volume and diameter of the daughter bubble.Noting that

4.CFD-PBE Modeling of a Bubble Column

4.1.Numerical details

Coupling the PBE with CFD, a solver was developed in the platform of Open FOAM-3.0.0.The simulated bubble column is the same as that in the experiment conducted by Bhole et al.[62].The diameter of the column is D=0.15m,and the static liquid height is h=0.9m.The column is operated at the superficial gas velocity of 0.02 m·s-1.The dispersed and continuous phases are air and water,respectively.Calculations for this essentially axisymmetric problem are performed on a 5°wedge geometry.The grid of30×300(radial by axial)hexahedral cells is adopted in this work to ensure grid-independent solutions.The bubble size range of 1 to 15 mm is divided into 15 classes.

Fig.3.Predicted results at different axial heights(h/D)at t=0.5 s with PBE solved both in and out the PIMPLE Loop.

The coupling between the pressure and the velocity is dealt with the transient PIMPLE algorithm,which is a merged PISO-SIMPLE algorithm.It is very expensive to solve a transient problem using the PISO algorithm for a long time period especially in complex geometries,as the time step in PISO is strictly limited by the Courant number which must be smaller than 1.In the PIMPLE algorithm,a PIMPLE loop in which SIMPLE is used is run for every time step.The PIMPLE algorithm has been employed in some recent PBE simulations based on OpenFOAM[16,63,64].Fig.2 illustrated the solution procedures for the CFD-PBE modeling in this work.It is straightforward to solve the PBE in the PIMPLE loop[Fig.2(a)],while it is possible to do this out the loop[Fig.2(b)].Under-relaxation is used for the PIMPLE loops except the last one to ensure stable solutions,for example,the value of 0.3 for the pressure andmomentumequations and 0.5 for other equations.No under-relaxation is used for the last loop,thus ensuring time consistency.

Fig.4.Predicted results at different axial heights(h/D)at t=3 s with PBE solved both in and out the PIMPLE Loop.

The simulated time period is 150 s. The implicit Euler scheme is used for the time integration.The adaptive time step is controlled according to the maximum Courant number,which is set to 0.4,1.0,1.5 and 3.0 for 0–10 s,10–15 s,15–40 s and after 40 s,respectively.The advective terms in the turbulence equations are interpolated with the Guass upwind scheme,whereas the flux limiter Gamma scheme is adopted for other equations.The diffusive terms are interpolated with the Guass linear orthogonal scheme,and for other explicit divergence terms the Guass linear scheme is used.The preconditioned conjugate gradient(PCG)solver combined with the pre-conditioner of diagonalincomplete-Cholesky(DIC)is used for solving the discretized pressure equation,and the smoother solver with the symmetric Guass-Seidel as the smoother is used for solving other sets of linear equations.All the simulations were performed in parallel mode(four nodes)on a Linux computer(Intel Core i7 CPU,2.8 GHz,4 GB RAM).

Table 1Time for solutions with PBE solved in and out PIMPLE loop

At the top outlet surface,the liquid axial velocity is set to zero while zero gradient is used for other variables.Wall functions are used for the turbulence equations.No-slip condition is enforced for the velocity fields at walls while zero gradient for other fields except for the pressure field.The fixed flux pressure boundary condition is applied to the pressure field at walls,which adjusts the pressure gradient so that the boundary flux matches the velocity boundary.

4.2.Results and discussion

Fig.5.Comparison of simulated results with experimental data[62].(Simulation 1:breakage rate of Lehr et al.[61];simulation 2:multiplying the breakage rate of Lehr et al.[61]by 10).

Simulations with PBE solved in and out the PIMPLE loop are both conducted.Figs.3 and 4 illustrate the results at different axial heights(h/D)from the two calculations for the time 0.5 and 3 s.It is shown that differences in results from these two solution procedures are very small and negligible in most of time.Differences in results from the two proceduresat other times(not shown)are similar to those shown in Figs.3 and 4.

Table 1 gave the executed time and the elapsed clock time for these two calculations.The clock time is somewhat larger than the executed time due to file reading and writing processes.When solving the PBE in the PIMPLE loop,the required computation time is almost four times that of solving out the loop.The studied problem in this work is two-dimensional and axisymmetric,and the geometry is very simple.However,many actual gas–liquid problems,such as the liquid phase methanol synthesis and Fischer-Tropsch synthesis etc.,belong to high turbulent flows with intrinsic transient and 3D characteristics[1].Computation time for CFD-PBE modeling of these problems becomes an important issue then.Moreover,for simulations on crystallization,the crystal size range is generally divided into more than 30 classes,and then almost an order of magnitude time is expected to be saved when solving the PBE out the PIMPLE Loop.

As reviewed above, simulated results are considerably dependent on adopted models for the interphase forces and the bubble-induced turbulence as well as combinations of these models.Besides these models,choice of aggregation and breakage models should also be considered.Chen et al.[65,66]studied bubble columns by CFD-PBE modeling.The predicted gas holdups and axial liquid velocities agreed well with the experimental results.They investigated the effect of different breakup and aggregation models,and concluded that the choice of these models did not have a significant impact on the simulated results.However,they usually increased the breakage rate by a factor of 10 in their simulations.Wang et al.[5,67]and Wang and Wang[46]also investigated the effect of different aggregation and breakage models,showing that the bubble size distributions were quite sensitive to the models used.Olmos et al.[68]multiplied the aggregation and breakage rates by a factor of 0.075 to fit their experimental data.Contradictory conclusions about the effect of aggregation and breakage models on the results have been drawn.This might be due to different other models adopted in their works.

Investigation of the aggregation and breakage models is not the focus of this work.However,we did a simulation with the breakage rate of Lehr et al.[61]multiplied by a factor of 10 following Chen et al.[65,66].Fig.5 shows the predicted results and experimental measurements by Bhole et al.[62].The aggregation model of Prince and Blanch[60]and the breakage model of Lehret al.[61]were used.The predicted results are in generally reasonable agreement with published experimental data whereas the modeling precision could be improved further with more sophisticated closure models, the bubble-induced turbulence models and aggregation and breakage models as well as combinations of these models.By multiplying the breakage rate by 10,the results are affected little except the bubble diameter.The predicted gas holdup is somewhat higher than the experimental data,meaning more suitable drag models should be used in the future.The radial distribution of the results is significantly affected by the adopted lift force and turbulent dispersion force models.And near the wall,the distribution is influenced by the lift force and wall lubrication force.The predicted turbulence dissipation rate is strongly dependent on adopted turbulence models and the bubble-induced turbulence models.This in turn influences the choice of the aggregation and breakage models as the aggregation and breakage rates are strongly dependent on the turbulent dissipation rate.

4.3.Further remarks

We just performed simulations on a simple geometry and one operating condition,but the developed solver will provide us a good platform for further investigations on the CFD-PBE modeling of bubble columns in OpenFOAM.It is obvious that more sophisticated closure models,bubble induced models and aggregation and breakage models as well as combinations of these models should be studied further to obtain better results.A suitable set of these models,which are at least applicable in the regime of bubbly flows,is hopefully provided in our solver in OpenFOAM.

5.Conclusions

A general CFD-PBE solver for gas–liquid poly-dispersed flows is developed in the framework of OpenFOAM in this work.This solver could be extended to poly-dispersed gas–liquid systems of different volume fractions(theoretically 0–100%with modifications to the mass transfer model when gas is not fully dispersed or the dispersed phase changes to the continuous one) or even stratified flows (using the single static pressure,its gradient will suffer great change across a stratified flow).Implementation of this solver in OpenFOAM is illustrated in detail.The PBE is solved by the cell average technique [26],which is an improved version of the widely used fixed pivot technique[22].The coupling between the pressure and the velocity is dealt with the transient PIMPLE algorithm.Conclusions could be drawn as follows:

(1)The PBE could be solved out the PIMPLE loop to save much computation time while still preserve the time information on variables.This is important for future CFD-PBE modeling of high turbulent flows with intrinsic transient and 3D characteristics.

(2)Results show generally reasonable agreement with the published experiment data whereas the modeling precision could be improved further with more sophisticated momentum transfer models,bubble-induced turbulence models and more suitable aggregation and breakage models as discussed by the authors.

Appendix A.Discretization of momentum transfer terms

The magnitude,linearity and uniformity of the interphase momentum transfer term in the momentum equation affect the robustness and stability characteristics of the two- fluid solution procedure.The treatment for the interphase transfer terms has been discussed by Rusche[11].The virtual mass is treated semi-implicitly:

where the explicit substantive derivatives are given by.

For the drag term, one part is implicitly discretized, and another part is not discretized in the momentum equations and is moved to the pressure equation.

The implicit treatment for the functional form of the lift force is extremely difficult.Thus,an explicit treatment is adopted though it will compromise convergence.The wall lubrication term is explicitly treated.Turbulent drag term is also not discretized in the momentum equations and is moved to the pressure equation.Following above procedures,the discretized form for the momentum transfer terms,Eq.(17),could be got.