High Order Moment Model for Polydisperse Evaporating Sprays Towards Interfacial Geometry Descriptionthanks: This work was funded by a PhD Grant for M. Essadki by IFPEn.

Mohamed Essadki IFP Energies nouvelles, 1-4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex and Laboratoire EM2C, CNRS, CentraleSupélec,Université Paris-Saclay, Grande Voie des Vignes, 92295 Châtenay-Malabry, Cedex - France (mohamed.essadki@centralesupelec.fr).    Stephane de Chaisemartin IFP Energies nouvelles, 1-4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex - France (stephane.de-chaisemartin@ifpen.fr).    Frédérique Laurent CNRS, Laboratoire EM2C, CNRS, CentraleSupélec,Université Paris-Saclay, Grande Voie des Vignes, 92295 Châtenay-Malabry, Cedex - France and Fédération de Mathématiques de l’Ecole Centrale Paris, FR CNRS 3487 (frederique.laurent@centralesupelec.fr).    Marc Massot Laboratoire EM2C, CNRS, CentraleSupélec,Université Paris-Saclay, Grande Voie des Vignes, 92295 Châtenay-Malabry, Cedex - France and Fédération de Mathématiques de l’Ecole Centrale Paris, FR CNRS 3487 (marc.massot@centralesupelec.fr).

In this paper we propose a new Eulerian modeling and related accurate and robust numerical methods, describing polydisperse evaporating sprays, based on high order moment methods in size. The main novelty of this model is its capacity to describe some geometrical variables of the droplet-gas interface, by analogy with the liquid-gas interface in interfacial flows [Drew_Geom]. For this purpose, we use fractional size-moments, where the size variable is taken as the droplet surface. In order to evaluate the evaporation of the polydisperse spray, we use a smooth reconstruction which maximizes the Shannon entropy [mead84]. However, the use of fractional moments introduces some theoretical and numerical difficulties, which need to be tackled. First, relying on a study of the moment space, we extend the Maximum Entropy (ME) reconstruction of the size distribution to the case of fractional moments. Then, we propose a new accurate and realizable algorithm to solve the moment evolution due to evaporation, which preserves the structure of the moment space. This algorithm is based on a mathematical analysis of the kinetic evolution due to evaporation, where it shown that the evolution of some negative order fractional moments have to be properly predicted, a peculiarity related to the use of fractional moments. The present model and numerical schemes yield an accurate and stable evaluation of the moment dynamics with minimal number of variables, as well as a minimal computational cost as with the EMSM model [massot2010, kah12], but with the very interesting additional capacity of coupling with diffuse interface model and transport equation of averaged geometrical interface variables, which are essential in oder to describe atomization.

Key words. High order moment method, moment space, realizable high order numerical scheme, polydisperse spray, evaporation, entropy maximization, interface geometry.

AMS subject classifications. 76T10, 35Q35, 65M08, 65M12, 65M99, 65D99.

1 Introduction

In the last decades, with the large change of the Earth's climate, the increase of energy demand and the high consumption of fossil resources, automotive industries are widely concerned to improve engine efficiency and reduce emissions. Actually, in automotive engines and particularly Diesel engine, the fuel is stored as a liquid phase and injected at high pressure in the combustion chamber. The flow of injected fuel in liquid form mixing with air in the combustion chamber is a complex two-phase flow, which has a direct impact on the combustion regime and pollutant emissions. The difficulties of getting accurate measurements inside an engine and the high cost of experiments [Bruneaux2005, Dumouchel_2001, itani2015, Malbec2013, pickett2016, ECN, pickett2015], make numerical simulations a promising and powerful tool. They aim at providing predictive simulations of the flow in automotive engines and understand the various physical mechanisms involved in this complex problem. The liquid phase is initially separated from the gaseous phase, in the neighborhood of the injector nozzle, we have to deal with what is called a separated phases two-phase flow. The interaction with the surrounding gas leads to a complex interface dynamics and eventually to the atomization of the dense liquid core to a polydisperse evaporating spray, in the downstream region, where evaporation and combustion are taking place.

Modeling the fuel injection flow faces major challenges, because of the multi-scale character of this two-phase flow problem. Indeed, it involves a large scale spectrum, varying from the large scales in the separated phases zone to the small scales in the disperse phase. Direct Numerical Simulations (DNS) have been widely used to solve the full flow dynamics in both phases with a full resolution of the interface. The flow in each phase is modeled through the monophasic Navier-Stokes equation and the interface is determined by using tracking methods (Lagrangian methods [Hirten74], Marker And Cell (MAC) methods [Harlow-Weltch-MAC]) or by interface capturing and reconstruction (VOF, Level-Set and hybrid method - see [menard07, lebas2009, LeChenadec201310, emre2014, Arienti14, desjardins2013, shinjo2010, shinjo2011] and references therein). These methods have shown a high accuracy in predicting a sharp interface for incompressible and compressible flows. But, they require a high mesh resolution to compute the solution and they fail to capture the full range of droplet sizes in realistic configuration of high Reynolds and Weber numbers. Therefore, reduced-order models, where averaging approaches are envisioned, have to be considered since a full resolution is out of reach for industrial applications.

The actual existing reduced-order models are not suitable to simulate the whole process of the fuel injection. Instead, two reduced-order model classes have been used depending on the flow region. In the dense core region, diffuse interface models can be used to simulate the flow with lower computational cost. In these models, the interface is considered as a mixing zone, such that the two phases coexist at the same macroscopic position, where each phase occupies a portion of the volume. Several strategies and equilibrium level can be used in order to derive such equations either following an averaging process [DrewPassman99], using fluid mechanics and thermodynamics of irreversible processes [B&N86] or using the principle of least action [drui2016b] (see also [letouze2015, lemartelot2014] and references therein). At the interface, the artificial fluid mixing leads to some level of diffusion of the interface, and as a consequence, we lose important information about the interface geometry. Indeed, the volume fraction is the only variable used to describe the interface geometry, whereas the details of the atomization process cannot be predicted so far with such averaged models. To gain more precision in capturing the atomization, some models add a transport equation of the interface area density [vallet2001, jay06, lebas2009]. However, two variables (volume fraction and interface area density) are not sufficient to describe the polydispersion in the disperse phase region.

In this region of the disperse phase, the dynamics and evaporation of droplets are conditioned by their size. Consequently, size distribution effects should be described in any attempt of modeling a polydisperse spray. The second reduced-order model class relies on a kinetic approach, which has been widely used to describe accurately a population of particle at mesoscopic level and has been shown to be efficient in a number of cases. In this approach, a cloud of droplets is modeled by a number density function (NDF), which satisfies a generalized population balance equation (GPBE), also called Williams-Boltzmann (WBE) equation [williams1958]. The internal variables of the NDF provide a statistical description of some relevant physical properties such as the droplet size, velocity and temperature. The numerical resolution of WBE can be achieved by the Lagrangian Monte-Carlo approach [bird94]. This method is considered to be the most accurate for solving WBE, but leads to a high computational cost for unsteady flows and requires complex load-balancing algorithms in parallel computation. Finally, we encounter some difficulties to couple this Lagrangian method for the dynamics of the spray with an Eulerian method used to solve the continuous phase. Alternative method consists in deriving an Eulerian moment model from the WBE. In this approach, a differential system of a finite set of moments of the NDF is closed through some closure assumptions on the velocity and size distributions, as long as the considered phase space variables are the velocity and the size. The unclosed velocity moment terms are closed by expressing the NDF as a function of the known velocity moments [kah12, boileau2010, boileau2015, vie2015] (monokinetic, Maxwell-Boltzmann, anisotropic Gaussian closures, ). For the modeling of the size distribution, three possible approaches can be used: 1-The first one consists in discretizing the size direction into sections and to use low order size moments in each of these sections. This approach is commonly known as Multi-fluid models [greenberg93, laurent01, laurent2015]. 2- The second one provides a closure of negative as well as fractional moments from the integer ones through a logarithm Lagrangian interpolation and is called the MOMIC method [FrenklachHarris86, Frenklach2002]; it has been essentially use for soot modeling and simulation. However, such an approach suffers from several issues; as Mueller et al [Mueller2009] have pointed that it has a hard time dealing with the boundary of the moment space, that is bimodal distributions as observed in experimental measurements [Zhao2003]. Besides, the interpolation procedure can not preserve the moment space in general and can suffer from important inaccuracy due to the way the negative order moments are approximated. 3- The third approach involves high order moments using either a quadrature of the distribution or a smooth reconstruction using a sum of kernel density functions and quadrature [nguyen2016, laurent-12] or Entropy maximization [kah12, emre2015, massot2010] on the whole size range. The Multi-fluid model suffers from numerical diffusion in the size direction, especially when we are dealing with evaporation, while high order moment approach do not encounter this limitation, instead it requires to solve complex algebra by using high order moments. Furthermore, the computation of the flux of the disappearance droplet through evaporation requires a pointwise evaluation of the size distribution at zero size, which leads to a singular flux in the case of quadrature closure (such as QMOM or DQMOM approaches). On the other hand, the continuous reconstruction through Entropy maximization, used the first time in the EMSM [kah12, massot2010] and CSVM [vieJCP2012] models, shows a great capacity in modeling the polydispersion and evaluating precisely the evaporation fluxes, with a minimal number of variables. Indeed, high order moments make polydisperse modeling possible with only one size-section. Eventually, we can also use an hybrid approach, where we couple high order moment model with the multi-fluid approach. In that case, few sections can be used to gain more accuracy in simulating the droplet dynamics [vieJCP2012].

So far, the existing models do not provide a unified description of the two regions of the flows (separated and disperse phases). In the present contribution, we propose a new high order moment model for the disperse phase with the capacity of describing the interface geometry by analogy with interfacial flows [Drew_Geom]. Our strategy is to resolve the polydispersion by using a set of variables, which can be identified as averages of the gas-liquid interface geometry. For this purpose, we show that some geometrical variables, which can give an accurate description of a complex interface, can be expressed as fractional size-moments of NDF in the disperse zone, when the size variable is given by the surface of the droplets. The present contribution aims at introducing the mathematics fundamentals of the model and showing that we can preserve all the advantages of the previously introduced high order moment methods in [kah12, massot2010] in terms of accuracy, realizability and robustness, but with a much higher potential in terms of coupling with a diffuse interface model. In fact, to evaluate the evaporation flux, we use the Maximum Entropy reconstruction as it was done in the EMSM model, but this time with fractional moments. Entropy maximization with fractional moments was used in another context [Novi2003, Gzyl2010] as a remedy to the ill-conditioning of the procedure when considering a high number of integer moments [Talenti1987]. The considered set of fractional moments is then recovered from the integer ones, and their orders are optimized to minimize the entropy difference with the real function. In the present contribution, a known and small set of fractional moments is used, deduced from physical considerations, in such a way that the problematic is different. The existence and uniqueness of this convex optimization problem under constraints is given in [mead84] in the case of integer moment. While some elements of proof are to be found in the fractional moment case in [Kapur1992], in the present contribution, we propose a generalization of the result in the case of a special set of fractional moments. Moreover, we need to generalize some useful properties of the fractional moment space such as canonical moment as well as lower principal representation [dette97]. These properties are relevant elements to design realizable schemes and algorithms to solve a high order moment system i.e. numerical scheme, which preserves the moment vector inside the moment space and yields accuracy and robustness of the numerical strategy. Finally, we propose a new realizable algorithm to solve the evolution of fractional moments due to the evaporation. The resolution of the evaporation is done by evaluating the disappearance flux, then we compute the internal size evolution by using a specific quadrature, which involves negative order of moments and requires an original strategy compared to the integer moment problem. The proposed strategy is then assessed by a careful investigation of the numerical errors as well as a detailed comparison with the original approach in 0D, 1D and 2D academic configurations. A companion paper [OGST] aims at implementing the proposed model into a massively parallel code with adaptive mesh refinement and showing the potential of the model and numerical method towards realistic engine simulations.

The paper is organized as follows. In a second section, the two-phase flow modeling of polydisperse evaporating sprays in a carrier gaseous flow field as well as the original high order moment modeling are introduced; we also recall the classical averaged geometrical description of interfacial flows in order to identify the relevant geometrical averaged variables. Section 3 is then devoted to the introduction of the new geometrical high order moments for a polydisperse spray as well as the resulting system of partial differential equations on moments derived from the WBE. The closure of the system through maximization of entropy is then presented and its mathematical properties detailed. Once the proposed system of equations is closed, Section 4 and 5 are dedicated to the numerical resolution of the obtained system. While Section 4 is devoted to the transport in phase space, Section 5 focuses on the transport in physical space. Section 6 is eventually concerned with the numerical verification and results in 0D, 1D and 2D, thus assessing the proposed modeling and numerical strategy, before concluding.

2 Two phase flows modeling

2.1 Kinetic modeling of polydisperse spray

The spray consists in a cloud of polydisperse droplets, which can be described statistically with the number density function (NDF) . This function represents the probable number of droplets located at position , travelling with velocity and having temperature and size . In general, the size of a spherical droplet can be given by its volume , its surface or its radius of the droplet. By considering a spherical form, these three geometrical variables are equivalent . In the following, we use the surface as the size variable. The NDF will be simply noted by . This function satisfies the Williams-Boltzmann Equation (WBE) [williams1958]:


where is the acceleration, , the evaporation rate, , the thermal transfer and , the source term, which includes collisions, secondary breakup and coalescence.

The WBE (LABEL:general-WB) gives a complete description of the spray dynamics in a general framework. However, the interaction terms (evaporation, acceleration, thermal transfer) with the gaseous phase and the source terms need to be modeled to close the WBE at the kinetic or mesoscopic level. The second step of modeling consists in reducing the large phase space dimension of this equation by using a moment method. Hence, some closures of the NDF have to be used and the set of assumptions clarified. For the sake of simplicity and the clarity of the presentation, we propose to work with a simplified WBE. We consider a dilute spray at high Knudsen number and small and spherical droplets at low Weber number, obtained after the atomization. Under these assumptions, secondary breakup, coalescence and collision can be neglected. We also assume that thermal transfer can be neglected and will mainly focus on one-way coupling. We refer the readers to the following articles and references therein, showing that such a mesoscopic approach is capable of describing coalescence [doisneau2013], break-up [doisneau13phd], heat transfer [laurent2004] and two-way coupling [emre2014, emre2015]. Furthermore, we consider that the drag force, due to the slip between the droplet and gas velocities, is the only force acting on the droplets. In the following, we model this force by the Stokes law:


where is the gas velocity and is the characteristic response time of the droplet. Finally, we use the law for evaporation:


where is the constant rate of evaporation.

The modeling and numerical strategy will be presented by considering these simplified assumptions. However, in Appendix LABEL:appendix:Realistic_closures, we discuss how to extend the proposed model and the numerical strategy to a general modeling level.

Considering these assumptions and by using non-dimensional variables, the dimensionless WBE reads:


where the superscript refers to dimensionless variables, and are defined as follows: if the concerned variable is , we note by the characteristic value of , then . The maximum size will be taken as the characteristic size, thus , and the characteristic gas time as the characteristic time of the flow. The Stokes number depends linearly on the droplet size :


where is the liquid mass density and is the gas dynamic viscosity. In the following, we consider only dimensionless variables and we omit the superscript .

2.2 High order size-moment modeling: related closure and moment space

The high dimension nature of the phase space of the WBE makes its discretization not convenient for complex industrial applications. Since the exact resolution of the NDF is not required and only macroscopic quantities of the flows are needed for such applications, an Eulerian moment method can be used to reduce the complexity of WBE. The size-velocity moments of the NDF are expressed as follows:


From the WBE, one can derive a system of equation of finite set of these moments. The obtained system is unclosed without any assumption on the NDF form. In the following, we consider that the velocity distribution has a prescribed form:


and does not depend on droplet size. The closure of the velocity distribution requires a specific treatment, especially when we are concerned with modeling particle trajectory crossing (PTC) at high Knudsen numbers. In such case, the collision operator has a very limited effect on the NDF and does not produce any kind of hydrodynamic equilibrium. For accurate modeling of PTC, one can use high order velocity-moment closed through an anisotropic Gaussian velocity-distribution [vie2015, sabat2016]. In the present work, we are not concerned with these modeling issues, and we will consider a monokinetic velocity distribution, which can be also analyzed as an hydrodynamic equilibrium distribution of Maxwell-Boltzmann type at zero temperature [dechaisemartin_th09]:


This closure does not take into account PTC, since only one velocity is defined per position and time. This assumption is valid for low inertia droplets, when the droplet velocities are rapidly relaxed to the local gas velocity. Thus, the droplets do not experience any crossing. It has been shown in [deChaisemartin2009] that for an evolving spray of initial distribution of the form (LABEL:eq:monokinetic-dis) and when the Stokes number is lower than the critical value , the monokinetic assumption is valid all the time. Nevertheless, the model can still be used in the limit of moderately inertial droplets , provided that we use appropriate numerical schemes to handle eventual delta-shocks creation when PTC occurs [dechaisemartin_th09, kah12, sabat2016]. Considering this simplified velocity distribution, we derive the following semi-kinetic system from equation (LABEL:W-B-dimensionless) by considering moments of order and in velocity:


In the present contribution, we adopt the high order size-moment method with a continuous reconstruction of the size distribution to model a polydisperse spray. This approach, based mainly on [massot2010, kah2010, kah12, kah2015], consists in deriving a dynamical system of finite set of size-moments of the NDF. The integer size-moments are defined as follows, with for a non dimensional size interval :


The system of the EMSM model obtained from an integration of the semi-kinetic system (LABEL:eq:semi-kinetic) over multiplied by :


Let us remark that, the term expresses the pointwise disappearance flux of droplets through evaporation. This term only appears in the first equation. But, it participates in the other moments evolution in the same way through the coupling terms , when we consider the integral formulation of this system [massot2010]. To close the system, we need to determine the size distribution from the known moments.

Before proposing a closure of this model, let us recall the definition of the moment space and some useful properties. We denote by the set of all probability density measures of the interval . Then the Nth ”normalized” moment space is define as follows:


Let us notice that , since we use probability density measures in this definition. In our case, we can normalize by to associate the moment vector to , where . The Nth ”normalized” moment space is a convex and bounded space.

Definition 1.

We define the Nth moment space, the set of vectors , whose normalized vector by the number density belongs to the normalized moment space.

Considering this definition and some results from [dette97]: if is in the interior of the moment space, there exists an infinity of size distributions, which represent this moment vector. In other word, there exists an infinity of size distributions , which are the solution of the following finite Hausdroff problem:


Massot et al [massot2010] proposed to use a continuous reconstruction of the size distribution through the maximization of Shannon entropy:


The existence and uniqueness of a size distribution which maximizes the Shannon entropy and is the solution of the finite Hausdroff moment problem (LABEL:moment-sys) was proved in [mead84] for the moment of integer order, when the moment vector belongs to the interior of the moment space and the solution is shown to have the following form:


where coefficients , are determined from system (LABEL:moment-sys). In the same article, the authors propose an algorithm to solve this constrained optimization problem, based on Newton-Raphson method. A discussion of the limitation of this algorithm when the moment vector is close to the boundary of the moment space, or equivalently when the ME reconstructed size-distribution approaches to a sum of delta Dirac functions, is given in [massot2010]. Vié et al [vieJCP2012] proposed some more advanced solutions to cope with this problem, by tabulating the coefficients depending on the moments and by using an adaptive support for the integral calculation, which enables an accurate computation of the integral moments when the NDFs are nearly singular.

This approach shows a high capacity in describing the dynamic of a polydisperse spray using only one size-section. Even though the high order moment formalism provides some key information about polydispersion, it is important to realize that it is restricted to the disperse phase zone. Coupling such an approach with a separated phases model requires some complementary information, which the usual approaches of diffuse interface models can not provide. Indeed, diffuse interface models, used to simulate a separated phases zone, consider the interface as a continuous band layer, where we have lost important information about the interface geometry. The first step, would be to enrich the diffuse interface models as in [drui2016b] in order to transport averaged geometrical variables to gain accuracy about the interface geometry. Nevertheless, even if the usual diffuse interface models would have some more information about the interfacial flow geometry, the coupling of two very different models is usually a rather cumbersome task and relies on some parameters, the described physics will depend on. Consequently we adopt a rather original strategy and we will build a model for the disperse phase, which involves variables that are describing the interface geometry in average, so that we end up with a set of variables that are common to the two zones and can potentially help in building a single unified model able to capture the proper physics in both zones.

In this way, our strategy consists in using averaged geometrical variables in the separated phases zone to model the atomization in future work, and to use the same variables in the disperse phase to describe the atomization by using the same concept as in the EMSM model. However, the integer moments used in the EMSM model do not all represent physical quantities of the flow and more precisely, these moments do not provide information about the interface geometry. Therefore, we need to use other geometrical variables in the disperse phase, while maintaining the attractiveness as well as the efficiency of the EMSM approach. In order to do so, let us introduce the natural geometrical variables in the separated phases zone, before extending this description to the disperse phase.

2.3 Geometrical description of interfacial flows

In many two-phase flow applications, the exact location of each phase is difficult to determine precisely because of the different unpredictable phenomena such as turbulence, interface instabilities and other small scale phenomena that cannot be simulated even with the most powerful supercomputers. Fortunately, in industrial applications, we are more concerned with the averaged features than to the small details of the flow. Therefore, we can use Diffuse Interface Models (DIM) [B&N86, Saurel_jcp09, Saurel_Abgrall_jcp09, lemartelot2014] to describe the interface location in terms of probability and averaged quantities based on averaged operators (ensemble averaged, time averaged or volume averaged). In the following, we define some averaged geometrical variables, which can be used to model the interface in separated phases for a complementary geometrical description. Their definitions are based here on the volume averaged operator following the derivation of Drew [Drew_Geom]. First, we define the phase function for a given phase by:

then, the volume-averaged operator is defined by:


where is a macroscopic space around the position , and is the occupied volume.

Let us emphasize that the DIM can be obtained by applying this operator to the monophasic Navier-Stokes equation. The obtained equations involve the volume fraction variable, which expresses the portion of the occupied volume by a given phase. Moreover, the volume fraction allows to locate the interface up to the averaging scale and is then a first piece of information about the interface geometry:


The second variable treated by Drew [Drew_Geom] and used also in other two-phase flow models [vallet2001, jay06, lebas2009] is the interfacial area density. The importance of this variable relies mainly on the modeling of exchange terms (evaporation, thermal transfer and drag force) as well as modeling the primary atomization. The interfacial area density is interpreted as the ratio of the surface area of an interface contained in a macroscopic volume and this volume.


So far, the interface modeling is still incomplete, since no information on the interface shape is being given. In fact, the small details of the interface can not be modeled accurately using only two geometrical variables. Drew proposed to introduce the two curvatures of the interface in his model. Indeed, these variables give a complementary description of the interface and are highly related to the atomization process, since they are involved in the jump relations at the interface.

The two local curvatures can be defined as follows: let be a point of the interface and is the normal vector at the point . Then, we take a plane that contains and parallel to . As the plane rotates around the normal vector, the intersection curve between the interface and the plane defines a curvature at point which corresponds to the curvature of this 2D curve. As the plane completes a full rotation, it can be shown that it has reached exactly two extremal curvature values: the two principal curvatures and .

Drew derived the dynamical equations of the mean and Gauss curvature respectively and from the differential definitions of and . These two variables read:


These variables are defined only at the interface. Therefore, we need a specific averaging for these interfacial variables. So, we introduce the interfacial averaging operator , defined as follows:


The interfacial averaged Gauss and mean curvature weighted by the interface are defined as follows :


These four geometrical variables can be transported and coupled with a diffuse interface model to gain the accuracy on the interface. Drew [Drew_Geom] derived conservative equations for these variables with source terms, which describe the stretch and the wrinkling of the interface. Its derivation is based on a cinematic evolution of an interface, when the interface velocity is given. But in real application, the interface velocity should be determined from the diffuse interface model. In the following, we would like to express these geometrical variables in the disperse phase as size-moments of the size distribution, and derive a new high order moment model using such variables.

3 Geometrical high order moment model

3.1 Interfacial geometrical variables for the disperse phase

Let us consider a population of spherical droplets represented by their size distribution . Then, by analogy with the separated phases, we express the averaged geometrical variables: volume fraction, interface area density, Gauss curvature and mean curvature in the disperse phase. The definition of these geometrical variables was based on the phase function . This function contains all the microscopic information about the interface. In the disperse phase, we use the statistical information about the droplet distribution, which is given by the size distribution . Considering this function, we define the different geometrical variables in the context of a polydisperse spray as follows:

  1. The volume fraction is the sum of the volume of each droplet divided by the contained volume at a given position:


    The droplet being spherical, .

  2. The interfacial area density is the sum of the surface of each droplet divided by the contained volume at a given position:

  3. The two local curvatures are equal for a spherical droplet . But since we use the mean and Gauss curvatures, we can define two different averaged quantities. Let us notice that, in the case of separated phases, the average mean and Gauss curvatures were defined as an average over a volume and weighed by the interfacial area. In the disperse phase case, this becomes:


These four geometrical variables are expressed as fractional moments of the size distribution :


These moments can be expressed as integer moment by simple variable substitution . However, we prefer to hold the droplet surface as the size variable, since we consider a evaporation law, where the evaporation rate is constant.

3.2 The governing moment equation

In this section, we derive from the kinetic equation (LABEL:W-B-dimensionless) a high order fractional moment model. This model gives the evolution of the mean geometrical interfacial variables due to transport, evaporation and drag force.


where represents the pointwise disappearance flux, and the moment of negative order, , naturally appears in the system after integrating by part the evaporation term in the WBE. In the following, these terms and the associated instantaneous fluxes will be closed by a smooth reconstruction of the size distribution through entropy maximisation.

The use of fractional moment introduces a new mathematical framework of modeling as well as some numerical difficulties, which require a specific treatment. Some useful mathematical properties of the fractional moments space are discussed in Appendix LABEL:appendix:fractional_moment_space. These results will be used to design realizable numerical schemes, i.e. schemes that ensure the preservation of the moment vector in the moment space.

3.3 Maximum Entropy reconstruction

NDF reconstruction through the maximum entropy provides a smooth reconstruction to close the moment system (LABEL:eq:EMSMG) as it was done in the EMSM model. The ME reconstruction consists in maximizing the following Shannon entropy:


under the condition that its first (in our case ) moments are equal to the computed moments


3.3.1 Existence and uniqueness of the solution

In this part, we give a proof of the existence and uniqueness of the ME distribution function. We mention that the case of the integer moments has been already treated in [mead84]. We have used some ideas of this work. But the present proof is completely different and simplified. Indeed, Mead et al [mead84] have used the monotonic properties of the moments, which is a characterisation of the integer moment space, to prove existence of the ME solution. In our case, we will only used the definition of the fractional moment space. The ME problem reads:


where is a moment vector in the interior of the fractional moment space (see Appendix LABEL:appendix:fractional_moment_space).

Lemma 2.

If the constrained optimization (LABEL:eq:ME_optimization_problem) problem admits a solution, then this solution is unique and can be written in the following form:


where .


The Lagrangian function associated to this standard constrained optimization problem is:


where is the vector of the Lagrange’s multipliers.
Let us suppose that, for a given moment vector , there exists a density function which is the solution of the ME problem (LABEL:eq:ME_optimization_problem). So, there exists a vector for which the differential of the Lagrange function at the point vanishes:


where is a positive distribution. Since the system (LABEL:eq:Existence-Lang) is valid for all , it yields:


The problem then consists in finding a vector in which satisfies the moment equations in the system (LABEL:SysME). This problem is equivalent to find an optimum of the potential function :


The Hessian matrix defined by is a positive definite matrix, which ensures uniqueness of an eventual existing solution.     

Lemma 3.

The function , defined in (LABEL:eq:potential_function), is a continuous function in , and goes to infinity when .


Let us suppose that the last assertion is wrong, so there exists a sequence such that when and .
Hence, there exists such that:


We write for each , , such that and .
Since the sequence is a bounded sequence, we can subtract a convergent subsequence , where is an increasing function and:


To simplify the notation, we can directly consider that when .
We note by and .

Since the vector is a moment vector, there exists a non-negative distribution function such that for , and


Since the first integral is positive


When tends to infinity, we get:


We have , and , and since is a continuous function, it follows from the inequality (LABEL:ineq-proof1) that there exists in which and . Since converges uniformly to in , then, for all and for the large enough values of :


Using these results in the inequality (LABEL:inequality), as well as , we get:


In the limit when goes to infinity, we get the contradiction , thus concluding the proof.     

Theorem 4.

If the vector belongs to the interior of the Nth fractional moment space, then the constrained optimization (LABEL:eq:ME_optimization_problem) problem admits a unique solution, which is in the following form:



The proof is straightforward by using the two last lemmas.     

3.4 Algorithm of the NDF reconstruction through the Entropy Maximisation

The reconstruction of the NDF through the maximization of Shannon entropy goes back to finding the Lagrange’s multipliers such that:


where . Solving this nonlinear system is equivalent to optimizing the convex function . We solve the problem by using the Newton iteration as proposed in [mead84]:

  Choose initial guess of the vector .
  while  do
     for  do
     end for
     for  do
     end for
  end while
Algorithm 1 EM algorithm

The integral computations are done by using Gauss-Legendre quadrature. In [mead84], it is shown that -point quadrature is very efficient to calculate accurately the different integral expressions involved in the algorithm.

4 Numerical resolution of the moment governing equations

The phenomena involved in the model can be classified in two main classes: the transport of the droplets in the physical space and the source terms, which induce the particle evolution in the phase space (velocity and size) through evaporation and drag. We use operator splitting techniques [doisneau14, descombes14] to separate the resolution of the different phenomena, we can then solve each operator separately.

4.1 Resolution of evaporation

Let us consider a pure evaporation without transport nor drag. The kinetic equation in this case reads:


The integration of this equation multiplied by the vector yields to the system of equation of the moment in the section


where the is obtained by ME algorithm.
Solving this system using classical integrator such as Euler or Runge-Kutta methods leads to serious stability problems. In fact, classical ODE integrators do not ensure the preservation of the moments in the moment space [massot2010]. This property is essential for robustness and accuracy to reconstruct a positive NDF.

4.1.1 Exact kinetic solution through the method of characteristics

The exact solution of the NDF evolution using a evaporation law, can be obtained easily by solving analytically the kinetic equation (LABEL:evap-kinet) using the method of characteristics:


For more general evaporation law, when the evaporation rate is a smooth function of the size, the kinetic equation is written as follows:


By multiplying the equation by , we obtain:


where . For a given initial time and size , we define the one variable function , such that is the characteristic curve verifying:


The derivative of vanishes, thus, we obtain the following expression:


Finally, we obtain the exact solution of the size distribution as follows:


4.1.2 Fully kinetic scheme

At each time step , the reconstructed NDF is determined using ME algorithm, then the exact kinetic solution of the NDF can be expressed analytically as a function of this initial NDF thanks to (LABEL:eq:exact_kinet_solution_d2). Finally, the updated moments are computed as follows:


The calculation of this integral can be achieved by using Gauss-Legendre quadrature as previously. The method is simple and accurate. However, its extension to complex evaporation law, where the evaporation rate depends on several parameters, could lead to heavy calculations. To clarify this point, we consider a smooth evaporation law, where depends on the droplet size. Considering the exact kinetic solution (LABEL:eq:exact-kinetic-sol2), we can obtain the following expression:


To compute this integral, we need to determine the size , for each abscissa of the -quadrature points, by solving the ODE system (LABEL:eq:Lagrang-size-evol).

4.1.3 Inefficiency of the original EMSM algorithm for evaporation

Massot et al [massot2010] proposed a realizable algorithm to solve correctly the evaporation moment system in the case of integer moment. The idea of this algorithm is also based on the known solution of the kinetic equation (LABEL:evap-kinet). In order to generalize the evaporation law for more complex laws, the calculation of the moments was done in three steps instead of computing directly the integral formulas (LABEL:eq:Kinetic-moment-sol) in the case of law or (LABEL:eq:Exact_kinet_mom_RS) in general case using a large number of Gauss-Legendre quadrature.
This algorithm reads in the case of fractional moments as follows:

  • From the moment vector , we reconstruct NDF by using the ME algorithm. Then, we compute the disappearance flux, which represents the droplets which will be totally evaporated at the next step:

  • The abscissas and the weights corresponding to the lower principal representation (LABEL:eq:low_prin_frac) of the moments111We can also use hybrid approach with multi-size sections, in this case we need to take into account the evaporation fluxes coming from the right section: , where indexes the section. are computed by using the PD algorithm [gordon1968].


    where is the moment vector in the support .

  • We calculate the moments at the next step:


In Figure LABEL:fig:emsm-algo, we compare the evolution of the reconstructed NDF from the moments computed with the fully kinetic algorithm, the reconstructed NDF from the moments computed with the above algorithm and the exact kinetic solution. The results show the inefficiency of this last algorithm to predict the right kinetic evolution of the NDF, when we use fractional moments. In the next section, we give a mathematical development of the exact kinetic solution of the fractional moment (LABEL:eq:Kinetic-moment-sol) in order to understand the limitation of this algorithm, then we propose a new original solution to adapt the algorithm.

Figure 1: Initial size distribution (dashed line) and the reconstructed size distribution at using: exact kinetic solution (solid line), EMSM algorithm (cross), fully kinetic scheme (circle)

4.1.4 Adapted evaporation scheme for fractional moments

In this section, we rewrite the updated moment as a function of an initial infinite set of moments on the support . This result will allow us to identify the missing point in the last algorithm. Then, a correction of the algorithm is proposed at the end of the section as well as its extension to more general evaporation laws.

Lemma 5.

For all positive integer and for all the function admits a power series which converges normally to the function .


The case of is an even integer is trivial. Let us consider the case where is an odd number. The function admits a power series and the convergence radius of the series is , and we for all :


we can show that the coefficients can be written as follows:


Now, by using the Stirling’s approximation:


We can prove the following equivalence relation when tends to infinity.


Therefore the series is convergent for any positive integer .     

We deduce that converges normally to for
. Thus, We can invert the sum and the integral in the moment expression:


Equation (LABEL:eq:series_algo) shows that the fractional moment at depends on an infinite set of the moments of support ( where ). In the case of the EMSM model, where only integer moments are used in the model, the same expansion of the exact kinetic solution of the integer moment involves only the four transported moments. For this reason, the evolution of the moment can be evaluated exactly by translating the abscissas (LABEL:eq:al_size_evol) of the lower principal representation (LABEL:eq:al_emsm_quad), which is not the case for the present model. Therefore, we need an adapted quadrature instead of the one used in (LABEL:eq:al_emsm_quad), such that the moment evolution in the equation (LABEL:eq:al_size_evol) approximates accurately the exact kinetic evolution (LABEL:eq:series_algo).

Definition 6.

We consider a density function of the support . The abscissas in and the weights are called Gauss quadrature rule of order corresponding to if and only if for all polynomial function of degree , we have:


The existence of such quadrature for each positive measure is proved in [Gautschi]. The Gauss quadrature of order is also the lower principal representation (see Appendix LABEL:appendix:fractional_moment_space) associated to the first integer moment of measure .

Our objective is to find an adequate Gauss quadrature with the lowest possible quadrature number , such that the following approximation


is accurate. More precisely, we would like to find an adequate Gauss quarature such that the difference


is at least to ensure the convergence of the numerical scheme. Unfortunately, it is a difficult task to prove this point and this problem goes back to seeking for an accurate estimation of the Gauss quadrature errors [Gautschi]. However, the task of providing a rigorous proof of the convergence is beyond of the objective of the present paper. In the following, we focus on practical issues and we will present a general strategy to decrease the error by cancelling a finite set of first terms in the infinite sum (LABEL:eq:algo-error). We mention that for even the error .

Let us consider the measure , where is a positive integer variable, defined in the support for the variable222This substitution makes the link between the fractional moment and its corresponding integer moment – see Appendix LABEL:appendix:fractional_moment_space. . Then, the Gauss quadrature of order of the measure , represented by the abscissas and weights , satisfies:


for (). The computation of the weights and the abscissas can be determined by the PD algorithm [gordon1968]. The coefficients and are defined as follows: