Implementation of a Moments Model in OpenFOAM

Implementation of a Moments Model in OpenFOAM for Polydispersed Multiphase Flows

Jo?a o N.E.Carneiro*,Volker Kaufmann,Wolfgang Polifke

Lehrstuhl f¨u r Thermodynamik,TU-M¨u nchen,Boltzmannstr.15,85748Garching,

Germany

*carneiro@td.mw.tum.de

Keywords:Polydispersed?ows,Moments model,Presumed number density function

Abstract

Accurate simulation of multiphase?ows with dispersed particles still represents

a great challenge and an open?eld for research.Examples can be found in several

applications such as bubbly?ows in aerated stirred vessels,atmospheric aerosols,

cloud physics,spray combustion,etc.One of the most important modeling require-

ments is the ability of capturing the polydispersity in these types of?ows.The scope

of this work is the further development of a novel approach to simulate multiphase

?ows with dispersed particles(either droplets or bubbles)based on the moments

of the size distribution function(therefore called Moments Model).The idea ex-

plored here is to use the Eulerian-Eulerian framework(using as a starting point

the twoPhaseEulerFoam solver),while taking into account the polydispersity by

transporting the moments of the particle size distribution function.Each moment

is transported with a di?erent velocity(the moment average velocities),determined

in an approximate manner through a concept using particle relaxation times and

force balances.Closure of the moments transport equations is achieved with help of

a presumed number density function(pNDF)approach.The momentum exchange

terms are evaluated using the local instantaneous Sauter Mean Diameter of the size

distribution function,providing coupling with the transported moments.

Results presented here–for a test case where buoyancy and(Stokes)drag are considered–show that the model implemented in OpenFOAM is indeed able to cap-

ture polydisperse e?ects e?ciently when compared to Multi-Fluid solutions obtained

with CFX.For more complex types of interaction(including lift and turbulent dis-

persion,for example),other closure relations are required,which will be the subject

of future work.

1Introduction

The simulation of polydispersed multiphase?ows normally requires a great compu-

tational e?ort with standard techniques.Recently,alternatives to the Lagrangian

or the Multi-Fluid approaches have been developed,which makes it possible to ac-

count for population dynamics in an e?cient manner.Most of these methods rely

upon a formulation involving moments(or some other information)of the parti-

cle size distribution function,and are divided into three categories:MOM(Method

of Moments),QMOM(Quadrature Method of Moments)and DQMOM(Direct

Quadrature Method of Moments)[1,7,8].

As an open source code,OpenFOAM is a very interesting alternative to commer-cial CFD solvers for the implementation of models for polydispersed?ows(see[10]

for example).The fundamental reason for that is the?exibility to incorporate addi-tional transport equations containing phase independent?uxes within a multiphase ?ow solver,and the possibility to test di?erent types of formulations for the mass, momentum and energy exchange terms with the continuous phase,as well as di?er-ent models to compute population dynamics due to break-up,coalescence or phase change,which are highly dependent on the speci?c application.

The objective of the present work is the implementation of a Moments Model to simulate polydispersed multiphase?ows in OpenFOAM,testing it for a simple con?guration comprising a channel?ow with a dilute mixture of water and air bubbles of di?erent sizes.Results obtained with a simpli?ed closure for the Moments Model are compared with a spectral solution obtained with the Multi-Fluid model in ANSYS-CFX11.0.Such a test case is particularly interesting once it separates the bubble transport from population dynamics.Hence,it allows to assess the ability of the model to correctly represent the segregation of bubbles with di?erent sizes due to buoyancy e?ects alone.

2Model Description

The idea explored here is to use the Eulerian-Eulerian framework of the Two-Fluid Model([5]),while considering polydispersity by transporting the moments of the particle size distribution function.The moments of the number distribution function f(D)are de?ned as:

M(k)=

Dmax

Dmin

D k f(D)dD.(1) A physical interpretation for each of the moments follows:

M(0):total number of particles(per unit cell-volume)

M(1):sum of the particle diameter( )

πM(2):total surface area of particles( )

π

6M(3):total volume of particles( );or the local volume fraction of the

dispersed phase(α)

The moments transport equations can be written in a generalized form as:

?M(k)

?t

+?(u(k)M(k))=S M(k),(2) with S M(k)the source terms resulting from the population balance dynamics.The moment average velocities u(k)’s(averaged over the k th moment)appearing in Eq.

(2)are in turn given by:

u(k)=

1

M(k)

Dmax

Dmin

uD k f(D)dD,(3)

which are the mean velocities at which the moments are convected1.There is no

reason,e.g.,for u(2)and u(3)to be the same,as they are surface area and mass 1Implicit here is the assumption that particles with the same size have at a given location the same velocity,i.e.,f=f(D,x,t),and not f=f(D,u,x,t)

averaged velocities,respectively.Moreover,except for u(3),which actually is deter-mined by the dispersed phase momentum equation,none of the u(k)’s is known a priori.Therefore,some modeling e?ort is needed at this point2.

The“Equilibrium Eulerian method”was introduced in[3,4].In the limit of very small particle relaxation times,particle terminal velocity deviates very little from the continuous phase velocity,scaling with the?rst order of the relaxation time.On the light of this idea,Bollweg et al.[2]propose a linear interpolation between the continuous phase and reference particle velocities(with size larger than minor groups)to determine minor particle size velocities.Here,we make use of this approximation in order to derive relations for the determination of moment average velocities from known velocities,thus allowing the development of a closure for Eq.

(2).

In the present study,the driving force which leads to the segregation of bub-bles with di?erent sizes results from a balance between drag and buoyancy.The Lagrangian equation of motion for each particle can be written in this case as:

du r,i dt =

1

τi

(u T,i?u r,i),(4)

with the relaxation timeτi(assuming Stokes Flow)de?ned as

τi=ρd D2i

18μc

,(5)

u r,i and u T,i are the relative and terminal velocity of bubble class i,respectively. They are in turn de?ned as:

u r,i=u i?u c(6)

and

u T,i=?ρgD2i

18μc

,(7)

with?ρ=ρd?ρc.By using the approximation proposed in[2],and integrating over the particle spectrum,it is possible to write:

u(k)=u c+u(k)

T ?

τ(k)

τ0

(u c?u0+u T,0),(8)

for the interpolation of moment average velocities with k≤2.In the above equation, the averaged particle relaxation timeτ(k)is de?ned as:

τ(k)=

1

M(k)

Dmax

Dmin

τD k f(D)dD,(9)

and the index0refers to reference properties used in the interpolation procedure. Since M(3)and u(3)are given by the solution of the continuity and momentum equations,they are primary candidates for reference properties.This approximation will show to be su?cient for the present purposes,but can not be expected to be valid in a broader context,e.g.if Stokes approximation is no longer valid or non-drag forces are considered.This will be subject of future investigation.

By introducing Eq.(5)into(9),it is possible to write:

τ(k)∝M(k+2)

M(k)

.(10)

2In fact,averaged equations can be formulated for each u(k),by integrating the dispersed phase momentum equations over the size spectrum,but this is out of the scope here

From the moments transport equation(2)–setting for simplicity the source term to null–and with help of relation(8),it follows that:

?M(k)?t +?(u c M(k))=?[(u c?u(k))

τ(k)

τ0

(u c?u0)

M(k)](11)

If Eq.(10)is substituted into(11),it follows that to determine a certain lower order moment M(k),M(k+2)must be known.This gives rise to a closure problem.

A simple way to overcome this issue is to assume the functional form of f(D) (a presumed Number Density Function–pNDF),calculating unknown moments by reconstruction of the distribution using lower order moments.Here,we choose the Gamma distribution function,as it is able to represent a wide variety of size distribution spectra.Furthermore,its reconstruction from lower order moments is straightforward,making it very suitable for the present purposes(although other functional forms can be easily incorporated in the model).Hence,if f(D)is given by a Gamma function:

f(D)=C0Dβ?1e?Dα

αβΓ(β)

(12)

with

C0=M(0)

α=M(2)M(0)?(M(1))2

M(0)M(1)

β=

(M(1))2

M(2)M(0)?(M(1))2

,(13)

any higher order moments can be explicitly calculated as:

M(k)=C0Γ(β+k)αk

Γ(β)

,k≥3.(14)

The complete model requires the solution of Eq.(2)for M(0)-M(3),with M(4) given by Eq.(14).The solution of the coupled moments transport equations enables the calculation of the Sauter Mean Diameter of the distribution,determined by

D32=M(3)

M(2).This mean diameter is in turn used to calculate the e?ective drag

used in the phase momentum equations,hence providing coupling to the moments equations.

3Implementation in OpenFOAM

OpenFOAM is a C++class library which can be easily used to develop CFD codes for a wide variety of problems,including multiphase?ows,either by starting from existing solvers or by developing completely new ones,using capabilities of the?nite volume method already available[6,12].

The objective here is to illustrate the incorporation of additional transport equa-tions for the moments within a two-phase?ow solver using a simple closure for u(k) through a concept using relaxation times and force balances,as illustrated above.

The solution of additional transport equations for the moments within the twoPhaseEulerFoam3is done following the pressure-velocity system,as it would be the case of the k? equations,for example.Regarding the implementation in the code itself,the moments transport equations can be de?ned within the tensor derivative class fvScalarMatrix,with all terms being treated implicitly through the class fvm,as depicted by the following code snippet:

surfaceScalarField C=tauk*(phib-phia)/taua;

fvScalarMatrix mkEqn

(

fvm::ddt(mk)

+fvm::div(phib,mk)

+fvm::div(C,mk)//convective correction

);

mkEqn.relax();

mkEqn.solve();

According to Weller’s notation[11],the discretised form of the moments trans-port equations is given by:

[?[M(k)]

]+[?·(φb[M(k)]f(φb,UD,))]+[?·(C[M(k)]f(C,UD,))]=0,(15)?t

where the Euler implicit scheme for the time derivative and upwind di?erencing face values of M(k)’s for the convective terms were used.The relaxation times taua and tauk for the correction coe?cient C are de?ned as surfaceScalarField types, given by Eqs.(5)–calculated using the Sauter Mean Diameter–and(10).Thus, segregation of the moments throughout the?ow?eld is achieved by considering di?erent convective corrections for each transport equation.Moreover,coupling with the phase equations is achieved through the drag force,where the dispersed phase diameter is simply the updated value of the Sauter Mean Diameter through use of the pNDF approach described earlier.

4Test Case

The channel?ow represents a simple con?guration which allows to assess the ability of the model to predict e?ects of polydispersity due to buoyancy through convection of the moments in the?ow?eld.The set up consists of a channel with2cm height and10cm length.A dilute mixture of water and air bubbles at standard conditions was considered.The dispersed phase volume fraction at the inlet corresponded to 10?5and the mixture Reynolds number was chosen to be approximately2000.A constant pressure was kept at the outlet,with zero gradient boundary conditions for all other variables.For simplicity,slip conditions at the walls were applied,in order to avoid accumulation of bubbles at the top wall.Boundary conditions for the moments were taken from the size distribution function applied to the Multi-Fluid case,which is shown in Fig.(1).The range of sizes considered was10-150μm,with

a total of15bubble size classes and maximum particle Reynolds number with order

of magnitude of1.

3A full description of the Two-Fluid Model implemented in OpenFOAM can be found in[9].

f (#/m 4)Figure 1:Initial Distribution used in the spectral case.

The simulations were carried out with a 2D hexahedral mesh of 50x150points

and a ?xed time step of

10?5s.For the solution of the discretised equations,the

pre-conditioned (bi-)conjugate gradient linear solvers were used,with a convergence

criteria of 10?10for both pressure and volume fraction.

Figure (2)shows the distribution of the Sauter Mean Diameter inside the channel

obtained with the Moments Model in OpenFOAM.An increase towards the upper

wall can be clearly observed.This re?ects the fact that the vertical velocity induced

by buoyancy tends to increase with the bubble diameter.Since the contribution of

bigger bubbles becomes more important for higher order moments,M (3)rises faster

than M (2),leading consequently to an increase in D 32.

Figure 2:Contours of the Sauter Mean Diameter obtained with the Moments Model in OpenFOAM.

Simulations were also performed with the Multi-Fluid Model available in AN-

SYS CFX 11.0and used as a reference solution to validate the Moments Model

implemented in OpenFOAM.Thus,sets of conservation equations for mass and

momentum were solved for the Eulerian phases (a total of 15dispersed phases and

1continuous phase),and the post-processed local particle size distribution func-

tion,as well as all its moments,were compared to the solution obtained with the

Moments Model.

The axial evolution of the moments (normalized by the inlet values)at both

walls is shown in Fig.(3)for both cases.The overall behavior is well represented

by the Moments Model,also comparing quantitatively well against the Multi-Fluid

approach:as bigger bubbles tend to migrate faster to the upper side of the channel,

higher moments tend to decrease faster at the bottom wall and increase faster at

the top wall.

(k )(k )(a)

(k )(k )(b)

Figure 3:Comparison of the evolution of the normalized moments along the channel at the bottom (a)and top (b)walls.

Figure (4)shows the evolution of the bubble size distribution function at several

positions inside the channel (corresponding to axial positions x =1,5;5,7cm and

vertical positions y =0,1;10;19,9mm).As bubbles tend to move towards the

upper side of the channel,the number density is higher for the vertical coordinate

corresponding to y =19,9mm,and smaller at y =0,1mm.Furthermore,a very

reasonable agreement can be noted between both methods,with the Moments Model

being able to reproduce not only the shape of the distribution function,but also its

variation along the channel.However,the distribution obtained with the Moments

Model near the upper wall at the most downstream axial position in the channel

deviates considerably from that of the Multi-Fluid solution.This occurs because

the shape of the spectral distribution at that location is dominated by the high

accumulation of big bubbles near the top wall,which is not appropriately reproduced

by a Gamma distribution.

Multi-Fluid (CFX)

Moments Model (OpenFOAM)f (#/m 4)f (#/m 4)f (#/m 4)Figure 4:Comparison of the evolution of the bubble size distribution function along the channel.

5Conclusions and Future Work

A new Moments Model for polydispersed multiphase ?ows with a simpli?ed closure

for the moment average velocities was presented and its implementation within the

twoPhaseEulerFoam solver of OpenFOAM was illustrated.The model was tested

for a simple con?guration comprising of a channel?ow with a two-phase mixture of water and air bubbles.Very reasonable agreement for the evolution of the bubble size distribution function and its moments inside the channel was obtained against simulations carried out with the Multi-Fluid Model available in ANSYS CFX11.0, showing that the model was able to capture e?ects of polydispersity due to buoyancy. It is clear,however,that a correct description of the polydispersity hinges on the determination of u(k),which was done here in an approximate manner through a concept using particle relaxation times and force balances.The further development of the model involves the implementation of extra averaged momentum equations to determine the moment average velocities,with special closure needed to treat other types of interactions(e.g.,lift and turbulent dispersion forces).OpenFOAM has proven to be a very appropriate tool for this task.

6Acknowledgements

The authors want to acknowledge the?nancial support given by CNPq(Brazilian National Council for Research and Technological Development)and GRS(Gesellschaft f¨u r Reaktorsicherheit).

References

[1]J.C.Beck and A.P.Watkins.On the development of a spray model based on

drop-size moments.Proc.R.Soc.Lond.A,459:1365–1394,2003.

[2]P.Bollweg,A.Kaufmann,and W.Polifke.Derivation and application of a

poly-celerid method for poly-dispersed two-phase?ows.6th International Con-ference on Multiphase Flow,Leipzig,Germany,2007.

[3]J.Ferry and S.Balachandar.A fast eulerian method for dispersed two-phase

?ow.Int.J.of Multiphase Flow,27:1199–1226,2001.

[4]J.Ferry and S.Balachandar.Equilibrium expansion for the Eulerian velocity

of small particles.Powder Technology,125:131–139,2002.

[5]M.Ishii.Thermo-Fluid Dynamic Theory of Two-Phase Flow.Eyrolles,Paris,

1975.

[6]H.Jasak.Error analysis and estimation for the?nite volume method with

applications to?uid?ows.PhD thesis,Imperial College of Science Technology and Medicine,London,UK,1996.

[7]V.John,I.Angelov,A.A.¨Onc¨u l,and D.Th′e venin.Techniques for the re-

construction of a distribution from a?nite number of moments.Chemical Engineering Science,62:2890–2904,2007.

[8] D.L.Marchisio and R.O.Fox.Multiphase Reacting Flows:Modelling and

Simulation.SpringerWienNewYork,2007.

[9]https://www.wendangku.net/doc/e410350372.html,putational Fluid Dynamics of Dispersed Two-Phase Flows at

High Phase Fractions.PhD thesis,Imperial College of Science Technology and Medicine,London,UK,2002.

[10]L.F.Silva and https://www.wendangku.net/doc/e410350372.html,ge.Implementation of an Eulerian multi-phase model

in OpenFOAM and its application to polydisperse two-phase?ow.OpenFOAM International Conference2007,London,UK,2007.

[11]H.G.Weller.A code independent notation for?nite volume algorithms.Tech-

nical report,TR/HGW/01,OpenCFD,2005.

[12]H.G.Weller,G.Tabor,H.Jasak,and C.Fureby.A tensorial approach to com-

putational continuum mechanics using object-oriented https://www.wendangku.net/doc/e410350372.html,puters in Physics,12(6):620–631,1998.