© M. Essadki et al., published by IFP Energies nouvelles, 2016

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 high-cost experiments, numerical simulation is nowadays considered as a good alternative in order to improve our understanding of two-phase flows and two-phase 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 two-phase 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 Navier-Stokes 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), level-set 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, reduced-order models have to be considered since a full resolution is out of reach. So far two reduced-order 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 reduced-order 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 William-Boltzmann (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 Monte-Carlo 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 load-balancing 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 size-velocity moments of this NDF. Considering some assumptions on the velocity distribution (monokinetic, Maxwell-Boltzmann distribution, etc.), we are able to close the equation in the velocity direction, and derive a semi-kinetic 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 multi-fluid 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 reduced-order 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 reduced-order 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 multi-scale 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 test-cases 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 reduced-order 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 non-conforming 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 break-up of droplets. Then, the dynamics of the droplets and their interaction with the surrounding gas is given by the William-Boltzmann (WB) equation, [18].$${\mathrm{\partial }}_{t}f+{\mathrm{\partial }}_{\vec{x}}\cdot \left(\vec{v}f\right)+{\mathrm{\partial }}_{\vec{v}}\cdot \left(\vec{F}f\right)+{\mathrm{\partial }}_S\left({\mathrm{R}}_{S}f\right)=0$$(1)

where is the evaporation rate and is the drag acceleration. The evaporation is modeled by a simple d ² 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:$$\vec{F}=\frac{{\vec{u}}_g-\vec{v}}{{\tau }_p\left(S\right)}$$(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:$${\tau }_p\left(S\right)=\frac{{\rho }_{l}S}{18\pi {\mu }_g}$$(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:$${\mathrm{\partial }}_{t}f+{\mathrm{\partial }}_{\vec{x}}\cdot \left(\vec{v}f\right)+{\mathrm{\partial }}_{\vec{v}}\cdot \left(\frac{{\vec{u}}_g-\vec{v}}{\mathrm{St}\left(S\right)}f\right)-\mathrm{K} {\mathrm{\partial }}_{S}f=0$$(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):$$M_{i,j,k,l}={\int }_0^1{\int }_{{\mathbb{R}}^3}^{}{S}^l{v}_x^i{v}_y^j{v}_z^{k}f\left(t,\vec{x},\vec{v},S\right)\mathrm{dS}{d}^3\vec{v}$$(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:$$f\left(t,\vec{x},\vec{v},S\right)=n\left(t,\vec{x},S\right)\varphi \left(\vec{v}-\vec{u}\left(t,\vec{x},S\right)\right), \int \varphi \left(\vec{v}\right)d^3\vec{v}=1$$(6)where is the mean velocity conditioned by size. In this paper, we even simplify our study to the case of size-independent macroscopic velocity . However, it is worth noticing that such a complete size-velocity correlation has already been considered and tackled in the CSVM (Coupled Size-Velocity 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 Maxwell-Boltzmann 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 delta-distribution:$$f\left(t,\vec{x},\vec{v},S\right)=n\left(t,\vec{x},S\right)\delta \left(\vec{v}-\vec{u}\left(t,\vec{x}\right)\right)$$(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 delta-shocks by what we call the Monokinetic-Equilibrium 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 semi-kinetic system for Equation (4):$$\begin{array}{c}{\mathrm{\partial }}_{t}n+{\mathrm{\partial }}_{\vec{x}}\cdot \left(n\vec{u}\right)=\mathrm{K} {\mathrm{\partial }}_{S}n\\ {\mathrm{\partial }}_{t}n\vec{u}+{\mathrm{\partial }}_{\vec{x}}\cdot \left(n\vec{u}\otimes \vec{u}\right)=\mathrm{K} {\mathrm{\partial }}_S\left(n\vec{u}\right)+n\frac{{\vec{u}}_g-\vec{u}}{\mathrm{St}\left(S\right)}\end{array}$$(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 multi-fluid model [23, 42], the number density is considered to be constant in each section and the computation then requires the calculation of one size-moment 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 size-moments and reconstruct a continuous size distribution function n(S) from this set of high order size-moments 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 two-size 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 multi-fluid 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:$$m_k={\int }_0^{S_{\max }}{S}^{k}n\left(S\right)\mathrm{dS}$$(9)

Using non-dimensional variables, S _max can be set to S _max  = 1.

Starting from the semi-kinetic system (8), we derive the system on N size-moments:$$\begin{array}{c}{\mathrm{\partial }}_t{m}_0+{\mathrm{\partial }}_{\vec{x}}.\left(m_0\vec{u}\right)={\psi }_{+}^0-{\psi }_{-}^0\\ {\mathrm{\partial }}_t{m}_1+{\mathrm{\partial }}_{\vec{x}}.\left(m_1\vec{u}\right)={\psi }_{+}^1-{\psi }_{-}^1-\mathrm{K} m_0\\ {\mathrm{\partial }}_t{m}_2+{\mathrm{\partial }}_{\vec{x}}.\left(m_2\vec{u}\right)={\psi }_{+}^2-{\psi }_{-}^2-2\mathrm{K} m_1\\ ⋮\\ {\mathrm{\partial }}_t{m}_N+{\mathrm{\partial }}_{\vec{x}}.\left(m_N\vec{u}\right)={\psi }_{+}^N-{\psi }_{-}^N-N\mathrm{K} m_{N-1}\\ {\mathrm{\partial }}_t\left(m_1\vec{u}\right)+{\mathrm{\partial }}_{\vec{x}}.\left(m_1\vec{u}\otimes \vec{u}\right)=-\mathrm{K} m_0\vec{u}+\left({\psi }_{+}^1-{\psi }_{-}^1\right)\vec{u}+m_0\frac{{\vec{u}}_g-\vec{u}}{\theta }\end{array}$$(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 Monokinetic-Equilibrium assumption, we only need one velocity moment to compute the velocity. In this case, we have considered:$$m_1\vec{u}\left(t,\vec{x}\right)={\int }_0^1{\int }_{{\mathbb{R}}^d}^{}S\vec{v}f\left(t,\vec{x},\vec{v},S\right)\mathrm{dSd}\vec{v}$$(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:$$H\left(n\right)=-{\int }_0^{1}n\left(s\right)\ln \left(n\left(s\right)\right)\mathrm{ds}$$(12)

The existence and uniqueness of a density function n which maximizes the Shannon entropy and satisfies:$$\begin{array}{c}{m}_0={\int }_0^{1}n\left(S\right)\mathrm{dS}\\ ⋮\\ m_N={\int }_0^1{S}^{N}n\left(S\right)\mathrm{dS}\\ \end{array}$$(13)was proved in [43], and the solution is shown to have the following form:$$n\left(S\right)=\exp \left(-\left({\lambda }_0+{\lambda }_{1}S+\dots +{\lambda }_N{S}^N\right)\right)$$(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 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 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

Two-phase 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 two-phase flow topologies, going from separated two-phase 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:$${\chi }_k(t,\vec{x})=\{\begin{array}{ll}1,& \mathrm{if} \vec{x} \mathrm{is} \mathrm{in} \mathrm{phase} k\\ 0,& \mathrm{otherwise}\end{array}$$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 delta-function 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:

$${\alpha }_k(t,\vec{x})=\frac{1}{\left|V\right|}{\int }_V^{}{\chi }_k(t,\vec{x})\mathrm{dV}\left(\vec{x}\right)$$(15)

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:

$$\mathrm{\Sigma }\left(t,\vec{x}\right)=\frac{1}{\left|V\right|}{\int }_V^{}\left|\left|\vec{\nabla }{\chi }_k\left(t,\vec{x}\right)\right|\right|\mathrm{dV}\left(\vec{x}\right)$$(16)

• 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 :

$$\begin{array}{ccc}H& =& \frac{1}{2}\left(k_1+k_2\right)\\ G& =& k_1{k}_2\end{array}$$(17)

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:$$\begin{array}{c}\mathrm{\Sigma }\tilde{H}={\int }_V^{}H\left|\left|\vec{\nabla }{\chi }_k\left(t,\vec{x}\right)\right|\right|\mathrm{dV}\left(\vec{x}\right)\\ \mathrm{\Sigma }\tilde{G}={\int }_V^{}G\left|\left|\vec{\nabla }{\chi }_k\left(t,\vec{x}\right)\right|\right|\mathrm{dV}\left(\vec{x}\right)\\ \end{array}$$(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:$$\begin{array}{c}{\mathrm{\Sigma }}_d{G}_d=4\pi m_0\\ {\mathrm{\Sigma }}_d{H}_d=2\sqrt{\pi }{m}_{1/2}\\ {\mathrm{\Sigma }}_d=m_1\\ {\alpha }_d=\frac{1}{6\sqrt{\pi }}{m}_{3/2}\end{array}$$(19)where the half-integer moments are defined by:$$m_{k/2}={\int }_0^1{S}^{k/2}n\left(S\right)\mathrm{dS}$$(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:$$\{\begin{array}{ll}{\mathrm{\partial }}_t{m}_0+{\mathrm{\partial }}_{\vec{x}}\cdot \left(m_0\vec{u}\right)=-{\psi }_{-}^0& \\ {\mathrm{\partial }}_t{m}_{1/2}+{\mathrm{\partial }}_{\vec{x}}\cdot \left(m_{1/2}\vec{u}\right)=-{\psi }_{-}^{1/2}-\frac{\mathrm{K}}{2}{m}_{-1/2}& \\ {\mathrm{\partial }}_t{m}_1+{\mathrm{\partial }}_{\vec{x}}\cdot \left(m_1\vec{u}\right)=-{\psi }_{-}^1-\mathrm{K} m_0& \\ {\mathrm{\partial }}_t{m}_{1/2}+{\mathrm{\partial }}_{\vec{x}}\cdot \left(m_{3/2}\vec{u}\right)=-{\psi }_{-}^{3/2}-\frac{3\mathrm{K}}{2}{m}_{1/2}& \\ {\mathrm{\partial }}_t\left(m_1\vec{u}\right)+{\mathrm{\partial }}_{\vec{x}}\cdot \left(m_1\vec{u}\otimes \vec{u}\right)=-{\psi }_{-}^1\vec{u}-\mathrm{K} m_0\vec{u}+m_0\frac{{\vec{u}}_g-\vec{u}}{\theta }& \\ & \\ & \end{array}$$(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:$$m_{-1/2}={\int }_0^1{S}^{-1/2}n\left(S\right)\mathrm{dS}$$(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:$$n^{\mathrm{EM}}\left(S\right)=\exp \left(-\left({\lambda }_0+{\lambda }_1{S}^{1/2}+{\lambda }_{2}S+{\lambda }_3{S}^{1/2}\right)\right)$$(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 operator-splitting 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 one-dimensional 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 .$$\begin{array}{c}{\mathrm{\partial }}_t{m}_0+{\mathrm{\partial }}_x\left(m_{0}u\right)=0\\ {\mathrm{\partial }}_t{m}_{1/2}+{\mathrm{\partial }}_x\left(m_{1/2}u\right)=0\\ {\mathrm{\partial }}_t{m}_1+{\mathrm{\partial }}_x\left(m_{1}u\right)=0\\ {\mathrm{\partial }}_t{m}_{3/2}+{\mathrm{\partial }}_x\left(m_{3/2}u\right)=0\\ {\mathrm{\partial }}_t\left(m_{1}u\right)+{\mathrm{\partial }}_x\left(m_1{u}^2\right)=0\\ {\mathrm{\partial }}_t\left(m_{1}v\right)+{\mathrm{\partial }}_x\left(m_1\mathrm{uv}\right)=0\\ \end{array}$$(24)

The above differential system is not strictly hyperbolic, and the semi-kinetic system from which the moment system is derived, is similar to the pressureless gas system:$$\{\begin{array}{c}{\mathrm{\partial }}_{t}n+{\mathrm{\partial }}_x\mathrm{nu}=0\\ {\mathrm{\partial }}_t\left(\mathrm{nu}\right)+{\mathrm{\partial }}_x\left(n{u}^2\right)=0\\ {\mathrm{\partial }}_t\left(\mathrm{nv}\right)+{\mathrm{\partial }}_x\left(\mathrm{nuv}\right)=0\end{array}$$(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:$$\{\begin{array}{l}{\mathrm{\partial }}_{t}f+{\mathrm{\partial }}_x\left(c_{x}f\right)=0\\ f\left(t,x,\vec{c},S\right)=n\left(t,x,S\right)\delta \left(c_x-u\right)\delta \left(c_y-v\right)\end{array}$$(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:$$\begin{array}{ccc}{m}_{k/2, i}^n& =& \frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{m}_{k/2}\left(t_n,x\right)\mathrm{dx}\\ & =& \frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{\int }_0^1{\int }_{{\mathbb{R}}^2}^{}{S}^{k/2}f\left(t_n,x,c_x,c_y,S\right)d{c}_{x}d{c}_y\mathrm{dSdx}\end{array}$$(27) $$\begin{array}{ccc}{\vec{p}}_i^n& =& m_{1,i}^n\left(\begin{array}{l}{u}_i^n\\ v_i^n\end{array}\right)\\ & =& \frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{m}_1\left(t_n,x\right)\left(\begin{array}{l}u\left(t_n,x\right)\\ v\left(t_n,x\right)\end{array}\right)\mathrm{dx}\\ & =& \frac{1}{\Delta x}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{\int }_0^1{\int }_{{\mathbb{R}}^2}^{}S\left(\begin{array}{l}{c}_x\\ c_y\end{array}\right)f\left(t_n,x,c_x,c_y,S\right)d{c}_{x}d{c}_y\mathrm{dS}\\ & & \end{array}$$(28)

Using the exact solution of the kinetic equation (26) for :$$f\left(t,x,\vec{c},S\right)=f\left(t_n,x-\left(t-t_n\right)c_x,S\right)$$(29)we establish the following scheme using the same steps as in [27]:$$\begin{array}{ccc}{M}_i^{n+1}& =& M_i^n-\frac{\Delta t}{\Delta x}\left(F_{i+1/2}-F_{i-1/2}\right)\\ p_i^{n+1}& =& p_i^n-\frac{\Delta t}{\Delta x}\left(G_{i+1/2}-G_{i-1/2}\right)\end{array}$$(30)where the fluxes can be decomposed into and such that:$$\left(\begin{array}{l}{F}_{i+1/2}^{+}\\ G_{i+1/2}^{+}\end{array}\right)=\frac{1}{\Delta t}{\int }_{x_{i-1/2}}^{x_{i+1/2}}\left(\begin{array}{l}{m}_0\left(t_n,x\right)\\ m_{1/2}\left(t_n,x\right)\\ m_1\left(t_n,x\right)\\ m_{3/2}\left(t_n,x\right)\\ m_{1}u\left(t_n,x\right)\\ m_{1}v\left(t_n,x\right)\end{array}\right)I_{\left\{x^{\mathrm{\prime }},x_{i+1/2}-\Delta \mathrm{tu}\left(t_n,x^{\mathrm{\prime }}\right)<x^{\mathrm{\prime }}\right\}}\left(x\right)\mathrm{dx}$$(31)and$$\left(\begin{array}{l}{F}_{i+1/2}^{-}\\ G_{i+1/2}^{-}\end{array}\right)=\frac{1}{\Delta t}{\int }_{x_{i+1/2}}^{x_{i+3/2}}\left(\begin{array}{l}{m}_0\left(t_n,x\right)\\ m_{1/2}\left(t_n,x\right)\\ m_1\left(t_n,x\right)\\ m_{3/2}\left(t_n,x\right)\\ m_{1}u\left(t_n,x\right)\\ m_{1}v\left(t_n,x\right)\end{array}\right)I_{\left\{x^{\mathrm{\prime }},x_{i+1/2}-\Delta \mathrm{tu}\left(t_n,x^{\mathrm{\prime }}\right)>x^{\mathrm{\prime }}\right\}}\left(x\right)\mathrm{dx}$$(32)

For a first order scheme, we consider a constant piecewise reconstruction for the moments and the velocities. Then the fluxes are:$$F_{i+1/2}=\left(\begin{array}{l}{m}_{0,i}^n\\ m_{1/2,i}^n\\ m_{1,i}^n\\ m_{3/2,i}^n\\ \end{array}\right)\max (u_i^n,0)+\left(\begin{array}{l}{m}_{0,i+1}^n\\ m_{1/2,i+1}^n\\ m_{1,i+1}^n\\ m_{3/2,i+1}^n\\ \end{array}\right)\min (u_{i+1}^n,0)$$(33)and$$G_{i+1/2}=\left(\begin{array}{l}{m}_{1,i}^n{u}_{0,i}^n\\ m_{1,i}^n{v}_i^n\\ \end{array}\right)\max \left(u_i^n,0\right)+\left(\begin{array}{l}{m}_{1,i+1}^n{u}_{0,i+1}^n\\ m_{1,i}^n{v}_{i+1}^n\\ \end{array}\right)\min \left(u_{i+1}^n,0\right)$$(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 $${\mathrm{\partial }}_{t}n-{\mathrm{\partial }}_S\left(\mathrm{K} n\right)=0.$$(35)

By integrating this equation, we get the following Ordinary Differential Equations (ODE) for the moments evolution:$$\begin{array}{ccc}{d}_t{m}_0& =& -n\left(S=0\right)\\ d_t{m}_{1/2}& =& -\frac{\mathrm{K}}{2}{m}_{-1/2}\\ d_t{m}_1& =& -\mathrm{K} m_0\\ d_t{m}_{3/2}& =& -\frac{3\mathrm{K}}{2}{m}_{1/2}\end{array}$$(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 Runge-Kutta 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 size-moments 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:

$${\mathrm{\Phi }}_{-}\left(t_n\right)={\int }_0^{\mathrm{K}\Delta t}{n}^{\mathrm{EM}}\left(t=t_n,S\right)\left(\begin{array}{l}1\\ s^{1/2}\\ s\\ s^{3/2}\end{array}\right)\mathrm{ds}$$(37)

• Calculate the negative order moments within the size-interval

$$m_{-\mathrm{a}/2}^{\left[\mathrm{K}\Delta t,1\right]}={\int }_{\mathrm{K}\Delta t}^1{s}^{-\mathrm{a}/2}{n}^{\mathrm{EM}}\left(s\right)\mathrm{ds}$$(38)

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:$$m_{\mathrm{k}/2}^{\left[\mathrm{K}\Delta t,1\right]}=m_{\mathrm{k}/2}\left(t_n\right)-{\mathrm{\Phi }}_{-,\mathrm{k}/2}\left(t_n\right)$$(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:

$$m_{\mathrm{k}/2}^{\left[\mathrm{K}\Delta t,1\right]}=\sum _{i=1}^{n_s}{w}_i{S}_i^{\mathrm{k}/2}$$(40)

where n _s  = 2 + M/2.

• Calculate the new set of evolved moments:

$$m_{\mathrm{k}/2}(t+\Delta t)=\sum _{i=1}^{n_s}{w}_i(S_i-\mathrm{K}\Delta t{)}^{\mathrm{k}/2}, k=-M,\dots ,3$$(41)

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:$$\begin{array}{ccc}{d}_t{m}_0& =& -n\left(S=0\right)\\ d_t{m}_{1/2}& =& -\frac{\mathrm{K}}{2}{m}_{-1/2}\\ d_t{m}_1& =& -\mathrm{K} m_0\\ d_t{m}_{3/2}& =& -\frac{3\mathrm{K}}{2}{m}_{1/2}\\ d_t\left(m_1\vec{u}\right)& =& -\mathrm{K} m_0\vec{u}+m_0\frac{{\vec{u}}_g-\vec{u}}{\theta }\end{array}$$(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 size-velocity, 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:$$\left(m_1\vec{u}{)}^{\left[\mathrm{K\Delta }t,1\right]}\right(t)=\sum _{i=1}^{n_s}{w}_i{S}_i\left(t\right){\vec{c}}_i\left(t\right)$$(43)where .

In the kinetic equation, the evolution of the size and the velocity for each abscissa is independent of the other abscissas.$$\{\begin{array}{lll}\frac{d{\vec{c}}_i}{\mathrm{dt}}& =& -\frac{{\vec{u}}_g-{\vec{c}}_i}{\theta S_i}\\ {\vec{c}}_i\left(t=t_n\right)& =& \vec{u}\left(t_n\right)\end{array}$$(44)and$$\frac{d{S}_i}{\mathrm{dt}}=-\mathrm{K}$$(45)

The solution of the velocity evolution after one time step is:$$\begin{array}{ccc}{S}_i\left(t_{\mathrm{n}+1}\right)& =& S_i\left(t_n\right)-\mathrm{K\Delta }t\\ {\vec{c}}_i\left(t_{n+1}\right)& =& {\vec{u}}_g+\left({\vec{c}}_i\left(t_n\right)-{\vec{u}}_g\right){\left(\frac{S_i\left(t_{n+1}\right)}{S_i\left(t_n\right)}\right)}^{\frac{-1}{\mathrm{K} \theta }}\end{array}$$(46)

Then the momentum and the final cell velocity are computed as follows:$$\begin{array}{ccc}(m_1\vec{u}{)}^{n+1}& =& \sum _{i=1}^{n_s}{w}_k\left(S_i\left(t_{n+1}\right){\vec{c}}_i\left(t_{n+1}\right)\right)\\ \vec{u}\left(t_{n+1}\right)& =& \frac{(m_1\vec{u}{)}^{n+1}}{m_1^{n+1}}\end{array}$$(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:

• Block-based method Figure 1a: the domain is composed of block-structures which are uniformly refined. This method is easier to implement;

  Figure 1. Illustration of AMR techniques.

• Cell-based 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 tree-structure, only the leaves of the tree need to be stored.

In the present paper, we use a Cell-based method, managing the quadrants in a tree-structure. The main advantage of Cell-based methods compared to Block-based methods, is the high rate of compression allowed. Indeed the tree structure used to store the mesh in the Cell-based 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 two-phase flow finite volume solver and the cell-based 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 cell-based 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 sub-domains, each of them being a conforming hexahedral macro-mesh (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 macro-meshes that cover the whole domain, constitute a forest of trees. The possibility of considering multi-trees, where each one covers a sub-domain, enables to represent complex geometries and not only simple square or cubic domain. In each tree, the corresponding sub-domain is defined by a macro-mesh (the tree root) which represents the coarsen mesh and its corresponding refinement level is . The macro-mesh is then recursively refined to multiple non conforming micro-meshes (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 z-order curve, see Fig. 2) to map the leaves of the tree.

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 call-back 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 “processor-decomposition” 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 two-phase model, separated two-phase model and astrophysics applications). These applications involve multi-scale 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.

In the following, we present the finite volume scheme used on the non-conforming 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:$$\begin{array}{ccc}{M}_q^{n+1}& =& M_q^n-\frac{\Delta t}{\Delta x_q}\left(F_{q+1/2}-F_{q-1/2}\right)\\ p_q^{n+1}& =& p_q^n-\frac{\Delta t}{\Delta x_q}\left(G_{q+1/2}-G_{q-1/2}\right)\end{array}$$(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:$$\begin{array}{ccc}{F}_{i+1/2}& =& F\left(\right(M,p{)}_q^n,(M,p{)}_{{\eta }_1^r}^n,(M,\mathrm{p}{)}_{{\eta }_2^r}^n)\\ G_{i+1/2}& =& G\left(\right(M,p{)}_q^n,(M,p{)}_{{\eta }_1^r}^n,(M,p{)}_{{\eta }_2^r}^n)\end{array}$$(49)

The exact expression of the fluxes, derived from the integral expression of the moments and the kinetic solution (29), is:$$\left(\begin{array}{l}{F}_{q+1/2}\\ G_{q+1/2}\end{array}\right)=\left(\begin{array}{l}{F}_{q+1/2}^{+}\\ G_{q+1/2}^{+}\end{array}\right)+\left(\begin{array}{l}{F}_{q+1/2}^{-}\\ G_{q+1/2}^{-}\end{array}\right)$$(50)such that$$\begin{array}{l}\left(\begin{array}{l}{F}_{q+1/2}^{+}\\ G_{q+1/2}^{+}\end{array}\right)=\frac{1}{\Delta y_q\Delta t}{\int }_{\left(x,y\right)\in q}^{}\left(\begin{array}{l}M\left(t_n,x,y\right)\\ p\left(t_n,x,y\right)\end{array}\right)I_{{\mathrm{\Sigma }}_q^{+}\left(\Delta t\right)}\left(x,y\right)\mathrm{dxdy}\end{array}$$(51)and$$\begin{array}{c}\left(\begin{array}{l}{F}_{q+1/2}^{-}\\ G_{q+1/2}^{-}\end{array}\right)=\frac{1}{\Delta y_q\Delta t}\left({\int }_{\left(x,y\right)\in {\eta }_1^r}^{}\left(\begin{array}{l}M\left(t_n,x,y\right)\\ p\left(t_n,x,y\right)\end{array}\right)I_{{\mathrm{\Sigma }}_{{\eta }_1^r}^{-}\left(\Delta t\right)}(x,y)\mathrm{dxdy}+{\int }_{\left(x,y\right)\in {\eta }_2^r}^{}\left(\begin{array}{l}M\left(t_n,x,y\right)\\ p\left(t_n,x,y\right)\end{array}\right)I_{{\mathrm{\Sigma }}_{{\eta }_2^r}^{-}\left(\Delta t\right)}(x,y)\mathrm{dxdy}\right)\end{array}$$(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 x-direction). 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:$$\left(\begin{array}{l}{F}_{q+1/2}\\ G_{q+1/2}\end{array}\right)=\frac{\Delta x_{{\eta }_1^r}}{\Delta x_q}\left(\begin{array}{l}F\left((M,p{)}_q,(M,p{)}_{{\eta }_1^r}\right)\\ G\left((M,p{)}_q,(M,p{)}_{{\eta }_1^r}\right)\end{array}\right)+\frac{\Delta x_{{\eta }_2^r}}{\Delta x_q}\left(\begin{array}{l}F\left((M,p{)}_q,(M,p{)}_{{\eta }_2^r}\right)\\ G\left((M,p{)}_q,(M,p{)}_{{\eta }_2^r}\right)\end{array}\right)$$(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 double-size neighbor on its left and two half-size 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):$$\begin{array}{l}{\tilde{m}}_0\left(x,y\right)=(\begin{array}{ll}{m}_{0,q}+D{m}_0^l\left(x-x_q\right)+D{m}_0^u\left(y-y_q\right),& \mathrm{if} \left(x,y\right)\in \mathrm{chil}{d}_1\\ m_{0,q}+D{m}_0^r\left(x-x_q\right)+D{m}_0^u\left(y-y_q\right),& \mathrm{if} \left(x,y\right)\in \mathrm{chil}{d}_2\\ m_{0,q}+D{m}_0^l\left(x-x_q\right)+D{m}_0^d\left(y-y_q\right),& \mathrm{if} \left(x,y\right)\in \mathrm{chil}{d}_3\\ m_{0,q}+D{m}_0^r\left(x-x_q\right)+D{m}_0^d\left(y-y_q\right),& \mathrm{if} \left(x,y\right)\in \mathrm{chil}{d}_4\\ & \end{array}\\ \end{array}$$(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 z-curve 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:$$\begin{array}{ccc}D{m}_0^l& =& \frac{m_{0,q}-m_{0,{\eta }^l}}{(\Delta x_q+\Delta x_{{\eta }^l})/2}\\ D{m}_0^r& =& \frac{(m_{0,{\eta }_1^r}+m_{0,{\eta }_2^r})/2-m_{0,q}}{(\Delta x_q+\Delta x_{{\eta }_1^r})/2}\end{array}$$(55)

Then we compute the L ¹ -norm difference between and the predicted .$$\mathrm{err}(m_0{)}_q={\int }_q^{}|m_{0,q}-{\tilde{m}}_0(x,y\left)\right|\mathrm{dxdy}$$(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 test-cases:

• the first case consists in a localized spray in the presence of Taylor-Green (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 Taylor-Green 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 non-favorable case since the spray is present in the whole domain.

For the three cases, we fix the CFL (Courant-Friedrichs-Lewy) at CFL = 0.9 and the constant used in the refinement criterion at .

4.1 Droplet Cloud in Taylor-Green Vortices

We simulate a non evaporating polydisperse spray in two dimensional space configuration in the presence of Taylor-Green vortices for the gas: a steady solution of the inviscid incompressible Euler equations. The non-dimensional velocity field of the gas is given by the following expression:$$\begin{array}{ccc}{u}_g(x,y)& =& \sin \left(2\pi x\right)\cos \left(2\pi y\right)\\ v_g(x,y)& =& -\cos \left(2\pi x\right)\sin \left(2\pi y\right)\end{array}$$(57)where .

Initially, we consider a motionless cloud located in the bottom-left 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:$$m_{k/2}(t=0,\vec{x})=\frac{2}{k+2}\left(S_{\mathrm{mid}}^{(k+2)/2}\right)\exp (-||\vec{x}-{\vec{x}}_c|{|}_2^2/r^2)$$(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 delta-shocks. Notice that these singularities are completely handled numerically, thanks to the ‘kinetic’ numerical scheme and the robustness is ensured.

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 .

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 Taylor-Green 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.

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 ¹-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 relative-error 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.

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 Taylor-Green 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.

4.2 Taylor-Green 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.

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 Taylor-Green 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.

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:$$\langle m_0\rangle \left(t\right)={\int }_V^{}{m}_0\left(t,\vec{x}\right)\mathrm{dV}$$(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.

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 three-dimensional 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 three-dimensional Navier-Stokes equations for the gas phase under the low-Mach assumption. The characteristics of this HIT is given in Table 5, where is the turbulent Reynolds number, is the velocity root-mean-square, 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.

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).

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 .

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 .

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 two-phase 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 Williams-Boltzmann 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 two-phase 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 ² 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 multi-scale 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 multi-scale 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 size-velocity 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 two-phase 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

All Tables

All Figures