[PPT] - Multiphase flows in complex geometries: a UQ perspective Matteo PowerPoint Presentation

SLIDE 1

Multiphase flows in complex geometries: a UQ perspective

Matteo Icardi

Stochastic Numerics group & SRI-UQ center (KAUST) ICES (UT Austin)

SRI-UQ Annual Workshop January 6-9, KAUST

SLIDE 2

Acknowledgements

KAUST

R. Tempone, H. Höel, B. Saad, Q. Long

SRI-UQ center and

Stochastic Numerics research group

POLITECNICO DI TORINO

R. Sethi, T. Tosco (Groundwater Eng.)
D. Marchisio,
G. Boccardo (Chem. Eng.)
C. Canuto, N. Quadrio (Math.)

P

. Asinari, D. Ceglia (Energy)

THE UNIVERSITY OF TEXAS AT AUSTIN

S. Prudhomme (now at EP Montreal)
M. Prodanovic (Petr. Eng.)
I. Babuska (ICES)
A. Passalacqua, R. Fox (Iowa State Univ.), D. Tartakovsky (UCSD)

2 / 34

SLIDE 3

Motivation I: Data-driven UQ

Forward and inverse UQ techniques are mature enough to be applied in real complex engineering problems Unfortunately in most applications, the main problem is defining/parametrizing the input uncertainties (parametric uncertainties, error models)

MOMENTS

Sometimes parameters are dependent and only known through finite order statistics (mixed moments)

Part I: Quadrature-based moment methods

SAMPLES

Many uncertainties can hardly be parametrized (random geometry)

Part II: Multi-level Monte Carlo

3 / 34

SLIDE 4

Motivation II: UQ as a physical closure problem

In multi-scale problems the concept of probability can be used for modeling

purposes. The forward UQ problem becomes a “closure” problem

REACTIONS: MICRO-MIXING PROBLEM

A+B k − → R SR = −kφAφB 〈SR〉 = −k〈φAφB〉 = −k〈φA〉〈φB〉 The chemical source term can be closed if we assume the existence of a joint probability density function of the concentrations 〈S(φ)〉 =

S(ψ)n(ψ;x,t)dψ

Part II: Quadrature-based moment methods

MULTIPHASE FLOWS: MOMENTUM TRANSFER

How to model the drag force Fd acting on fluid-fluid or fluid-solid suspensions when shape, position, velocity and size of interphase vary in time and space ? Given a configuration space (probability space), statistics (mean, variance, ...) can be computed

Part II: Multi-level Monte Carlo

4 / 34

SLIDE 5

Multiscale models

Microscale model

Fully resolved direct numerical simulation

Mesoscale model

Population balance equation Boltzmann equation Kinetic equation (Lagrangian methods)

Macroscale model

Moment equations Mixture model Two-fluid model (Eulerian methods) Density function closure Moment transform Moment closure Volume average Reconstruction

5 / 34

SLIDE 6

PART I: Quadrature-based Moment Methods

5 / 34

SLIDE 7

Method of Moments

Important advantage is that the moments correspond to macroscopic

quantities that have meaningful physical interpretations and are therefore directly measurable

In many applications the NDF is not directly measured, but is inferred from

measurements of integral quantities

Two main issues arise when using MOM: number of moments to be tracked

and closure problem

The closure problem is the impossibility of writing a term as a function of a

finite set of moments

When the PDF equation contains transport in the phase space, we have a

cascade of equations (Mi depends on Mi+1)

All non-polynomial integral source terms also need closure

6 / 34

SLIDE 8

Closure problem

When applying the moment transform (i.e. integral)

THE CLOSURE PROBLEMS APPEAR IN THE FOLLOWING FORM:

I =

Ωξ

g(ξ)n(ξ)dξ, (1) where n(ξ) is the NDF to be approximated and g a generic function to be closed

PRESUMED-PDF APPROACHES

A possible approach is to assume a functional form for the PDF/NDF This is reasonable when there is a fast relaxation towards a known

equilibrium

How to find the best closure for a general case?

7 / 34

SLIDE 9

Quadrature-based Moment Method

The closure problem can be overcome by using the following quadrature

formula:

Ωξ

n(ξ)g(ξ)dξ ≈

N

α=1

wαg(ξα), (2) where wα and ξα are the weights and the nodes/abscissas of the quadrature formula, and N is the number of nodes

GAUSSIAN QUADRATURE IS THE “BEST” CHOICE

The NDF is the weight function or measure; the moments Mk = 〈ξk〉 =

Ωξ

n(ξ)ξk dξ, k = 0,1,2,... (3) are used to compute a Gaussian quadrature rule with a degree of accuracy of 2N −1,.

8 / 34

SLIDE 10

Quadrature-based moment methods

For univariate problems things are simple: transport equations for the first

2N moments are solved, the closure problem is overcome by using a quadrature approx.

The quadrature is very accurate in approximating integrals of smooth

functions in which the NDF appears (such as high-order moments and other integrals)

The method implicitly assumes that the NDF is a summation of Dirac delta

functions centered on the nodes and weighted by the weights of the quadrature approximation

The quadrature approximation of order N is calculated from the first 2N

moments with direct inversion algorithms: product-difference or Wheeler algorithms

9 / 34

SLIDE 11

Quadrature method of moments

Let us consider the following example (pure growth process):

Dn Dt = S = − ∂ ∂ξ

Gpn
,

(4)

Moment equations

DMk Dt = k ∞ 〈Gp|ξ〉ξk−1n(ξ)dξ (5)

The closure problem is overcome as follows:

DMk Dt = k

N

α=1

〈Gp|ξα〉(ξα)k−1wα (6)

weights wα and nodes ξα are calculated from the PD or Wheeler algorithm

from the moments

This method is called Quadrature Method of Moments (the variables are the

moments)

10 / 34

SLIDE 12

Example: Boltzmann equation

Particle Velocity Distribution f (Up;x,t) =

R

n(L,Up;x,t)dL

MODEL EQUATION

scalar molecular speed ξ distribution of velocities f (ξ) no spatial derivatives

nly collisions Q(f ,f )

Hard-sphere model x = cos(ξ1∠n),

y = cos(ξ2∠n)

collision frequency |ξ2y −ξ1x|

11 / 34

SLIDE 13

Homogeneous Isotropic Boltzmann equation

HIBE evolution equation: df dt = Q(f ,f )

Q = 2π2 a2 +∞ ξ2

2

+1

−1

+1

−1

f (ξ′

1)f (ξ′ 2)−f (ξ1)f (ξ2)

|ξ2y −ξ1x|dxdydξ2

Let us introduce the even moments M2p = Φp (density Φ0, energy Φ1, etc.). After change of variable E = ξ2

2 and a few manipulations:

dΦp dt = 4π

2

∞ Q(f ,f )Ep+1/2dE =

16π3 2 ∞ ∞ +1

−1

+1

−1

|q|

(C+

p )−(C− p )

f (E)f (E∗)(EE∗)1/2dxdydE∗dE

q = yE1/2

∗

−xE1/2, C+

p =

E(1−x2)+E∗y2p

, C−

p = Ep

12 / 34

SLIDE 14

Quadrature approximation for HIBE

COLLISIONAL INTEGRAL

∞ ∞ +1

−1

+1

−1

|q|

(C+

p )−(C− p )

f (E)f (E∗)(EE∗)1/2dxdydE∗dE

M/2

i=0

M/2

j=0

+1

−1

+1

−1

|q|

(C+

p )−(C− p )

wiwjdxdydE∗dE

Written in terms of pre-collisional quantities only → quadrature approximation ˜ f (E) ≈

M/2

i=1

wi 4π

2E

δ(E −Ei) ⇒ dΦp dt ≈ π

2

M/2

i=1

M/2

j=1

wiwjΛij,p at equilibrium Q = 0 → dΦp

dt = 0 BUT this is not guaranteed with the quadrature

approximation Collisions drive the flows towards the wrong steady state

13 / 34

SLIDE 15

Quadrature-Based Moment Method (QBMM)

dΦp dt = π

2

M/2

i=1

M/2

j=1

wiwjΛij,p = Sp wi and Ei calculated from Φp with inversion algorithms1

QMOM fails to approach the equilibrium → source term must be corrected Static correction (QMOM+SC):

dΦp dt = Sp −Seq

p

Dynamic correction (QMOM+DC) weighted with a distance from

equilibrium: dΦp dt = Sp −

Φp(t)−Φp(0)

Φeq

p (t)−Φp(0)

h

Seq

p ,

h = 1.5

1Product-Difference or Wheeler algorithms

14 / 34

SLIDE 16

QMOM approximation vs DVM, Grad, oLBM

Let us consider a closed box of particles far from equilibrium (initial velocity distribution NOT Gaussian): QMOM approximation selecting M moments (M/2 nodes and weights)

COMPARISON WITH:

Discrete Velocities Method (DVM): reference results with 400 discrete

velocities - HomIsBoltz open-source Matlab code (Asinari, 2010)

Grad expansion method (GM) of order M with generalized Laguerre

polynomials

Quadrature approximation with M fixed Laguerre nodes → Lattice

Boltzmann Method (LBM) with off-lattice prescribed velocities

15 / 34

SLIDE 17

QMOM approximation vs DVM, Grad, oLBM – 1

M = 4 Initial condition

16 / 34

SLIDE 18

QMOM approximation vs DVM, Grad, oLBM – 2

M = 4 Relaxation to equilibrium Rp =

Φp−Φeq

p

Φeq

p

2nd energy moment 3rd energy moment

17 / 34

SLIDE 19

QMOM approximation vs DVM, Grad, oLBM – 3

M = 4 Relative error on Φp 2nd energy moment 3rd energy moment

18 / 34

SLIDE 20

HIBE steady state

Comparison of steady-state moments

M = 4 M = 6 M = 8 ×108 Φ2 % Φ2 % Φ2 % Grad (exact) 1.0034 0.0 1.0034 0.0 1.0034 0.0

LBM

1.0033 0.01 1.0034 0.003 1.0034 0.001 QMOM 0.9704 3.29 0.9979 0.55 1.0022 0.12 QMOM+SC/DC 1.0034 0.0 1.0034 0.0 1.0034 0.0 ×1011 Φ3 % Φ3 % Φ3 % Grad (exact) 1.4921 0.0 1.4921 0.0 1.4921 0.0

LBM

1.4919 0.02 1.4921 0.004 1.4921 0.002 QMOM 1.4707 1.44 1.4874 0.31 1.4910 0.08 QMOM+SC/DC 1.4921 0.0 1.4921 0.0 1.4921 0.0 ×1014 Φ4 % Φ4 % Grad (exact)

2.8529

0.0 2.8529 0.0

LBM
2.8527

0.005 2.8528 0.002 QMOM

2.8302

0.80 2.8478 0.18 QMOM+SC/DC

2.8529

0.0 2.8529 0.0 ×1017 Φ5 % Φ5 % Grad (exact)

6.6666

0.0 6.6666 0.0

LBM
6.6662

0.006 6.6665 0.002 QMOM

6.6244

6.34 6.6548 0.18 QMOM+SC/DC

6.6666

0.0 6.6666 0.0

19 / 34

SLIDE 21

HIBE steady state

Quadrature error at the true equilibrium

20 / 34

SLIDE 22

QBMM - multivariate case

The Generalized Population Balance Equation for M-dimensional phase space: ∂n ∂t + ∂ ∂x ·

Upn
= ∂

∂x ·

D∂n

∂x

− ∂

∂ξ ·(Gn)+Q, (7) After having defined a generic moments : Mk ≡

Ωξ

ξk1

1 ···ξkM M n(t,ξ)dξ

(8) equivalent notations Mk1,...,kM = Mk = M(k1,...,kM) = M(k), we solve the resulting transport equations: ∂Mk ∂t + ∂ ∂x ·

Uk

pMk

− ∂

∂x ·

Dk ∂Mk

∂x

=
Ωξ

ξk

− ∂

∂ξ ·

〈Gp|ξ〉n
+Q
dξ

(9) and the closure problem is overcome by using the quadrature approximation.

21 / 34

SLIDE 23

Multivariate case

N(M +1) moments ( optimal moment set) are used to determine the quadrature approximation Brute-force: direct solution of the following non-linear system (ill-conditioned) Mki1,ki2,...,kiM = M(ki) =

N

α=1

wα

M

β=1

ξ

kiβ β,α,

1 ≤ i ≤ N(M +1) by employing the Newton-Raphson iterative scheme Alternatively the Tensor-Product QMOM can be used using marginals or the Conditional QMOM introducing an order in the variables and using conditional moments

22 / 34

SLIDE 24

QBMM - multivariate case

Bivariate Gaussian distr. with M = 9 Direct (Newton) (diamond), Tensor product (circle) and Conditional (square)

23 / 34

SLIDE 25

Conclusions

The idea of QBMM is similar to the “arbitrary Polynomial Chaos Expansion” (aPC, Oladyshkin and Nowak, 2012)

Input random variables ξi known only by moments Moments induce orthogonal polynomials and quadrature rule Stochastic spectral and collocation approaches for computing the response

surface

Multidimensional correlated variables treated directly Forward propagation of uncertainty from moments of input to moments of

response, without passing through the PDF Good alternative for correlated variables in moderate dimensions

24 / 34

SLIDE 26

Future works

The main advantage in using only moments is the possibility of adaptively updating the quadrature rule, when the underlying PDF changes in time/space

NON-LINEAR FOKKER-PLANCK

Attar and Vedula, 2008; Otten and Vedula, 2011 Drift-diffusion equations

˙ xi = hi(x,t)+gi(x,t)Wi, Wij

Applications to mixing in porous media

NON-LINEAR FILTERING AND BAYESIAN UPDATES

Xu and Vedula, 2009, 2010 Propagation step through Fokker-Planck Bayesian update based on quadrature or on EnKF Mixture of gaussians that preserves moments

25 / 34

SLIDE 27

PART II: MLMC for dispersed multiphase flows

25 / 34

SLIDE 28

Suspensions/dispersed flows

DNS of fluid-solid interface (particle-resolved simulations) are often unfeasible

WHY

Limits in resolution Physical models too complex

(e.g. particulate processes)

We want a steady solution or

2D/1D models

Unknown particle

shape/properties

EXAMPLES

Fluidized beds, separators Environmental flows Pore-scale processes Long piping systems

We need effective closures for momentum transfer (e.g. drag force correlations) that could lump together some of the complexities of the system.

26 / 34

SLIDE 29

Drag correlations

Many different closures has been proposed:

Isolated sphere: Stokes, Schiller-Neumann, Clift... High density/porous media: Ergun (Kozeny Carman), ... Periodic arrangements: Sangani and Acrivos, ... Granular flows: Gidaspow, Beetstra et al (2007), Tenneti et al (2011)

We generally seek a relation Fd = Fd(Re,φ) where Fd is the normalized equivalent mean drag force over a particle ensemble and φ volume fraction. Many systems (e.g., fluidized beds) suffer from high sensitivity to the drag law

27 / 34

SLIDE 30

Stochastic interpretation

Is the problem deterministic or stochastic?

Two-fluid (or mixture) models are derived with volume/ensemble averaging Deterministic averaging/homogenization limit exist for drag force in high

density systems ∆p ∝ Fd ∝ µ K 〈U〉 (Darcy Eq with constant permeability)

The ensemble/volume size however tends to infinity for low porosity Stochastic homogenization: the actual averaging 〈·〉 is not deterministic

1. What is the effect of particle random position and velocity on the drag?
2. What is the role/magnitude of drag force fluctuations over (and inside) the

averaging volume? We seek Fd = Fd(Re,φ,θ) and F′

d = F′ d(Re,φ,θ)

θ is a “granular temperature” and F′

d are the random drag fluctuations

28 / 34

SLIDE 31

Statistical estimator

Even if we are interested only in the average drag on a fixed ensemble/volume size, the actual computations requires thousands 3D expensive CFD realizations

MLMC ALGORITHM

1. Each estimator, defined with prescribed tolerance and confidence, is

initialized with L levels

2. Each level is initialized with Mℓ random realizations (samples)
3. For each sample, generate the geometry (packing algorithms), discretize

(meshing), solve (finite volume) and extract the QoI (integrate) at level ℓ and level ℓ−1

4. Re-assemble the multilevel estimator. If numerical or statistical error above

threshold add respectively levels or samples and go back to 3

5. Estimate independently more times to double check accuracy criteria

ONGOING EXTENSIONS

Combination of grid-based multilevel estimator with other discretization

parameters (domain size, tolerances, solvers) in the MIMC spirit

High-order extrapolation of mean and variances using weighted sums of

multiple levels

29 / 34

SLIDE 32

Problem formulation

We developed fully automatic code to solve different stochastic PDE problems and compute effective properties:

ADIMENSIONAL STEADY INCOMPRESSIBLE NAVIER-STOKES

u(x;ω)·∇u(x;ω)+∇p(x;ω)− 1 Re ∇·(∇u(x;ω)) = 0, x ∈ Γ(ω) ⊂ [0,1]d, ∇·u = 0 no-slip: u = 0, ∇np = 0

n

Ig inlet: ∇nu = 0, p = pi(x)

n

Ii = {x : x1 = 0}

utlet:

∇nu = 0, p = po(x) = 0

n

Io = {x : x1 = 1} symmetry: ∇nut = 0, un = 0, ∇np = 0

n

Is = {x : x2,x3 = 0} The quantity of interest Q(u) is, in this case, the mean flux 〈ux〉 or, equivalently, the drag force Fd = ∆p(1−φ)d2

18〈ux〉

The algorithm is tested for creeping flow through random array of mono dispersed

spheres. Extended to non-spherical, polydispersed, and moving2 particles.

2in local equilibrium

30 / 34

SLIDE 33

Example: mean flux on mono dispersed spheres (φ = 5%)

Unstructured anisotropic mesh with non-uniform localized refinements Computational savings: 158x faster than MC, (variance 132x) Mean: 0.113 Variance: 0.0136 CoV: 103% Error estimation: Numerical bias 0.4% (0.6%), Statistical 1.4% (3.1%)

31 / 34

SLIDE 34

Example: mean flux on mono dispersed spheres (φ = 20%)

Unstructured anisotropic mesh with non-uniform localized refinements Computational savings: 131x faster than MC, (variance 129x) Mean: 0.0124 Variance: 0.000165 CoV: 33% Error estimation: Numerical bias 0.5% (0.9%), Statistical 1.4% (3.0%)

32 / 34

SLIDE 35

Why does the weak error stagnate?

This can be explained by the non-uniform refinement and the non-monotone non-asymptotic behavior total discretization error Deterministic convergence rate for simple regular packings Convergence rate for pure diffusion (left) and Navier-Stokes with porosity 0.47 (center) and 0.366 (right). Voxelized mesh, Locally refined, Uniform

33 / 34

SLIDE 36

Conclusions

ACHIEVEMENTS

Developed automatic (and parallel) workflow for setup, meshing, simulation

and post-processing for flow through random granular materials

Assessed accuracy (convergence rates) and computed uncertainty

(variations) for low-Re systems

A-posteriori error control: both statistical and numerical Preliminary results for mono disperse spheres in Stokes regime

OPEN PROBLEMS AND ONGOING WORK – High numerical accuracy (<1%) still unobtainable for large samples

For new drag laws, huge amount of data for different Re and φ needed Influence and representativity of the packing algorithm Limit the effect of the boundary conditions (periodicity and lift estimation) How to use high-order moments computed by MLMC in Eulerian models?

34 / 34

SLIDE 37

References (for Part I)

D. Marchisio and R. Fox, “Computational Models for Polydisperse Particulate

and Multiphase Systems” (Cambridge, 2012)

Rodney Fox, “Computational Models for turbulent reacting flows”

(Cambridge, 2003)

H. Struchtrup, “Macroscopic Transport Equations for Rarefied Gas Flows”

(Springer, 2005)

Cercignani, “The Boltzmann equation and its applications” (Springer, 1988)

W. Gautschi, “Orthoghonal Polynomials” (Oxford Sci., 2004)
D. Ramkrishna, “Population Balances” (Academic Press, 2000)

my PhD thesis (2012)

O. Le Maitre and O. Knio, “Spectral Methods for Uncertainty Quantification”

(Springer, 2010)

POSTERS: Flow, transport and diffusion in random geometries I and II

34 / 34

SLIDE 38

BACKUP SLIDES

34 / 34

SLIDE 39

Smoluchowski coagulation equation

Describe particle aggregation/coagulation ∂n(L,t) ∂t = 1 2 L K(L− ˜ L, ˜ L)n(L− ˜ L,t)n(˜ L,t)d˜ L− ∞ K(L, ˜ L)n(L,t)n(˜ L,t)d˜ L

Not to be confused with the Smoluchowski equation (drift-diffusion

equation)

Integro-differential equation Number density n of particles with size L No spatial dependence RHS (Q) has a gain and a loss term written in terms of the aggregation

frequency K

K is usually very complex and non-linear

Very interesting mathematical topic (stability, linearization, large-time asymptotics), see works of Klemens Fellner (Graz)

35 / 34

SLIDE 40

Population Balance Equation (PBE)

Extension of the Smoluchowski coagulation equation, very popular in Chemical Engineering, Biology and Social Models ∂n ∂t + ∂ ∂x ·(Un)+ ∂ ∂ξ ·(Gn) = Q

n = n(ξ;x,t) number of particles per unit volume ξi generic internal variables (volume, size, area, composition, other

properties)

U given advection velocity field Q describe the generic particulate process with Birth (nucleation), Breakage

(Fragmentation) and Aggregation (Coagulation)3 Q = Q0 +Q1(n)+Q2(n,n)

G describe growth of particles (velocity in phase space ξ) Usually written for spatially inhomogeneous systems and coupled with CFD

to simulate dispersed bubbles and particles

3with frequency kernels that can be seen as upscaled jump processes

36 / 34

SLIDE 41

The physical basis of Boltzmann equation

Molecular encounters for hard-spheres potential

d is the molecular diameter g = U−U∗ is the pre-collision velocity

difference

g′ = U′ −U′

∗ is the post-collision velocity

difference Number and momentum are conserved during a collision m+m∗ = m′ +m′

∗;

U+U∗ = U′ +U′

∗

as well as kinetic energy (in the case of elastic collisions) U·U+U∗ ·U∗ = U′ ·U′ +U′

∗ ·U′ ∗

37 / 34

SLIDE 42

Boltzmann Equation (BE)

From the Liouville n-particles equation, assuming indistinguishable particle and the "Stosszahlansatz" (molecular chaos4 hypothesis), the Boltzmann equation is

btained

∂n(U) ∂t + ∂ ∂x ·(Un(U))+ ∂ ∂U ·(An(U)) =

R3
S+
n(U′)n(U′

∗)−n(U)n(U∗)

β
g,x
dsdU∗

A is the acceleration, S+ is the solid angle β

g,x
is the collision kernel

Equations with similar structure: Vlasov (no collisions), Williams-Boltzmann spray equation, ...

4particles are uncorrelated, two-particles PDF is product of one-particle PDF

38 / 34

SLIDE 43

Boltzmann Equation and Equilibrium

Most of the studies and methods for BE are based on the concept of Equilibrium Distribution

The equilibrium distribution feq is such that Q(feq,feq) = 0 Also called "Maxwell-Boltzmann" or "Maxwellian" distribution, a Gaussian in

the standard case

Fixing the first 2 moments (3 in case of NDF), only one Gaussian exist If one wants to study non-equilibrium phenomena, more moments must be

studied

Grad’s moment method is based on small deviations from equilibrium

f ≈ feq(1+a1H1(v)+a2H2(v)+...)

Other methods (such as Lattice Boltzmann) rely on a linearized collision term

(BGK, Bhatnagar-Gross-Krook approximation) Q ≈ 1 τr (f −feq) All these methods require the prior knowledge of the equilibrium distribution!

39 / 34

SLIDE 44

The Method of Moments

MOM originates in kinetic theory of gases5: n(t,x,U) The moment of order zero is related to the density:

ρ(t,x) ≡ mM0,0,0(t,x) = m +∞

−∞

n(t,x,U)dU

The moment of order one is used to define the average velocity:

〈U〉 = M1,0,0 M0,0,0 , M0,1,0 M0,0,0 , M0,0,1 M0,0,0

By applying the MOM to the Boltzmann equation, an infinite system of

equations (cascade) appears

Euler, Navier-Stokes and Burnett equations are obtained by with the

Chapman-Enskog method (asymptotic expansions at different orders)

5The same name however may refer to many different methods in different fields

40 / 34

SLIDE 45

Flow regimes

The Boltzmann Equation is not only for rarefied gases!

The flow regime of a granular gas depends on Knudsen number Kn = λ

L0

The hydrodynamic description based on the Navier-Stokes-Fourier (NSF)

equation is valid only for low Kn:

Continuous regime (Kn < 0.01): NSF with no-slip BC. Slip regime (0.01 < Kn < 0.1): NSF and partial slip BC at walls. For Kn > 0.1: Full Boltzmann equation.

Kn

Euler N-S Boltzmann equation Collisionless Boltzmann eq. ∞ 0.01 0.1 1 100 10

Adapted from: G. A. Bird (1994)

41 / 34

SLIDE 46

The GPBE fluid-particle systems

n(t,x,Up,ξp,Uf,ξf), represents the number of particles (per unit volume)

with velocity equal to Up, internal coordinate ξp and see a continuous phase with velocity Uf and internal coordinate ξf.

The evolution of the NDF is dictated by the GPBE:

∂n ∂t + ∂ ∂x ·

Upn
+

∂ ∂Up ·

Afp +Ap
n
+

∂ ∂ξp ·

Gpn
+

∂ ∂Uf ·

Apf +Af
n
+ ∂

∂ξf ·

Gfn
= Q

(10)

phase-space velocity for particle velocity: Afp +Ap ( average particle

acceleration → acceleration model (drag, lift, ...)

phase-space velocity for particle internal coordinate: Gp (e.g. particle growth

rate)

Discontinuous jump term: Q (e.g. particle collision, aggregation, breakage,

nucleation, etc.)

42 / 34

SLIDE 47

History of quadrature-based moment methods

McGraw (1997) introduced the Quadrature Method Of Moments (QMOM) in

the context of Population Balance Equation

Marchisio and Fox (2005) developed the Direct Quadrature Method of

Moments (DQMOM) that became very popular in the Chemical Engineering community because very easy to implement in CFD codes

Different extensions to multivariate cases have been proposed. The more

general one is the Conditional QMOM (CQMOM) (Yuan and Fox, 2011)

Kernel density type reconstruction have been incorporated in QMOM with

the Extended QMOM (EQMOM) (Yuan, Laurent, Fox, 2012)

Different applications have been studied: particulate processes, dispersed

flows (gas-liquid, gas-solid, fluid-solid), turbulent micro-mixing, rarefied gases, generic kinetic/Fokker-Planck equations, non-linear filtering, ...

Recently proposed for Uncertainty Quantification (Passalacqua and Fox,

2012; Attar and Vedula, 2012)

A large number of variants have been proposed (FCMOM, TBMM,

DuQMoGeM, SQMOM , ...) All the closures based on Gaussian quadratures computed directly from moments took the name of QUADRATURE-BASED MOMENT METHODS (QBMM)

43 / 34

SLIDE 48

Quadrature-based moment methods

HOW DO WE COMPUTE THE GAUSSIAN QUADRATURE?

The N weights and N abscissas can be determined by solving the following non-linear system: M0 =

N

α=1

wα, . . . M2N−1 =

N

α=1

wαξ2N−1

α

. (11) using the Newton-Raphson method, or any other non-linear equation solver (very expensive, ill-posed and very good initial guess needed) Much more efficient are the product-difference or Wheeler algorithms

44 / 34

SLIDE 49

Gaussian quadrature

A set of polynomials {P0(ξ),P1(ξ),...,Pα(ξ),...} with Pα(ξ) = kα,0xα +kα,1xα−1 +···+kα,α, is said to be orthogonal in the integration interval Ωξ, with respect to the weight function, if

Ωξ

n(ξ)Pα(ξ)Pβ(ξ)dξ

= 0

for α = β, > 0 for α = β, (12) and, of course, is said to be orthonormal if

Ωξ

n(ξ)Pα(ξ)Pβ(ξ)dξ =

for α = β,

1 for α = β. (13)

45 / 34

SLIDE 50

Gaussian quadrature

Any set of orthogonal polynomials {Pα(ξ)} has a recurrence formula relating any three consecutive polynomials in the following sequence: Pα+1(ξ) = (ξ−aα)Pα(ξ)−bαPα−1(ξ), α = 0,1,2,... (14) with P−1(ξ) ≡ 0 and P0(ξ) ≡ 1 and where aα =

Ωξ n(ξ)ξPα(ξ)Pα(ξ)dξ
Ωξ n(ξ)Pα(ξ)Pα(ξ)dξ ,

α = 0,1,2,... (15) bα =

Ωξ n(ξ)Pα(ξ)Pα(ξ)dξ
Ωξ n(ξ)Pα−1(ξ)Pα−1(ξ)dξ,

α = 1,2,.... (16) One can calculate a0, then P1(ξ), then a1 and b1 and so on...

46 / 34

SLIDE 51

Gaussian quadrature

The coefficients aα and bα can be written in terms of the moments The coefficients necessary for the construction of a polynomial of order N

can be calculated from the first 2N −1 moments of the NDF

For example with M0, M1, M2 and M3, it is possible to calculate the following

coefficients: a0 = M1 M0 , a1 = M3M2

0 +M3 1 −2M2M1M0

M2M0 +M2

1 −2M2 1M0

, b1 = M2M0 +M2

1 −2M2 1M0

M2 , (17) which suffice for the calculation of the polynomial P2(ξ)

47 / 34

SLIDE 52

Gaussian quadrature

GAUSSIAN QUADRATURE

The necessary and sufficient condition for the following formula:

Ωξ

n(ξ)g(ξ)dξ =

N

α=1

g(ξα)wα +RN(g), (18) to be a Gaussian quadrature approximation is that its nodes {ξα} coincide with the N roots of the polynomial PN(ξ) of order N orthogonal in Ωξ with respect to the weight function n(ξ).

48 / 34

SLIDE 53

Direct quadrature method of moments

The fact that the closure problem is overcome with the quadrature

approximation:

Ωξ

n(ξ)g(ξ)dξ ≈

N

α=1

wαg(ξα), (19)

Is equivalent to the assumption that the NDF is as follows:

n(ξ) =

N

α=1

wαδ(ξ−ξα), (20)

Instead of tracking the evolution for the moments, the evolution of the

weights and nodes in the quadrature approximation could be directly tracked: Direct quadrature method of moments

49 / 34

SLIDE 54

Direct quadrature method of moments

Assuming that the weights and nodes are differentiable in space/time the

following transport equation is obtained:

N

α=1

δ(ξ−ξα) Dwα Dt

−

N

α=1

δ′(ξ−ξα)

wα

Dξα Dt

= S(ξ),

(21)

If the weighted nodes (or weighted abscissas) ςα = wαξα are introduced:

N

α=1

δ(ξ−ξα) Dwα Dt

−

N

α=1

δ′(ξ−ξα)

−ξα

Dwα Dt + Dςα Dt

= S(ξ).

(22)

50 / 34

SLIDE 55

Direct quadrature method of moments

We now define aα and bα to be the source terms:

Dwα Dt = aα, Dςα Dt = bα. (23)

Using these definitions Eq. (22) can be rewritten in a simpler form:

N

α=1
δ(ξ−ξα)+δ′(ξ−ξα)ξα
aα −

N

α=1

δ′(ξ−ξα)bα = S(ξ). (24)

This equation can now be used to determine the unknown functions aα and

bα by applying the moment transformation.

51 / 34

SLIDE 56

Direct quadrature method of moments

DQMOM can be applied for any independent set of moments (number of

moments MUST be equal the number of unknown functions)

Knowing that:

+∞

−∞

ξkδ(ξ−ξα)dξ = (ξα)k, +∞

−∞

ξkδ′(ξ−ξα)dξ = −k(ξα)k−1, (25)

The moment transform of Eq. (24) yields

(1−k)

N

α=1

ξk

αaα +k N

α=1

ξk−1

α

bα = Sk, (26) with k = k1,k2,...,k2N.

52 / 34

SLIDE 57

Direct quadrature method of moments

The linear system in Eq. (26) can be written in matrix form:

Aα = d. (27) where α =

a1

··· aN b1 ··· bN T = a b

,

(28) d =

Sk1

··· Sk2N−1 T , (29)

The components of the matrix A are

aij =   

1−ki
ξki

j

if 1 ≤ j ≤ N, kiξki−1

j

if N +1 ≤ j ≤ 2N. (30)

53 / 34

SLIDE 58

Direct quadrature method of moments

If (as in QMOM) the first 2N integer moments are chosen (i.e., k = 0,...,2N −1), the

matrix of the linear system is A =         1 ··· 1 ··· ··· 1 ··· 1 −ξ2

1

··· −ξ2

N

2ξ1 ··· 2ξN . . . . . . . . . . . . 2(1−N)ξ2N−1

1

··· 2(1−N)ξ2N−1

N

(2N −1)ξ2N−2

1

··· (2N −1)ξ2N−2

N

        . (31)

A does not depend on the weights wα and if the abscissas ξα are unique, then A will be

full rank.

54 / 34

SLIDE 59

Direct quadrature method of moments

This method is called Direct quadrature method of moments and follows this

procedure:

The evolution equations for weights and nodes of the quadrature

approximation are solved: Dwα Dt = aα, Dwαξα Dt = bα. (32)

The source terms are calculated by inverting the linear system and by using

the following initial condition: wα(0) = w0

α,

ςα(0) = w0

αξ0 α

for k = 1,...,N. (33) in turn calculated from the initial moments

55 / 34

SLIDE 60

QMOM & DQMOM

QMOM & DQMOM are very accurate in tracking the evolution of the

moments of the NDF: 4-8 moments do the same job of many (e.g. 100) classes

r sections (see for example the work of Marchisio et al., 2003 and Vanni,2000)

The Wheeler algorithm is very robust (if the moments are realizable) and for

particular cases the Wheeler algorithm is successful when PD fails

QMOM & DQMOM are identical for spatially homogeneous systems (if the

nodes are distinct and if the problem is continuous in time)

Important differences arise when treating spatially inhomogeneous systems

(discussed next)

In general increasing the number of nodes of the quadrature approximation

and of moments to be tracked increases the accuracy

Problems can appear when kernels and NDFs are discontinuous or when

they are localized in the phase-space (e.g. fine dissolution)

56 / 34

SLIDE 61

Product-difference algorithm

A smarter way is to employ the recursive relationship for the orthogonal polynomials: ξ         P0(ξ) P1(ξ) . . . PN−2(ξ) PN−1(ξ)         =          a0 1 b1 a1 1 ... ... 1 aN−1                  P0(ξ) P1(ξ) . . . PN−2(ξ) PN−1(ξ)         +         . . . PN(ξ)         . (34) The nodes of the quadrature approximation {ξα}, are the eigenvalues of the tridiagonal matrix appearing in the equation. This matrix is often re-written in terms of an equivalent tridiagonal symmetric matrix!

57 / 34

SLIDE 62

Product-difference algorithm

In fact the matrix can be made symmetric (preserving the eigenvalues) by a diagonal similarity transformation to give a Jacobi matrix: J =                  a0

b1
b1

a1

b2
b2

a2

b3
b3

a3 ... ... ... ... ... aN−2

bN−1
bN−1

aN−1                  (35) This procedure transforms the ill-conditioned problem of finding the roots of a polynomial into the well-conditioned problem of finding the eigenvalues and eigenvectors of a tridiagonal symmetric matrix. The N weights can then be calculated as wα = M0ϕ2

α1 where ϕα1 is the first component of

the αth eigenvector ϕα of the Jacobi matrix.

58 / 34

SLIDE 63

Product-difference algorithm

1. Construct the matrix P with components Pα,β:

Pα,β = P1,β−1Pα+1,β−2 −P1,β−2Pα+1,β−1 β ∈ 3,...,2N +1 and α ∈ 1,...,2N +2−j. (36)

2. where the first row of the matrix is:

Pα,1 = δα1 α ∈ 1,...,2N +1, (37)

3. where δα1 is the Kronecker delta and where the components in the second

column of P are Pα,2 = (−1)α−1Mα−1 α ∈ 1,...,2N. (38)

4. Calculate the coefficients of the continued fraction {ζα}:

ζα = P1,α+1 P1,αP1,α−1 α ∈ 2,...,2N. (39)

59 / 34

SLIDE 64

Product-difference algorithm

1. The coefficients of the symmetric tridiagonal Jacobi matrix are then obtained

from sums and products of ζα: aα = ζ2α +ζ2α−1 α ∈ 1,...,N (40) bα = −

ζ2α+1ζ2α

α ∈ 1,...,N −1. (41)

2. For example for N = 2 the P matrix is

       1 M0 M1 M0M2 −(M1)2 M0

M3M1 −(M2)2

−M1 −M2 −(M0M3 −M2M1) M2 M3 −M3        . (42)

60 / 34

SLIDE 65

Extended QMOM

nly integral representation, no localized processes in phase space,

smoothness of kernels β, G,...

It fails in the computation of the entropy

Ωξ f logf dξ

However smooth basis functions can be used instead of delta functions

(Extended Quadrature Method of Moments, EQMOM)

PDF reconstructed as a mixture of Gaussians (or other distributions) like the

Kernel Density Estimation

The choice of basis functions can be done arbitrarily or with additional

unknown parameters (e.g., mean and variance of normal distribution)

Additional parameters must be found by additional equations (i.e., more

moments must be solved)

The inversion problems can become difficult to solve

An alternative is the reconstruction through Maximum Entropy Method or the reconstruction through the orthogonal polynomials

61 / 34

SLIDE 66

Extended Quadrature-based moment methods

EQMOM for a dissolution problem GAUSSIAN DISTR. SOLUTION WITH N = 4 BETA DISTR.

−2 −1 1 2 3 4 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

ξ n(ξ)

0.125 0.25 0.5 1

a

aYuan, C., Laurent, F

., Fox, R.O. An extended quadrature method of moments for population balance equations (2012) Journal of Aerosol Science, 51, pp. 1-23.

δσ = 1

2πσ

exp

− (ξ−ξα)2

2σ2

;δσ = ξ

ξα σ −1 (1−ξ) (1−ξα) σ −1

B(ξα,σ)

62 / 34

SLIDE 67

Optimal moment set

OPTIMAL MOMENT SET

1. An optimal moment set consists of N(M +1) distinct moments.
2. An optimal moment set will result in a full-rank square matrix A for all

possible sets of N distinct, non-degenerate abscissas.

3. An optimal moment set includes all linearly independent moments of a

particular order γi before adding moments of higher order.

4. An optimal moment set must result in a perfectly symmetric treatment of the

internal coordinates.

63 / 34

SLIDE 68

Optimal moment set

Table: Moments used to build a bivariate quadrature approximation (M = 2) for N = 2. In this case M0,3 is chosen as the third-order moment to saturate the degrees of freedom.

M(2,0) M(1,0) M(0,0) M(0,1) M(0,2) M(0,3)

Table: Moments used to build a bivariate quadrature approximation (M = 2) for N = 3. In this case M2,1, M1,2 and M0,3 are chosen among the third-order moments to saturate the degrees

f freedom.

M(2,0) M(2,1) M(1,0) M(1,1) M(1,2) M(0,0) M(0,1) M(0,2) M(0,3)

64 / 34

SLIDE 69

Optimal moment set

Table: Optimal moment set used to build a bivariate quadrature approximation (M = 2) for N = 4. Only when N1/M is an integer, there exists an optimal moment set (that fulfills the symmetry requirement). M(3,0) M(3,1) M(2,0) M(2,1) M(1,0) M(1,1) M(1,2) M(1,3) M(0,0) M(0,1) M(0,2) M(0,3) Table: Optimal moment set used to build a bivariate quadrature approximation (M = 2) for N = 9. M(5,0) M(5,1) M(5,2) M(4,0) M(4,1) M(4,2) M(3,0) M(3,1) M(3,2) M(2,0) M(2,1) M(2,2) M(2,3) M(2,4) M(2,5) M(1,0) M(1,1) M(1,2) M(1,3) M(1,4) M(1,5) M(0,0) M(0,1) M(0,2) M(0,3) M(0,4) M(0,5)

65 / 34

SLIDE 70

Multivariate Conditional QMOM

Methods based on conditional density functions:

n(ξ1,ξ2) = n1(ξ1)f21(ξ2|ξ1) = n2(ξ2)f12(ξ1|ξ2)

Univariate quadrature (N1) calculated from the first 2N1 −1:

   M0,0,...,0,0 . . . M2N1−1,0,...,0,0    PD/Wheeler − − − − − − − − − →    w1 . . . wN1       ξ1;1 . . . ξ1;N1   . resulting for example in: n(ξ1,ξ2) = N1

α1=1 wα1δ

ξ1 −ξ1;α1
f21(ξ2|ξ1)

The generic moment becomes:

Mk1,k2 =

n(ξ1,ξ2)ξk1

1 ξk2 2 dξ1 dξ2

=

N1

α1=1

wα1ξk1

1;α1

f21(ξ2|ξ1;α)ξk2

2 dξ2

Conditional moment:

ξk2

2

α1

=

f12(ξ2|ξ1;α1)ξk2

2 dξ2

66 / 34

SLIDE 71

Multivariate Conditional QMOM

For each of these N1 nodes, 2N2 −1 conditional moments are calculated, and univariate quadratures N2 are determined (in direction ξ2): Conditional QMOM

r CQMOM

1 ; 1

ξ ξ ξ ξ

1

ξ ξ ξ ξ

2

ξ ξ ξ ξ

1 , 1 ; 2

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1

w

( ( ( ( ) ) ) )

1 ; 1

w

2 ; 1

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2

w

1 1

; 1

− − − − N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1

1 −

− − − N

w

1

; 1 N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1

N

w

2 , 1 ; 2

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 ; 1

w

1 , 1 ; 2

2 −

− − − N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1 ; 1

1 −

− − − N

w

2

, 1 ; 2 N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1

; 1 N

w

1 , 2 ; 2

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1 ; 2

w

2 , 2 ; 2

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 ; 2

w

( ( ( ( ) ) ) )

1 ; 2

2 −

− − − N

w

2

, 2 ; 2 N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2

; 2 N

w

1 , 2 ; 2

2 −

− − − N

ξ ξ ξ ξ

1 , ; 2

1

N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1 ;

1

N

w

2 , ; 2

1

N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 ;

1

N

w

( ( ( ( ) ) ) )

1 ;

2 1

− − − − N N

w

2 1 ,

; 2 N N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 1;N

N

w

1 , ; 2

2 1

− − − − N N

ξ ξ ξ ξ

1 , 1 ; 2

1 −

− − − N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

1 ; 1

1 −

− − − N

w

2 , 1 ; 2

1 −

− − − N

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 ; 1

1 −

− − − N

w

( ( ( ( ) ) ) )

1 ; 1

2 1

− − − − − − − − N N

w

2 1

, 1 ; 2 N N − − − −

ξ ξ ξ ξ

( ( ( ( ) ) ) )

2 1

; 1 N N

w

− − − −

1 , 1 ; 2

2 1

− − − − − − − − N N

ξ ξ ξ ξ

67 / 34

SLIDE 72

Multivariate Conditional QMOM

Table: Moments used to build a bivariate quadrature approximation (M = 2) for N1 = N2 = 3 using CQMOM with ξ2 conditioned on ξ1 (top) and ξ1 conditioned on ξ2 (bottom). M(5,0) M(4,0) M(3,0) M(2,0) M(2,1) M(2,2) M(2,3) M(2,4) M(2,5) M(1,0) M(1,1) M(1,2) M(1,3) M(1,4) M(1,5) M(0,0) M(0,1) M(0,2) M(0,3) M(0,4) M(0,5) M(5,0) M(5,1) M(5,2) M(4,0) M(4,1) M(4,2) M(3,0) M(3,1) M(3,2) M(2,0) M(2,1) M(2,2) M(1,0) M(1,1) M(1,2) M(0,0) M(0,1) M(0,2) M(0,3) M(0,4) M(0,5)

68 / 34

SLIDE 73

Inhomogeneous Quadrature Method of Moments

The spatially inhomogeneous GPBE operating on n(ξ) reads as follows:

∂n ∂t + ∂ ∂x ·

〈Up|ξ〉n
= ∂

∂x ·

D(ξ)∂n

∂x

− ∂

∂ξ

〈Gp|ξ〉n
+S

(43)

The application of the moment transform, Mk =

n(ξ)ξk dξ, generates several

closure problems: ∂Mk ∂t + ∂ ∂x ·

Uk

pMk

= ∂

∂x ·

Dk ∂Mk

∂x

+Sk

(44) where of course Uk

p =

〈Up|ξ〉n(ξ)ξk dξ

Mk

(similarly Dk) and S k =

∂

∂ξ

〈Gp|ξ〉n
+S
ξk dξ.

The solution with QMOM is as usual ...

69 / 34

SLIDE 74

Quadrature Method of Moments

From {M0,M1,...,M2N−1} the Gaussian quadrature with N nodes is

constructed resulting in the following approximation Mk(t,x) =

Ωξ

n(t,x,ξ)ξk dξ ≈

N

α=1

wα(t,x)(ξα(t,x))k (45)

that can be used to overcome the different closure problems, for example:

Uk

p(t,x) =

Ωξ〈Up|ξ〉n(t,x,ξ)ξk dξ

Mk ≈ N

α=1〈Up|ξα〉wαξk α

Mk (46)

By using different velocities, Uk

p(t,x), for the different moments, we can

describe mixing and segregation patterns

70 / 34

SLIDE 75

Direct quadrature method of moments

Moreover if wα and ξα are continuous functions of space and time, DQMOM can also be applied (slightly different now due to the diffusion term): ∂wα ∂t + ∂ ∂x ·

〈Up|ξα〉wα
= ∂

∂x ·

D(ξα)∂wα

∂x

+aα

(47) ∂wαξα ∂t + ∂ ∂x ·

〈Up|ξα〉wαξα
= ∂

∂x ·

D(ξα)∂wαξα

∂x

+bα

(48) with α ∈ 1,...,N and where the source terms are as usual determined by solving the following linear system: (1−k)

N

α=1

ξk

αaα +k N

α=1

ξk−1

α

bα = Sk +Ck, (49) where Ck = k(k −1)

N

α=1

ξk−2

α

Cα, Cα = wαD ∂ξα ∂x · ∂ξα ∂x

,

(50)

71 / 34

SLIDE 76

Spatial discretization and moment corruption

A moment set is said to be valid or realizable, if the Hankel-Hadamard determinants

are all non-negative: ∆k,l =

Mk

Mk+1 ... Mk+l Mk+1 Mk+2 ... Mk+l+1 . . . . . . . . . . . . Mk+l Mk+l+1 ... Mk+l+l

≥ 0,

k = 0,1, l ≥ 0 (51) for k = 0,1 and l ≥ 0.

For l = 1 we get:

MkMk+2 −M2

k+1 ≥ 0

k = 0,1,2,.... (52) which for k = 0 becomes the constraint of positive variance

Equivalently this becomes:

ln(Mk)+ln(Mk+2) 2 ≥ ln(Mk+1) k = 0,1,2,...; (53)

r in other words convexity of the function ln(Mk) with respect to k

72 / 34

SLIDE 77

Spatial discretization and moment corruption

Less stringent condition: convexity of function ln(Mk) with respect to k

Ad-hoc correction algorithms are needed!

73 / 34

SLIDE 78

Spatial discretization in QMOM

An alternative is to solve the transport equation for a generic moment Mk with discretization schemes that PRESERVE the moments ("Realizable schemes"): ∂Mk ∂t +

✚✚✚ ✚

Up ∂Mk ∂x = Sk −Up ∂Mk ∂x After spatial discretization (finite-volume):

W P E Ub x

dMP

k

dt = SP

k −

Up ∆x

Me

k −Mw k

With first-order upwind Me

k = MP k and Mw k = MW k

Spatial discretization schemes based on first-order upwind always result in VALID moments. Higher-order schemes (CDS, second-order upwind, QUICK, MUSCL) can result in INVALID moments.

74 / 34

SLIDE 79

Spatial discretization in QMOM

One solution would be to evaluate the moments at the faces Me

k and Mw k

through the quadrature approximation

We know the value of the moments at the center of the cells MW

k , MP k , ME k W P E Ub x

From these moments we can evaluate the corresponding weights wP

α and

abscissas ξP

α

If weights at the center of the face are interpolated with pth-order spatial

reconstruction and the abscissas are interpolated 1st-order spatial the resulting moments will be valid

This allows to improve the numerical accuracy preserving the moments! Alternatively one can use semi-Lagrangian schemes (Attili and Bisetti, 2013)

75 / 34

SLIDE 80

Spatial discretization in DQMOM

In QMOM the governing equations are for the moments that are ‘conserved’

variables

In DQMOM the governing equations are for weights and weighted abscissas

that are ‘primitive’ variables

When QMOM is used (if proper discretization schemes are used) only the

transported 2N moments are preserved and conserved (and their linear combination)

When DQMOM is used only weights and weighted abscissas are conserved

and their linear combination: M0 and M1

One disadvantage of DQMOM is that only two moments (of the 2N selected)

are conserved and saved from numerical errors!

Another way to look at the problem is to consider that the equations in

finite-volume codes are solved with an error called numerical diffusion whose coefficient is unknown!

76 / 34

SLIDE 81

Spatial discretization in DQMOM

The original DQMOM requires the solution of these equations:

∂wα ∂t +Up ∂wα ∂x = aα ∂wαξα ∂t +Up ∂wαξα ∂x = bα

∂Mk

∂t +Up ∂Mk ∂x = Sk

When the finite-volume discretization is applied:

dwP

α

dt + Up ∆x

we

α −Mw k

=

aP

α

d(wαξα)P dt + Up ∆x

(wαξα)e −(wαξα)w

= bP

α

dMP

k

dt + Up ∆x

Me

k −Mw k

=

SP

k

failing in conserving higher-order moments (k ≥ 2)

77 / 34

SLIDE 82

Spatial discretization in DQMOM

Solution: fully conservative version of DQMOM ( DQMOM-FC):

dMP

k

dt = − Up ∆x

Me

k −Mw k

+SP

k

dwP

α

dt = −aP

c,α +aP α

d(wαξα)P dt = −bP

c,α +bP α

now successfully conserving all the moments of the NDF

Alternatively one can use semi-Lagrangian schemes (Attili and Bisetti, 2013) Summarizing if the equations for the moments have large physical diffusion terms

DQMOM can be safely used (numerical diffusion will be smaller than real diffusion and the solution will contain no discontinuities and no shocks)

When only the equation for the moments do not contain any physical diffusion term

than DQMOM-FC should be used

78 / 34

SLIDE 83

Spatial discretization and moment corruption

The convexity of the function ln(Mk) with respect to k can be easily verified by building a difference table of ln(Mk). Example: VALID SET; moment of a Gaussian distribution (M0 = 1, M1 = 5, M2 = 26, M3 = 140, M4 = 778, M5 = 4450, M6 = 26140, M7 = 157400 ) k d0 = ln(Mk) d1 d2 d3 1.609 0.039

0.0043

1 1.609 1.648 0.034

0.0033

2 3.258 1.683 0.031

0.0027

3 4.941 1.715 0.028

0.0022

4 6.656 1.743 0.026

0.0019

5 8.400 1.770 0.024 6 10.171 1.795 7 11.966

79 / 34

SLIDE 84

Spatial discretization and moment corruption

The convexity of the function ln(Mk) with respect to k can be easily verified by building a difference table of ln(Mk). Example: INVALID SET; moment of a Gaussian distribution (M0 = 1, M1 = 5, M2 = 25, M3 = 140, M4 = 778, M5 = 4450, M6 = 26140, M7 = 157400) k d0 = ln(Mk) d1 d2 d3 1.609 0.113 1 1.609 1.609 0.113

0.121

2 3.218 1.722

0.007

0.036 3 4.941 1.715 0.028

0.002

4 6.656 1.743 0.026

0.001

5 8.400 1.770 0.024 6 10.171 1.795 7 11.966

80 / 34

SLIDE 85

Spatial discretization and moment corruption

When using first-order upwind spatial discretization schemes and first-order

explicit time discretization schemes the validity of the moment set should be preserved

In all the other cases it is very easy to CORRUPT the moment set and it is

anyway safe to have algorithms that DETECT CORRUPTION AND CORRECT invalid moment set

If we transform the moment set so that d2 is positive, we are almost sure that

the moment set is valid

But how positive? The moments of a log-normal distribution have the smallest d3

n(ξ) = NT σ

2πξ

exp

−
ln(ξ)−µ

2 2σ2

,

(54) Mk = NT exp

kµ+ k2σ2

2

,

(55)

The log-normal distribution is the smoothest distribution!

81 / 34

SLIDE 86

Spatial discretization and moment corruption

CORRECTION ALGORITHM BY MCGRAW

1. Build the difference table and check if d2 is negative
2. Identify the moment order k that causes the biggest change in d3
3. Change the moment (by multiplying it for a constant) in order to MINIMIZE

d3

4. Go back to point 1

CORRECTION ALGORITHM BY WRIGHT

1. Build the difference table and check if d2 is negative
2. Replace the moments with those of a log-normal distribution with

µ = j ij −i2 ln Mi M0

+

i ij −j2 ln Mj M0

(56)

σ2 = 1 1−i/j 2 j2 ln Mj M0

− 2

ij ln Mi M0

(57)

82 / 34

SLIDE 87

QBMM approximation vs DVM, Grad, oLBM – 4

M = 4 Normalized relaxation rate 2nd energy moment 3rd energy moment

83 / 34

SLIDE 88

QMOM approximation vs DVM, Grad, oLBM – 5

M = 4 BGK-equivalent relaxation time νp(t) = dΦp

dt 1 Φeq

p −Φp

2nd energy moment 3rd energy moment

84 / 34

SLIDE 89

Analytical equations for QBMM

dΦp dt = π

2

M/2

i=1

M/2

j=1

wiwjΛ+

ij,p − M/2

i=1

M/2

j=1

wiwjΛ−

ij,p

δij =
Ej/Ei;

qij = q(x,y,Ei,Ej); C+

ij,p = C+ p (x,y,Ei,Ej);

C−

i,p = C− p (Ei);

Λ+

ij,p =

+1

−1

+1

−1 |qij|C+ ij,p dxdy;

Λ−

ij,p =

+1

−1

+1

−1 |qij|C− i,p dxdy

Λ+

ij,p = 2Ep i

E+

p

α=0
p

α p−α

β=0
p−α

β

(−1)γij E−β−α

i

Eα

+ Eβ −

  2 2β+1   1−r2α+2

ij

2α+2 + r2α+2

ij

2α+2β+3  + 1 β+1   r2α+2

ij

2α+1 − r2α+2

ij

2α+2β+3     Λ−

ij,p = 2Ep i

E+

 1+ r2

ij

3  

γij =

α

ifEi ≥ Ej β ifEi < Ej E+ = max

Ei,Ej
; E− = min
Ei,Ej
; rij = min
Ei

Ej ,

Ej

Ei

85 / 34

SLIDE 90

Turbulent reactive flows

Ns species denoted by capital letters A,B,C,...,R,S Vector of concentrations φ = (φ1,...,φNs) Consider a simple reaction

A+B k − → R

Transport equations for concentrations with source terms

SA = −kφAφB = SB

In turbulent flows a RANS average (in time) or LES filter (in space) is

performed 〈SA〉 = −k〈φAφB〉 = −k〈φA〉〈φB〉

In general the term 〈S(φ)〉 = S(〈φ〉)

86 / 34

SLIDE 91

Turbulent reactive flows

The chemical source term can be closed if we assume the existence of a joint

probability density function of the concentrations 〈S(φ)〉 =

S(ψ)P(ψ;x,t)dψ

Rigorously, when using LES we are dealing with a filtered-density function

P(ψ;x,t) =

δ(ψ−φ(x,t))G(r−x)dr

where G is the LES filter

87 / 34

SLIDE 92

Turbulent reactive flows

For the box filter this could be no more a PDF so we simply assume a PDF P

that represents the spatial dishomogeneity in the cells such that

ψP(ψ;x,t)dψ = 〈φ〉

is the filtered scalar and

(ψ2

α −〈ψα〉2)P(ψα;x,t)dψα = 〈φ′2 α 〉

the scalar fluctuation

88 / 34

SLIDE 93

PDF Transport Equation

P can be solved in a advection-reaction-diffusion equation with Ns +4

independent variables

The turbulent diffusion term comes from the assumption that the

conditional fluctuations 〈U′

i|ψ〉P = −DT ∂P ∂xi

Molecular mixing term needs to be modeled

89 / 34

SLIDE 94

Mixing models

The most simple and popular model for Eulerian simulations is the IEM

〈D∇2φα|ψ〉 = − Cφ τ (ψ−〈φ〉)

It is a linear relaxation of the scalars to the mean value with time scale τ and a

parameter Cφ

τ is chosen to be a turbulence time scale, 2k

ǫ for RANS and 2∆2 D+DT where ∆ is

the filter width

Cφ is the scalar-to-mechanical time-scale ratio and it depends on the local

Schmidt and Reynolds numbers (for gases at high Re Cφ ≈ 2)

90 / 34

SLIDE 95

PDF discretization

In the Quadrature Method of Moments (QMOM) the integrals are

approximated using a quadrature rule

g(ψ)P(ψ;x,t)dψ ≈

M

i=1

g(φi)wi where φi are the abscissas and wi the weights of the quadrature rule

For a given set of 2M moments, the M abscissas and weights can be

calculated using inversion algorithms (Wheeler or Product-Difference)

This means that we are approximating the exact PDF P with a

multi-environment PDF f fφ(ψ;x,t) =

M

i=1

wi(x,t)δ[ψ−φi(x,t)] where δ is a multi-dimensional delta function

91 / 34

SLIDE 96

DQMOM-IEM

In the DQMOM transport equations for wi and wiψi are solved instead of

equations for µi

M(1+N) equations with some constraints (e.g.M

i=1 wi = 1)

The source term is such that the first M(1+N) moments are coherent with

the transported ones

Let us consider a competitive reaction scheme simplified using a mixture

fraction ξ and a reaction progress Y (linear combination of species concentration). This results in N = M = 2 and φ = (ξ,Y )

92 / 34

SLIDE 97

DQMOM-IEM

∂w1 ∂t +Ui ∂w1 ∂xi − ∂ ∂xi

DT

∂w1 ∂xi

= 0,

w2 = 1−w1 ∂w1ξ1 ∂t +Ui ∂w1ξ1 ∂xi − ∂ ∂xi

DT

∂w1ξ1 ∂xi

=

Cφ τ w1w2[ξ2 −ξ1]+ DT ξ1 −ξ2

w1

∂ξ1 ∂xi ∂ξ1 ∂xi +w2 ∂ξ2 ∂xi ∂ξ2 ∂xi

,

∂w2ξ2 ∂t +Ui ∂w2ξ2 ∂xi − ∂ ∂xi

DT

∂w2ξ2 ∂xi

=

Cφ τ w1w2[ξ1 −ξ2]+ DT ξ2 −ξ1

w1

∂ξ1 ∂xi ∂ξ1 ∂xi +w2 ∂ξ2 ∂xi ∂ξ2 ∂xi

,

93 / 34

SLIDE 98

DQMOM-IEM

∂w1Y1 ∂t +Ui ∂w1Y1 ∂xi − ∂ ∂xi

DT

∂w1Y1 ∂xi

= w1S(ξ1,Y1)+

Cφ τ w1w2[Y2 −Y1] + DT Y1 −Y2

w1

∂Y1 ∂xi ∂Y1 ∂xi +w2 ∂Y2 ∂xi ∂Y2 ∂xi

,

∂w1Y2 ∂t +Ui ∂w1Y2 ∂xi − ∂ ∂xi

DT

∂w1Y2 ∂xi

= w2S(ξ2,Y2)+

Cφ τ w1w2[Y1 −Y2] + DT Y2 −Y1

w1

∂Y1 ∂xi ∂Y1 ∂xi +w2 ∂Y2 ∂xi ∂Y2 ∂xi

.

Only the progress variable Y has a source term S

94 / 34

SLIDE 99

Poly-dispersed particles

Particle Size Distribution n(L;x,t) =

R3 n(L,Up;x,t)dUp

St = τp τf ≈ ρpD2 18µτf Particles of different size behave very differently with a non-linear relation between size and relative velocity

95 / 34

SLIDE 100

Poly-dispersed particles

Quadrature-based Moment Method (QBMM) nodes and weights to approximate f (L) 2 moments → 1 node, 1 weight Mean size

96 / 34

SLIDE 101

Real-space advection

Population balance equation: ∂tn+∂x〈Up|L〉n+✘✘

✘ ✿ no growth term

∂LnG =✘✘✘

✘ ✿ no aggregation,breakage,nucleation

Q(n,n) Models for conditional velocity 〈Up|L〉 = Up(L)

PSEUDO-HOMOGENEOUS MODEL, St ≪ 1

Particles flow with fluid velocity Up(L) = U

ALGEBRAIC MODELS (MIXTURE), St ≤ 1

Up(L) is calculated with algebraic relations from the knowledge of U

DIFFERENTIAL MODELS (MULTI-FLUID), St > 1

Up(L) is solved with a momentum balance equation

97 / 34

SLIDE 102

Equilibrium Model (Ferry/Balachandar)

similar to the classical Algebraic Slip Model expansion of Maxey-Riley equation for small particle response time τp

Up ≈ U+O(τp) ⇒ Up −U = −τp DU Dt −g

1st order correction

Up −U = −τp

I+τp∇UT−1 DU

Dt −g

High Rep effects considered using a modified τ∗

p

can be extended to take into account other forces and two-way coupling

98 / 34

SLIDE 103

SGS model for particles

In Large Eddy Simulation (LES) we have only filtered velocity field U

APPROXIMATE DECONVOLUTION METHOD (ADM)

To estimate the “real” unfiltered velocity we use an approximate inverse filter U∗

i = Nc

k=0

(I−G)k Ui

G is a test filter Nc is the deconvolution order if Nc = 1 it is like an inverse dynamic model considered Nc = 1,5 with Gaussian weights

99 / 34

SLIDE 104

PBE approximation

Single momentum equation for the mixture Equilibrium model to calculate velocity Up(L) Population balance discretized with Quadrature-Based Moment Method

(QBMM)

Solved directly in terms of nodes Li and weights wi (DQMOM formulation)

SINGLE VELOCITY MODEL

All particles flow with an overall mean velocity Up(L) = Up(L)

MULTIPLE VELOCITIES MODEL

Each node (particle class) considered as a separate phase with its own relative velocity

100 / 34

SLIDE 105

Turbulent poly-dispersed channel flow

101 / 34