Closure of non-integer moments arising in multiphase flow phenomena

Chemical Engineering Science 75 (2012) 424–434

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Closure of non-integer moments arising in multiphase flow phenomena S. Haeri, J.S. Shrimpton n Tizard Building, School of Engineering Sciences, Highfield Campus, University of Southampton, Southampton S017 1BJ, UK

a r t i c l e i n f o

a b s t r a c t

Article history: Received 26 October 2011 Received in revised form 26 March 2012 Accepted 30 March 2012 Available online 9 April 2012

In this paper we consider closure problem for non-integer moments of form /xp S, p A R which are encountered in the moment evolution equations of many multiphase flow systems involving polydispersed particles, heat and mass transfer, chemical reactions, agglomeration or break-up. Two different categories of methods are considered: the first method is based on the reconstruction of the underlying PDF using Laguerre polynomials and the other is based on the direct calculation of noninteger moments using the fractional derivatives of moment generating function (MGF). By applying the results of fractional calculus an explicit equation is derived to express non-integer moments as a function of any arbitrary number of integer moments. The proposed methods are tested on several highly non-Gaussian analytical PDFs in addition to experimental agglomeration data and direct numerical simulation of fluid–particle turbulent multiphase flows. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Fractional moments Moment generating function Laguerre polynomial Closure problem Multiphase flow Eulerian–Eulerian turbulent models

1. Introduction The governing equations for many physical phenomena are precisely known, however numerical simulation by direct discretization of the governing equations is not generally feasible when complicated phenomena such as turbulence, chemical reaction, agglomeration and break-up are involved. The starting point to tackle such complicated phenomena is usually to write the averaged equations. The process starts by defining an arbitrary probability density function (PDF) for the process and deriving an equation for the evolution of the PDF. The PDF evolution equation is not an ordinary PDE with derivatives with respect to time and space but is a PDE in a phase space of several dimension which is at least 7 (i.e. time, space and velocity) and direct discretization of this equation is still not feasible. To reduce the PDF evolution equations, they can be integrated over the defined phase space resulting in a set of moment evolution equation. For a review on the general method and application to polydispersed turbulent flows see Haeri and Shrimpton (2011), Minier and Peirano (2001), and Minier et al. (2004) and the references therein. Rigopoulos (2010) also discusses similar methods from an engineering point of view using the population balance equation (PBE). In any case the equations contain unclosed terms that may also contain fractional moments for which closures are not trivial.

n

Corresponding author. Tel.: þ44 23 8059 4894. E-mail address: [email protected] (J.S. Shrimpton).

0009-2509/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2012.03.052

Non-integer moments can be introduced in the equations if the size distribution of particles retained in Eulerian field equations (moment evolution) (Archambault et al., 2003a,b; Haeri and Shrimpton, 2011; Minier and Peirano, 2001). Phenomena such as coagulation, nucleation, coalescence and breakup are present in the PBE (Diemer and Olson, 2002a,b,c, 2006; Frenklach and Harris, 1987) and heat and mass transfer considered (Beck and Watkins, 2002, 2003). There are some ad hoc solutions to the problem such as interpolation between the moments. Beck and Watkins (2002, 2003) used geometric interpolation and Frenklach and Harris (1987) used Laguerre interpolation between the logarithms of the integer moments. Note that interpolation between the moments is not possible for negative fractional moments since usually only equations for positive moments are available. Fractional negative moments can appear in the moment evolution equations and will be discussed in Section 4.5. PDF reconstruction methods are also attempted (Diemer and Olson, 2002a,b,c) which will be discussed in more detail below. In this paper we examine the problem from a mathematical point of view and consider the moments mp ¼ /xp S of a general univariate PDF given by Z mp ¼ xp PðxÞ dx, p A R, ð1Þ O

where O is the domain on which x is defined. Also note that in this paper we only consider the univariate PDFs. It is obvious that having the PDF any moment can readily be calculated, however in many physical phenomena it is easy to determine the moments but it is extremely difficult to determine the distributions themselves (John et al., 2007). In addition during a numerical

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

simulation the only information available are the solved variables and in this case the solved variables are the first few integer moments not the evolution equation of the PDF itself (Haeri and Shrimpton, 2011). Thus we require a method to enable us to estimate the real order moments mp , p A R, using limited number of integer moments mi , i ¼ 1; 2, . . . ,n. To accomplish this we will consider two general methods: a PDF reconstruction method based on Laguerre polynomials and a direct fractional method of moments (DFMMs) based on derivatives of the moment generating function (MGF). The first obvious approach is to try to reconstruct the PDF using integer moments mi which leads to the well-known finite moment problem (John et al., 2007) and can be regarded as a finite dimensional version of the Hausdorff moment problem (Tagliani, 1999; Inverardi et al., 2003, 2005; Pintarellia and Vericat, 2003) which is in general ill-posed lacking one or more conditions of a well posed problem (i.e. existence, uniqueness and stability) (Talenti, 1987). The reconstruction methods can be classified under linear and non-linear methods as discussed by Volpe and Baganoff (2003). They classified the maximum entropy method (MEM) as a non-linear method and reconstruction methods based on expansion around some parent distribution as linear methods. The MEM was first utilized by Koopman (1969) based on the idea of information theory of Shannon (1948). The MEM (Koopman) method is attractive for several reasons (Paris and Vencovska, 1997): (i) sound conceptual foundations; (ii) interdependence between even and odd order moments; (iii) nonnegative probabilities; (iv) produces the most unbiased distribution possible, i.e. if a distribution with less entropy (uncertainty) was used that would imply the existence of additional knowledge (Archambault et al., 2003a). The last property of MEM (i.e. most unbiased distribution) might be interesting in some applications but not for the current problem since MEM will strictly produce Gaussian distributions if for example two moments were used (Haeri and Shrimpton, 2011). PDF reconstruction based on some priori simple shape is the other option. However this method has a global realizability issue, in that an assumed PDF, valid at one location in space and time may evolve into a form that violates the assumption of the PDF at another position (Beck and Watkins, 2003). A more advanced method in this category is expansion using orthogonal polynomials. Two basic methods in this category which are extensively used in the literature are Gram–Charlier and Edgeworth series expansion (Blinnikov and Moessner, 1998). In both methods the PDF P is evaluated using a truncated expansion in terms of Hermite’s polynomials Hn(x) PðxÞ ¼

N X

C n Hn ðxÞN ðm, sÞ,

ð2Þ

n¼1

where N ðm, sÞ is a Gaussian distribution with parameters m and s, Cn are the coefficients containing the higher order moments (Gram–Charlier) or higher order cumulants (Edgeworth) and Hn is the Hermite polynomial which can explicitly be written (Abramowitz and Stegun, 1972) by Hn ðxÞ ¼ n!

bn=2c X

ð1Þk xn2k

k k ¼ 0 k!ðn2kÞ!2

:

ð3Þ

Majumdar et al. (2007) among others (Gaztanaga et al., 2000; Kendall et al., 1991) discussed that both methods can be divergent due to the fact that the series expansion is sensitive to the behavior of PðxÞ at infinity. For the series to be convergent PðxÞ 2 should fall faster to zero than ex =4 . Generalized Laguerre polynomial expansion is the more promising approach in this category which does not have any oscillatory or divergence problems of

425

Edgeworth or Gram–Charlier expansions (Majumdar et al., 2007). Although Laguerre polynomial expansion has been used in the literature to reconstruct different types of PDFs, to our knowledge its accuracy in estimating non-integer moments has not been tested and will be discussed in this paper. There are other methods which cannot be considered general such as spline reconstruction method proposed by John et al. (2007) which fits cubic splines to the available moments however for the method to be effective usually a large number of moments are required. The methods discussed above are all indirect methods i.e. we first reconstruct the PDF using one of the discussed methods then numerically integrate Eq. (1) to calculate the required moments. We will also formulate a direct method using fractional derivatives of the MGF. The property of the moment generating function is that its kth derivative calculated at s¼0 is equal to the kth moment of the distribution PðxÞ (Mood et al., 1974)  k d G mk ¼ k  : ð4Þ ds  s¼0

Now if we can extend this equation to non-positive non-integer values of k then we have a way of expressing non-integer moments as a function of integer moments. This extension is formalized in the field of fractional calculus and will be discussed in Section 3. Then we will discuss several test cases using moments of analytical PDFs, moments of agglomeration phenomena measured experimentally and moments extracted from direct numerical simulation of isotropic fluid–particle turbulent systems.

2. Generalized Laguerre polynomial expansion The Laguerre expansion of a PDF P can be written by (Majumdar et al., 2007; Mustapha and Dimitrakopoulos, 2010) PðxÞ ¼

1 X

r n Lan ðbxÞGðx; a þ 1,bÞ,

ð5Þ

n¼0

where Lan is the generalized Laguerre polynomial with parameter a and Gðx; a,bÞ is a Gamma distribution with parameters a and b given by a

Gðx; a,bÞ ¼

b xa1 exp ðbxÞ, GðaÞ

x Z0, a,b 40,

where GðaÞ is the gamma function given by Z 1 GðaÞ ¼ t a1 exp ðtÞ dt:

ð6Þ

ð7Þ

0

If the argument is an integer then GðaÞ ¼ ða1Þ!. The function can be approximated by directly calculating the integral which is not a stable method and a better approach would be the approximation method proposed by Cody (1976) which is adopted for the current study. Generalized Laguerre polynomials are orthogonal with respect to the measure G therefore using the Rodriguez formula (Lebedev, 1972), Ln can be defined by  n 1 d ½xn Gðx; a þ 1; 1Þ, ð8Þ Lan ðxÞGðx; a þ1; 1Þ ¼ n! dx and using Leibniz’s theorem for differentiation of product it can explicitly be written by (Lebedev, 1972)  n  X n þ a ðxÞm , ð9Þ Lan ðxÞ ¼ m! nm m¼0 where the binomial coefficients are generalized using the gamma function. The coefficients rn can be found using the orthogonality of the Laguerre polynomials. A sequence of polynomials Q 1 ðxÞ, . . . ,Q n ðxÞ are said to be orthogonal on the interval ½x1 ,x2 

426

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

with respect to a weight function W(x) if Z x2 Q n Q m WðxÞ dx ¼ 0 8m an:

3.2. RL fractional derivative ð10Þ

x1

For Laguerre polynomials it can easily be shown that Z 1 Gðn þa þ1Þ Lan ðbxÞLam ðbxÞGðx; a þ 1,bÞ dx ¼ dnm , n!Gða þ1Þ 0

ð11Þ

Introducing Eq. (9) into Eq. (12) and using the definition of G function and Eq. (1) we have  Z 1 n  X n þ a ðbxÞi n!Gða þ 1Þ dx PðxÞ rn ¼ Gðn þ a þ 1Þ 0 i! ni i¼0 ¼ n!Gða þ 1Þ

i

ðbÞ mi , i!ðniÞ! Gða þ iþ 1Þ i¼0

ð13Þ

where mi are the available moments of order i. The only remaining issue is the calculation of the parameters a and b which are defined such that the expansion (5) has the same mean and variance as the original distribution PðxÞ. A simple derivation of these parameters can be found in Mustapha and Dimitrakopoulos (2010) and are given by 2m21 m2 a¼ m2 m21

m1 and b ¼ : m2 m21

ð16Þ

a

where dnm ¼ 1 if m ¼n and is 0 otherwise. Now multiplying the right and left hand sides of Eq. (5) by Lam ðxÞ and using Eq. (11) the result is Z 1 n!Gða þ 1Þ rn ¼ PðxÞLan ðbxÞ dx: ð12Þ Gðn þ a þ 1Þ 0

n X

Consider the integral Z t 1 f ðtÞ ¼ f ðtÞ dt,

ð14Þ

and define Z t Z t1 2 dt1 f ðtÞ dt: f ðtÞ ¼

ð17Þ

a

a

It can be shown that, see Appendix A, Z t 2 ðttÞf ðtÞ dt, f ðtÞ ¼

ð18Þ

a

and by induction Z t 1 n ðttÞn1 f ðtÞ dt, f ðtÞ ¼ ðn1Þ! a

n Z1:

ð19Þ

This can be extended to non-integer values of n using the Gamma function Z t 1 p D f ðtÞ ¼ ðttÞp1 f ðtÞ dt, p 40: ð20Þ a t GðpÞ a The fractional derivative can then be defined by (Cottone and Di Paola, 2009; Schiavone and Lamb, 1990)  k Z t 1 d p D f ðtÞ ¼ ðttÞkp1 f ðtÞ dt, k1 rp r k: ð21Þ a t GðkpÞ dt a Note that this definition is not arbitrary and it is defined to be equivalent to the GL definition which is defined as an extension to backward differences. Having the differentiation extended to arbitrary order, it can be shown (Appendix B)

mp ¼ 1 Dp0 Gð0Þ, 3. Direct fractional method of moments In this section we formulate an approach to estimate fractional moments as a series of the integer moments. The MGF of a positive valued distribution can be defined by (Mood et al., 1974) Z 1 GðsÞ ¼ PðuÞesu du: ð15Þ 0

The main property of the MGF is that the moments of the original PDF PðxÞ can be derived from this function by Eq. (4). Extension of this equation to non-integer values is now discussed.

where mp is the moment of the density PðxÞ of order p. Although the RL definition of fractional derivatives and integrals is more convenient to mathematically link the fractional moments to fractional derivatives it does not provide a numerical method for calculating the moments. Therefore we revert to the GL definition to write an expression for fractional moments. 3.3. GL fractional derivatives GL definition of fractional derivatives is more intuitive and starts by observing the series of backward differences 0

f ðtÞ ¼ 3.1. Fractional derivatives and integrals

df f ðtÞf ðthÞ ¼ lim , dt h-0 h 2

The derivative of arbitrary real order p can be considered as an interpolation to the operators of a sequence of n-fold integration and nth order derivative. Fractional derivatives are presented by p a Dt ðtÞ where a and t are the limits related to the operation of the fractional differentiation and are commonly called terminals of the fractional differentiation. These are essential to avoid ambiguities in application of fractional derivative to real numbers (Hilfer, 2000; Podlubny, 1999; Ross, 1977). There are two equivalent approaches to the definition of fractional differentiation namely the Grunwald–Letnikov (GL) approach and the Rienmann–Liouville (RL) approach. It is customary (Podlubny, 1999) to use the RL formulation for problem setup and use GL approach to obtain a numerical solution. Following the same approach we use RL definition to establish the relation between the evaluation of the moments and fractional calculus and GL approach to provide an explicit equation for fractional moments as a function of integer moments.

ð22Þ

00

f ðtÞ ¼

d f dt

2

¼ lim

f ðtÞ2f ðthÞ þ f ðt2hÞ

3

000

f ðtÞ ¼

d f dt

3

2

h

h-0

¼ lim

,

f ðtÞ3f ðthÞ þ3f ðt2hÞf ðt3hÞ 3

h

h-0

,

and by induction f

ðnÞ

ðtÞ ¼

  n n d f 1 X r n ¼ lim ð1Þ f ðtrhÞ: n n h-0 h r dt r¼0

Now consider the following generalization:   n p 1 X ðpÞ ð1Þr f ðtrhÞ f h ðtÞ ¼ p r h r¼0

ð23Þ

ð24Þ

for arbitrary natural numbers p and n, such that p r n, we have p

lim f h ðtÞ ¼ f

h-0

p

ðpÞ

ðtÞ ¼

d f p, dt

ð25Þ

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

because all the coefficients in the numerator after ðppÞ are identically zero. Eq. (24) can be inverted (see Appendix C) to get   i X i r ðrÞ h f ðtÞ: f ðtihÞ ¼ ð1Þr ð26Þ r r¼0 Eq. (26) is the key for writing the non-integer moments as functions of integer moments and it will become clear in the course of this section. To extend the definition to negative values a careful definition of the upper bound of the summation is required otherwise the limit in Eq. (24) would be strictly zero for any n. Therefore by taking h ¼ ðtaÞ=n, a being a real number, and writing ðpÞ

lim f h

h-0 nh ¼ ta

ðtÞ ¼ a Dp t f ðtÞ:

ð27Þ

Note that Eq. (27) is actually the definition of an integral, for example writing the series for p ¼1 the limit is simply the Rt definition of a f ðtÞ dt. By the method of induction it is possible to show (Hilfer, 2000; Podlubny, 1999)   Z t n X p 1 p p ð1Þr ðttÞp1 f ðtÞ dt: f ðtrhÞ ¼ a Dt f ðtÞ ¼ lim h ðp1Þ! a r r¼0 h-0 nh ¼ ta

ð28Þ It can easily be shown (Podlubny, 1999) Z t Z tp1 Z t1 p dt1 dt2    f ðtp Þ dtp : a Dt f ðtÞ ¼ a a a |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð29Þ

p times

Therefore p a Dt f ðtÞ ¼

p

lim h

h-0 nh ¼ ta

n X

ð1Þr

  p

r¼0

r

f ðtrhÞ

ð30Þ

is indeed a general expression for p-fold integration and derivatives of order n. Eq. (30) can be extended to non-integer p values and by calculating the limits directly it can be shown that (Cottone and Di Paola, 2009; Podlubny, 1999) Z t 1 p ðttÞp1 f ðtÞ dt ð31Þ a Dt f ðtÞ ¼ GðpÞ a and p a Dt f ðtÞ ¼

ðrÞ m X f ðaÞðtaÞp þ r r¼0 Z t



Gðp þ r þ1Þ ðttÞmp f

ðm þ 1Þ

þ

1 Gðp þ m þ 1Þ

ðtÞ dt

8m 4 p1:

ð32Þ

a

Note the equivalence between Eqs. (20) and (31). The equivalence between Eqs. (32) and (21) is harder to note but direct differentiation and integration by parts shows that both definitions are indeed equivalent (Hilfer, 2000). Having the framework established we will derive an equation to explicitly write the noninteger moments as functions of integer moments in the next section. 3.4. Estimating the non-integer moments A first order approximation to p-order derivative using the GL definition can be written as   bðtaÞ=hc X p p p p p ð1Þr f ðtrhÞ: ð33Þ a Dt f ðtÞ ¼ lim h a Dt f ðtÞ  h h-0 r r¼0 The number of addends in Eq. (33) becomes very large for t ba. This series can be truncated using the short memory principle; taking into account the behavior of f(t) only in the recent past. This means that the truncated series (33) is in particular a very good approximation for calculating fractions

427

that are approximately equal to the number of integer moments retained in the expansion and the error can be quantified as suggested by Deng (2007). In Eq. (33), f ðtrhÞ can be considered as the MGF but since the function is not explicitly available the equation is not usable in this form. Eq. (26) provides an equation for this term, now inserting Eq. (26) into Eq. (33) noting the definitions (22) and (4) we get !  X N r X r p jp mp  ð1Þr ð1Þj ð34Þ h mj , j r r¼0

j¼0

where N is the number of moments retained in the expansion. We also acknowledge that a similar equation is recently derived by Gzyl and Tagliani (2010) using Taylor expansions. However they applied the series to very large number of integer moments Oð100Þ which is not applicable to the problems considered in this paper. The only remaining issues are to provide an equation to   p calculate the coefficient wpr ¼ ð1Þr and propose a step size h. r p The coefficients, wr , are generalizations of binomial coefficients to non-integer values and can be calculated using the definition of the gamma function. However another possible approach is to use the recursive relations (Podlubny, 1999)   p þ1 wpr1 , r ¼ 1; 2, . . . , wp0 ¼ 1, wpr ¼ 1 ð35Þ r which eliminates the need for the evaluation of gamma functions and provides better computational efficiency and stability. In addition it provides a unique opportunity for higher order estimations of the fractional derivative which is discussed next. Considering the coefficients of power series for the function ð1zÞp ð1zÞp ¼

1 X

wpr zr :

ð36Þ

r¼0

Substituting z ¼ eiy coefficients can be expressed in terms of Fourier transforms Z 2p 1 ð1eiy Þp eiky dy: ð37Þ wpr ¼ 2pi 0 Since we always use only a finite number of moments any FFT library can be used to calculate the coefficients using Eq. (37). Eqs. (35) and (37) are only first order approximations however higher order approximations can be constructed using higher order polynomials. Lubich (1986) suggests polynomials up to sixth order which can be used in conjunction with Eq. (37) and an FFT library to find coefficients for the higher order approximation to the fractional derivatives. However we only use Eq. (35) in this study due to its simplicity and efficiency. The step size h should ideally be very small. If we consider 1 H ¼ h in Eq. (34) prescribing a value for H is equivalent to the width of the PDF considered in the calculations. Gzyl and Tagliani (2010) showed that the series is always convergent for 1=2 oh r 1 or 1 r H o 2 however as will be discussed in the next section this value causes a severe under-estimation when only a few moments are used. It seems to be reasonable to relate H to the mean and the variance of the available data H ¼ m þ ls,

ð38Þ

where m is the mean of the data, s is the standard deviation and l is an adjustable parameter. It seems to be valid that if the fractional moment of interest is not much larger than the number of moments retained in the series setting l ¼ 1=ðp þ1Þ provides an accurate estimate. This will be elaborated further in the next section when analyzing the results where we explore the hypothesis that l ¼ 1=ðp þ1Þ is reasonable.

428

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

4. Results and discussion 4.1. Log-normal distribution We test the results on analytic log-normal distribution given

Laguerre polynomial approach. Intentionally only the first two moments (normalized on [0,12]) are used to test the ability of the two proposed methods in estimating the non-integer moments using only limited data. In this paper all errors are calculated using

by 2 2 1 PðxÞ ¼ pffiffiffiffiffiffi eðlog xm Þ =2s , 2psx

Err ¼

where s and m are scale and location parameters respectively. Higher order moments are analytically given by 2

mj ¼ eðð1=2Þj

s2 þ jm Þ :

9ma me 9

ð39Þ

ð40Þ

Fig. 1 shows the simulation of the log-normal distribution using the Laguerre polynomials where the PDF is reconstructed for 200 equidistant points in the range ½0; 12. Although our goal is to simulate the non-integer moments the accuracy of the calculated moments are directly proportional to the accuracy of the fit. Generally a Laguerre polynomial is able to produce a very good fit to this distribution with a limited number of moments because of the similarity between the shape of gamma and log-normal distributions. Fig. 2 shows the fractional moments mp , p ¼ f1:2,2:3,3:4,4:5g calculated using the fractional moments approach and the

ma

,

ð41Þ

where subscripts a and e stand for the analytical and estimated values respectively and also %Err ¼ 100  Err. In Fig. 2 error bars, for 7100  Err, are also provided with maximum error of 27:5% in calculating the m4:5 using the fractional moments method with all other errors being less than 5%. The reason for this large error will be discussed in the next section and a remedy will be suggested.

4.2. Mixture of normal distributions A mixture of normal distributions can be created using PðxÞ ¼

N X i¼1

wi N ðmi , si Þ,

N X

wi ¼ 1,

ð42Þ

i1

where N ðm, sÞ is the normal distribution with mean m and variance s. We consider the mixture distribution with two terms PðxÞ ¼ að0:8N ð0:5,1Þ þ 0:2N ð4; 1ÞÞ1½0;8 ,

Fig. 1. A log-normal distribution with scale and position parameters 0.5 and 1 can accurately be reproduced using small number of moments.

Fig. 2. Estimation of non-integer moments using LPM and DFMM. Error bars are presented for each calculation.

ð43Þ

where 1A is the identity function which is equal to 1 if x A A and is zero otherwise. The parameters are chosen such that the final mixture be a bimodal distribution, see Behboodian (1970) and Schilling et al. (2002) for the sufficient conditions. The normalR8 ization constant, a¼1.327743884718793, ensures 0 PðxÞ dx ¼ 1. First 12 moments of this function are given in Table 1 which are calculated by direct integration using adaptive quadrature method (Shampine, 2008) with relative and absolute tolerance of 1e8 and 1e15 respectively. Note that the PDF is normalized and therefore m0 ¼ 1. Fig. 3 shows the simulation of the mixture distribution using the Laguerre polynomials for 200 discrete points. More moments are needed in this example to correctly capture the tail of the distribution. Table 2 shows the fractional moments of the mixture Gaussian distribution using Laguerre polynomials method (LPM) and DFMM. DFMM estimates the fractional moments particularly well as long as the number of integer moments is near the value of the fractional moments. For example using first three integer moments the estimated values for m1:2 , m2:3 , m3:4 are particularly precise with maximum error of 5:2% for m3:4 which are better than the estimates provided by LPM. However for m4:5 a large error is detected using the DFMM. Same calculations are performed using the first 5 integer moments and the results are listed in Table 3. In this example very precise estimates are provided using DFMM which are all better than those calculated using LPM. In the LPM case increasing the number of moments to five actually increases the error. We performed the test with seven and nine moments and detect a mild oscillatory convergence which is not reported in other studies using Laguerre series. On the other hand DFMM is based on a firm mathematical ground with predictable behavior which is a direct consequence of the short memory principle discussed in Section 3.4. Parameter l can easily be adjusted to provide better results, for example setting l ¼ 2=ðp þ1Þ results in m4:5 ¼ 209:9239 with %Err ¼ 1:57, see Table 3.

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

429

Table 1 First 12 moments of mixture Gaussian distribution. Moment

Value

Moment

Value

m1 m3 m5 m7 m9 m11

1.803323915192423 22.21194221994151 467.7931304312276 12 010.60922300377 350 667.3348484823 11 345 464.98539287

m2 m4 m6 m8 m10 m12

5.618825039917679 98.56381003061070 2326.755938195312 63 991.84370378382 1 971 324.861855061 66 731 974.28785729

Fig. 3. A mixture normal distribution (Eq. (43)) is reconstructed using 8 and 12 moments.

Table 2 Fractional moments of mixture Gaussian distribution estimated using first three integer moments (m0 , m1 , m2 ). Moment

Value

LPM

DFMM

LPM %Err

DFMM %Err

m1:2 m2:3 m3:4 m4:5

2.1974 8.3418 39.9035 213.2800

2.1093 7.5806 37.1028 219.6070

2.1333 8.5486 37.8250 140.6065

4.0093 9.1249 7.0187 2.9665

2.9171 2.4793 5.2088 34.0742

Table 3 Fractional moments of mixture Gaussian distribution estimated using first five integer moments (m0 , m1 , m2 , m3 , m4 ). Moment

Value

LPM

DFMM

LPM %Err

DFMM %Err

m1:2 m2:3 m3:4 m4:5

2.1974 8.3418 39.9035 213.2800

2.2203 8.7174 45.4930 271.8476

2.1804 8.3543 39.9454 209.9239

1.0431 4.5029 14.0075 27.4604

0.7741 0.1501 0.1050 1.5736

4.3. Rice-Nakagami distribution A Rice-Nakagami distribution can be written by (Majumdar and Gamo, 1982) Pðx; Ic , sÞ ¼

1

s2

exp

 

x þ Ic

s2

  pffiffiffiffiffiffi xIc I0 2 2 ,

s

ð44Þ

where s is the standard deviation and Ic ¼ ms2 , with m being the mean. I0 is the zeroth order modified Bessel function of the first kind. Higher order moments are analytically given by

Fig. 4. Rice-Nakagami distribution reconstructed on 200 points using three and five moments.

(Majumdar, 1984)     I I mj ¼ s2j exp  c2 Gðj þ 1ÞM j þ 1; 1, c2 ,

s

s

ð45Þ

where M is the confluent hyper-geometric function (Kummar function). Parameters Ic and s2 are set to 1 and 1.9 respectively. Fig. 4 shows the reconstructed PDF (Eq. (44)) in ½0; 12 for 200 discrete points with the normalization constant calculated to be 0.99994848708433 on this interval. Evidently a Laguerre expansion can capture the features of this distribution with very high accuracy even with very limited number of moments. Note that here we are only using the first three and five moments to produce Fig. 4 whereas in Figs. 1 and 3 we used a larger set of moments. Fig. 5 shows the exact values of the fractional moments and the values calculated using LPM and DFMM methods. Since this distribution can accurately be reconstructed using LPM, very accurate estimations for the fractional moments can be achieved with errors never exceeding 2.5%. Despite this DFMM still produces acceptable results with the maximum error of 22% in calculating m2:3 and error for m1:2 , m3:4 , m4:5 being in the same range as those calculated by LP method. We also attempted the estimation of the fractional moments using a constant value as suggested by Gzyl and Tagliani (2010). The results are plotted in Fig. 5 and the only reliable results are those calculated between the available integer moments, i.e. m1:2 , m2:3 . The values for m3:4 , m4:5 are severely under-estimated and practically unusable. Note also that in this example we opt for the maximum of the suggested interval, i.e. H¼2, using smaller values causes even larger under-estimations. However it should be stated that the adaptive value for H is selected to work with very small number of integer moments, 3 r N r5, which has practical applications, for example, when solving Eulerian field equations. Therefore if a larger number of integer moments are available in the model one should use a more conservative value

430

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

Fig. 5. Fractional moments of the Rice-Nakagami distribution calculated using LPM and DFMM methods using both an adaptive and constant value for H.

Table 4 First 11 moments of particle mass distribution (PMD) measured by Oncul et al. (2005).

Fig. 6. Reconstruction of PMD by Laguerre polynomials using 11 moments using experimental measurements of Oncul et al. (2005).

Table 5 Fractional moments of experimental PMD estimated using first three integer moments (m0 , m1 , m2 ).

Moment

Value

Moment

Value

Moment Value

LPM

DFMM

m0 m2 m4 m6 m8 m10

1.028571810801850e  02 3.731662761082810e  08 2.961317428763996e  13 3.106813835529885e  18 3.973765578013985e  23 5.901352427828747e  28

m1 m3 m5 m7 m9 m11

1.627619528771249e  05 1.005399673163290e  10 9.324232740417577e  16 1.087507001932166e  20 1.507036647558311e  25 NA

m1:2 m2:3 m3:4 m4:5

4.4901E  06 5.6016E  09 8.2381E  12 1.0276E  14

4.6738E  06 3.5979 6.3647E  09 8.0759 1.0012E  11 11.5492 1.8520E  14 33.9737

for H or include the number of integer moments, N, in the definition of l. 4.4. Crystallization process: experimental data In this section experimental data of Oncul et al. (2005), see also John et al. (2007), for a crystallization process are used to further test the two proposed methods. The moments of the particle mass distribution (PMD) are provided in Table 4. Fig. 6 shows the reconstructed PMD using 11 moments by Laguerre polynomials which cannot capture the first sharp peak accurately. Only the experimental moments are provided in John et al. (2007) and Oncul et al. (2005) and not an equation for the PMD. We need the PDM to analytically calculate the fractional moments and compare them to the values estimated either by the LPM or DFMM. The exact values of the fractional moments are calculated using the following graphical procedure: the graph of the data is sampled using 103 points, these digitized points are then resampled using cubic spline interpolation to produce a fine uniform sample of 5000 points and then integrated to get the fractional moments. Comparison of the integer moments available in the paper with recalculated integer moments from the abovementioned process shows that the error is always less than 1% for the first five moments and therefore we use the same process to calculate the fractional moments m1:2 , . . . , m4:5 assuming the calculated values from the abovementioned integration procedure to be exact. Table 5 shows the results of fractional moments calculation using only first three moments. The LP method produces large errors for m4:5 since the method is not able to capture the sharp peak in the distribution. DFMM still produces acceptable results which are better than the estimates produced by LP method. Note that l ¼ 1=ð1 þ pÞ is used for p A 1:2,2:3,3:4 and l ¼ 2=ð1 þpÞ for p ¼4.5.

4.6577E  06 6.0937E  09 9.3137E  12 1.5563E  14

LPM %Err DFMM %Err 0.3461 4.4469 7.4976 18.9984

4.5. Polydispersed fluid–particle system: DNS simulation A pseudo-spectral DNS code capable of two-way coupled fluid–particle, homogeneous turbulence simulations, ‘PANDORA’ (Scott, 2006; Scott et al., 2009), is used in this section to provide the particle distribution data. Two sample simulations were used with 323 and 1283 grid points equivalent to a Taylor Reynolds number, Rel , of 24.24 and 83.40 respectively. More details on these simulation can be found in Scott (2006) and Scott et al. (2009). Es is the enstrophy sampled along the particle path and /eS is the mean fluid enstrophy. We will not discuss the problem and related simulations in more detail here and the reader can simply consider the PDFs discussed in this section as some form of conditional distribution of a globally uniform particle size distribution. The interested reader can consult Aliseda et al. (2002), Eaton and Fessler (1994), Squires and Eaton (1990), and Sundaram and Collins (1997) for further detail on the physical phenomena that causes a deviation in the conditional PDFs from globally uniform distribution. Despite the very simple shape of these PDFs using traditional methods such as MEM to reconstruct them is almost hopeless due to the fact that they do not tend to zero near the limits of the phase space. Fig. 7 shows the reconstruction results of the diameter distribution of a globally uniform diameter size distribution conditioned on the samples enstrophy for 323 and 1283 runs. Eleven moments are used to generate these figures. Although LPM can capture the general behavior of the conditional particle size distribution, it does not provide a good fit specially near both ends of the distribution even with a large set of integer moments m1    m11 . Another issue is that even if we had very larger number of moments we were still restricted by the truncation errors since usually higher order moments are extremely small (large) which are multiplied by large (small) coefficients in the series. This results in error that can be as large as all the significant digits in the calculation. Double precision

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

431

Fig. 8. Fractional moments ðm1:2 , m2:3 , m3:4 , m4:5 Þ calculated using first three and five moments using LPM and DFMM for 323 simulation with 2/eS o Es o 3/eS. In all cases better accuracy for DFMM is evident. Fig. 7. Conditional PDF reconstructed using 11 moments, (a) 323 and (b) 1283 simulation.

arithmetic was sufficient for all the calculations in this paper however if larger number of moments are required one can consider using arbitrary precision libraries such as the GMP (Granlund, 2011). Figs. 8 and 9 summarize the fractional moment calculations using both methods for 323 simulation. In this section exact values are directly calculated from the DNS data and relative errors are calculated using these values. Similarly l ¼ 1=ðp þ 1Þ is used for p A f1:2,2:3,3:4g and 2=ð1 þpÞ used for p ¼4.5. LPM produces errors that can be as large as 45% when only first three moments ðm0 , m1 , m2 Þ are used, which is not surprising since the PDF fit is not accurate. DFMM produces accurate results even when the first three moments are used with errors that never exceeds 20%. Adding two more moments and DFMM produces extremely accurate results with errors all under 1% but LPM still produces errors as large as 30%. Fig. 10 shows the same calculation using the 1283 simulation. We only present the results for 4/eSo Es o5/eS which shows the same trend as the 323 simulation with the errors in DFMM results never exceeding 10% in this case, even for estimation with three moments. Another important case is the estimation of non-integer negative moments that can be encountered in polydispersed Eulerian–Eulerian models (Haeri and Shrimpton, 2011). To test the ability of the proposed methods in calculation of

non-integer negative moments we use p A f0:5,1:5,2:5g. We use the same functional form for the l but since p is negative in this case it is more convenient to use l ¼ 1=ð9p9 þ 1Þ. It is also possible to write l as a function of both the number of retained integer moments and fractional order to guarantee both convergence properties and accuracy of the series as discussed in Section 4.3. The results are summarized in Table 6. DFMM generally produces better estimates with results that are in particular accurate for p 4 2. However one of the advantages of the DFMM is the possibility of adjusting the value of the free parameter l for specific applications to yield exceptionally accurate values. It is worth mentioning that having an adjustable parameter is usually not considered as an advantage since it is not possible in practice to compare the calculated results to a known solution. However since this series is derived analytically the parameter H (or l) only depends on the truncation of the series. This is directly related to the number of integer moments retained in the calculation but not on a specific form of the PDF or application. Therefore one can tune the parameter on standard test cases for a specific range of available moments and safely use it for any problem. For example l ¼ 1=ðp þ 1Þ suggested in this paper can safely be used for any problem as long as p rN þ 1 and for slightly large values up to p rN þ 2, l ¼ 2=ðp þ 1Þ is a good choice. This range of available integer moments to our knowledge covers a very broad range of engineering applications.

432

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

Fig. 9. Fractional moments ðm1:2 , m2:3 , m3:4 , m4:5 Þ calculated using first three and five moments using LPM and DFMM for 323 simulation with 4/eSo Es o 5/eS. In all cases better accuracy for DFMM is evident.

5. Conclusion In this paper the non-integer moments closure problem which is encountered in many simulations concerning multiphase flows is considered. We considered two methods, a direct method and a method based on reconstruction of the underlying PDF. Different PDF reconstruction methods were considered and LPM is suggested due to its non-oscillatory behavior and better convergence properties compared to other orthogonal polynomial fitting methods. We also formulate a direct method based on the results of fractional calculus and the fact that the MGF can also generate the non-integer moments. Several problems are considered, including analytical, experimental and numerical simulations. All the test cases involve highly non-Gaussian distributions and DFMM method produced better results in almost all the test cases. DFMM method involves a free parameter that should be adjusted for different problems. This parameter is directly related to the assumptions made in the mathematical derivation of DFMM equation and has simple mathematical explanation similar to the step size used in the finite difference approximations. A simple general expression is provided for the free parameter l which produces accurate results in our test cases, however it can be adjusted for different problems to yield higher accuracies. Generally both methods can be considered general with DFMM possessing more desirable properties, which in addition

Fig. 10. Fractional moments ðm1:2 , m2:3 , m3:4 , m4:5 Þ calculated using first three and five moments using LPM and DFMM for 1283 simulation. Maximum error never exceeds 10% for DFMM.

Table 6 Fractional moments of DNS simulation estimated using first three and five integer moments including m0 . No. integer moments

Moment Value

3

m0:5 m1:5 m2:5 m0:5 m1:5 m2:5

5

LPM

DFMM

LPM %Err

4.3310 4.0025 4.6390 7.58 118.4843 90.0237 105.9396 24.20 5260.7300 2643.1820 3139.4362 49.76 4.3310 4.1165 4.7404 4.95 118.4843 107.2110 119.72330 9.51 5260.7300 4108.2148 4677.2128 21.90

DFMM %Err 7.11 10.58 40.32 9.45 1.05 11.09

to its higher accuracy observed in our test cases, can be summarized as follows:

 No intermediate PDF is involved. This is very important



property since one should always consider the possibility of producing non-realizable PDFs (negative probabilities) when reconstructing a PDF. There is no guarantee that LPM always produces realizable PDFs specially near the tails however it is easy to check for this issue and a simple spline fit to the small negative region can correct the behavior. DFMM is computationally much more efficient than any indirect method involving an intermediate PDF. If an indirect method is used in addition to reconstructing the intermediate PDF one also needs to perform numerical integration to

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434



calculate the value of the non-integer moment whereas in DFMM method a simple explicit equation can be derived which can be implemented in less than 100 lines of code. Coordinate shifts are needed if we let the random variable assume negative values in LPM whereas no such manipulations are needed for DFMM.

However it should be noted that since in the LPM method we are merely using series expansions, extension to higher dimensions is straightforward but in the case of the DFMM method a rigorous extension of the theory is required.

In this section we prove Eq. (18) which is used to extend the nfold integration to arbitrary orders. We start by setting Rt hðt1 Þ ¼ a 1 f ðtÞ dt and g 0 ðt1 Þ ¼ 1 then Z t Z t Z t1 dt1 f ðtÞ dt ¼ hðt1 Þg 0 ðt1 Þ dt1 a a a Z t Z t Z t t ¼ hðt1 Þgðt1 Þ9a  t1 f ðt1 Þ dt1 ¼ t f ðtÞ dt tf ðtÞ dt a

Z

a

ðttÞf ðtÞ dt:

ðA:1Þ

a

It is now easy to use Eq. (A.1) n times to show by induction that Eq. (19) holds for any integer n.

Appendix B. Generation of arbitrary order moments from the MGF First note that using Eq. (21) with n ¼ kp we have  k Z t 1 d p ðttÞn1 f ðtÞ dt: 1 Dt f ðtÞ ¼ GðnÞ dt 1

ðB:1Þ

Using Z ¼ tt  k Z 1 1 d p D f ðtÞ ¼ Zn1 f ðtZÞ dZ 1 t GðnÞ dt 0  k Z 1 1 d ¼ Zn1 f ðtZÞ dZ dt GðnÞ 0  k Z t 1 d ¼ ðttÞn1 f ðtÞ dt: dt GðnÞ 1

ðB:2Þ

The last line of Eq. (B.2) is the Caputo (1967) definition of the fractional derivative with a-1, i.e. C1 Dpt Df ðtÞ. The RL and Caputa definitions are not generally equivalent however Eq. (B.2) shows that both definitions are equivalent for the limit a-1. Now we can use (Cressie and Borkent, 1986) p C ct 1 Dt De

¼ cp ect

8t, 1o t o1, 8c 40, p ct 1 Dt De

ðB:3Þ

p ct

¼ c e . Assuming G(s) is analytic in and consequently ð1,0 and defining n ¼ kp, one can use Eq. (21) and change the order of integration to get

p 1 1 Dt GðsÞ9s ¼ 0 ¼ GðnÞ

Z

1

GðnÞ

¼ Z

"

0 1

¼ 0

Z ¼

0

1

d

k k

ds   p su  D ðe Þ 1 s 

Z

d

k

ds

p

Z

s

ðstÞn1

Z 0

1

1

   etu PðuÞ du dt 

#   ðstÞn1 etu dt  PðuÞ du  1 s¼0  Z 1  PðuÞ du ¼ up esu  PðuÞ du

s¼0

s

s¼0

1

k

u PðuÞ du ¼ mp :

0

Eq. (24) can be inverted by expanding the series for the first few terms f

ð2Þ

ðtÞ ¼ h

2

f ðtÞ2h

2

f

ð1Þ

ðtÞ ¼ h

1

f ðtÞh

f

ð0Þ

ðtÞ ¼ f ðtÞ:

1

2

f ðthÞ þ h

f ðt2hÞ,

ðC:1Þ

f ðthÞ,

ðC:2Þ ðC:3Þ 1

Introducing Eq. (C.3) into Eq. (C.2) to eliminate h ð1Þ

ðtÞ,

f ðtÞ we get ðC:4Þ

similarly by introducing Eqs. (C.3) and (C.2) into Eq. (C.1) to eliminate first two terms on the RHS we have f ðt2hÞ ¼ f ðtÞ2hf

ð1Þ

2 ð2Þ

ðtÞ þ h f

ðtÞ:

ðC:5Þ

Eq. (26) easily follows by induction. References

a

t

¼

Appendix C. Inversion of the GL series

f ðthÞ ¼ f ðtÞhf

Appendix A. Proof of n-fold integration formula

433

s¼0

ðB:4Þ

Abramowitz, M., Stegun, I., 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications. Aliseda, A., Cartellier, A., Hainaux, F., Lasheras, J., 2002. Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence. J. Fluid Mech. 468, 77–105. Archambault, M.R., Edwards, C.F., McCormack, R.W., 2003a. Computation of spray dynamics by moment transport equations. I. Theory and development. Atomization Sprays 13, 63–87. Archambault, M.R., Edwards, C.F., McCormack, R.W., 2003b. Computation of spray dynamics by moment transport equations. II. Application to calculation of a quasi-one-dimensional spray. Atomization Sprays 13, 89–102. Beck, J., Watkins, A., 2002. On the development of spray submodels based on droplet size moments. J. Comput. Phys. 182, 586–621. Beck, J.C., Watkins, A.P., 2003. The droplet number moments approach to spray modelling: the development of heat and mass transfer sub-models. Int. J. Heat Fluid Flow 24, 242–259. Behboodian, J., 1970. On the modes of a mixture of two normal distributions on the modes of a mixture of two normal distributions. Technometrics 12, 131–139. Blinnikov, S., Moessner, R., 1998. Expansions for nearly Gaussian distributions. Astron. Astrophys. Suppl. Ser. 130, 193–205. Caputo, M., 1967. Linear models of dissipation whose q is almost frequency independent—II. Geophys. J. R. Astron. Soc. 13, 529–539. Cody, W., 1976. An overview of software development for special functions. Lecture Notes in Mathematics, vol. 506; 1976, pp. 38–48. Cottone, G., Di Paola, M., 2009. On the use of fractional calculus for the probabilistic characterization of random variables. Probab. Eng. Mech. 24, 321–334. Cressie, N., Borkent, M., 1986. The moment generating function has its moments. J. Stat. Plann. Inference 13, 337–344. Deng, W., 2007. Short memory principle and a predictor corrector approach for fractional differential equations. J. Comput. Appl. Math. 206, 174–188. Diemer, R., Olson, J., 2006. Bivariate moment methods for simultaneous coagulation, coalescence and breakup. J. Aerosol Sci. 37, 363–385. Diemer, R.B., Olson, J.H., 2002a. A moment methodology for coagulation and breakage problems. Part 1. Analytical solution of the steady-state population balance. Chem. Eng. Sci. 57, 2193–2209. Diemer, R.B., Olson, J.H., 2002b. A moment methodology for coagulation and breakage problems. Part 2. Moment models and distribution reconstruction. Chem. Eng. Sci. 57, 2211–2228. Diemer, R.B., Olson, J.H., 2002c. A moment methodology for coagulation and breakage problems. Part 3. Generalized daughter distribution functions. Chem. Eng. Sci. 57, 4187–4198. Eaton, J., Fessler, J., 1994. Preferential concentration of particles by turbulence. Int. J. Multiphase Flow 20, 169–209. Frenklach, M., Harris, S.J., 1987. Aerosol dynamics modeling using the method of moments. J. Colloid Interface Sci. 118, 252–261. Gaztanaga, E., Fosalba, P., Elizalde, E., 2000. Gravitational evolution of the largescale probability density. Astrophys. J. 539, 522–531. Granlund, T., 2011. The GNU Multiple Precision Arithmetic Library, 5th edition. Gzyl, H., Tagliani, A., 2010. Hausdorff moment problem and fractional moments. Appl. Math. Comput. 216, 3319–3328. Haeri, S., Shrimpton, J., 2011. A mesoscopic description of polydispersed particle laden turbulent flows. Prog. Energy Combust. Sci 37, 716–740. Hilfer, R., 2000. Fractional Calculus in Physics. World Scientific. Inverardi, P., Pontuale, G., Petri, A., Tagliani, A., 2003. Hausdorff moment problem via fractional moments. Appl. Math. Comput. 144, 61–74.

434

S. Haeri, J.S. Shrimpton / Chemical Engineering Science 75 (2012) 424–434

Inverardi, P., Pontuale, G., Petri, A., Tagliani, A., 2005. Stieltjes moment problem via fractional moments. Appl. Math. Comput. 166, 664–677. John, V., Angelov, I., Oncul, A.A., Thevenin, D., 2007. Techniques for the reconstruction of a distribution from a finite number of its moments. Chem. Eng. Sci. 62, 2890–2904. Kendall, M., Stuart, A., Ord, J., 1991. Kendall’s advanced theory of statistics. Distribution Theory, vol. 1. John Wiley & Sons Inc. Koopman, B.O., 1969. Relaxed motion in irreversible molecular statistics. Stochastic Processes Chem. Phys. 15, 37–63. Lebedev, N., 1972. Special Functions and Their Applications. Dover Publications. Lubich, C., 1986. Discretized fractional calculus. SIAM J. Math. Anal. 17, 704–719. Majumdar, A., 1984. Uniqueness of statistics derived from moments of irradiance fluctuations in atmospheric optical propagation. Opt. Commun. 50, 1–7. Majumdar, A.K., Gamo, H., 1982. Statistical measurements of irradiance fluctuations of a multipass laser beam propagated through laboratory-simulated atmospheric turbulence. Appl. Opt. 21, 2229–2235. Majumdar, A.K., Luna, C., Idell, P., 2007. Reconstruction of probability density function of intensity fluctuations relevant to free-space laser communications through atmospheric turbulence. In: Proc. SPIE., http://dx.doi.org/10.1117/12. 728699. Minier, J.-P., Peirano, E., 2001. The PDF approach to turbulent polydispersed twophase flows. Phys. Rep. 352, 1–214. Minier, J.-P., Peirano, E., Chibbaro, S., 2004. PDF model based on Langevin equation for polydispersed two-phase flows applied to a bluff-body gas–solid flow. Phys. Fluids A 16, 2419–2431. Mood, A., Graybill, F.A., Boes, D., 1974. Introduction to the Theory of Statistics. McGraw-Hill. Mustapha, H., Dimitrakopoulos, R., 2010. Generalized laguerre expansions of multivariate probability densities with moments. Comput. Math. Appl. 60, 2178–2189. Oncul, A., Thevenin, D., Elsner, M., Seidel-Morgenstern, A., 2005. Numerical analysis of the preferential crystallization of enantiomers. In: 11th Workshop on Two-Phase Flow Predictions.

Paris, J., Vencovska, A., 1997. In defense of the maximum entropy inference process. Int. J. Approx. Reasoning 17, 77–103. Pintarellia, M., Vericat, F., 2003. Generalized Hausdorff inverse moment problem. Phys. A: Stat. Mech. Appl. 324, 568–588. Podlubny, I., 1999. Fractional Differential Equations. Academic Press. Rigopoulos, S., 2010. Population balance modelling of polydispersed particles in reactive flows. Prog. Energy Combust. Sci. 36, 412–443. Ross, B., 1977. Fractional calculus: an historical apologia for the development of a calculus using differentiation and antidifferentiation of non integral orders. Math. Mag. 50, 115–122. Schiavone, S., Lamb, W., 1990. A fractional power approach to fractional calculus. J. Math. Anal. Appl. 149, 377–401. Schilling, M., Watkins, A., Watkins, W., 2002. Is human height bimodal? Am. Stat. 56, 223–229. Scott, S.J., 2006. A PDF Based Method for Modelling Polysized Particle Laden Turbulent Flows Without Size Class Discretisation. Ph.D. Thesis. Imperial College London. Scott, S.J., Karnik, A.U., Shrimpton, J.S., 2009. On the quantification of preferential accumulation. Int. J. Heat Fluid Flow 30, 189–795. Shampine, L., 2008. Vectorized adaptive quadrature in matlab. J. Comput. Appl. Math. 211, 131–140. Shannon, C.E., 1948. A mathematical theory of communication. Bell Syst. Tech. J. 27, 379–623. Squires, K., Eaton, J., 1990. Preferential concentration of particles by turbulence. Phys. Fluids 3, 1169–1178. Sundaram, S., Collins, L., 1997. Collision statistics in an isotropic particle-laden turbulent suspension. Part 1. Direct numerical simulations. J. Fluid Mech. 335, 75–109. Tagliani, A., 1999. Hausdorff moment problem and maximum entropy: a unified approach. Appl. Math. Comput. 105, 291–305. Talenti, G., 1987. Recovering a function from a finite number of moments. Inverse Probl. 3, 501–517. Volpe, E., Baganoff, D., 2003. Maximum entropy PDFs and the moment problem under near-Gaussian conditions. Probab. Eng. Mech. 18, 17–29.