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 oLBM 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 oLBM 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 oLBM - - 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 oLBM - - 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
HIBE steady state Quadrature error at the true equilibrium 20 / 34
QBMM - multivariate case The Generalized Population Balance Equation for M -dimensional phase space: ∂ n ∂ t + ∂ = ∂ � D ∂ n � − ∂ � � ∂ x · U p n ∂ x · ∂ ξ · ( G n ) + Q , (7) ∂ x After having defined a generic moments : � ξ k 1 1 ··· ξ k M M k ≡ M n ( t , ξ )d ξ (8) Ω ξ equivalent notations M k 1 ,..., k M = M k = M ( k 1 ,..., k M ) = M ( k ), we solve the resulting transport equations: ∂ M k + ∂ − ∂ � D k ∂ M k � � − ∂ � � � � U k ξ k � � ∂ x · p M k ∂ x · = ∂ ξ · 〈 G p | ξ 〉 n + Q d ξ (9) ∂ t ∂ x Ω ξ and the closure problem is overcome by using the quadrature approximation. 21 / 34
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) N M k i β � � M k i 1 , k i 2 ,..., k iM = M ( k i ) = w α ξ β , α , 1 ≤ i ≤ N ( M + 1) α = 1 β = 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
QBMM - multivariate case Bivariate Gaussian distr. with M = 9 Direct (Newton) (diamond), Tensor product (circle) and Conditional (square) 23 / 34
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
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 N ON - LINEAR F OKKER -P LANCK � Attar and Vedula, 2008; Otten and Vedula, 2011 � Drift-diffusion equations x i = h i ( x , t ) + g i ( x , t ) W i , ˙ W ij � Applications to mixing in porous media N ON - LINEAR FILTERING AND B AYESIAN 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
PART II: MLMC for dispersed multiphase flows 25 / 34
Suspensions/dispersed flows DNS of fluid-solid interface (particle-resolved simulations) are often unfeasible W HY � Limits in resolution E XAMPLES � Physical models too complex � Fluidized beds, separators (e.g. particulate processes) � Environmental flows � We want a steady solution or � Pore-scale processes 2D/1D models � Long piping systems � Unknown particle shape/properties 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
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 F d = F d ( Re , φ ) where F d 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
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 ∝ F d ∝ µ 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 F d = F d ( Re , φ , θ ) and F ′ d = F ′ d ( Re , φ , θ ) θ is a “granular temperature” and F ′ d are the random drag fluctuations 28 / 34
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 A LGORITHM 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 O NGOING 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
Problem formulation We developed fully automatic code to solve different stochastic PDE problems and compute effective properties: A DIMENSIONAL STEADY INCOMPRESSIBLE N AVIER -S TOKES u ( x ; ω ) ·∇ u ( x ; ω ) +∇ p ( x ; ω ) − 1 x ∈ Γ ( ω ) ⊂ [0,1] d , Re ∇· ( ∇ u ( x ; ω )) = 0, ∇· u = 0 u = 0, ∇ n p = 0 no-slip: on I g ∇ n u = 0, p = p i ( x ) I i = { x : x 1 = 0} inlet: on ∇ n u = 0, p = p o ( x ) = 0 I o = { x : x 1 = 1} outlet: on symmetry: ∇ n u t = 0, u n = 0, ∇ n p = 0 on I s = { x : x 2 , x 3 = 0} The quantity of interest Q ( u ) is, in this case, the mean flux 〈 u x 〉 or, equivalently, the drag force F d = ∆ p (1 − φ ) d 2 18 〈 u x 〉 The algorithm is tested for creeping flow through random array of mono dispersed spheres. Extended to non-spherical, polydispersed, and moving 2 particles . 2 in local equilibrium 30 / 34
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
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
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
Conclusions A CHIEVEMENTS � 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 O PEN 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
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
BACKUP SLIDES 34 / 34
Smoluchowski coagulation equation Describe particle aggregation/coagulation � L � ∞ ∂ n ( L , t ) = 1 K ( L − ˜ L , ˜ L ) n ( L − ˜ L , t ) n (˜ L , t ) d ˜ K ( L , ˜ L ) n ( L , t ) n (˜ L , t ) d ˜ L − L ∂ t 2 0 0 � 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
Population Balance Equation (PBE) Extension of the Smoluchowski coagulation equation, very popular in Chemical Engineering, Biology and Social Models ∂ n ∂ t + ∂ ∂ x · ( U n ) + ∂ ∂ ξ · ( G n ) = 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 = Q 0 + Q 1 ( n ) + Q 2 ( 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 3 with frequency kernels that can be seen as upscaled jump processes 36 / 34
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
Boltzmann Equation (BE) From the Liouville n-particles equation, assuming indistinguishable particle and the "Stosszahlansatz" (molecular chaos 4 hypothesis), the Boltzmann equation is obtained ∂ n ( U ) + ∂ ∂ x · ( U n ( U )) + ∂ ∂ U · ( A n ( U )) = ∂ t � � n ( U ′ ) n ( U ′ � � � � ∗ ) − n ( U ) n ( U ∗ ) β g , x d s d U ∗ R 3 S + � 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, ... 4 particles are uncorrelated, two-particles PDF is product of one-particle PDF 38 / 34
Boltzmann Equation and Equilibrium Most of the studies and methods for BE are based on the concept of Equilibrium Distribution � The equilibrium distribution f eq is such that Q ( f eq , f eq ) = 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 ≈ f eq (1 + a 1 H 1 ( v ) + a 2 H 2 ( v ) + ...) � Other methods (such as Lattice Boltzmann) rely on a linearized collision term (BGK, Bhatnagar-Gross-Krook approximation) Q ≈ 1 ( f − f eq ) τ r All these methods require the prior knowledge of the equilibrium distribution! 39 / 34
The Method of Moments � MOM originates in kinetic theory of gases 5 : n ( t , x , U ) � The moment of order zero is related to the density: � +∞ ρ ( t , x ) ≡ mM 0,0,0 ( t , x ) = m n ( t , x , U )d U −∞ � The moment of order one is used to define the average velocity: � M 1,0,0 , M 0,1,0 , M 0,0,1 � 〈 U 〉 = M 0,0,0 M 0,0,0 M 0,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) 5 The same name however may refer to many different methods in different fields 40 / 34
Flow regimes The Boltzmann Equation is not only for rarefied gases! � The flow regime of a granular gas depends on Knudsen number Kn = λ L 0 � 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. Collisionless Boltzmann equation Boltzmann eq. Euler N-S Kn 0.01 0.1 1 10 ∞ 0 100 Adapted from: G. A. Bird (1994) 41 / 34
The GPBE fluid-particle systems � n ( t , x , U p , ξ p , U f , ξ f ), represents the number of particles (per unit volume) with velocity equal to U p , internal coordinate ξ p and see a continuous phase with velocity U f and internal coordinate ξ f . � The evolution of the NDF is dictated by the GPBE: ∂ n ∂ t + ∂ ∂ ∂ �� � � � � � � ∂ x · U p n + · A fp + A p n + · G p n ∂ U p ∂ ξ p ∂ + ∂ �� � � � � + · A pf + A f n · G f n = Q (10) ∂ U f ∂ ξ f � phase-space velocity for particle velocity: A fp + A p ( average particle acceleration → acceleration model (drag, lift, ...) � phase-space velocity for particle internal coordinate: G p (e.g. particle growth rate) � Discontinuous jump term: Q (e.g. particle collision, aggregation, breakage, nucleation, etc.) 42 / 34
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
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: N � M 0 = w α , α = 1 . . (11) . N w α ξ 2 N − 1 � M 2 N − 1 = . α α = 1 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
Gaussian quadrature A set of polynomials { P 0 ( ξ ), P 1 ( ξ ),..., P α ( ξ ),...} with P α ( ξ ) = k α ,0 x α + k α ,1 x α − 1 +···+ k α , α , is said to be orthogonal in the integration interval Ω ξ , with respect to the weight function, if � = 0 for α �= β , � n ( ξ ) P α ( ξ ) P β ( ξ )d ξ (12) > 0 for α = β , Ω ξ and, of course, is said to be orthonormal if � 0 for α �= β , � n ( ξ ) P α ( ξ ) P β ( ξ )d ξ = (13) 1 for α = β . Ω ξ 45 / 34
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 P 0 ( ξ ) ≡ 1 and where � Ω ξ n ( ξ ) ξ P α ( ξ ) P α ( ξ )d ξ a α = α = 0,1,2,... Ω ξ n ( ξ ) P α ( ξ ) P α ( ξ )d ξ , (15) � � Ω ξ n ( ξ ) P α ( ξ ) P α ( ξ )d ξ b α = Ω ξ n ( ξ ) P α − 1 ( ξ ) P α − 1 ( ξ )d ξ , α = 1,2,.... (16) � One can calculate a 0 , then P 1 ( ξ ), then a 1 and b 1 and so on... 46 / 34
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 2 N − 1 moments of the NDF � For example with M 0 , M 1 , M 2 and M 3 , it is possible to calculate the following coefficients: a 0 = M 1 , M 0 M 3 M 2 0 + M 3 1 − 2 M 2 M 1 M 0 a 1 = , (17) M 2 M 0 + M 2 1 − 2 M 2 1 M 0 M 2 M 0 + M 2 1 − 2 M 2 1 M 0 b 1 = , M 2 0 which suffice for the calculation of the polynomial P 2 ( ξ ) 47 / 34
Gaussian quadrature GAUSSIAN QUADRATURE The necessary and sufficient condition for the following formula: N � � n ( ξ ) g ( ξ )d ξ = g ( ξ α ) w α + R N ( g ), (18) Ω ξ α = 1 to be a Gaussian quadrature approximation is that its nodes { ξ α } coincide with the N roots of the polynomial P N ( ξ ) of order N orthogonal in Ω ξ with respect to the weight function n ( ξ ). 48 / 34
Direct quadrature method of moments � The fact that the closure problem is overcome with the quadrature approximation: N � � n ( ξ ) g ( ξ )d ξ ≈ w α g ( ξ α ), (19) Ω ξ α = 1 � Is equivalent to the assumption that the NDF is as follows: N � n ( ξ ) = w α δ ( ξ − ξ α ), (20) α = 1 � 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
Direct quadrature method of moments � Assuming that the weights and nodes are differentiable in space/time the following transport equation is obtained: N N � D w α � � D ξ α � δ ′ ( ξ − ξ α ) � � δ ( ξ − ξ α ) − = S ( ξ ), w α (21) D t D t α = 1 α = 1 � If the weighted nodes (or weighted abscissas) ς α = w α ξ α are introduced: N N � D w α � � D w α + D ς α � δ ′ ( ξ − ξ α ) � � δ ( ξ − ξ α ) − − ξ α D t D t D t α = 1 α = 1 = S ( ξ ). (22) 50 / 34
Direct quadrature method of moments � We now define a α and b α to be the source terms: D w α D ς α = a α , D t = b α . (23) D t � Using these definitions Eq. (22) can be rewritten in a simpler form: N N δ ( ξ − ξ α ) + δ ′ ( ξ − ξ α ) ξ α δ ′ ( ξ − ξ α ) b α = S ( ξ ). � � � � a α − (24) α = 1 α = 1 � This equation can now be used to determine the unknown functions a α and b α by applying the moment transformation. 51 / 34
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 , −∞ (25) � +∞ ξ k δ ′ ( ξ − ξ α )d ξ = − k ( ξ α ) k − 1 , −∞ � The moment transform of Eq. (24) yields N N ξ k ξ k − 1 � � (1 − k ) α a α + k b α = S k , (26) α α = 1 α = 1 with k = k 1 , k 2 ,..., k 2 N . 52 / 34
Direct quadrature method of moments � The linear system in Eq. (26) can be written in matrix form: A α = d . (27) where � a � � T = � α = a 1 ··· a N b 1 ··· b N , (28) b � T , d = � S k 1 ··· S k 2 N − 1 (29) � The components of the matrix A are ξ k i � � 1 − k i if 1 ≤ j ≤ N , j a ij = (30) k i ξ k i − 1 if N + 1 ≤ j ≤ 2 N . j 53 / 34
Direct quadrature method of moments � If (as in QMOM) the first 2 N integer moments are chosen (i.e., k = 0,...,2 N − 1), the matrix of the linear system is A = 1 ··· 1 0 ··· 0 0 ··· 0 1 ··· 1 − ξ 2 − ξ 2 ··· 2 ξ 1 ··· 2 ξ N 1 N . . . . . . . . . . . . . 2(1 − N ) ξ 2 N − 1 2(1 − N ) ξ 2 N − 1 (2 N − 1) ξ 2 N − 2 (2 N − 1) ξ 2 N − 2 ··· ··· 1 N 1 N (31) � A does not depend on the weights w α and if the abscissas ξ α are unique, then A will be full rank. 54 / 34
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: D w α ξ α D w α = a α , = b α . (32) D t D t � The source terms are calculated by inverting the linear system and by using the following initial condition: w α (0) = w 0 ς α (0) = w 0 α ξ 0 for k = 1,..., N . α , (33) α in turn calculated from the initial moments 55 / 34
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 or 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
Product-difference algorithm A smarter way is to employ the recursive relationship for the orthogonal polynomials: a 0 1 P 0 ( ξ ) P 0 ( ξ ) 0 b 1 a 1 1 P 1 ( ξ ) P 1 ( ξ ) 0 ... . . . . = . + . ξ . (34) . . . ... P N − 2 ( ξ ) P N − 2 ( ξ ) 0 1 P N − 1 ( ξ ) P N − 1 ( ξ ) P N ( ξ ) a N − 1 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
Product-difference algorithm In fact the matrix can be made symmetric (preserving the eigenvalues) by a diagonal similarity transformation to give a Jacobi matrix: � a 0 b 1 � � b 1 a 1 b 2 � � b 2 a 2 b 3 ... � b 3 a 3 J = (35) ... ... ... ... � a N − 2 b N − 1 � b N − 1 a N − 1 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 α = M 0 ϕ 2 α 1 where ϕ α 1 is the first component of the α th eigenvector ϕ α of the Jacobi matrix. 58 / 34
Product-difference algorithm 1. Construct the matrix P with components P α , β : P α , β = P 1, β − 1 P α + 1, β − 2 − P 1, β − 2 P α + 1, β − 1 β ∈ 3,...,2 N + 1 and α ∈ 1,...,2 N + 2 − j . (36) 2. where the first row of the matrix is: P α ,1 = δ α 1 α ∈ 1,...,2 N + 1, (37) 3. where δ α 1 is the Kronecker delta and where the components in the second column of P are P α ,2 = ( − 1) α − 1 M α − 1 α ∈ 1,...,2 N . (38) 4. Calculate the coefficients of the continued fraction { ζ α }: P 1, α + 1 ζ α = α ∈ 2,...,2 N . (39) P 1, α P 1, α − 1 59 / 34
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 M 0 M 2 − ( M 1 ) 2 M 3 M 1 − ( M 2 ) 2 � � 1 M 0 M 1 M 0 0 − M 1 − M 2 − ( M 0 M 3 − M 2 M 1 ) 0 M 2 M 3 . (42) − M 3 0 0 60 / 34
Extended QMOM � only integral representation, no localized processes in phase space, smoothness of kernels β , G ,... � It fails in the computation of the entropy � Ω ξ f log f 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
Extended Quadrature-based moment methods EQMOM for a dissolution problem G AUSSIAN DISTR . S OLUTION WITH N = 4 BETA DISTR . 1.6 0.125 0.25 1.4 0.5 1 1.2 1 n( ξ ) 0.8 0.6 0.4 0.2 0 −2 −1 0 1 2 3 4 ξ a a Yuan, 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 − ξα ) σ − 1 (1 − ξ ) − 1 � − ( ξ − ξ α ) 2 � σ 1 ; δ σ = ξ δ σ = � exp 2 σ 2 B ( ξ α , σ ) 2 πσ 62 / 34
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
Optimal moment set Table: Moments used to build a bivariate quadrature approximation ( M = 2) for N = 2. In this case M 0,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 M 2,1 , M 1,2 and M 0,3 are chosen among the third-order moments to saturate the degrees of 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
Optimal moment set Table: Optimal moment set used to build a bivariate quadrature approximation ( M = 2) for N = 4. Only when N 1/ 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
Multivariate Conditional QMOM � Methods based on conditional density functions: n ( ξ 1 , ξ 2 ) = n 1 ( ξ 1 ) f 21 ( ξ 2 | ξ 1 ) = n 2 ( ξ 2 ) f 12 ( ξ 1 | ξ 2 ) � Univariate quadrature ( N 1 ) calculated from the first 2 N 1 − 1: M 0,0,...,0,0 w 1 ξ 1;1 . PD/Wheeler . . − − − − − − − − − → . . . . . . . M 2 N 1 − 1,0,...,0,0 w N 1 ξ 1; N 1 resulting for example in: n ( ξ 1 , ξ 2 ) = � N 1 � � α 1 = 1 w α 1 δ ξ 1 − ξ 1; α 1 f 21 ( ξ 2 | ξ 1 ) � The generic moment becomes: � n ( ξ 1 , ξ 2 ) ξ k 1 1 ξ k 2 M k 1 , k 2 = 2 d ξ 1 d ξ 2 N 1 � w α 1 ξ k 1 f 21 ( ξ 2 | ξ 1; α ) ξ k 2 � = 2 d ξ 2 1; α 1 α 1 = 1 � ξ k 2 � f 12 ( ξ 2 | ξ 1; α 1 ) ξ k 2 � Conditional moment: = � 2 d ξ 2 2 α 1 66 / 34
Multivariate Conditional QMOM For each of these N 1 nodes, 2 N 2 − 1 conditional moments are calculated, and univariate quadratures N 2 are determined (in direction ξ 2 ): Conditional QMOM or CQMOM ( ( ( ( ) ) ) ) ξ ξ ξ ξ ξ ξ ξ ξ w 2 ; N − − − − 1 , N N − − − − 1 N ; 2 1 2 1 2 ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) ξ ξ ξ ξ ξ ξ ξ ξ w w 2 ; 2 , N 2 ; N − − − − 1 , N − − − − 1 2 N ; 1 ; 1 2 N − − − − N − − − − 1 2 2 1 2 ( ( ( ( ) ) ) ) ξ ξ ξ ξ ( ( ( ( ) ) ) ) ξ ξ ξ ξ ξ ξ w ξ ξ w 2 ; N 1 , N 2 ; 1 , N 1 ; N N 2 1 N ; 2 2 ; 2 , N 2 − − − − 1 1 2 ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) ξ ξ ξ ξ ξ ξ ξ ξ w w 2 ; , 1 2 ; 1 , N 2 − − − − 1 2 ; N 2 − − − − 1 N N − − − − 1 ; N 1 − − − − 1 ( ( ( ( 1 2 ) ) ) ) ξ ξ ξ ξ w 2 ; N 1 − − − − 1 , 2 ; 1 N N − − − − ( ( ( ( ) ) ) ) 1 2 ( ( ( ( ) ) ) ) ξ ξ ξ ξ w ( ( ( ( ) ) ) ) ξ ξ ξ ξ w N 1 − − − − 1 ; 2 w ( ( ( ( ) ) ) ) 2 ; 2 , 2 ξ ξ ξ ξ 2 ; 2 2 ; 1 , 2 1 ; 2 ξ ξ ξ ξ w 2 ; , 2 N N ; 2 ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) 2 ; 1 − − − − 1 , 1 1 ξ ξ ξ ξ ξ ξ N 1 ξ ξ ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) w w ξ ξ ξ ξ 2 ; 1 , 1 2 ; 2 , 1 1 ; 1 2 ; 1 w w 2 ; , 1 N N 1 − − − − 1 ; 1 N ; 1 1 1 ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ � 1 ; 1 1 ; 2 1 ; 1 N ; 1 N − − − − 1 1 1 ( ( ) ) ( ( ( ( ) ) ) ) ( ( ) ) ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) � w w w w N 1 − − − − 1 N 1 2 1 67 / 34
Multivariate Conditional QMOM Table: Moments used to build a bivariate quadrature approximation ( M = 2) for N 1 = N 2 = 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
Inhomogeneous Quadrature Method of Moments � The spatially inhomogeneous GPBE operating on n ( ξ ) reads as follows: ∂ n ∂ t + ∂ = ∂ � D ( ξ ) ∂ n � − ∂ � � � � ∂ x · 〈 U p | ξ 〉 n ∂ x · 〈 G p | ξ 〉 n + S (43) ∂ x ∂ξ n ( ξ ) ξ k d ξ , generates several � The application of the moment transform, M k = � closure problems: ∂ M k + ∂ = ∂ � D k ∂ M k � � � U k ∂ x · p M k ∂ x · + S k (44) ∂ t ∂ x 〈 U p | ξ 〉 n ( ξ ) ξ k d ξ � where of course U k (similarly D k ) and p = M k �� � ξ k d ξ . ∂ � � S k = 〈 G p | ξ 〉 n + S ∂ξ � The solution with QMOM is as usual ... 69 / 34
Quadrature Method of Moments � From { M 0 , M 1 ,..., M 2 N − 1 } the Gaussian quadrature with N nodes is constructed resulting in the following approximation N � n ( t , x , ξ ) ξ k d ξ ≈ w α ( t , x )( ξ α ( t , x )) k � M k ( t , x ) = (45) Ω ξ α = 1 � that can be used to overcome the different closure problems, for example: Ω ξ 〈 U p | ξ 〉 n ( t , x , ξ ) ξ k d ξ � � N α = 1 〈 U p | ξ α 〉 w α ξ k α U k p ( t , x ) = ≈ (46) M k M k � By using different velocities, U k p ( t , x ), for the different moments, we can describe mixing and segregation patterns 70 / 34
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 α + ∂ = ∂ � D ( ξ α ) ∂ w α � � � ∂ x · 〈 U p | ξ α 〉 w α ∂ x · + a α (47) ∂ t ∂ x ∂ w α ξ α + ∂ = ∂ � D ( ξ α ) ∂ w α ξ α � � � ∂ x · 〈 U p | ξ α 〉 w α ξ α ∂ x · + b α (48) ∂ t ∂ x with α ∈ 1,..., N and where the source terms are as usual determined by solving the following linear system: N N ξ k ξ k − 1 � � (1 − k ) α a α + k b α = S k + C k , (49) α α = 1 α = 1 where N � ∂ξ α ∂ x · ∂ξ α � ξ k − 2 � C k = k ( k − 1) C α , C α = w α D , (50) α ∂ x α = 1 71 / 34
Spatial discretization and moment corruption � A moment set is said to be valid or realizable, if the Hankel-Hadamard determinants are all non-negative: � M k M k + 1 ... M k + l � � � � M k + 1 M k + 2 ... M k + l + 1 � � � � � ∆ k , l = ≥ 0, k = 0,1, l ≥ 0 (51) . . . . � � . . . . � � . . . . � � � � M k + l M k + l + 1 ... M k + l + l � � for k = 0,1 and l ≥ 0. � For l = 1 we get: M k M k + 2 − M 2 k + 1 ≥ 0 k = 0,1,2,.... (52) which for k = 0 becomes the constraint of positive variance � Equivalently this becomes: ln( M k ) + ln( M k + 2 ) ≥ ln( M k + 1 ) k = 0,1,2,...; (53) 2 or in other words convexity of the function ln( M k ) with respect to k 72 / 34
Spatial discretization and moment corruption Less stringent condition: convexity of function ln( M k ) with respect to k Ad-hoc correction algorithms are needed! 73 / 34
Spatial discretization in QMOM An alternative is to solve the transport equation for a generic moment M k with discretization schemes that PRESERVE the moments ("Realizable schemes"): ✚ ✚✚✚ ∂ M k ∂ M k ∂ M k + U p = S k − U p ∂ t ∂ x ∂ x After spatial discretization (finite-volume): U b d M P U p = S P k M e k − M w � � k − k d t ∆ x x W P E With first-order upwind M e k = M P k and M w k = M W 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
Spatial discretization in QMOM � One solution would be to evaluate the moments at the faces M e k and M w k through the quadrature approximation � We know the value of the moments at the center of the cells M W k , M P k , M E k U b x W P E � From these moments we can evaluate the corresponding weights w P α and abscissas ξ P α � If weights at the center of the face are interpolated with p th -order spatial reconstruction and the abscissas are interpolated 1 st -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
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 2 N moments are preserved and conserved (and their linear combination) � When DQMOM is used only weights and weighted abscissas are conserved and their linear combination: M 0 and M 1 � One disadvantage of DQMOM is that only two moments (of the 2 N 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
Spatial discretization in DQMOM � The original DQMOM requires the solution of these equations: ∂ w α ∂ w α + U p = a α ∂ t ∂ x ∂ w α ξ α ∂ w α ξ α + U p = b α ∂ t ∂ x � ∂ M k ∂ M k + U p = S k ∂ t ∂ x � When the finite-volume discretization is applied: d w P U p α w e α − M w a P + � � = α k d t ∆ x d( w α ξ α ) P U p ( w α ξ α ) e − ( w α ξ α ) w � b P + � = α d t ∆ x � d M P U p S P k � M e k − M w � + = k k d t ∆ x failing in conserving higher-order moments ( k ≥ 2) 77 / 34
Spatial discretization in DQMOM � Solution: fully conservative version of DQMOM ( DQMOM-FC): d M P U p + S P k � � M e k − M w = − k k d t ∆ x � d w P α − a P c , α + a P = α d t d( w α ξ α ) P − b P c , α + b P = α d t 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
Spatial discretization and moment corruption The convexity of the function ln( M k ) with respect to k can be easily verified by building a difference table of ln( M k ). Example: VALID SET; moment of a Gaussian distribution ( M 0 = 1, M 1 = 5, M 2 = 26, M 3 = 140, M 4 = 778, M 5 = 4450, M 6 = 26140, M 7 = 157400 ) d 0 = ln( M k ) k d 1 d 2 d 3 0 0 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 0 6 10.171 1.795 0 0 7 11.966 0 0 0 79 / 34
Spatial discretization and moment corruption The convexity of the function ln( M k ) with respect to k can be easily verified by building a difference table of ln( M k ). Example: INVALID SET; moment of a Gaussian distribution ( M 0 = 1, M 1 = 5, M 2 = 25, M 3 = 140, M 4 = 778, M 5 = 4450, M 6 = 26140, M 7 = 157400) d 0 = ln( M k ) k d 1 d 2 d 3 0 0 1.609 0 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 0 6 10.171 1.795 0 0 7 11.966 0 0 0 80 / 34
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 d 2 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 d 3 � 2 � � � − ln( ξ ) − µ N T n ( ξ ) = � exp , (54) 2 σ 2 σ 2 πξ � � k µ + k 2 σ 2 M k = N T exp , (55) 2 � The log-normal distribution is the smoothest distribution! 81 / 34
Spatial discretization and moment corruption CORRECTION ALGORITHM BY MCGRAW 1. Build the difference table and check if d 2 is negative 2. Identify the moment order k that causes the biggest change in d 3 3. Change the moment (by multiplying it for a constant) in order to MINIMIZE d 3 4. Go back to point 1 CORRECTION ALGORITHM BY WRIGHT 1. Build the difference table and check if d 2 is negative 2. Replace the moments with those of a log-normal distribution with � M i � M j j � i � µ = ij − i 2 ln + ij − j 2 ln (56) M 0 M 0 � 2 � M j � M i 1 � − 2 �� σ 2 = j 2 ln ij ln (57) 1 − i / j M 0 M 0 82 / 34
QBMM approximation vs DVM, Grad, oLBM – 4 M = 4 Normalized relaxation rate 2nd energy moment 3rd energy moment 83 / 34
QMOM approximation vs DVM, Grad, oLBM – 5 M = 4 BGK-equivalent relaxation time ν p ( t ) = d Φ p 1 Φ eq dt p − Φ p 2nd energy moment 3rd energy moment 84 / 34
Analytical equations for QBMM � M /2 � M /2 M /2 M /2 d Φ p = π w i w j Λ + w i w j Λ − � � � � � ij , p − ij , p dt 2 i = 1 j = 1 i = 1 j = 1 � C + ij , p = C + δ ij = E j / E i ; q ij = q ( x , y , E i , E j ); p ( x , y , E i , E j ); C − i , p = C − Λ + � + 1 � + 1 − 1 | q ij | C + Λ − � + 1 � + 1 − 1 | q ij | C − ij , p = ij , p = p ( E i ); ij , p dxdy ; i , p dxdy − 1 − 1 p � � p − α � � p − α p ( − 1) γ ij E − β − α Λ + ij , p = 2 E p E α + E β � � � E + − i i α β α = 0 β = 0 1 − r 2 α + 2 r 2 α + 2 r 2 α + 2 r 2 α + 2 2 1 ij ij ij ij + + 2 α + 1 − 2 β + 1 2 α + 2 2 α + 2 β + 3 β + 1 2 α + 2 β + 3 r 2 ij ij , p = 2 E p Λ − � 1 + E + i 3 � �� � � α if E i ≥ E j Ej Ei � � � � γ ij = E + = max ; E − = min ; r ij = min E i , E j E i , E j Ej , Ei if E i < E j β 85 / 34
Turbulent reactive flows � N s species denoted by capital letters A , B , C ,..., R , S � Vector of concentrations φ = ( φ 1 ,..., φ N s ) � Consider a simple reaction A + B k − → R � Transport equations for concentrations with source terms S A = − k φ A φ B = S B � In turbulent flows a RANS average (in time) or LES filter (in space) is performed 〈 S A 〉 = − k 〈 φ A φ B 〉 �= − k 〈 φ A 〉〈 φ B 〉 � In general the term 〈 S ( φ ) 〉 �= S ( 〈 φ 〉 ) 86 / 34
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 ) d r where G is the LES filter 87 / 34
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 ) P ( ψ α ; x , t ) d ψ α = 〈 φ ′ 2 ( ψ 2 α 〉 the scalar fluctuation 88 / 34
PDF Transport Equation � P can be solved in a advection-reaction-diffusion equation with N s + 4 independent variables � The turbulent diffusion term comes from the assumption that the conditional fluctuations 〈 U ′ i | ψ 〉 P = − D T ∂ P ∂ x i � Molecular mixing term needs to be modeled 89 / 34
Mixing models � The most simple and popular model for Eulerian simulations is the IEM C φ 〈 D ∇ 2 φ α | ψ 〉 = − τ ( ψ −〈 φ 〉 ) � It is a linear relaxation of the scalars to the mean value with time scale τ and a parameter C φ 2 ∆ 2 � τ is chosen to be a turbulence time scale, 2 k ǫ for RANS and D + D T 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
PDF discretization � In the Quadrature Method of Moments (QMOM) the integrals are approximated using a quadrature rule M � � g ( ψ ) P ( ψ ; x , t ) d ψ ≈ g ( φ i ) w i i = 1 where φ i are the abscissas and w i the weights of the quadrature rule � For a given set of 2 M 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 M � f φ ( ψ ; x , t ) = w i ( x , t ) δ [ ψ − φ i ( x , t )] i = 1 where δ is a multi-dimensional delta function 91 / 34
DQMOM-IEM � In the DQMOM transport equations for w i and w i ψ i are solved instead of equations for µ i i = 1 w i = 1) � M (1 + N ) equations with some constraints (e.g. � M � 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
DQMOM-IEM ∂ w 1 ∂ w 1 − ∂ � ∂ w 1 � + U i = 0, D T ∂ t ∂ x i ∂ x i ∂ x i w 2 = 1 − w 1 ∂ w 1 ξ 1 ∂ w 1 ξ 1 − ∂ � ∂ w 1 ξ 1 � + U i D T ∂ t ∂ x i ∂ x i ∂ x i C φ D T � ∂ξ 1 ∂ξ 1 ∂ξ 2 ∂ξ 2 � = τ w 1 w 2 [ ξ 2 − ξ 1 ] + w 1 + w 2 , ξ 1 − ξ 2 ∂ x i ∂ x i ∂ x i ∂ x i ∂ w 2 ξ 2 ∂ w 2 ξ 2 − ∂ � ∂ w 2 ξ 2 � + U i D T ∂ t ∂ x i ∂ x i ∂ x i C φ D T � ∂ξ 1 ∂ξ 1 ∂ξ 2 ∂ξ 2 � = τ w 1 w 2 [ ξ 1 − ξ 2 ] + w 1 + w 2 , ξ 2 − ξ 1 ∂ x i ∂ x i ∂ x i ∂ x i 93 / 34
DQMOM-IEM ∂ w 1 Y 1 ∂ w 1 Y 1 − ∂ � ∂ w 1 Y 1 � + U i D T ∂ t ∂ x i ∂ x i ∂ x i C φ = w 1 S ( ξ 1 , Y 1 ) + τ w 1 w 2 [ Y 2 − Y 1 ] � � D T ∂ Y 1 ∂ Y 1 ∂ Y 2 ∂ Y 2 + w 1 + w 2 , Y 1 − Y 2 ∂ x i ∂ x i ∂ x i ∂ x i ∂ w 1 Y 2 ∂ w 1 Y 2 − ∂ � ∂ w 1 Y 2 � + U i D T ∂ t ∂ x i ∂ x i ∂ x i C φ = w 2 S ( ξ 2 , Y 2 ) + τ w 1 w 2 [ Y 1 − Y 2 ] D T � ∂ Y 1 ∂ Y 1 ∂ Y 2 ∂ Y 2 � + w 1 + w 2 . Y 2 − Y 1 ∂ x i ∂ x i ∂ x i ∂ x i Only the progress variable Y has a source term S 94 / 34
Poly-dispersed particles � Particle Size Distribution n ( L ; x , t ) = R 3 n ( L , U p ; x , t ) d U p ρ p D 2 τ p St = ≈ τ f 18 µτ f Particles of different size behave very differently with a non-linear relation between size and relative velocity 95 / 34
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
Recommend
More recommend
