Eulerian–Lagrangian simulation of bubble coalescence in bubbly flow using the spring-dashpot model

Jing Xue ,Feiguo Chen ,Ning Yang ,*,Wei Ge

1 State Key Laboratory of Multiphase Complex Systems,Institute of Process Engineering,Chinese Academy of Sciences,P.O.Box 353,Beijing 100190,China

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

1.Introduction

In bubbly flows the bubble size distribution(BSD)is of critical importance to characterize the fluid flow and mass transfer and reflect the bubble coalescence and breakage.In the traditional Eulerian–Eulerian methods,both the liquid and the gas phase are treated as interpenetrating continuum,and the behavior of individual bubbles such as the position,velocity and size of each bubble is not directly simulated.To reconstruct the BSD,a population balance equation(PBE)is required to be integrated into the Eulerian–Eulerian method.The PBE also needs the coalescence kernel function and breakage kernel function,both of which are based on statistical models.Various PBE kernel functions have been proposed in literature,such as Coulaloglou and Tavlarides[1],Prince and Blanch[2],Luo and Svendsen[3],Alopaeuset al.[4],Lehret al.[5].

In Eulerian–Lagrangian method,the bubble size distribution could be explicitly modeled since each bubble is tracked individually.When two bubbles collide with each other,they may bounce off or coalesce,depending on the relative velocity and the radii.The criteria to determine coalescence and breakage events need to be established.Based on the hard-sphere method of Delnoijet al.[6],a few researchers have applied some coalescence and break-up models in Eulerian–Lagrangian framework.Sommerfeldet al.[7]proposed a coalescence model based on the stochastic model in their simulation without considering the bubble breakage.Darmanaet al.[8]reported a deterministic coalescence model.Lauet al.[9]established a complete coalescence and breakage model and predicted the BSD for the experiments of the Deen case.In the work of Sungkornet al.[10],a stirred reactor was simulated with bubble breakage and coalescence.Jainet al.[11]simulated a micro-structured bubble column.Gruberet al.[12]established a stochastic coalescence and break-up model within the Eulerian–Lagrangian framework.Lauet al.[13]presented a review of break-up and coalescence models and studied the influence of several parameters on the simulation of BSD.

To our knowledge there is no report about the spring-dashpot collisional model in the simulation of bubble coalescence,though the model is more physical than the hard-sphere model to model the soft bubbles and calculate the bubble contact time.This may be partly attributed to the lack of bubble interaction parameters for spring-dashpot model since few literatures concerned on the performance of the spring-dashpot model in gas–liquid systems.We have already verified the applicability of the spring-dashpot model for gas–liquid systems in our previous work[14],and in this paper we further establish a preliminary modeling framework to incorporate the bubble collisions by using the spring-dashpot model and the bubble coalescence and breakage into the Eulerian–Lagrangian simulation.

2.Model Formation

2.1.Bubble dynamics

Each bubble is tracked individually in the Eulerian–Lagrangian simulation.All the forces exerted on a bubble is summed up,and then the velocityvand positionris updated within a fixed time step by:

where fbrepresents the buoyancy force,fgthe gravity,fdthe drag force,flthe lift force and fvmthe virtual mass force.The buoyancy and gravity force are written together as:

2.2.Liquid phase hydrodynamics

The averaged Navier–Stokes equation is employed for liquid phase.The gas volume fraction and the momentum exchange rate are added to reflect the existence of bubbles:

The item Φ represents the interfacial forces exerted by fluid on a bubble like buoyancy,drag,lift and virtual mass,then averaged on each cell:

The final equation becomes:

2.3.Collision model

The simplified linear model by Cundall and Strack[16]is used to calculate the contact force between bubbles.The normal contact force between two bubblesaandbis calculated from:

where δnrepresents the overlap between two bubbles,κnthe normal spring stiffness coefficient,ηnthe normal damping coefficient,nabthe normal unit vector and uabthe relative velocity.The tangential component is given in a similar manner:

δtis the tangential displacement,κtthe tangential spring stiffness coefficient, ηnthe tangential damping coefficient,and μ the friction coefficient.

In our previous work,the experimental systems of Becker case[17]and Deen case[18]were simulated without considering bubble coalescence or breakage.We found that the simulation of Becker case was insensitive to the normal stiffness coefficient,and the normal stiffness coefficient 1 N·m?1gave slightly more accurate results than 0.01 N·m?1or 50 N·m?1.The simulation of Deen case indicated that 1 N·m?1of the normal stiffness coefficient predicted the closest fluid velocity to the experimental data.In this paper,we will further incorporate the bubble coalescence and breakage models and investigate the influence of spring-dashpot model parameters on bubble contact time.

2.4.Coalescence model

The film drainage model of Prince and Blanch[2]is implemented in the simulation.Coalescence occurs when two bubbles collide and the time for the liquid film between the bubbles to drain is smaller than the bubble contact time.The film drainage time is given as:

where θ0is the initial film thickness and θfthe final thickness.For airwater system θ0is taken as 10?4m and θfis 10?8m[19].req.is the equivalent radius as defined by Chesters and Hofman[20]:

In hard-sphere models the contact time should be zero since the collision is assumed to be instantaneous;whereas for modeling the coalescence the contact time should be a considerable positive value.Sommerfeldet al.[7]supplemented an additional approximation which assumes that the contact time is proportional to the deformation distance divided by the normal relative velocity.Darmanaet al.[8]and Lauet al.[13]adopted this method in their simulation.However,this approximation can be eliminated in the spring-dashpot model as the contact time is detected and recorded during the modeling of collision process.Hence the coalescence criterion in our model is given by:

2.5.Break-up model

There are several mechanisms and corresponding models in literatures accounting for break-up,including turbulent fluctuation,shear stress,interfacial instability.In this study we only utilize a simple break-up model for simplicity.The deformation stress arises from the velocity fluctuations and bubbles will break when the Weber number is greater than a critical value:

when a bubble is judged to break,a uniform random numberxbetween 0 and 1 is generated by the computer,then the inverse transform method is used to obtain the random number γ for the size distribution function.

The mother bubble is resized to become the larger daughter bubble and the other new smaller bubble is randomly placed around the larger one.It should be pointed out that the U-shape distribution is selected in this paper because its cumulative distribution function is easily obtained in our method.However,an M-shape distribution[5,22]can be more physical since it incorporates the influence of both additional pressure and surface energy.

2.6.Code implementation

The open source CFD package OpenFOAM is utilized in our work to establish a complete program for flow field calculation,bubble tracking,coalescence and break-up.The default gas–solid solver DPMFoam is rebuilt to make it consistent with the governing equations above.And we add the spring-dashpot collision model and the bubble coalescence and breakage models mentioned above.The break-up events are also handled in the code.For simplicity we assumed that the momentum of bubbles is conserved during the breakage or coalescence process and the velocity of each daughter bubble is equal to that of the mother bubble.A brief flowsheet of the process is shown in Fig.1:

Fig.1.Flowchart of the solver.

Fig.2 illustrates coalescence and break-up process.The new diameter,velocity and position after coalescence are calculated according to:

where Δr is a random vector with a magnitude of 0.6 × (da+db).

Fig.2.Illustration of the coalescence and break-up process.

3.Simulation Setup

3.1.The Deen case

A square bubble column(Fig.3)of the Deen case is simulated in the study.The column has a height of 0.45 m and a cross section of 0.15 m×0.15 m.The gas distributor is fixed at the center of the bottom with 49 pores.The bubble diameter is initially 4 mm according to Deenet al.[18].The superficial gas velocity is 4.9 mm·s?1,thus about 3290 bubbles are injected into the system per second.The no-slip condition for the liquid is imposed on the walls except the top where the freeslip condition is used.The collision between bubbles and sidewall are handled by the hard-sphere model with a restitution coefficient of 1.The large eddy simulation(LES)is selected for turbulence.The simulation time is set to 180 s,of which the first 20s are discarded due to start-up effects.More detailed parameters are given in Table 1.

3.2.Selection of the normal stiffness coefficient

In the spring-dashpot model for bubbles,the normal stiffness coefficient plays the important role.It determines the magnitude of the overlap between bubbles,thus the contact time and coalescence rate.However there is hitherto no experimental data for bubble stiffness coefficients,and we refer to the work of Chesters[23]in which the coalescence processes were modeled.

Fig.3.Schematic view of the bubble column.

Table 1Detailed simulation settings

Chesters estimated the interaction force between bubbles of the same radiusRas:

wherearepresents the intersection radius.If the normal overlap δnis small,this radius can be expressed as:

Sincef～ κnδnin the spring-dashpot model,the normal stiffness coefficient can be derived as:

For the air-water system in this paper,the surface tension is about 0.075 N·m?1,thus the normal stiffness coefficient is in the order of 0.5 N·m?1.In our previous work 1 N·m?1is found most suitable in the prediction of fluid velocity for the coalescence-free cases.In this paper,we further evaluated the different normal stiffness coefficient(0.5 N·m?1,0.75 N·m?1and 1.0 N·m?1)in the calculation.

4.Results and Discussion

Lauet al.[13]reported that different daughter size distributions and critical Weber numbers resulted in similar multimodal bubble size distribution in their deterministic coalescence and break-up models.The critical Weber number ranging from 1 to 12 leads to different heights of the peaks of bubble size distribution.Actually the critical number should be evaluated from experiments or direct numerical simulations.Here we only use the U-shape daughter bubble size distribution and a critical Weber number of 1 in the three cases.A snapshot of the simulation and the time-averaged fluid field are shown in Fig.4.Fig.5 shows the time-averaged vertical liquid velocity pro file at the height of 0.25 m.The solid line represents the experimental data of Deenet al.[24],in which kitchen salt was added to inhibit the bubble coalescence and break-up.The maximum velocity magnitude is higher than the same case without coalescence or break-up.The simulation without considering the coalescence and break-up is in better agreement with the experiment.The simulations of different stiffness coefficient show similar trends,and the use of models of coalescence and breakup plays more important roles in the simulation of vertical velocity.

Fig.4.(a)Snapshot of the bubbles in the column.(b)time-averaged liquid velocity on the mid-depth plane.

Fig.5.Mean vertical velocity at the height of 0.25 m.

4.1.Coalescence rates

The coalescence rates on the mid-depth plane are illustrated in Fig.6.Since the number density of bubbles is quite high at the inlet,collisions occur frequently and results in high coalescence rate there.Besides,coalescence also prefers to occur in the left-top and right-top corners.This can be attributed to the vortices near the upper corners,which make the bubbles to recirculate and increase the possibility to collide(Fig.7).The coalescence rate is sensitive to the normal stiffness coefficient.Higher normal stiffness coefficient leads to shorter contact time between bubbles and lower probability of coalescence.The normal stiffness coefficient of 1 N·m?1obtains a maximum coalescence rate of only 50·s?1;whereas for 0.5 N·m?1the maximum coalescence rate is 180·s?1.Also the coalescence rate at inlet seems to be underestimated at 1 N·m?1.

4.2.Bubble size distribution

The predicted BSDs(Fig.8)show the same trend with the particle imaging technique(PIV)measurement by Hansen[25],except the peak at 4 mm.As mentioned before,in this study only the shear stress for break-up is considered,the break-up rates are likely to be underestimated.Thus a large number of bubbles with initial diameters 4 mm go through the column without breakage.This may also be related to the critical Weber number or the method to calculate the mean square velocity difference in Eq.(16).

The simulation of κn=0.5 N·m?1seems better than that of κn=1 N·m?1and κn=0.75 N·m?1,which contains less bubbles with diameters of 4 mm and more bubbles with diameters＜2 mm.Since the stiffness coefficients are not directly relevant to the break-up processes,the difference can be only attributed to the coalescence processes.For κn=0.5 N·m?1,more bubbles with initial diameters 4 mm can first coalesce and then break up.Fig.9 indicates that at 0.5 N·m?1the break-up rate is also higher.This may imply that a lower normal stiffness coefficient achieves better balance of coalescence and break-up.

All the tails of the three BSDs terminate at about 6 mm,which is consistent with the results of stochastic approach of Gruberet al.[12],but in Hansen's experiment there are a significant number of large bubbles with diameters above 6 mm.This accordance might indicate that more suitable parameters for the film drainage model is required.

5.Conclusions

Fig.6.Coalescence rates on mid-depth plane.(a)–(c)Normal stiffness coefficient 1.0 N·m?1,0.75 N·m?1,0.5 N·m?1.

We have simulated a 3D bubble column with consideration of bubble coalescence and break-up under Eulerian–Lagrangian framework.We utilized the spring-dashpot model in bubbly flow to handle the collisions between bubbles and record the contact time.The latter is then incorporated into the film-drainage coalescence model to determine the occurrence of bubble coalescence.A simple model related to the critical Weber number is used for bubble break-up.The simulations with different normal stiffness coefficients were carried out to study the applicability of the spring-dashpot in bubble coalescence.It indicates that the coalescence rates and BSDs are both pertinent to the normal stiffness coefficient.A slight increase in normal stiffness coefficient results in a remarkable lower coalescence rate.A normal stiffness coefficient of 0.5 N·m?1seems give the highest coalescence rate and the best BSD.It seems the film drainage model underestimates the coalescence rate of large bubbles since both our simulation and the stochastic approach of Gruberet al.[12]contain less bubbles with diameters above 6 mm than the experimental observation of Hansen[25].This paper here presents a preliminary study and analysis of the Eulerian–Lagrangian simulation with the spring-dashpot model.It should be pointed out that the simulation depends on the model parameters in the spring-dashpot model,and the bubble coalescence and breakage models.These parameters may have to be determined by more micro-scale experiments or direct numerical simulations.Also the simulation results are influenced by numerical algorithms in collision detection and contact time record method,all of which needs further in-depth study.

Fig.7.Collision frequency on mid-depth plane.(a)–(c)Normal stiffness coefficient 1.0 N·m?1,0.75 N·m?1,0.5 N·m?1.

Fig.8.Overall bubble size distribution.(a)–(c)normal stiffness coefficient 1.0 N·m?1,0.75 N·m?1,0.5 N·m?1.

Fig.9.Break-up rates on mid-depth plane.(a)–(c)normal stiffness coefficient 1.0 N·m?1,0.75 N·m?1,0.5 N·m?1.

Nomenclature

Cddrag coefficient(dimensionless)

dbbubble diameter,m

EoEotvos number(dimensionless)

g gravitational acceleration,m2·s?1

ReReynolds number(dimensionless)

rbbubble location,m

Δttime step,s

Ufliquid velocity,m·s?1

Vbbubble volume,m3

vbbubble velocity,m·s?1

WeWeber number(dimensionless)

δnnormal displacement,m

δttangential displacement,m

εfliquid volume fraction(dimensionless)

ηnnormal damping coefficient,kg·s?1

ηttangential damping coefficient,kg·s?1

κnnormal stiffness coefficient,N·m?1

κttangential stiffness coefficient,N·m?1

ρbbubble density,kg·m?3

ρfliquid density,kg·m?3

σ surface tension,N·m?1

Acknowledgments

The authors wish to acknowledge the long term support from the National Natural Science Foundation of China(Grant No.91434121),Ministry of Science and Technology of China(Grant No.2013BAC12B01)and State Key Laboratory of Multiphase complex systems(Grant No.MPCS-2015-A-03)and Chinese Academy of Sciences(Grant No.XDA07080301).

[1]C.A.Coulaloglou,L.L.Tavlarides,Description of interaction processes in agitated liquid–liquid dispersions,Chem.Eng.Sci.32(1977)1289–1297.

[2]M.J.Prince,H.W.Blanch,Bubble coalescence and break-up in air-sparged bubble columns,AIChE J.36(1990)1485–1499.

[3]H.Luo,H.F.Svendsen,Theoretical model for drop and bubble breakup in turbulent dispersions,AIChE J.42(1996)1225–1233.

[4]V.Alopaeus,J.Koskinen,K.I.Keskinen,J.Majander,Simulation of the population balances for liquid–liquid systems in a non-ideal stirred tank.Part 2 — Parameter fitting and the use of the multiblock model for dense dispersions,Chem.Eng.Sci.57(2002)1815–1825.

[5]F.Lehr,M.Millies,D.Mewes,Bubble-size distributions and flow fields in bubble columns,AIChE J.48(2002)2426–2443.

[6]E.Delnoij,F.A.Lammers,J.A.M.Kuipers,W.P.M.van Swaaij,Dynamic simulation of dispersed gas–liquid two-phase flow using a discrete bubble model,Chem.Eng.Sci.52(9)(1997)1429–1458.

[7]M.Sommerfeld,E.Bourloutski,D.Br?der,Euler/Lagrange calculations of bubbly flows with consideration of bubble coalescence,Can.J.Chem.Eng.81(2003)508–518.

[8]D.Darmana,N.G.Deen,J.A.M.Kuipers,Parallelization of an Euler–Lagrange model using mixed domain decomposition and a mirror domain technique:Application to dispersed gas–liquid two-phase flow,J.Comput.Phys.220(1)(2006)216–248.

[9]Y.M.Lau,N.G.Deen,J.A.M.Kuipers,Bubble breakup in Euler–Lagrange simulations of bubbly flow,7th Int.Conf.on Multi-phase Flow(ICMF),Tampa,FL,2010.

[10]R.Sungkorn,J.J.Derksen,J.G.Khinast,Euler–Lagrange modeling of a gas–liquid stirred reactor with consideration of bubble breakage and coalescence,AIChE J.58(2012)1356–1370.

[11]D.Jain,Y.M.Lau,J.A.M.Kuipers,N.G.Deen,Discrete bubble modeling for a microstructured bubble column,Chem.Eng.Sci.100(2013)496–505.

[12]M.C.Gruber,S.Radl,J.G.Khinast,Coalescence and break-up in bubble columns:Euler–Lagrange simulations using a stochastic approach,Chem.Ing.Tech.7(2013)1118–1130.

[13]Y.M.Lau,W.Bai,N.G.Deen,J.A.M.Kuipers,Numerical study of bubble break-up in bubbly flows using a deterministic Euler–Lagrange framework,Chem.Eng.Sci.108(2014)9–22.

[14]J.Xue,Applicability of the soft-sphere model in bubbly flowMaster Thesis University of Chinese Academy of Sciences,China,2016.

[15]A.Tomiyama,Struggle with computational bubble dynamics,Third International Conference on Multiphase Flow,ICMF'98,Lyon,France,1998.

[16]P.A.Cundall,O.D.L.Strack,A discrete numerical model for granular assemblies,Geotechnique29(1979)47–65.

[17]S.Becker,A.Sokolichin,G.Eigenberger,Gas–liquid flow in bubble columns and loop reactors:Part II.Comparison of detailed experiments and flow simulations,Chem.Eng.Sci.49(1994)5747–5762.

[18]N.G.Deen,B.H.Hjertager,T.Solberg,Comparison of PIV and LDA measurement methods applied to the gas–liquid flow in a bubble column,Presented at the 10th international symposium on applications of laser techniques to fluid mechanics,Lisbon,Portugal,2000.

[19]R.Kirkpatrick,M.Lockett,Influence of approach velocity on bubble coalescence,Chem.Eng.Sci.29(1974)2362–2373.

[20]A.K.Chesters,G.Hofman,Bubble coalescence in pure liquids,Appl.Sci.Res.38(1982)353–361.

[21]C.Tsouris,L.L.Tavlarides,Breakage and coalescence models for drops in turbulent dispersions,AIChE J.40(1994)395–406.

[22]T.Wang,J.Wang,Y.Jin,A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow,Chem.Eng.Sci.58(2003)4629–4637.

[23]A.K.Chesters,The modeling of coalescence processes in fluid–liquid dispersions:A review of current understanding,Chem.Eng.Res.Des.Trans.Inst.Chem.Eng.Part A69(1991)259–270.

[24]N.G.Deen,T.Solberg,B.H.Hjertager,Large Eddy simulation of the gas–liquid flow in a square cross-sectioned bubble column,Chem.Eng.Sci.56(2001)6341–6349.

[25]R.Hansen,Computational and experimental study of bubble size in bubble columnsPh.D.Thesis Aalborg University,Denmark,2009.