Adaptive Mesh Refinement and High Order Geometrical Moment Method for the Simulation of Polydisperse Evaporating Sprays
Raffinemnt de maillage adaptatif et méthods de moments géométriques d’odre élevé pour la simulation des sprays en évaporation
^{1}
Laboratoire EM2C, CNRS, CentraleSupélec, Université ParisSaclay, Grande Voie des Vignes, 92295
ChâtenayMalabry Cedex  France
^{2}
IFP Energies nouvelles, 14 avenue de BoisPréau, 92852
RueilMalmaison Cedex  France
^{3}
Fédération de Mathématiques de l’École Centrale Paris, CNRS FR 3487, Grande Voie des Vignes, 92295
ChâtenayMalabry Cedex  France
^{4}
Institut Carnot IFPEN Transports Energie, 14 avenue de BoisPréau, 92852
RueilMalmaison Cedex  France
email: stephane.dechaisemartin@ifpen.fr
^{*} Coressponding author
Received:
19
February
2016
Accepted:
7
June
2016
Predictive simulation of liquid fuel injection in automotive engines has become a major challenge for science and applications. The key issue in order to properly predict various combustion regimes and pollutant formation is to accurately describe the interaction between the carrier gaseous phase and the polydisperse evaporating spray produced through atomization. For this purpose, we rely on the EMSM (Eulerian MultiSize Moment) Eulerian polydisperse model. It is based on a high order moment method in size, with a maximization of entropy technique in order to provide a smooth reconstruction of the distribution, derived from a WilliamsBoltzmann mesoscopic model under the monokinetic assumption [O. Emre (2014) PhD Thesis, École Centrale Paris; O. Emre, R.O. Fox, M. Massot, S. Chaisemartin, S. Jay, F. Laurent (2014) Flow, Turbulence and Combustion 93, 689722; O. Emre, D. Kah, S. Jay, Q.H. Tran, A. Velghe, S. de Chaisemartin, F. Laurent, M. Massot (2015) Atomization Sprays 25, 189254; D. Kah, F. Laurent, M. Massot, S. Jay (2012) J. Comput. Phys. 231, 394422; D. Kah, O. Emre, Q.H. Tran, S. de Chaisemartin, S. Jay, F. Laurent, M. Massot (2015) Int. J. Multiphase Flows 71, 3865; A. Vié, F. Laurent, M. Massot (2013) J. Comp. Phys. 237, 277310]. The present contribution relies on a major extension of this model [M. Essadki, S. de Chaisemartin, F. Laurent, A. Larat, M. Massot (2016) Submitted to SIAM J. Appl. Math.], with the aim of building a unified approach and coupling with a separated phases model describing the dynamics and atomization of the interface near the injector. The novelty is to be found in terms of modeling, numerical schemes and implementation. A new high order moment approach is introduced using fractional moments in surface, which can be related to geometrical quantities of the gasliquid interface. We also provide a novel algorithm for an accurate resolution of the evaporation. Adaptive mesh refinement properly scaling on massively parallel architectures yields a precise integration of transport in physical space limiting both numerical dissipation as well as the memory trace of the solver. A series of testcases is presented and analyzed, thus assessing the proposed approach and its parallel computational efficiency while evaluating its potential for complex applications.
Résumé
La simulation prédictive de l’injection diphasique dans les chambres de combustion automobiles représente un enjeu majeur scientifique et applicatif. La description détaillée de l’interaction entre le brouillard de gouttes polydispersé produit par atomisation et l’écoulement gazeux est fondamentale pour prédire les régimes de combustion et la formation de polluant. Pour décrire la phase liquide, le modèle Eulérien polydisperse EMSM (Eulerian MultiSize Moment) est choisi. Cette approche de type moments d’ordre élevé en taille avec reconstruction continue par maximisation d’entropie est construite à partir d’un modèle mésoscopique de WilliamsBoltzmann sous l’hypothèse monocinétique [O. Emre (2014) PhD Thesis, École Centrale Paris; O. Emre, R.O. Fox, M. Massot, S. Chaisemartin, S. Jay, F. Laurent (2014) Flow, Turbulence and Combustion 93, 689722; O. Emre, D. Kah, S. Jay, Q.H. Tran, A. Velghe, S. de Chaisemartin, F. Laurent, M. Massot (2015) Atomization Sprays 25, 189254; D. Kah, F. Laurent, M. Massot, S. Jay (2012) J. Comput. Phys. 231, 394422; D. Kah, O. Emre, Q.H. Tran, S. de Chaisemartin, S. Jay, F. Laurent, M. Massot (2015) Int. J. Multiphase Flows 71, 3865; A. Vié, F. Laurent, M. Massot (2013) J. Comp. Phys. 237, 277310]. La présente contribution propose une extension majeure de ce modèle [M. Essadki, S. de Chaisemartin, F. Laurent, A. Larat, M. Massot (2016) Submitted to SIAM J. Appl. Math.] dans l’optique de le coupler à un modèle de type phases séparées avec résolution d’interface. La nouveauté se situe en termes de modélisation, de schémas numériques et d’implémentation. La nouvelle approche de moments d’ordre élevé en taille repose sur des moments fractionnaires de la surface que l’on peut relier à des quantités géométriques de l’interface gazliquide. Un nouvel algorithme permet une intégration précise de l’évaporation. L’utilisation d’un maillage adaptatif pour le transport dans l’espace physique passant à l’échelle sur architectures parallèles offre un contrôle de la dissipation numérique tout en limitant la trace mémoire de l’application. Une série de castests permettant de vérifier la qualité de l’approche, son efficacité parallèle et son potentiel pour des applications complexes est présentée et analysée.
© M. Essadki et al., published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Introduction
In the last decades, automotive industries have been widely concerned by the efficiency of combustion devices and the reduction of pollutant emissions. Considering the difficulties of getting accurate measurements inside an engine, implying highcost experiments, numerical simulation is nowadays considered as a good alternative in order to improve our understanding of twophase flows and twophase combustion and to design new combustion chambers. Actually, in automotive engines, the fuel is stored as a liquid phase and injected at high pressure in the combustion chamber. Close after the outlet nozzle, the fuel is still in liquid form and separated from the gaseous phase by an interface; we have to deal in this region with what is called a separated phases twophase flow. The atomization process is then going to produce a complex dynamics of the interface and eventually yields a polydisperse evaporating spray, which will control the combustion regime as well as pollutant formation. It has been shown that the polydisperse character of the spray has a key influence and should be described in any attempt of modeling such flows. The challenge is to provide a model as well as a numerical strategy, which is able to capture the large scale spectrum of such physics and provide predictive numerical simulations. In fact, the challenge is as much a scientific as an application one.
The first possibility is to consider a DNS (Direct Numerical Simulation) approach and to envision solving the full dynamics of both phases and of the interface from inside the injector until the polydisperse spray has been formed. The interface between the two fluids is considered as a sharp discontinuity. In this context, its exact location is needed to determine precisely an existing region for each phase and then apply the monophasic NavierStokes equations separately in each dedicated region. The dynamics of the interface can be traced by some tracking methods (Lagrangian methods [1], Marker And Cell (MAC) methods [2]) or by interface capturing and reconstruction (VOF (Volume Of Fluid), levelset and hybrid method  see [3–6] and references therein). Even if such approaches are powerful and have been used on parallel architecture, they are mainly focused on mechanically incompressible flows and can not resolve the full spectrum of scales present at representative Reynolds and Weber numbers of realistic fuel injection. More specifically, the resolution of droplet size distribution is far from being reached in such approaches, even though they are essential in order to analyze and model atomization processes.
Consequently, reducedorder models have to be considered since a full resolution is out of reach. So far two reducedorder approaches have been used depending on the distance from the injector mouth. In the dense core region, one possibility is to use diffuse interface models. The interface is considered as a mixing zone, such that the two phases coexist at the same macroscopic position and each phase occupies a portion of the volume. In this type of model, a continuous color function is introduced: the volume fraction. This variable takes values close to 1 in a region of a pure given phase and values close to 0 in the regions where this phase nearly disappears. Crossing the interface, the volume fraction varies continuously from 0 to 1 or the other way around. At the interface, the artificial fluid mixing may lead to some thermodynamic difficulties [7] when too much numerical diffusion is involved. Therefore, high resolution is recommended near the interface to limit this diffusion. Several strategies and equilibrium level can be used in order to derived such equations either following an averaging process [8], using fluid mechanics and thermodynamics of irreversible processes [9] or using the principle of least action [10] (see also [11, 12] and references therein). The modeling choices have a crucial impact in terms of numerical schemes [10, 13–15]. Furthermore, some mean geometrical quantities such as the interface area density can be transported to gain a more accurate description of the interface [5, 16, 17] and decrease the lack of information on the interface shape due to the initial averaging process. This way, an accurate modeling of the dynamics of the interface and a good prediction of the interaction between the two phases can be reached at large scale [13], whereas the details of the atomization process cannot be predicted so far with such averaged models.
Another reducedorder model is related to the description of the polydisperse spray far from the injection point, once the atomization process has been conducted, and relies on a statistical approach for the droplet population by using a Number Density Function (NDF) which satisfies a Population Balance Equation (PBE), also called WilliamBoltzmann (WB) equation [18]. The resulting coupled system of equations involving a model for the gaseous carrier phase has to be resolved using either a Lagrangian or an Eulerian approach as far as the liquid phase is concerned. Solving this equation can be achieved by using Lagrangian Direct Simulation MonteCarlo method (DSMC [19]). This method shows a good efficiency in numerous cases. Nevertheless, the method suffers from several drawbacks. Indeed, this Lagrangian method requires a large number of particle samples to reach the convergence which can lead to high Central Processing Unit (CPU) time and memory consumption, especially in polydisperse unsteady cases. Furthermore, we encounter some difficulties to couple this Lagrangian method for the dynamics of the spray with the Eulerian method usually used to solve the gas flow [20]. Finally, complex loadbalancing algorithms are needed to reach high performance computing. Therefore, we prefer to use an Eulerian model for the description of the spray, at least as an alternative strategy. More precisely, such an Eulerian model is based on a moment method which, instead of solving for the whole NDF in a large parameter space, solves for a finite set of sizevelocity moments of this NDF. Considering some assumptions on the velocity distribution (monokinetic, MaxwellBoltzmann distribution, etc.), we are able to close the equation in the velocity direction, and derive a semikinetic system for the size distribution function [21]. Then, three possible approaches can be used to solve this size distribution. 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 multifluid models [22, 23]. The second approach involves higher order moments using either a quadrature of the distribution or a smooth reconstruction using kernels and quadrature or entropy maximization (see [24–26] and references therein) on the whole size range. In the third approach, we use a hybrid method [27–31] using high order moment method on discretized size intervals. This last approach has been used in recent contribution of some of the authors in [4, 24, 25, 27, 31, 32] and seems to be a very good compromise in terms of accuracy versus computational cost for automotive engine simulations.
The main drawback of the two types of reducedorder models is due to their very different nature and the difficulty to couple both models in order to reach a predictive simulation of the whole process. Recent contributions in [12, 33, 34] have proposed a way to couple these two reducedorder models, but several issues are still wide open in order to increase the level of prediction one can reach with such a model coupling. We insist on the fact that the key issue is the ability of predicting the polydisperse character of the spray and the purpose of the present work is to use a recently introduced new hybrid high order size moment model [35] with the same level of accuracy as what has been done in [4, 24, 25, 27, 31, 32], but with a much higher potential in terms of coupling with a diffuse interface model. The key ingredient in this new model is the use of fractional moments, which describe the interfacial topology of the droplet population in analogy with interface geometry used for separated phases modeling [36]. Since geometrical variables are now transported in the spray region, it can be coupled with a separated phases flow model, which would also incorporate information on the average interface geometry. In this lies the first main novelty of the current paper. Besides, a new realizable numerical scheme to solve for the evaporation process in the case of fractional moments is developed and presented. Eventually, even though this article only presents simulations of disperse flows, we are still concerned with the coupling of our new model with a separated phases simulation. Yet, such separated phases flows need an accurate description of the dynamics of the interface between the two phases, since small scales are known to be of great influence on the whole multiscale phenomena. Moreover, we would limit the numerical diffusion at the interface to reduce artificial mixing and inherent numerical difficulties. This justifies our will to integrate mesh adaptation in the numerical simulation of this model. This has been made possible by the use of the Adaptive Mesh Refinement (AMR) library p4est [37]. Nonetheless, the use of AMR techniques in the spray simulation is also very helpful to capture the spray segregation [38] by using a finer mesh in the high concentrated regions and keeping a coarser mesh in the low concentrated or vacuum regions. This way we obtain an accurate result with a relatively low CPU time and memory usage. A series of testcases is presented and analyzed, thus assessing the proposed approach and its parallel computational efficiency while evaluating its potential for complex applications.
This paper is organized as follows: in Section 1 we present the framework for the Eulerian spray modeling, including the moment method used to derive a reducedorder modeling from the mesoscopic level and an introduction to our new geometrical moment model. The new numerical schemes to solve evaporation for the new model are presented in Section 2. Section 3 contains a brief presentation of the different AMR techniques and an overview of the main functionalities of the p4est library used to manage the adaptive mesh refinement. We also present our code, named CanoP, which assists a “simple” use of the AMR library p4est with finite volume schemes for different data models. Furthermore, we present a first order numerical scheme for the discretization of the whole procedure in space on nonconforming meshes, thus suitable for AMR. An extension of this spatial discretization to higher order, still ensuring a realizability condition, is also given. Finally, Section 4 is concerned with numerical results and related discussions.
1 Eulerian Spray Modeling
We adopt a statistical approach to describe the droplet population. This description relies on the evolution of the NDF , which represents the probable number of droplets located at position , traveling with velocity , having size S and temperature T. In the following, we neglect the thermal inertia of the droplets, therefore the NDF is reduced to . By considering a very dilute spray of small droplets, we can also neglect collisions and secondary breakup of droplets. Then, the dynamics of the droplets and their interaction with the surrounding gas is given by the WilliamBoltzmann (WB) equation, [18].(1)
where is the evaporation rate and is the drag acceleration. The evaporation is modeled by a simple d ^{2} law [39], so that R_{ S }(S) = −K < 0 is constant. On the other hand, the drag force term is given by the Stokes law. It is valid for low Reynolds number (calculated with the droplet size S). The force then reads:(2)
where is the gas velocity (not disturbed by the presence of the droplet) and τ_{ p }(S) is a characteristic response time of the droplet to the Stokes drag force:(3)
where ρ _{ l } is the liquid mass density and μ _{ g } is the dynamic gas viscosity.
In the following, we use dimensionless variables, therefore the WB Equation (1) can be rewritten as follows:(4)
where is the Stokes number, which characterizes the response of the droplet to the gas dynamic and τ _{ g } is the time scale of the gaseous flow. Besides, K is the ratio of τ _{ g } to an evaporation time. Let us underline that we focus here on simplified droplet models for the sake of simplicity of the presentation. The present approach can be extended to complex models as presented in [4].
1.1 Eulerian Moment Method
Given the high dimension nature of the parameter space of the WB equation, its discretization for industrial applications is not attainable. Since a high accuracy on the resolution of the distribution of the droplets is not necessary for the applications we are concerned with, and only macroscopic quantities are needed, an Eulerian moment method can be a promising alternative to Lagrangian approaches for such problems. In this approach, we consider the moments of the NDF over the phase space (velocity and size):(5)
In this section, we present the general derivation of these methods, based on the work of Emre [4], Emre et al. [24, 32], Kah [21], Kah et al. [25, 27]. To simplify the modeling, we considered that the velocity distribution profile ϕ is independent of the size variable, so that the NDF can be written:(6)where is the mean velocity conditioned by size. In this paper, we even simplify our study to the case of sizeindependent macroscopic velocity . However, it is worth noticing that such a complete sizevelocity correlation has already been considered and tackled in the CSVM (Coupled SizeVelocity Moment) model [31].
In the EMSM model developed by Kah [21], we consider a monokinetic velocity distribution, which can be also analyzed as an hydrodynamic equilibrium distribution of MaxwellBoltzmann type at zero temperature. It implies that the droplets at a given position and time have the same velocity. Therefore, the distribution ϕ is given by a Dirac deltadistribution:(7)
Then we have a uniquely defined velocity per position and can not take into account Particles Trajectories Crossing (PTC) under this assumption. Therefore, this model is valid only for low inertia droplets, the critical Stokes number for the appearance of PTC having been defined in [40]. Otherwise, ballistic trajectories of higher inertia droplets will generate crossings between particles, what would be captured as deltashocks by what we call the MonokineticEquilibrium model. Extension to PTC through the anisotropic Gaussian model can be considered and we refer to [38, 41] for more details.
Considering this simplified velocity distribution, we derive the following semikinetic system for Equation (4):(8)
Two possible approaches can be used to solve the size distribution. On the one hand, the size phase can be discretized into small intervals, wherein the NDF n(S) is considered to have a simple form. In the multifluid model [23, 42], the number density is considered to be constant in each section and the computation then requires the calculation of one sizemoment per section. On the other hand, one can use high order moment methods either using quadrature methods or preserving the continuity of the size distribution, by computing the evolution of high order sizemoments and reconstruct a continuous size distribution function n(S) from this set of high order sizemoments in order to close the model. We use a hybrid method where a size discretization can be considered, but using high order moments in the potential sections and relying on a smooth reconstruction. Two moments per section with affine reconstruction lead to a twosize moment method developed in [28, 29], whereas we rather choose to use four moments in each section and maximization of entropy in order to reconstruct the distribution function [27, 30, 31]. In this last approach, we can use one or very few sections (compared to the multifluid model, which requires a rather fine discretization to properly predict the evolution of the global moments). Both EMSM and CSVM models make use of this last approach.
For a given density repartition in size n(S), the moment of order k is defined by:(9)
Using nondimensional variables, S _{ max } can be set to S _{ max } = 1.
Starting from the semikinetic system (8), we derive the system on N sizemoments:(10)where (resp. ) is the instantaneous disappearance flux (resp. appearance flux) coming from other sections. Since we consider only one section in the EMSM model and S _{ min } = 0. Also, the Stokes number has been supposed to depend linearly on the droplet surface: .
The system of moments consists of equations, where d is the space dimension. The first equations represent the evolution of the size moments through free transport and evaporation. The last vectorial equation is the momentum equation, used to determine the velocity. Under the MonokineticEquilibrium assumption, we only need one velocity moment to compute the velocity. In this case, we have considered:(11)
Yet, the system is not fully closed because the boundary value is unknown, since we solve for the moments. It can be interpreted as the instantaneous droplets disappearance flux by evaporation. Even though this term appears only in the first size moment equation, it influences recursively the other equations through the last terms .
In [30], the integration of the system in the case of pure evaporation (no transport) is solved in two steps for a given time step: a) compute the disappearance flux and subtract this amount from the set of moments, b) shift the size distribution of the remaining droplets using the link between moments and abscissas and weights of the lower principal representation. The link with a quadrature method of moments approach [30] is explained in the paper. To close the system, we need a smooth reconstruction of the NDF and this is done by the maximization of the following Shannon entropy:(12)
The existence and uniqueness of a density function n which maximizes the Shannon entropy and satisfies:(13)was proved in [43], and the solution is shown to have the following form:(14)where coefficients , are determined from system (13). In the same article, the authors propose an algorithm to solve this constrained optimization problem, based on NewtonRaphson 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 distribution is close to a sum of delta Dirac functions, is given in [30]. Vié et al. [31] proposed some solutions to cope with this problem, by tabulating the lambdas 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 NDF are singulars.
1.2 Geometrical Moment Model
Twophase flows involve a variety of phenomena, which depend on the flow topology and the way the two phases interact across the interface. Then, the dynamic of the interface is a crucial feature of an accurate modeling. Indeed, in the fuel injection process, we encounter different classes of twophase flow topologies, going from separated twophase flow to a dilute disperse phase. In each distinct region, different mechanisms occur depending on the topology of the interface. In the separated phases region, the interface deformation leads to the atomization process. The accuracy of the resolution of the interface will then highly influence the size distribution of the spray in the forthcoming regions. In the disperse phase, the liquid has a spherical and more stable shape. Thus, we do not need to capture precisely the interface of each droplet better than just considering the distribution of an ensemble of spherical droplets. In this particular region, the leading phenomena are evaporation and drag and they depend mainly on the droplet size. This is why taking into account the polydisperse character of the spray is a key feature of Eulerian models for the disperse phase.
In the following, we briefly present some important geometrical variables used to model the interface in the case of separated phases. Then, we express the same variables in the context of a disperse phase. Finally, from the kinetic model (4), we derive a new moment model using these geometrical variables [35].
We consider the phase function:to define the volume averaged variables that we use in diffuse interface models. We also use the gradient of this function, which can be seen as a Dirac deltafunction at the interface, to define some interfacial averaged variables.
In [36], Drew proposed a description of the interface based on the following variables:

the volume fraction variables express the volume portion occupied by a given phase:
where is a macroscopic space around the position and is the total volume.

the interfacial area density is defined as the mean quantity of interface per unit volume:

the average Gauss and mean curvatures are defined as an interfacial averaging of the Gauss and mean curvatures. The local Gauss G (resp. mean H) curvature is the product (resp. mean) of the two interfacial principal curvatures and :
The two local curvatures are defined as follows: we consider a point P at the interface and its normal vector . Then, we take a plane that contains P and its normal and let it rotate around . In each position the intersection curve with the interface defines a curvature at point P. As the plane completes a full rotation, it can be shown it has reached exactly two extremal curvature values: the two principal curvatures.
Then the average Gauss and mean curvature are:(18)
In the case of a disperse phase, we express the same variables for a droplet population described by the NDF .
In fact, these four geometrical variables have been expressed in [35] as fractional moments:(19)where the halfinteger moments are defined by:(20)
These moments can indeed be expressed as integer moments by simple variable substitution . However, we prefer to hold with the droplet surface S as the size variable, since the evaporation rate K is constant for the law.
Following the same derivation as the EMSM model, we derive the following system for the new moments:(21)where K is the constant evaporation rate for the evaporation law and is the instantaneous disappearance flux. Since , then when .
Now, the negative order moment appears in the system after integrating by part the evaporation term in the WB equation:(22)
This term and the associated instantaneous flux need to be closed and we choose to do that with an adequate reconstructed NDF. We adopt a reconstruction by entropy maximization as it has already been done in the previous section. The existence and uniqueness of such a reconstruction are proved in [35]. The reconstructed NDF has the following form:(23)
To determine the coefficients, we use the same inversion algorithm as in [43].
2 Numerical Method
In this section, we present the discretization method used to solve numerically the moment system (21). Three different phenomena are involved in this system: the spatial evolution of the moments by the spatial transport, the size variation by the evaporation and the velocity evolution by the drag force. It is then interesting to separate these phenomena by using an operatorsplitting technique [44]. This method simplifies the complexity of the system by solving alternatively and independently the different dynamics. Moreover, the use of Cartesian grid allows us to use a dimensional splitting. Therefore, the onedimensional numerical scheme can be easily extended to multidimensional space.
2.1 Transport Scheme
In the dimensional splitting framework, we consider a free transport in one direction (we chose the direction) of the droplets without the evaporation nor the drag force. We choose to present the scheme in a two dimensional space to lighten the notations and we note .(24)
The above differential system is not strictly hyperbolic, and the semikinetic system from which the moment system is derived, is similar to the pressureless gas system:(25)
Following the idea of Bouchut et al. [45], de Chaisemartin et al. [42], Kah et al. [27] used this approach to develop a finite volume kinetic scheme for the EMSM model. We use the same scheme to solve numerically the system (24). This scheme is established by considering the equivalence between the pressureless system and the following kinetic system:(26)where is the microscopic velocity of the droplets. The proof of the equivalence between the two last systems can be found in [45].
In this section, we consider a uniform discretization of the space. In Section 3.2.1, we present the discretization of the moment system on non conforming meshes. We note x _{ i } = (i + 1/2)Δx the center of cell i, and (resp. ) the left (resp. right) extremity of cell . Finally, we define the averaged moment and the momentum attributed to cell by:(27) (28)
Using the exact solution of the kinetic equation (26) for :(29)we establish the following scheme using the same steps as in [27]:(30)where the fluxes can be decomposed into and such that:(31)and(32)
For a first order scheme, we consider a constant piecewise reconstruction for the moments and the velocities. Then the fluxes are:(33)and(34)
For higher order scheme, we refer the reader to the original article [27].
2.2 Evaporation Solver
Now, we consider the evaporation of a polydisperse spray without transport. The kinetic evolution of the spray is obtained by considering a homogeneous number density (35)
By integrating this equation, we get the following Ordinary Differential Equations (ODE) for the moments evolution:(36)
We use Entropy Maximization (EM), see Massot et al. [30], to reconstruct the most probable underlying number density fonction , in order to close the system and estimate the instantaneous disappearance flux of the droplets. The NDF is computed at each time step. Simple ODE integrators like Euler or RungeKutta method cannot be used to solve this system because they do not ensure the preservation of the moments within the moment space. By moment space, one has to understand the set of moments which are the sizemoments of a positive distribution . When this condition is not satisfied, the NDF reconstruction algorithm does not converge. In order to ensure the numerical preservation of the moment space, a kinetic scheme was developed in [30]. In that article, the evolution over of the integer moment system (EMSM) due to evaporation was done in two main steps:

droplet with size smaller than will completely evaporate. The disappearance flux is then evaluated and subtracted from the moments;

other droplets will see their size reduced by . It is shown that this operation can be done numerically by simply translating the two abscissas of the lower principal representation of the resulting moments.
For our new model, the equivalence between the evolution of the evaporation through the kinetic equation (35) and the moments system (36) is still valid, but it requires a much more involved algebra [35]. In this paper, we admit it and rely on the following algorithm.
Algorithm

Reconstruct which corresponds to the moment vector by entropy maximization, then calculate the disappearance flux:

Calculate the negative order moments within the sizeinterval
for , where is the number of negative moments used to correct the algorithm presented in [30] with fractional moment. The origin of these negative moments is discussed in [35]. The singularity of the negative moment integral when becomes very small limits the use of high value of . In practice, and the choice or renders accurate enough solutions. The other moments of positive order are computed using the disappearance flux:(39)where .

The abscissas and the weights corresponding to the lower principal representation of the moments for are computed using the PDA [46], such that:
where n _{ s } = 2 + M/2.

Calculate the new set of evolved moments:
This algorithm is an original extension of the one used in [30] since it has been adapted to fractional moments and copes with the inherent difficulties, while preserving a high accuracy. The different integral terms involved in the algorithm are computed by using Gauss quadrature approximation.
2.3 Evaporation and Drag Force
In this paragraph, we present a coupled solver for the spray evolution under evaporation and drag force. The corresponding system of equations is:(42)
Since the first four equations do not depend on the last equation, the moments are computed using the algorithm presented in the last section. We use the method developed in [31] to solve the evolution of the velocity by the drag force coupled with evaporation.
The momentum evolution is done in two steps: first, we remove the droplets which will disappear through evaporation during the next time step. Then the size distribution of the remaining droplets is approximated by the lower principal representation. In this last step, we consider a correlated sizevelocity, such that, to each droplet of size S _{ i }, we attribute the velocity .
Thus, after substracting the disappearance flux from the moments, the momentum vector writes as follows:(43)where .
In the kinetic equation, the evolution of the size and the velocity for each abscissa is independent of the other abscissas.(44)and(45)
The solution of the velocity evolution after one time step is:(46)
Then the momentum and the final cell velocity are computed as follows:(47)
This method can be generalized for more complex evaporation laws by replacing in Equation (45) by general evaporation rate .
Even though it is possible to solve drag and evaporation separately and thus fully split the treatment of the three operators (transport, evaporation and drag), gathering the treatment of the two source terms is a relevant choice in terms of splitting error decrease, especially since this coupling is after all not a complex task.
3 AMR Techniques
In numerical simulations, especially for DNS, fine meshes can be needed in some locations in order to capture small scale effects. In fact, the presence of multiscale phenomena that interact with each other at variable positions over time, imposes to use a very small cell size for the whole domain to capture the smallest scale. In this case, the numerical simulation may require high computational resources. On the other hand, the use of adapted grid instead of uniform grid can greatly save CPU time and memory. In some applications, we can be interested to track the front of the flame, a plasma, or the interface in multiphase flows. Therefore mesh refinement is a useful tool to track these discontinuities/fronts and use fine meshes in these regions while keeping a coarser mesh elsewhere. For dilute sprays, the spatial concentration of the droplets is not uniform in the combustion chamber, and some locations are also not reached by the spray. Furthermore, in turbulent jets, the droplets are ejected from the core of high vortices and concentrate in low vorticity regions, creating large vacuum regions inside the chamber. Therefore, the mesh can remain coarse in these regions without affecting the accuracy of the simulation. However, the choice of an adequate refinement criterion is crucial to perform a good refinement/coarsening job, while limiting the numerical diffusion arising around shocks and high gradients.
We can briefly classify AMR approaches into two main categories:

Blockbased method Figure 1a: the domain is composed of blockstructures which are uniformly refined. This method is easier to implement;
Figure 1. Illustration of AMR techniques.

Cellbased method Figure 1b: in this method, a coarse mesh cell called the parent cell, is refined to 4 (2D) or 8 (3D) children cells, and recursively each quadrant (2D) or octant (3D) child can be refined in the same way. In this treestructure, only the leaves of the tree need to be stored.
In the present paper, we use a Cellbased method, managing the quadrants in a treestructure. The main advantage of Cellbased methods compared to Blockbased methods, is the high rate of compression allowed. Indeed the tree structure used to store the mesh in the Cellbased method, allows us to refine accurately only at the location of high variations (depending on the refinement criterion). But this heavy use of the tree structure to store and modify meshes continuously, complicates the software design and the implementation of numerical methods. Indeed, the need of flexible access to the information in each mesh, balancing load between computing processors, and optimum use of cache memory during a parallel AMR simulation should be taken into account carefully to design “powerful” AMR software. Recently, many new software propose AMR framework as Chombo [47], p4est [37], etc. In a recent project, Drui et al. [13] developed a new interface between a twophase flow finite volume solver and the cellbased library: p4est. The code shows a high performance computing, efficient parallel simulation and strong scaling. These first results were an encouragement to continue the development of a more general interface between p4est and finite volume schemes for several applications: separated phases, spray and astrophysics. This code was named CanoP, and the results presented in this paper, Section 4, are obtained using this code.
3.1 P4est Library
p4est is a cellbased AMR library written in C code. It provides several functions to manage adaptive meshes in parallel. In p4est, the spatial domain can first be divided into one or several subdomains, each of them being a conforming hexahedral macromesh (the cells are sharing only one edge (2D) or face (3D) with each other). In the tree structure, these macro meshes are represented by the root of a tree. Thus, the different macromeshes that cover the whole domain, constitute a forest of trees. The possibility of considering multitrees, where each one covers a subdomain, enables to represent complex geometries and not only simple square or cubic domain. In each tree, the corresponding subdomain is defined by a macromesh (the tree root) which represents the coarsen mesh and its corresponding refinement level is . The macromesh is then recursively refined to multiple non conforming micromeshes (where a cell can share a face or an edge with more than one neighboring cell) of high level , such that the mesh size of level is in p4est grid (the real size can be computed by a simple transformation defined initially by the user). Each micro mesh is associated to a unique reference in , where is the space dimension and is the maximum level of refinement. The generated tree is stored in a linear array, using the Morton space filling [48] (also called zorder curve, see Fig. 2) to map the leaves of the tree.
Figure 2.
zorder traversal of the quadrants in a tree and load partition into four processes. Dashed line: zorder curve. Quadrant label: zorder index. Color: MPI processes [13]. 
In the following, we present the main functionalities offered by p4est [37]:

create new refined forest and equipartition load between MPI processes;

iterate among each quadrant and call a callback function programmed by the user to mark the quadrants to be refined (creation of new children) or coarsened (remove the quadrants from the tree);

ensure the 2:1 balance (the size ratio between two neighbouring quadrant does not exceed two);

partition the load between the processes after each modification in the trees;

communicate the data of quadrants located at the boundaries of each “processordecomposition” domain to the neighboring processors, by using ghost quadrants.
3.2 CanoP Code
CanoP is C/C++ code dedicated to finite volume numerical AMR simulation for various applications (dispersed twophase model, separated twophase model and astrophysics applications). These applications involve multiscale phenomena. Therefore, CanoP uses the AMR p4est library to manage the meshes and benefits from its high performance.
CanoP can be presented at three stages:

the user: implements easily a new application by defining appropriate data model (for example: the density, the pressure, etc.), finite volume scheme and refinement criterion. CanoP provides a flexible and simple framework to develop new applications;

the CanoP core: adapts the user data model to its data structure, manages the time loop, calls p4est function to iterate over the meshes, then apply the numerical scheme and update the data and finally call the p4est functions for mesh adaptation and processor load balancing;

the libraries: CanoP integrates different libraries to benefit from their performance: p4est (AMR), HDF5 (Parallel output) and Lua (Input files).
Figure 3 illustrates the CanoP structure and the different stages.
Figure 3.
CanoP code structure. 
In the following, we present the finite volume scheme used on the nonconforming meshes and the refinement criterion developed in CanoP for the spray applications.
3.2.1 Kinetic Scheme in Non Conforming Meshes
The numerical scheme presented in Section 2.1 dedicated to the transport part, needs to be adapted to non conforming meshes. We consider a quadrant denoted by in a 2D non conforming mesh verifying the 2:1 size balance, which means that can have a common edge with at most two neighboring quadrants. We denote by its moment vector and by its momentum vector.
The finite volume scheme in writes as follows:(48)
To illustrate how we adapt the scheme in a 2:1 balanced meshe, we consider the case where the quadrant has two neighbour quadrants on his right, and as is illustrated in Figure 4. Therefore, the fluxes across the common edge are:(49)
Figure 4.
Example of non conforming mesh. 
The exact expression of the fluxes, derived from the integral expression of the moments and the kinetic solution (29), is:(50)such that(51)and(52)where (resp. ) is the set of droplets of quadrant that will reach its right (resp. left) side during the time .
Note that when we use a dimensional splitting operator, we consider the variation of the concerned variables only in the concerned direction (the xdirection). But in the case of non conforming meshes and especially for high order scheme, we need to consider the variation in the other perpendicular directions.
Since the scheme is first order, the solution is constant per quadrant and it is easy to derive the expression of the fluxes:(53)
The fluxes and are computed as in the case of a uniform meshe, using the expressions (33), (34). is the length of the quadrant .
For a higher order scheme, we use a polynomial reconstruction of the moments and an affine reconstruction for the velocity in each quadrant depending on the state of this quadrant and his neighbors. In order to keep the reconstructed moments within the space of realizability, this moment reconstruction is based on an affine reconstruction of the canonical moments [27]. Finally the reconstruction should ensure:

the conservation of the conservative variables,

a maximum principle for the canonical moments and the velocities.
In the case of non conforming meshes, there are different ways to reconstruct the solution based on the quadrant neighbor states and their positions. In CanoP, the solution is reconstructed a simple way. It is illustrated this in the context of Figure 4: quadrant has a doublesize neighbor on its left and two halfsize neighbors and on its right. The reconstructed variables (canonical moments and the velocity) are computed first using quadrants , then again, using this time . The final reconstructed state is simply the minmod limiter of the two computed slopes. Since quadrants of different size are used and the limitation of the slopes is rather restrictive, second order is not guaranteed. However the robustness and the realizability of the method are ensured. Moreover the results show a good accuracy improvement compared to the first order scheme.
3.2.2 Refinement Criterion
The determination of efficient criterion remains a complex subject, since the simulation depends on the physical phenomena and on the numerical methods. Multiresolution analysis [49–51] provides some robust tools to estimate and control rigorously the error at each refinement level. In CanoP, these methods are difficult to implement because the tree structure does not respect some necessary properties (like treaded tree, see [50]) for multiresolution methods. Therefore, we use a heuristic criterion based on some consideration of the main results cited in [50].
Harten [50] introduced a method to control the error, based on an interpolating function to approximate the variables at refined level from the coarsen level , then he computes the details which is the difference of the predictive variables and the true variables. The details are considered as the errors at the level . Following Harten’s idea, we implement a refinement criterion adapted to CanoP for spray applications. The method consists of predicting the moment for the children quadrants of a given parent quadrant . The prediction is based on the slope computation in each direction for each face depending on the neighbor quadrants (not necessary of the same size as the quadrant . But for Harten’s predictor operator, he used neighbors of the same level to predict the variable in the next high level):(54)
The subscripts refer respectively to the left, right, down and up edge. is the slope variation at the direction if or and otherwise. Finally, the children index follows the zcurve order.
Let us emphasize that, this reconstruction do not respect some physical and mathematical properties (the conservation, maximum principle and realizability). In fact, the reconstruction is only used to estimate the error and is not used in the numerical scheme.
The slopes are computed from the neighbors quadrant. To illustrate this calculation for the slopes at direction , let consider the case in the Figure 4:(55)
Then we compute the L ^{ 1 }norm difference between and the predicted .(56)where is the total surface of the quadrant .
This estimated error is used as refinement criterion, such that the error is compared to a given threshold and if then the quadrant is refined. If all siblings of a given quadrant verify then the quadrant is marked for coarsening. The parameter c ≤ 1 is chosen by the user.
4 Results
Once a model has been obtained with a clear relationship with geometrical quantities in separated phases averaged models and once a realizable, and thus robust, and accurate numerical scheme coupled to a mesh adaptation strategy has been designed, several steps have to be conducted in order to assess the proposed strategy. Beyond providing a verification of the proposed new geometrical moment and related algorithm compared to the EMSM model, we aim at proving that: 1 the AMR solver reaches a very nice level of accuracy once we have chosen a proper refinement criterion; 2 the proposed strategy is valid for complex droplet models, that is for source term with a much higher level of complexity leading to much higher arithmetic intensity; 3 that we have a high level of scalability and efficiency of the parallel implementation of the numerical strategy relying on realizable schemes and splitting.
To do so, we have chosen three testcases:

the first case consists in a localized spray in the presence of TaylorGreen (TG) vortices, Figure 5. These cases though simple, do illustrate a realistic configuration occurring in automotive engine, where the droplets are concentrated in some local regions and where large vacuum regions can be found. The first case is a non evaporating spray. For this case, we assess the performance of AMR calculation compared to a uniform grid calculation and study scalability for a given level of accuracy;
Figure 5. Stationary gaseous velocity vector field of the TaylorGreen vortices and spray initial number density.

the second case is the same as the first one, but with a major change. We consider an evaporating spray. The objective of this simulation is twofold; first we show the robustness and the accuracy of the evaporation algorithm coupled with transport and drag. Second, the evaporation algorithm is representative of the high arithmetic intensity within an embarrassingly parallel configuration using operator splitting we will encounter with more complex droplet models. We aim at illustrating the impact of the AMR techniques coupled to operator splitting on the computational cost for intensive source terms such as the inversion algorithm, used to reconstruct the approximated spray density function;

the third case is a non evaporating Homogeneous Isotropic Turbulence (HIT) in 2D and 3D. A comparison of the AMR Eulerian solution with a uniform Eulerian solution and a Lagrangian reference solution based on the segregation is performed to show the efficiency of the AMR Eulerian solution to predict correctly the main physical features with low computational resources in the framework of a much richer configuration in terms of physics, while in a really nonfavorable case since the spray is present in the whole domain.
For the three cases, we fix the CFL (CourantFriedrichsLewy) at CFL = 0.9 and the constant used in the refinement criterion at .
4.1 Droplet Cloud in TaylorGreen Vortices
We simulate a non evaporating polydisperse spray in two dimensional space configuration in the presence of TaylorGreen vortices for the gas: a steady solution of the inviscid incompressible Euler equations. The nondimensional velocity field of the gas is given by the following expression:(57)where .
Initially, we consider a motionless cloud located in the bottomleft vortex in Figure 5.
The initial spatial distribution of the spray is Gaussian inside a small disk of a radius and equals to zero outside. The initial size distribution, uniform inside the disk, is equal to for and otherwise. The moments are thus given by:(58)where , and .
The initial exact size distribution as well as its EM reconstructed density are presented in Figure 6. We consider a low inertia droplets such that the Stokes number is less than the critical value . It has been shown in de Chaisemartin [40], that for St < St_{ c } the droplets stay in their origin vortex and cannot travel from one to another, while, for St ≥ St_{ c } they are ejected out their original vortices. Thus, we encounter particle trajectory crossings, that are not taken into account in the monokinetic assumption and lead to deltashocks. Notice that these singularities are completely handled numerically, thanks to the ‘kinetic’ numerical scheme and the robustness is ensured.
Figure 6.
The initial NDF in the dashed line and its reconstruction through EM in the solid line. 
In Figure 7, we display the spatial evolution of the number density using a second order scheme and AMR grid at different times, as well as the mesh grid given for two instants in Figure 8. The maximum refinement level is and the minimum level is . We use the refinement criterion presented in Section 3.2.2 and choose a small threshold .
Figure 7.
TaylorGreen simulation using “second order scheme” in adaptive refinement grid, the maximum level is and the minimum level is . 
Figure 8.
The AMR grid in the case of TaylorGreen no evaporating spray, with , and threshold . 
4.1.1 Refinement Criterion and Solution Quality
In this section, we study the influence of the refinement criterion the accuracy of the solution and its impact on the computational time. The objective is to show that we can keep an accurate solution while reaching a high compression rate, thus we save time and memory. It also allows to set up a fair framework in order to conduct a parallel efficiency study in the next subsection. After investigating the influence of the refinement criterion depending on the order of the scheme used, we focus on the intrinsic accuracy of the solution.
In Tables 1 and 2, we present respectively the following results of the first order and the second order of TaylorGreen simulation using 36 processors: the error computed relatively to the solution on a refined uniform grid , considered as reference solution, using respectively the first order and second order. The compression rate and the CPU time are also provided. The purpose is clearly to evaluate the impact of the compression rate on both accuracy and computational efficiency.
Relative L ^{1}error, compression rate and CPU time depending on the threshold using the first order scheme
Relative L ^{1}error, compression rate and CPU time depending on the threshold using the second order scheme
The two tables show a significant CPU time saving compared to the reference solution, the CPU time of which is 29.15 s. Using the smallest threshold , we reduce the computation time by a factor of five with a very good level of precision. We can save more CPU time and memory (see the compression rate) by increasing the threshold but we lose the accuracy of the solution. For , the relative error is more than . Since the computing time of the second order volume finite solver is much larger than the communication and mesh adaptation time, we save more CPU time in this case, whereas the CPU time used in the computation of the second order reference solution is much larger: 231.63 s. By using a mesh adaptation with the threshold , we reduce this CPU time by a factor of eight, maintaining a very good precision. We can save more CPU time and memory by choosing with a good accuracy (relative L ^{1}error ∼ 5%). Let us underline that these two tables do not compare the solution accuracy depending on the scheme order because the relative error is computed for each order as the difference between the solution on a uniformly refined grid and the solution for an adaptive refined grid using the same order.
In a second study, we consider a common reference solution for the first and second order simulations. The reference solution is computed using a second order scheme on a uniform grid , which corresponds to the level . Figures 9 and 10 show that the relative error does not evolve much when we vary the minimum level except when we reach a uniform grid. However, for a given maximal level of refinement, there is a strong impact of on the ability of the mesh adapted solution (l _{ min } < l _{ max }) to capture an accurate solution. In fact for each maximal level of refinement, there is a threshold in terms of refinement criterion above which compressing the solution deteriorates its quality. Its is clear in Figures 9 and 10 that is below this threshold and is above it. For the relativeerror remains almost constant and we can use mesh refinement while preserving an accurate solution, whereas for , we can observe that mesh refinement for results in a clear increase of the error. Once we have found a well chosen threshold, for which the accuracy of the uniform and AMR solutions are equivalent, we can then study the computational efficiency and parallel performance of the numerical strategy.
Figure 9.
error for the first order scheme in logarithm scale versus the minimum level of compression plotted for different maximum refined levels . 
Figure 10.
error for second order scheme in logarithm scale versus the minimum level of compression plotted for different maximum refined levels . 
4.1.2 Parallel Performance
We assess the parallel performance by measuring the computational time of the various operations:

solver time corresponds to the time needed in the finite volume scheme and source evaluation in the splitting strategy;

adaptation time corresponds to the time of refinement, coarsening, partition and 2:1 size balance;

and finally the total time which is the sum of these two last times without I/O.
We use the second order TaylorGreen simulation in an adaptive refinement grid such that the maximum level l = 9 corresponds to the finest meshes and the minimum level l = 4 corresponds to the coarsest meshes. We choose the previously obtained threshold to ensure the quality of the compressed solution.
In Figure 11, we display the strong scaling. The simulation is performed in Icepar156q CentraleSuplec cluster (memory per node Go, total core number ) using MPI processors number up to . The solver time (transport and drag) shows an efficiency more than . The adaptation time shows a weak efficiency due to the large time spent in the communication between processors essentially for the balance algorithm. Nevertheless, the global time of computation shows a good parallel performance and the efficiency reaches for processors. In fact, the solver time represents of the total CPU time. Let us underline that the number of effective cells in this simulation with a compression rate of about is close to 30,000 so that an efficiency of 81% on 60 cores is very satisfactory.
Figure 11.
Strong scaling of the second order scheme in AMR grid, the maximum refinement level is and the minimum refinement level is . 
4.2 TaylorGreen Evaporating Spray
Using the exact same configuration with Stokes number , we switch to an evaporating spray with K = 1/2. Since we use an inversion algorithm in order to reconstruct the NDF from the moments by EM in each cell at each time step for the evaluation of the disappearance flux of droplet at zero size, we expect a significant computational cost related to source terms, which is representative of the complexity we will have to face by switching to more complex droplet models in realistic configurations.
Table 3 presents the computational time ratio spent in solving each part of the solver, transport and source, for two uniform grids (256 × 256 and 512 × 512) on processors. As expected, solving for the source term requires around of the total time and it is interesting to investigate how mesh adaptation is going to deal with such a configuration in the framework of operator splitting, that is in a numerical strategy where the source term in various cells are resolved completely independently leading to an embarrassingly parallel configuration.
Time ratio spent in transport solver and source solver (including evaporation and drag) using 96 MPI processes on a uniform mesh
In the previous study without evaporation, was shown to be an adequate threshold for a high compression and good solution quality, therefore we maintain this value in the present case. Figure 12 shows the spatial number density at four different times of the evaporating TaylorGreen simulation on AMR grid ( and ) and Figure 13 presents the corresponding meshes. Compared to the simulation running on the uniform grid , the relative error of the AMR solution computed at is , the compression rate is and we have reduced the computational time by a factor of more than . In fact, the CPU time used in the resolution of evaporation is highly larger than the CPU time that we spend in the resolution of the transport and the mesh adaptation. More precisely, in this simulation, the CPU time spent in the evaporation resolution is 76.5 s, the transport resolution 2.4 s and the mesh adaptation 5.2 s, which means that the evaporation alone takes more than of the total CPU time, Table 4 gives further details of the time ratio spending at each operation using different MPI process number. In Figure 14, we present the strong scaling of the various step of the solving process (mesh adaptation, transport solver, term source solver and the total time). Up to 60 processes, we observe a very good scalability, which is also to be seen in Table 4. Once again, let us insist on the fact that for 60 processes and considering the level of compression reached, we end up with about cells per core, which is very few, and still have a very nice level of efficiency. The transport part still has about 80% of efficiency and the source term scales vary properly considering that load balancing could be also a problem (p4est also provides the ability to integrate computational complexity in the the load balance algorithm, which has not yet been used in the present simulations). Increasing the number of processes to , we have a remarkable decrease of the scalability, which reaches , and that can be explained by the very small number of cells per process, which is going to kill the efficiency in terms of communications. However, the total efficiency is still about , which shows both the efficiency of CanoP and p4est, as well as the proper choice of numerical strategy.
Figure 12.
Evaporating TaylorGreen simulation using second order scheme for transport in AMR grid with , and . 
Figure 13.
The AMR grid in the case of TaylorGreen evaporating spray, with , and threshold . 
Figure 14.
Strong scaling of the TaylorGreen evaporated case in AMR grid, the maximum refinement level is and the minimum refinement level is . 
The time ratio spending in each operation: Transport solver and source term solver (Evaporation+drag) for different MPI process number
4.3 Verification of the Proposed Numerical Strategy
In order to verify the evaporation algorithm, we calculate the residual number density in the whole domain:(59)
We compute this quantity by integrating the number density obtained from the simulation using the evaporation algorithm presented in Section 2.2 and we compare it with the one obtained by a kinetic solution. The kinetic solution is based on the exact solution of the NDF. In Figure 15, the blue dashed line represents the number density of the kinetic solution using as initial NDF . The black line represents the kinetic solution using an EM reconstructed NDF from the initial moments. The red line represents the simulated solution using the evaporation algorithm. The figure shows a very good precision of the evaporation algorithm to predict the total remainder droplets compared to the exact kinetic solution and provide us with a verification of the proposed algorithm.
Figure 15.
Evolution of the number density over the domain compared to a exact solution using a kinetic solution of the evaporation. 
4.4 Homogeneous Isotropic Turbulence in 2D/3D
In order to assess the code ability in a more complex gaseous flow field representative of mode complex applications and in a challenging case, where spray is to be found everywhere in the domain, we consider a non evaporating spray in a two and threedimensional configurations with a frozen HIT gas. The HIT gas field was generated independently by Sabat et al. [52] with the ASPHODELE code of CORIA [53], which solves the threedimensional NavierStokes equations for the gas phase under the lowMach assumption. The characteristics of this HIT is given in Table 5, where is the turbulent Reynolds number, is the velocity rootmeansquare, is the mean dissipation rate, is the smallest structures scale, is the integral scale of the turbulence, is the Kolmogorov time scale of the turbulence, and is the eddy turnover time. We consider a no evaporating spray with Stokes number of . The initial spatial distribution is uniform.
Turbulence properties of the frozen HIT gaseous flow field
In order to limit the PTC and to have comparable results between the monokinetic Eulerian simulation and the Lagrangian simulation, we consider low inertia droplets such that the Stokes number is taken based on . Numerical dissipation in the simulation of Eulerian models for spray dynamics is a key issue, since we have to use very robust numerical schemes due to the presence of singularities and asymptotic limits (zero density), but still want an accurate resolution in the smooth regions. In particular, numerical dissipation has a negative effect in predicting properly segregation of the spray in turbulent flows. In order to limit this diffusion, we use the second order scheme and AMR in order to reach an accurate resolution in the high concentrated regions, while limiting the computational cost and memory trace. Figure 16 shows the spatial distribution of the number density at time (this dimensionless time corresponds to t = 1.8 s in the real time, such that the characteristic time in this case is τ _{ k } = 0.36 s).
Figure 16.
Number density at t = 1.8 s with AMR ( and ) and with threshold . 
The droplets are ejected from vortices and concentrated in zone of low vorticity.
In Figure 17, we compare the time evolution of the segregation for a Lagrangian simulation [52] and two Eulerian simulations:

on uniform grid 1024 × 1024;

using AMR with and .
Figure 17.
Evolution of the segregation with time for the Lagrangian (black line), Eulerian on uniform grid (blue line) and Eulerian with AMR (red line). 
In the two Eulerian simulations, we make a small error on the segregation compared to the Lagrangian simulation. However, the three results stay closely comparable. Furthermore, the segregation evolution of the two Eulerian simulations is closely the same. Therefore, we preserve the same solution quality while using less cells in AMR simulation. We have reduced the computational time by a factor of two compared to the uniform grid computation and we obtain a compression rate of . Let us underline that the chosen configuration is probably the worst configuration for such a strategy. However, we wanted to investigate the ability of the proposed strategy to cope with such a case, knowing that in realistic configurations, the flow is going to be turbulent, whereas the spray is going to occupy only a fraction of the computational domain.
A 3D simulation of the spray in the presence of a frozen HIT was performed to test the AMR capacity in compressing the solution. The characteristic of the HIT is given in Table 6. The spray simulation was performed in MPI processes using AMR grid, where and with the threshold . The compression rate is . The result was compared with a simulation obtained in uniform grid of cells. The relative error between the two solutions is . The computational time is reduced by a factor of compared to the uniform grid computation. In order to decrease the relative error, a second simulation was performed using . The relative error obtained in this case is , the compression rate is and the computational time is reduced by a factor of .
Turbulence properties of the frozen HIT gaseous flow field in 3D
Conclusion
The main contribution of this paper is threefold:

a new Eulerian moment model for the spray with a high potential for coupling with a separated phases twophase flow model has been proposed;

a robust and accurate numerical scheme has been derived for this new model including evaporation, transport and drag;

Adaptive Mesh Refinement (AMR) techniques have been introduced in this context to save both CPU time and memory.
In each stage of these developments, we have aimed at bringing solutions to predict an accurate physical simulation using low computational resources. Firstly, at the modeling stage, the kinetic model, which consists of the kinetic WilliamsBoltzmann equation, is reduced to an Eulerian moment model. Thus, the dimension of the phase space is reduced. The novelty in this article is the new geometrical moment model, which gives the dynamical evolution of some physically identified interfacial variables using fractional moments. This new model has been developed in order to consider a coupling with a diffuse interface model for the separated phases twophase flow models. Secondly, we have used operator splitting to separate the resolution of transport in physical space from the one of the source terms, which represents transport in phase space. For transport, we have extended the numerical scheme developed by Kah et al. [27]. For the source terms, a new evaporation solver has been developed for fractional moments and d ^{2} evaporation law. However, this algorithm can be adapted for more complex droplet models such as in [4]. This solver has been optimized by using fractional and negative quadrature moments to evaluate the size variation of the droplets due to evaporation. Finally, we have used AMR techniques to reach a fine mesh resolution in high variation regions while keeping low resolution in other regions. The results achieved by the CanoP code show that we can reach the same level of accuracy for AMR solution as uniform mesh solution with low computational time and memory. Furthermore, the use of AMR in a unified simulation of the fuel injection (from the separated phases to the evaporating spray) should allow to have an accurate simulation of this multiscale problem with reasonable computational need. Indeed, for diffuse interface model, we need fine meshes to limit the interface diffusion. Moreover, in the fuel injection, we encounter some vacuum region and some high concentrated droplets region. Therefore, a large spectrum of size is involved in this multiscale problem and the use of AMR grid should be very beneficial.
To reach the ultimate goal of designing adequate model, numerical scheme and HPC application for a fully Eulerian simulation of fuel injection, we still need a deep study to understand this kind of flow, then, to improve the present model for more realistic case. First, we aim at extending the actual spray model to take into account the sizevelocity correlation. Furthermore, to tackle the PTC, which will certainly occur in transition zone (because of inertial drops), we would like to consider other closing assumption for the velocity distribution. These types of model have already been studied in [31, 41]. Secondly, the high order schemes used in the present work, still require a careful treatment and study in non conforming meshes. Finally, we need to focus on understanding the mechanism of the interface evolution in separated phases twophase flow, which is the key mechanism of generating polydisperse spray.
Acknowledgments
This research was supported by Ph.D. grant for Mohamed Essadki from IFP Energies nouvelles and EM2C. We would like to thank our colleagues from La Maison de Simulation (Samuel Kokh, Florence Drui and Pierre Kestener) for their collaboration in developing CanoP code, Florence Drui for giving helpful comments on writing this paper and Aymeric Vié for several interesting discussions. The help of Macole Sabat in providing the Lagrangian numerical results of HIT is also gratefully acknowledged to validate the numerical results.
References
 Hirt C.W., Amsden A.A., Cook J.L. (1974) An arbitrary lagrangianeulerian computing method for all flow speeds, Journal of Computational Physics 14, 3, 227–253. [Google Scholar]
 Harlow F., Welch J.E. (1965) Numerical calculation of timedependent viscous incompressible flow of fluid with free surface, Physics of Fluid 8, 12, 2182–2189. [Google Scholar]
 Le Chenadec V., Pitsch H. (2013) A 3D unsplit forward/backward volumeoffluid approach and coupling to the level set method, Journal of Computational Physics 233, 10–33. [CrossRef] [Google Scholar]
 Emre O. (2014) Modélisation de la polydispersion des brouillards de gouttes sous l’effet des interactions twoway turbulentes pour l’injection directe à haute pression dans les moteurs. PhD Thesis, École Centrale Paris, Available at https://tel.archivesouvertes.fr/tel01089937. [Google Scholar]
 Lebas R., Menard T., Beau P.A., Berlemont A., Demoulin F.X. (2009) Numerical simulation of primary breakup and atomization: DNS and modelling study, Int. J. Multiphase Flows 35, 3, 247–260. [CrossRef] [Google Scholar]
 Menard T., Tanguy S., Berlemont A. (2007) Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary breakup of a liquid jet, International Journal of Multiphase Flow 33, 5, 510–524. [CrossRef] [Google Scholar]
 Saurel R., Lemetayer O. (2001) A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation, Journal of Fluid Mechanics 431, 2, 239–271. [CrossRef] [Google Scholar]
 Drew D.A., Passman S.L. (1999) Theory of multicomponent fluids, Appl. Math. Sci. 135. [Google Scholar]
 Baer M.R., Nunziato J.W. (1986) A twophase mixture theory for the deflagrationtodetonation transition (DDT) in reactive granular materials, International Journal of Multiphase Flow 12, 6, 861–889. [Google Scholar]
 Drui F., Kokh S., Larat A., Massot M. (2016) A hierarchy of simple hyperbolic twofluid models for bubbly flows, Submitted to Physics of Fluids, pp. 1–40. [Google Scholar]
 Le Martelot S., Saurel R., Nkonga B. (2014) Towards the direct numerical simulation of nucleate boiling flows, International Journal of Multiphase Flow 66, 62–78. [CrossRef] [MathSciNet] [Google Scholar]
 Le Touze C. (2015) Coupling between separated and dispersed twophase flow models for the simulation of primary atomization in cryogenic combustion, Theses, Universite Nice Sophia Antipolis, December. URL https://tel.archivesouvertes.fr/tel01250527 [Google Scholar]
 Drui F., Fikl A., Kestener P., Kokh S., Larat A., Le Chenadec V., Massot M. (2016) Experimenting with the p4est library for AMR simulations of twophase flows. to appear in ESAIM: Proceedings and Surveys, pp. 1–16. [Google Scholar]
 Saurel R., Abgrall R. (1999) A multiphase godunov method for compressible multifluid and multiphase flows, Journal of Computational Physics 150, 2, 425–467. [CrossRef] [MathSciNet] [Google Scholar]
 Saurel R., Petitpas F., Berry R.A. (2009) Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures, Journal of Computational Physics 228, 1678–1712. [Google Scholar]
 Jay S., Lacas F., Candel S. (2006) Combined surface density concepts for dense spray combustion, Combustion and Flame 144, 3, 558–577. [CrossRef] [Google Scholar]
 Vallet A., Burluka A., Borghi R. (2001) Development of a eulerian model for the “atomization” of a liquid jet, Atomization and Sprays 11, 619–642. [CrossRef] [Google Scholar]
 Williams F.A. (1958) Spray combustion and atomization, Physics of Fluids 1, 541–545. [CrossRef] [Google Scholar]
 Bird G.A. (1994) Molecular gas dynamics and the direct simulation of gas flows, Oxford Science Publications, 42. [Google Scholar]
 Garcia M. (2009) Development and validation of the Euler Lagrange formulation on a parallel and unstructured solver for largeeddy simulation, PhD Thesis, Université Toulouse III, Available online at http://ethesis.inptoulouse.fr/archive/00000715/. [Google Scholar]
 Kah D. (2010) Taking into account polydispersity in the framework of a coupled EulerLagrange approach for the modeling of liquid fuel injection in internal combustion engines, PhD Thesis, École Centrale de Paris, Available online at http://tel.archivesouvertes.fr/tel00618786/en/. [Google Scholar]
 Greenberg J.B., Silverman I., Tambour Y. (1993) On the origin of spray sectional conservation equations, Combustion and Flame 93, 90–96. [CrossRef] [Google Scholar]
 Laurent F., Massot M. (2001) Multifluid modeling of laminar polydispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods, Combust, Theory and Modelling 5, 537–572. [CrossRef] [Google Scholar]
 Emre O., Kah D., Jay S., Tran Q.H., Velghe A., De Chaisemartin S., Laurent F., Massot M. (2015) Eulerian Moment Methods for Automotive Sprays, Atomization and Sprays 25, 189–254. [CrossRef] [Google Scholar]
 Kah D., Emre O., Tran Q.H., de Chaisemartin S., Jay S., Laurent F., Massot M. (2015) High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines, International Journal of Multiphase Flows 71, 38–65. [CrossRef] [Google Scholar]
 Nguyen T.T., Laurent F., Fox R.O., Massot M. Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes, Journal of Computational Physics, pp. 1–42. URL https://hal.archivesouvertes.fr/hal01247390. Submitted. [Google Scholar]
 Kah D., Laurent F., Massot M., Jay S. (2012) A high order moment method simulating evaporation and advection of a polydisperse spray, J. Comput. Phys. 231, 2, 394–422. [CrossRef] [Google Scholar]
 Laurent F. (2006) Numerical analysis of eulerian multifluid models in the context of kinetic formulations for dilute evaporating sprays, Mathematical Modeling and Numerical Analysis 3, 431–468. [CrossRef] [EDP Sciences] [Google Scholar]
 Laurent F., Sibra A., Doisneau F. (2016) Twosize moment Eulerian multifluid model: a flexible and realizable highfidelity description of polydisperse moderately dense evaporating sprays, Communications in Computational Physics, pp. 1–42, URL https://hal.archivesouvertes.fr/hal01169730. Accepted. [Google Scholar]
 Massot M., Laurent F., Kah D., de Chaisemartin S. (2010) A robust moment method for evaluation of the disappearance rate of evaporating sprays, SIAM J. Appl. Math. 70, 3203–3234. [CrossRef] [Google Scholar]
 Vié A., Laurent F., Massot M. (2013) Sizevelocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues, J. Comp. Phys. 237, 277–310. [Google Scholar]
 Emre O., Fox R.O., Massot M., de Chaisemartin S., Jay S., Laurent F. (2014) Towards eulerian modeling of a polydisperse evaporating spray under realistic internalcombustionengine conditions, Flow, Turbulence and Combustion 93, 689–722. [CrossRef] [Google Scholar]
 Devassy M.B., Habchi C., Danoem E. (2015) Atomization modelling of liquid jets using a twosurfacedensity approach, Journal of Atomization and Sprays 25, 1, 47–80. [CrossRef] [Google Scholar]
 Sibra A. (2015) Eulerian MultiFluid modeling and simulation of evaporation and combustion of polydisperse sprays in solid rocket motors, Theses, Université ParisSaclay, November. URL https://tel.archivesouvertes.fr/tel01260314. [Google Scholar]
 Essadki M., de Chaisemartin S., Laurent F., Larat A., Massot M. (2016) Topological moment model for polydisperse evaporating sprays, Submitted to SIAM Journal on Applied Mathematics. [Google Scholar]
 Drew D.A. (1990) Evolution of geometric statistics, SIAM J. Appl. Math. 50, 3, 649–666. [CrossRef] [Google Scholar]
 Burstedde C., Wilcox L.C., Ghattas O. (2011) p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM Journal on Scientific Computing 33, 3, 1103–1133. [Google Scholar]
 Sabat M. (2016) Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays, PhD Thesis, Université ParisSaclay, CentraleSupélec. [Google Scholar]
 Godsave G.A.E. (1953) Studies of the combustion of drops in a fuel spray: the burning of single drops of fuel, Proceedings of the 4th Symp. (International) on Combustion, The Comb. Institute, Baltimore, pp. 818–830. [Google Scholar]
 de Chaisemartin S. (2009) Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporation sprays, PhD Thesis, École Centrale Paris, France, Available on TEL: http://tel.archivesouvertes.fr/tel00443982/en/. [Google Scholar]
 Vié A., Doisneau F., Massot M. (2015) On the Anisotropic Gaussian closure for the prediction of inertialparticle laden flows, Communication in Computational Physics 17, 1, 1–46. [CrossRef] [Google Scholar]
 de Chaisemartin S., Fréret L., Kah D., Laurent F., Fox R.O., Reveillon J., Massot M. (2009) Eulerian models for turbulent spray combustion with polydispersity and droplet crossing, Comptes Rendus Mécanique 337, 438–448, Special Issue ‘Combustion for Aerospace Propulsion’. [CrossRef] [Google Scholar]
 Mead L.R., Papanicolaou N. (1984) Maximum entropy in the problem of moments, J. Math. Phys. 25, 8, 2404–2417. ISSN 00222488. [NASA ADS] [CrossRef] [Google Scholar]
 Descombes S., Massot M. (2004) Operator splitting for nonlinear reactiondiffusion systems with an entropic structure: singular perturbation and order reduction, Numer. Math. 97, 4, 667–698. [CrossRef] [MathSciNet] [Google Scholar]
 Bouchut F., Jin S., Li X. (2003) Numerical approximations of pressureless and isothermal gas dynamics, SIAM J. Num. Anal. 41, 135–158. [CrossRef] [Google Scholar]
 Gordon R.G. (1968) Error bounds in equilibrium statistical mechanics, Journal of Mathematical Physics 9, 655–663. [CrossRef] [Google Scholar]
 Adams M., Colella P., Graves D.T., Johnson J.N., Keen N.D., Ligocki T.J., Martin D.F., McCorquodale P.W., Modiano D., Schwartz P.O., Sternberg T.D., Van Straalen B. (2013) Chombo software package for amr applications  design document. Technical Report LBNL6616E, Lawrence Berkeley National Laboratory. [Google Scholar]
 Morton G.M. (1966) A computer oriented geodetic data base and a new technique in file sequencing. Technical report, Ottawa, Ontario, Canada. [Google Scholar]
 Müller S. (2003) Adaptive Multiscale Schemes for Conservation Laws, Springer. [CrossRef] [Google Scholar]
 Harten A. (1995) Multiresolution algorithms for the numerical solution of hyperbolic conservation laws, Communication on Pure Applied Mathematics 48, 1305–1342. [CrossRef] [MathSciNet] [Google Scholar]
 Duarte M. (2011) Adaptive numerical methods in time and space for the simulation ofmultiscale reaction fronts, PhD Thesis, École Centrale de Paris. Available online at https://tel.archivesouvertes.fr/tel00667857. [Google Scholar]
 Sabat M., Larat A., Vié A., Massot M. (2014) On the development of high order realizable schemes for the Eulerian simulation of disperse phase flows: a convexstate preserving Discontinuous Galerkin method, Journal of Computational Multiphase Flows 6, 3, 247–270. [CrossRef] [MathSciNet] [Google Scholar]
 Reveillon J., Demoulin F.X. (2007) Effects of the preferential segregation of droplets on evaporation and turbulent mixing, Journal of Fluid Mechanics 583, 273–302. [CrossRef] [Google Scholar]
Cite this article as: M. Essadki, S. de Chaisemartin, M. Massot, F. Laurent, A. Larat and S. Jay (2015). Adaptive Mesh Refinement and High Order Geometrical Moment Method for the Simulation of Polydisperse Evaporating Sprays, Oil Gas Sci. Technol 71, 61.
All Tables
Relative L ^{1}error, compression rate and CPU time depending on the threshold using the first order scheme
Relative L ^{1}error, compression rate and CPU time depending on the threshold using the second order scheme
Time ratio spent in transport solver and source solver (including evaporation and drag) using 96 MPI processes on a uniform mesh
The time ratio spending in each operation: Transport solver and source term solver (Evaporation+drag) for different MPI process number
All Figures
Figure 1.
Illustration of AMR techniques. 

In the text 
Figure 2.
zorder traversal of the quadrants in a tree and load partition into four processes. Dashed line: zorder curve. Quadrant label: zorder index. Color: MPI processes [13]. 

In the text 
Figure 3.
CanoP code structure. 

In the text 
Figure 4.
Example of non conforming mesh. 

In the text 
Figure 5.
Stationary gaseous velocity vector field of the TaylorGreen vortices and spray initial number density. 

In the text 
Figure 6.
The initial NDF in the dashed line and its reconstruction through EM in the solid line. 

In the text 
Figure 7.
TaylorGreen simulation using “second order scheme” in adaptive refinement grid, the maximum level is and the minimum level is . 

In the text 
Figure 8.
The AMR grid in the case of TaylorGreen no evaporating spray, with , and threshold . 

In the text 
Figure 9.
error for the first order scheme in logarithm scale versus the minimum level of compression plotted for different maximum refined levels . 

In the text 
Figure 10.
error for second order scheme in logarithm scale versus the minimum level of compression plotted for different maximum refined levels . 

In the text 
Figure 11.
Strong scaling of the second order scheme in AMR grid, the maximum refinement level is and the minimum refinement level is . 

In the text 
Figure 12.
Evaporating TaylorGreen simulation using second order scheme for transport in AMR grid with , and . 

In the text 
Figure 13.
The AMR grid in the case of TaylorGreen evaporating spray, with , and threshold . 

In the text 
Figure 14.
Strong scaling of the TaylorGreen evaporated case in AMR grid, the maximum refinement level is and the minimum refinement level is . 

In the text 
Figure 15.
Evolution of the number density over the domain compared to a exact solution using a kinetic solution of the evaporation. 

In the text 
Figure 16.
Number density at t = 1.8 s with AMR ( and ) and with threshold . 

In the text 
Figure 17.
Evolution of the segregation with time for the Lagrangian (black line), Eulerian on uniform grid (blue line) and Eulerian with AMR (red line). 

In the text 