Konstantinos Feroukas
email: konstantinos.feroukas@etu.univ-amu.fr

Alexandre Chiapolino
email: alexandre.chiapolino@rs2n.eu

Richard Saurel
email: richard.saurel@univ-amu.fr

Simplified interfacial area modeling in polydisperse two-phase flows under explosion situations

Keywords: two-phase flows, explosions, interfacial area, polydisperse particles, radius distribution laws

The aim of the present work is to account for polydisperse effects in a two-phase flow with a simple and fast method. Polydisperse two-phase flows arise in numerous applications. Fire sprinkler systems are relevant examples as they release clouds of polydisperse droplets. Another relevant example is the polydisperse two-phase flow created by the detonation of an explosive charge surrounded by a liquid layer. In such a situation, material interfaces are initially present and the created two-phase flow consists of a carrier gas phase and a liquid phase involving many droplets of various sizes. Spherical particles or droplets are usually assumed in two-phase flow computations. When dealing with explosion situations involving both dense and dilute flow regimes, multiple particle diameters can be addressed but at the price of introducing as many additional equations that describe mass, momentum and energy balance of the various particle classes. Consequently, the computation time needed to address numerical resolution increases tremendously. Under explosion situations involving many particle diameters, the method becomes intractable and is usually reduced to a single diameter, which is often insufficient. A simplified approach is developed in the present work to account for a substantial number of particles of different sizes with few extra computational cost. The approach is said simplified as a single velocity and a single temperature are considered for all the spherical particles, regardless of their diameters. This type of modeling seems adapted for the target explosion situations. The focus is placed on the interfacial area, which is the main parameter involved in the coupling of the two phases. In the present work, Gamma-like continuous probability distributions are considered to address the various sizes of particles. The effects of the size distribution are only summarized in the specific interfacial area, yielding consequently few code modifications while taking into account the polydisperse aspect of the two-phase flow.

Introduction

Many two-phase flows involve a very large number of liquid or solid particles per unit volume (typically 10 12 particles • m -3 ) having a significant impact on the behavior of the flow, as a result of interactions. Such two-phase flows appear in various contexts, like the transport of solid particles by air [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF]), the collection of ice on buildings and aircraft structures (Lewis and Brun, 1956 [2], Gelder et al., 1956 [3]), or fluidized beds and other two-phase flow phenomena of interest in chemical processing (Torobin and Gauvin, 1959 [4]). Other relevant examples are fire sprinkler systems that release clouds of polydisperse droplets, fluid-droplet sprays that are of particular interest for applications involving combustion (Williams, 1958 [5]), and fuel suspensions resulting from explosion situations. This last case involves combustion as well and is of particular interest in the present contribution.

Whether in liquid sprays or when a liquid layer is set in motion by the detonation of an explosive charge, many liquid droplets of various sizes are created and consequently involve a large exchange surface with the gas phase. Determination of the interfacial area is a key problem in combustion and two-phase flow modeling (Drew and Passman, 2006 [6]). An illustrative example of a liquid layer set in motion by an explosive charge is presented in Fig. 1.

Ambient air

Liquid water Dense gas at elevated initial pressure Figure 1: A cylindrical explosive charge is initially surrounded by a liquid layer. When the charge explodes the liquid layer transforms to a cloud of droplets forming highly dynamical particle jets. Experimental results are presented on the left. Same jetting effects appear when the liquid is replaced by a granular layer. These jets are present in cylindrical and spherical dispersal explosions. On the right, a schematic representation of the initial cylindrical gas-liquid explosive system is depicted. The internal cylinder is initially filled with a dense gas at high pressure. The external cylinder is initially filled with liquid water at atmospheric pressure. Atmospheric air surrounds both cylinders. Material interfaces are initially present.

The liquid phase is said to be polydisperse as it contains a substantial number of droplets of different sizes, making major effects on the two-phase flow. Various methods are available to address multi-dimensional computations of polydisperse flows, see for instance [START_REF] Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporation sprays[END_REF], [START_REF] Fréret | Eulerian models and three-dimensional numerical simulation of polydisperse sprays[END_REF], [START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF], [START_REF] Kah | Taking into account polydispersity for the modeling of liquid fuel injection in internal combustion engines[END_REF], [START_REF] Fréret | Direct numerical simulation of polydisperse evaporating sprays in 3D jet configuration using Euler-Euler and Euler-Lagrange formalisms[END_REF], [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], [START_REF] Essadki | Adaptive mesh refinement and high order geometrical moment method for the simulation of polydisperse evaporating sprays[END_REF] in the context of polydisperse sprays where the droplets are diluted in a carrier phase. When dealing with explosion situations, such as the one of Fig. 1, the liquid layer transforms to a cloud of polydisperse droplets as well. However, material interfaces are initially present and bring out additional numerical difficulties.

One way to account for multiple sizes of particles, or droplets in the present explosion context where both dense and dilute flow regimes occur, is to address each class of particles with its own set of balance equations (mass, momentum, total energy, and specific number of particles). In the present work, a class of particles represents a population of spherical particles having a distinct size (radius). Each class of particles can be described by its own radius, velocity, temperature, specific number and specific interfacial area that can be easily determined, provided that the particles are supposed to be spherical. This approach appears ideal and appealing but becomes intractable when dealing with real explosion situations involving a very large number of classes of particles and consequently an even larger number of additional equations, tremendously increasing the computation time needed to address numerical simulations. In such a case, the method is usually reduced to a single or few classes of particles, which is often insufficient.

The present paper attempts to reduce substantially the simulation time while increasing the accuracy of the solution. Multiple approaches have been developed to address polydisperse effects without the burden of introducing many additional balance equations. For example, [START_REF] Fan | Application of the direct quadrature method of moments to polydisperse gassolid fluidized beds[END_REF] [START_REF] Fan | Application of the direct quadrature method of moments to polydisperse gassolid fluidized beds[END_REF] consider a population balance equation (PBE), coupled to the continuity and momentum balance equations with the help of the direct quadrature method of moments (DQMOM). This method seems to be the most popular to deal with polydisperse effects. However it is restricted to dilute flows, where the volume fraction of the dispersed phase as well as related terms are neglected. The flow model is in the same sense an extended version of Marble's model (1963) [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF]. This method then appears unable to account for initial material interfaces, like the initial liquid-gas interface depicted in Fig. 1. In such a situation, a material interface is present and the flow ranges from dense to dilute concentration of particles. The DQMOM method is consequently unsuitable to address the present target applications involving initial material interfaces. However, it appears well-suited to deal with spray flows, [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: Direct quadrature method of moments and multi-fluid method[END_REF] [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: Direct quadrature method of moments and multi-fluid method[END_REF].

Baer-and-Nunziato-type (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] two-phase flow models are then preferred. Various variants are available, such as [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF]Saurel et al.'s ( , 2017) ) [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF][START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF]. In the present work, unlike the DQMOM method, the liquid phase is described by a single velocity and a single temperature but multiple classes of particles. This simplification has been used for instance in [START_REF] Olmos | Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up[END_REF] [START_REF] Olmos | Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up[END_REF] in the context of bubble-column reactors, where the different classes of particles are convected with the same mean algebraic velocity. It may appear contradictory at first glance but implies that the essential of the polydisperse effects is addressed through the specific interfacial area which accounts for a whole spectrum of radius distributions. In the present explosion context, this assumption appears appropriate as timescales related to velocity relaxation are small. Moreover, such a modeling appears realistic when the liquid phase reaches the saturation temperature which is independent of the sizes of the droplets. The polydisperse character of the liquid phase is only summarized in the specific interfacial area, and only two phases are needed, i.e. the gas phase and the liquid phase that is made of polydisperse droplets. Consequently, only two sets of balance equations are required and the computation time requested to perform a simulation is comparable to the one needed by a conventional computation, considering a unique droplet size in a control volume.

A droplet size distribution is needed nonetheless. Many experimental studies aim to describe the size distribution in various applications. Recent works include for instance [START_REF] Chandrakar | Aerosol indirect effect from turbulence-induced broadening of cloud-droplet size distributions[END_REF] [START_REF] Chandrakar | Aerosol indirect effect from turbulence-induced broadening of cloud-droplet size distributions[END_REF] where the influence of aerosol concentration on a cloud-droplet size distribution is investigated in a laboratory chamber, or Rousseau et al. (2021) [START_REF] Rousseau | Droplet spatial distribution in a spray under evaporating and reacting conditions[END_REF] where spray combustion is studied with the help of an experimental test rig. Some theoretical studies attempt to reproduce favorably experimental data. For instance, [START_REF] Carrica | A polydisperse model for bubbly two-phase flow around a surface ship[END_REF] [START_REF] Carrica | A polydisperse model for bubbly two-phase flow around a surface ship[END_REF] used a statistical description of two-phase flows based on the Boltzmann theory of dispersed gases, and described a bubble distribution function with the help of the bubble mass, position and time. [START_REF] Li | Droplet size distribution in sprays based on maximization of entropy generation[END_REF] [START_REF] Li | Droplet size distribution in sprays based on maximization of entropy generation[END_REF] proposed a droplet size distribution model based on the concept of the maximization of entropy generation during the liquid atomization process. Zhang et al. (2019) [START_REF] Zhang | Size distribution theory for jumping-droplet condensation[END_REF] developed a theoretical framework based on the contact and the coalescence of droplets.

In the present paper, the droplet size distribution is based on a continuous probability distribution. Many mathematical functions are available and used to address particle size distribution (Yoon, 2005 [25], Igel and Van Den Heever, 2017 [START_REF] Igel | The importance of the shape of cloud droplet size distributions in shallow cumulus clouds. Part II: Bulk microphysics simulations[END_REF], Urbán and Józsa, 2018 [START_REF] Urbán | Investigation of fuel atomization with density functions[END_REF], Hareli et al., 2021 [28] to cite a few). The commonly used size distribution functions in fluid dynamics include the Normal, Log-Normal, Nukiyama-Tanasawa, Rosin-Rammler, Beta, modified Beta and Gamma-type distributions. Discussions about these size distribution functions can be found in Ge (2006) [START_REF] Ge | Probability density function modeling of turbulent non-reactive and reactive spray flows[END_REF] for instance. The proposed method is presented with the help of Gamma-type distributions but may be used with various functions.

The paper is organized as follows. Section 2 presents the two-phase model. Viscous drag interaction effects, having a major impact on the two-phase flow, are introduced in Section 3. Section 4 presents the proposed method and attempts to address, through a simple size distribution function, the polydisperse effects of the liquid phase. Computational examples are provided in Section 5.

Flow model

The model of [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] is based upon a mixture evolving in total disequilibrium. It is based upon the inviscid Euler equations for each pure phase. The balance equations for phases 1 and 2 are,

                                                       ∂α 1 ∂t + u I ∂α 1 ∂x = µ(p 1 -p 2 ), ∂ (α 1 ρ 1 ) ∂t + ∂ (α 1 ρ 1 u 1 ) ∂x = 0, ∂ (α 1 ρ 1 u 1 ) ∂t + ∂ α 1 ρ 1 u 2 1 + p 1 ∂x = p I ∂α 1 ∂x + λ (u 2 -u 1 ) , ∂ (α 1 ρ 1 E 1 ) ∂t + ∂ (α 1 (ρ 1 E 1 + p 1 ) u 1 ) ∂x = p I u I ∂α 1 ∂x -µp ′ I (p 1 -p 2 ) + λu ′ I (u 2 -u 1 ) , ∂ (α 2 ρ 2 ) ∂t + ∂ (α 2 ρ 2 u 2 ) ∂x = 0, ∂ (α 2 ρ 2 u 2 ) ∂t + ∂ α 2 ρ 2 u 2 2 + p 2 ∂x = p I ∂α 2 ∂x -λ (u 2 -u 1 ) , ∂ (α 2 ρ 2 E 2 ) ∂t + ∂ (α 2 (ρ 2 E 2 + p 2 ) u 2 ) ∂x = p I u I ∂α 2 ∂x + µp ′ I (p 1 -p 2 ) -λu ′ I (u 2 -u 1 ) , ∂N 2 ∂t + ∂ (N 2 u 2 ) ∂x = 0.
(2.1)

The notations are conventional in the two-phase flow literature. A frame of reference (x) is chosen and the time variable is denoted by t. α k , ρ k , p k , E k = e k + 1 2 u 2 k denote respectively the volume fraction, density, pressure and total energy of phase k. e k represents the internal energy and u k represents the center of mass velocity of phase k. N 2 represents the specific number of particles, i.e. the total number of particles per unit volume. In the rest of the paper, the carrier gas phase will then be indexed by 1 and the liquid phase by 2. The mixture internal energy is defined as e = Y k e k where Y k = (α k ρ k )/ρ denotes the mass fraction of phase k. The mixture density and pressure are defined as ρ = α k ρ k and p = α k p k . Equation system (2.1) is a two-phase model for mixture flows evolving in pressure, velocity and temperature disequilibria. The choice of interfacial average velocities u I and pressures p I was originally expressed with the relations u I = u 2 and p I = p 1 , the symmetric choice u I = u 1 and p I = p 2 , being possible as well. More general and symmetric estimates have been proposed by [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF] [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF],

       u I = u ′ I + sgn ∂α 1 ∂x p 2 -p 1 Z 1 + Z 2 , with u ′ I = Z 1 u 1 + Z 2 u 2 Z 1 + Z 2 , p I = p ′ I + sgn ∂α 1 ∂x (u 2 -u 1 ) Z 1 Z 2 Z 1 + Z 2 with p ′ I = Z 1 p 2 + Z 2 p 1 Z 1 + Z 2 , (2.2)
where Z k = ρ k c k is the acoustic impedance and c k is the speed of sound of fluid k. This latter is provided by a convex equation of state for each phase. The analysis which has led to these estimates is based upon a homogenization method developed by [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF]. Equations (2.1) have been extended to 3D in Franquet and Perrier (2012) [START_REF] Franquet | Runge-Kutta discontinuous Galerkin method for the approximation of Baer and Nunziato type multiphase models[END_REF].

The first equation of (2.1) is non-conservative and represents the transport of the first volume fraction α 1 at interfacial velocity u I . During the advection stage, volume variations caused by pressure differences between the phases appear through the relaxation term µ(p 1 -p 2 ), where µ controls the rate at which pressure equilibrium is reached. The above-mentioned analysis provided this coefficient as well µ = AI Z1+Z2 where A I represents the specific interfacial area of the mixture. For instance, when dealing with a cloud of liquid droplets of a single radius R 2 , the specific interfacial area is A I = 3α2 R2 . The volume variations of the phases are then directly proportional to the pressure difference between the phases and the speed at which the equilibrium is reached is controlled by the µ coefficient. This latter depends only upon the acoustic impedance of the phases and upon the specific interfacial area A I . The second and fifth equations of (2.1) describe mass balance of the corresponding phase while the third and sixth equations are related to their momentum balance. These last two relations are non-conservative. The velocity relaxation terms on the right-hand side of the momentum equations read ±λ (u 2 -u 1 ), where λ is the product of the specific interfacial area with the drag coefficient. λ is a positive function (or tensor if there are more than two fluids). It is involved in the viscous drag force between the two phases and controls the rate at which velocities tend towards equilibrium. The non-conservative term p I ∂α1 ∂x represents the pressure force acting at the liquid droplet cloud boundaries with p I denoting the interfacial pressure given by Eq. (2.2). This non-conservative term represents a "differential drag force" as its amplitude is high in zones of high volume gradients and vanishes when the volume fraction is uniform. It has been shown in [START_REF] Chiapolino | Numerical investigations of two-phase finger-like instabilities[END_REF] [START_REF] Chiapolino | Numerical investigations of two-phase finger-like instabilities[END_REF] that this term is of main importance in the formation of particle jets in the explosion situation depicted in Fig. 1. The fourth and seventh equations of (2.1) describe the energy balance of phase k. Those latter ones are non-conservative as well due to the presence of the term p I u I ∂α k ∂x and the relaxation terms on the right-hand side. Finally, the last equation describes the conservation of the specific number N 2 of liquid particles (droplets). This equation is conservative as fragmentation effects are not considered in the present work. Thermal exchange effects and mass transfer have been omitted as well for the sake simplicity.

Equation system (2.1) considers mixtures in pressure, velocity and temperature disequilibria. Its extension to more than two phases is possible [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF]. It is able to deal with material interfaces encountered for instance at the early times of the explosion situation depicted in Fig. 1, as well as two-phase suspensions occurring at later times. For the target explosion situations, the two-phase equation system (2.1) of [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF] [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF], variant of Baer and Nunziato's (BN) model (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], is consequently preferred over the model of [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] [START_REF] Marble | Dynamics of a gas containing small solid particles[END_REF] or the DQMOM formulation [START_REF] Fan | Application of the direct quadrature method of moments to polydisperse gassolid fluidized beds[END_REF], [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: Direct quadrature method of moments and multi-fluid method[END_REF], [START_REF] Marchisio | Computational models for polydisperse particulate and multiphase systems[END_REF], suitable only in the limit of disperse flows. Such formulations are indeed suitable for dilute spray flows but cannot address resolved-interface situations where material interfaces separate two pure or nearly pure media.

The two-phase equation system (2.1) ensures the satisfaction of the interface conditions through the non-conservative terms and the associated interfacial variables (Saurel and Pantano, 2018 [35]). Equation system (2.1) is indeed able to fulfill the expected interface condition of mechanical equilibrium (continuity of pressures and normal velocities) in the two limits, λ, µ = 0 and λ, µ → ∞. In the first option, the interface conditions are ensured by the non-conservative terms as the interfacial variables in Eq. (2.2) model contact interface conditions as general solutions of local Riemann problems. This method has been used for example by Layes and Le [START_REF] Layes | Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[END_REF] [START_REF] Layes | Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[END_REF] to study shock interaction with a gas bubble. The two-phase formulation is also able to deal with permeable interfaces (boundaries of bubbles or droplets clouds). For instance, it has been used to address permeable granular interfaces by Saurel et al. (2014) [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF], who extended (2.1) to account for granular effects. The second option relies on stiff mechanical relaxation. In this limit the mixture evolves with a single pressure and a single velocity. This approach was proposed by [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] in a splitting formulation where the hyperbolic part of (2.1) is solved during a time step in the absence of source terms, followed by pressure and velocity relaxation steps with sources µ(p 1 -p 2 ) and λ (u 2 -u 1 ) respectively, with both λ, µ → ∞.

Indeed, relaxation phenomena can be added depending upon the flow condition of the multiphase medium and may yield total or partial equilibrium depending upon the rate at which the corresponding equilibrium is supposed to be reached. For example, instantaneous pressure equilibrium may be found whereas velocities remain in disequilibrium. These circumstances are typical of explosion situations, as typical timescales associated with the pressure equilibrating process are small [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF], [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF], and are of particular interest of the present work. Note that such an instantaneous pressure relaxation is not equivalent to strict pressure equilibrium that yields non-hyperbolic or conditionally hyperbolic models. Details of stiff relaxation solvers can be found in Le [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF] [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF] for instance.

Under the form (2.1), the formulation is restricted to two phases, i.e. a carrier gas phase (indexed 1) and a liquid phase (indexed 2) made of a single class of N 2 droplets of radius R 2 . Equation system (2.1) can be extended to account for multiple classes of droplets with the help of as many additional sets of balance equations. In such a case, each class of droplet is described by its own radius, velocity, temperature, specific number and specific interfacial area. However, the present paper attempts to account for the polydisperse character of the liquid phase with a simplified method. As only flow situations involving pressure equilibrium between the gas and liquid phases are of interest in the present work, the polydisperse aspect of the liquid droplets impacts only the viscous drag force between the two phases through the specific interfacial area A I . The viscous drag effects are addressed hereafter.

Viscous drag

For the sake of clarity, computation of the viscous drag force is presented in the context of a single class of particles. In a control volume, the liquid phase is then made of N 2 droplets of a single radius R 2 . Polydisperse effects will be accounted for later through the specific interfacial area A I . The viscous drag parameter λ present in (2.1) controls the rate at which velocity equilibrium is reached between the gas and liquid phases. The λ (u 2 -u 1 ) term represents the viscous drag force and λu ′ I (u 2 -u 1 ) the power of this force (per unit volume). In the following, this force is denoted as N 2 F 1→2 where N 2 represents the specific number of droplets. For the sake of simplicity the droplets are considered spherical in this work and viscous drag effects are treated via the following Stokes relation,

F 1→2 = 6π µ 1 R 2 (u 1 -u 2 ) , (3.1)
where R 2 is the radius of the droplets and µ 1 the dynamic viscosity of the carrier phase. The particle Reynolds number is now introduced,

Re 2 = 2R 2 ρ 1 |u 2 -u 1 | µ 1 . (3.2)
It is important to note that such viscous drag representation is only valid for low Reynolds numbers. In such conditions the viscous drag coefficient reads C d = 24

Re2 . With the help of the previous relations, the viscous drag force can be written concisely as,

F 1→2 = C d R 2 2 π ρ 1 2 |u 2 -u 1 | (u 1 -u 2 ) . (3.3)
However, in order to extend the present viscous drag law to higher Reynolds numbers, the viscous drag coefficient is reconsidered to account for turbulent effects following Naumann and Schiller (1935) [START_REF] Naumann | A drag coefficient correlation[END_REF], (3.4)

C d =    24
As the droplets are considered spherical with a radius R 2 , the specific number of droplets reads,

N 2 = α 2 4 3 π R 3 2 , (3.5)
and the total viscous drag force in a control volume becomes,

N 2 F 1→2 = 3 8R 2 α 2 C d ρ 1 |u 2 -u 1 | (u 1 -u 2 ) , (3.6)
that is to say,

N 2 F 1→2 = λ (u 1 -u 2 ) with λ = 3 8R 2 α 2 C d ρ 1 |u 2 -u 1 |. (3.7)
It is however more convenient to express Eq. (3.7) in terms of specific interfacial area A I . The droplets being spherical, the specific interfacial area reads,

A I = 4πR 2 2 N 2 . (3.8)
The combination of Eqs. (3.8) and (3.5) yields,

A I = 3α 2 R 2 . ( 3.9)
The viscous drag parameter λ is consequently expressed as,

λ = 1 8 C d ρ 1 |u 2 -u 1 | × A I . (3.10)
Relation (3.10) is used to compute the viscous drag effects between the carrier gas phase 1 and the liquid phase 2. Turbulent effects are summarized through the C d coefficient, present in the first term of Eq. (3.10), and is determined with the help of Eq. (3.4). It depends on the particle Reynolds number (3.2) which itself depends on the radius R 2 of the droplets. This latter is determined by Eq. (3.5), as the number of droplets N 2 as well as the volume fraction of the liquid phase α 2 are known from the balance equations of (2.1). The size of the droplets is addressed through the specific interfacial area A I , present as the second term of Eq. (3.10), and determined via Eq. (3.9).

In the following section, multiple sizes of droplets are addressed via a simplified method that accounts for a whole spectrum of particle radii through the specific interfacial area A I , that is reconsidered. Computation of A I relies on a continuous probability distribution and yields only few code modifications. Depending on the context multiple distribution functions can be used as long as their distribution moments are available, as will be seen hereafter.

Polydisperse particles

Computation of the viscous drag force, with polydisperse droplets, is based on the previous monodisperse relations. However the specific interfacial area A I is adjusted to account for the multiple droplet sizes. The drag force is computed via Eq. (3.7) with the help of the λ coefficient determined by Eq. (3.10). This later depends on the specific interfacial area A I that will now be reconsidered.

General relations

Two additional variables are introduced to deal with the multiple classes of particles. Those are α 2,k and N 2,k that denote respectively the volume fraction and the specific number (per unit volume) of particles of the k th class. Recall that index 2 denotes the liquid phase composed of spherical particles, or droplets in the present context.

Liquid volume fraction and number of particles

The particles being spherical, the following relation appears,

α 2,k = 4 3 πR 3 2,k N 2,k . (4.1)
Moreover, the volume fraction of the liquid phase α 2 is defined as the sum of the volume fractions α 2,k constituting the liquid phase,

α 2 = N classes k=1 α 2,k , (4.2)
where N classes is the number of classes of particles. In the present context, the volume fractions α 2,k of the various classes of droplets depend on the various radii R 2,k . As the liquid volume fraction α 2 is the sum of the α 2,k volume fractions, the following relation arises:

α 2 = ∞ 0 α 2,k (R 2 )dR 2 . (4.3)
The volume fraction α 2,k , describing the volume occupied by the class k of particles in the liquid phase [START_REF] Lewis | Impingment of Water Droplets on a Rectangular Half Body in a Two-Dimensional Incompressible Flow Field[END_REF], is defined by Eq. (4.1) involving the specific number N 2,k of particles of the k th class. In order to account for multiple sizes of particles, a continuous probability distribution is addressed in the following. To this end, the probability

P (R 2 ∈ [a, b]
) is introduced and corresponds to the probability for the radius of the particles to belong in the interval [a, b]. This probability reads,

P (R 2 ∈ [a, b]) = b a f (R 2 ) dR 2 , (4.4)
where f is the probability density function (PDF). It is important to emphasize that the PDF specifies the probability of radius R 2 falling within a specific range of values, as opposed to taking on any one value. With the help of Eq. (4.4) the specific number of particles having a radius in the interval [a, b] is expressed as,

N (R 2 ∈ [a, b]) = N 2 P (R 2 ∈ [a, b]) = N 2 b a f (R 2 ) dR 2 , (4.5)
with N 2 the specific total number of liquid particles, all sizes included. Note that a probability density function satisfies the properties:

f (R 2 ) ≥ 0, ∀R 2 and +∞ -∞ f (R 2 ) dR 2 = 1. (4.6)
However, in the present two-phase context, the interval is reduced to R 2 ∈ [0, +∞[. As a result, the total number of particles N 2 is recovered in the [0, +∞[ interval,

N (R 2 ∈ [0, +∞[) = N 2 +∞ 0 f (R 2 ) dR 2 = N 2 . (4.7)
The specific number N 2,k of particles of the k th class, present in Eq. (4.1), is then to be expressed with a range of values for R 2 . Consequently, an infinitesimal interval [R, R + dR] is considered. Equation (4.1) is then used under the form,

α 2,k (R 2 ) = 4 3 πR 3 2 N (R 2 ∈ [R, R + dR]) . (4.8) Equation (4.8) is now introduced in (4.3), α 2 = ∞ 0 4 3 πR 3 2 N (R 2 ∈ [R, R + dR]) dR 2 . (4.9)
Then, Expression (4.5) is embedded in this last relation and yields,

α 2 = ∞ 0 4 3 πR 3 2 N 2 R+dR R f (R 2 ) dR 2 dR 2 . (4.10)
The second integral term is now analyzed,

R+dR R f (R 2 ) dR 2 = F (R + dR) -F (R) = dF (R 2 ) = f (R 2 ) , (4.11)
where F denotes a primitive of the probability density function f (R 2 ) and dF (R 2 ) its derivative, i.e. the probability density function f (R 2 ). Introducing Eq. (4.11) into (4.10), the following relation appears,

α 2 = ∞ 0 4 3 πR 3 2 N 2 f (R 2 ) dR 2 , (4.12)
that is to say,

α 2 = 4 3 πN 2 ∞ 0 R 3 2 f (R 2 ) dR 2 . (4.13)
Recall that N 2 is provided by the last equation of (2.1).

Interfacial area

The same reasoning is now repeated for the interfacial area A I , representing the specific exchange surface between the gas phase and the liquid phase containing many droplets of various sizes. The specific interfacial area A I is defined as the sum of the specific interfacial areas A I,k of the various droplet classes. The droplets being spherical, A I is defined as,

A I = N classes k=1 A I,k = N classes k=1 4πR 2 2,k N 2,k . (4.14)
As a continuous distribution of particles is considered, the previous relation becomes,

A I = ∞ 0 A 2,k (R 2 ) dR 2 = ∞ 0 4πR 2 2 N (R 2 ∈ [R, R + dR]) dR 2 = 4πN 2 ∞ 0 R 2 2 R+dR R f (R 2 ) dR 2 dR 2 , (4.15)
where Eq. (4.5) has been introduced with an infinitesimal interval [R, R + dR] for the above-mentioned reason. With the help of Eq. (4.11), this last relation becomes,

A I = 4πN 2 ∞ 0 R 2 2 f (R 2 ) dR 2 . (4.16)

Moments of the probability distribution

Analyzing Eq. (4.7) for the total number of particles N 2 , Eq. (4.13) for the volume fraction α 2 , and Eq. (4.16) for the specific interfacial area A I , it appears that all relations involve an integral term, related to the PDF, that can be written under the generic form,

m n = ∞ 0 R n 2 f (R 2 ) dR 2 , (4.17)
with n = 0, 2, 3. Equation (4.17) consists of the formulation of the distribution moments, which are well-known for most PDFs. The last two moments (n = 2 and n = 3) represent respectively the variance and the skewness of the continuous probability distribution. Moreover, the first moment m 1 is defined as the mean value of the distribution.

In the present context, it consists of the mean radius R 2 of the particles constituting the liquid phase. The various relations can then be expressed under the following form,

    N 2 R 2 A I α 2     =     N 2 1 4πN 2 4 3 πN 2         m 0 m 1 m 2 m 3     (4.18)
and may be used with various PDFs as long as their m 0 , m 1 , m 2 and m 3 moments are available. The choice of the PDF, in accordance with the present gas-liquid context, is discussed in the following section.

Distribution law

As mentioned in the Introduction, the commonly used size distribution functions in fluid dynamics include the Normal, Log-Normal, Nukiyama-Tanasawa [START_REF] Nukiyama | Experiments on the atomization of liquids in an air stream, report 3, on the dropletsize distribution in a atomized jet[END_REF], Rosin-Rammler [START_REF] Rosin | Laws governing the fineness of powdered coal[END_REF], Beta, modified Beta and Gamma-type distributions. The present work attempts to account for clouds of droplets encountered in explosion situations where the size distribution displays an asymmetric bell curve. Asymmetric continuous probability distributions are then considered. The size distribution laws are meant to be simple and made of as few adjustable parameters as possible. Moreover, the droplet size distribution is to be defined in the R 2 ∈ [0, +∞[ interval and is to satisfy the conditions, lim

R2→Rmin (f (R 2 )) = 0, lim R2→Rmax (f (R 2 )) = 0, (4.19)
where R min and R max are finite values representing respectively the minimum radius and the maximum radius of the particles.

The Log-Normal, Rosin-Rammler, modified Beta and Gamma-type distributions are well-suited in the present context, as they display an asymmetric behavior, involve only two parameters and satisfy the conditions (4.19). However, the Normal, Nukiyama-Tanasawa and Beta distributions appear unsuitable for the target application. Indeed, the Normal distribution presents a symmetric behavior, the Nukiyama-Tanasawa distribution involves four parameters, involving consequently a lot of adjustments, and the Beta distribution is defined only in the [0, 1] interval, which appears more suitable to describe mass fractions of multi-species flows.

The primary focus of the present work is to address the polydisperse aspect of a gas-liquid flow with a simple and fast method. The Gamma and Inverse-Gamma distributions are considered in the following. However, the method is not restricted to these functions and can be extended to the aforementioned distributions (the Log-Normal and Rosin-Rammler distributions are depicted in Appendix A). The Gamma distributions consist of two-parameter families of continuous probability distributions. Their probability density functions read:

f Gamma (R 2 ) = β κ Γ(κ) R κ-1 2 e -βR2 , (4.20)
f Inverse Gamma (R 2 ) = β κ Γ(κ) 1 R 2 κ+1 e -β R 2 , (4.21)
where β and κ are real positive coefficients used to adjust the desired particle size distribution depending on the studied situation, and Γ(κ) is the Gamma function [START_REF] Gradshteyn | Table of integrals, series, and products[END_REF] defined for complex numbers with a positive real part. The Gamma function is defined via a convergent improper integral,

Γ (κ) = ∞ 0 ξ κ-1 e -ξ dξ, (4.22)
with ξ the integration variable. This integral is a non elementary function commonly used as an extension of the factorial function to complex numbers. Note that computation of the Gamma function Γ(κ) is already included in some computer languages as an intrinsic function. It is the case with the Fortran language that is used in the present work. The Gamma function interpolates the factorial function (4.22). Note that for positive integer values of κ, the Gamma function simplifies to Γ (κ) = (κ -1)!. The combination of the Gamma (4.20) and Inverse-Gamma (4.21) PDFs with Eq. (4.17) yields, after some algebraic manipulations, the moments of the probability distributions,

m n =        Γ(κ + n) Γ(κ) 1 β n , Gamma, Γ(κ -n) Γ(κ) β n , Inverse Gamma. (4.23)
As the Gamma function Γ (κ) is defined for all complex numbers except the non-positive integers, the following conditions arise: κ Gamma + n > 0 and κ Inverse Gamma -n > 0. The κ coefficient is necessarily positive. Besides, according to (4.18), the moments are only of interest up to order 3 (n = 0, n = 1, n = 2, n = 3). Consequently, the previous restrictions become κ Gamma > 0 and κ Inverse Gamma > 3. The zero moment is trivial and reads: m 0 = 1 for both Gamma and Inverse Gamma PDF. Moreover, an interesting property of the Gamma function is Γ(κ + 1) = κΓ (κ) , ∀κ ∈ R * + . As a result, Eqs. (4.23) can be written as follows for n ∈ N + ,

m n≥1 =          l=n l=1 (k + (l -1)) β n , Gamma, β n l=n l=1 (κ -l) , Inverse Gamma. (4.24)
The moments of the Gamma and Inverse-Gamma probability distributions are then expressed only in terms of the β and κ coefficients. In the following, for the sake of convenience, only κ is considered as a free choice, and becomes a constant parameter. The β coefficient is computed according to the two-phase flow variables. After some algebraic manipulations, the combination of the fourth relation of (4.18) and Eq. (4.24) leads to,

β =        3 4πN 2 3α 2 (κ + 2)(κ + 1)κ, Gamma, 3 3α 2 4πN 2 (κ -3)(κ -2)(κ -1), Inverse Gamma. (4.25)
Note that Eq. (4.25), related to the Gamma distribution, demands α 2 > 0 for the β coefficient to be defined. Equation (4.25) becomes irrelevant in the absence of liquid phase. However, for numerical reasons, the two-phase model (2.1) considers ǫ ≤ α 2 ≤ 1 -ǫ (with ǫ of the order of 10 -6 ), see Saurel and Pantano (2018) [START_REF] Saurel | Diffuse Interfaces and Capturing Methods in Compressible Two-Phase Flows[END_REF] in the context of the diffuse interface method (DIM). Note also that Eq. (4.25), related to the Inverse-Gamma distribution, demands N 2 > 0 for the β coefficient to be defined. Yet, as α 2 ≥ ǫ, traces of liquid are present and β remains unambiguously defined. With the help of Relation (4.25), the specific interfacial A I , summarizing the polydisperse effects, is expressed with a single parameter related to the PDF, i.e. the κ parameter. The specific interfacial area is computed via the third relation of (4.18) and the second moment of the probability distribution defined by Eq. (4.17). In the present paper, the moments result from the Gamma and Inverse-Gamma distributions and are provided by Eq. (4.24). The specific interfacial area consequently reads,

A I =          4πN 2 (κ + 1)κ 4πN 2 3α 2 (κ + 2)(κ + 1)κ -2 3 , Gamma, 4πN 2 (κ -2)(κ -1) 3α 2 4πN 2 (κ -3)(κ -2)(κ -1) 2 3
, Inverse Gamma,

(4.26)
where κ is a constant parameter. α 2 = 1 -α 1 and N 2 are known from the corresponding balance equations of the twophase equation system (2.1). Nevertheless, the initial value of N 2 requires specific attention. This point is addressed hereafter.

Determination of the initial conditions

The initial flow composition is usually considered as initial data. The volume fraction α t=0 2 is then initially known. For a conventional computation, dealing with monodisperse particles, the initial radius of the droplets R t=0 2 is also usually given as an input data. The initial specific number of droplets N t=0 2 is consequently determined via Eq. (3.5). Nevertheless, in the present polydisperse case, the initial specific number N t=0 2 has to satisfy the fourth relation of (4.18), that includes the third moment of the probability distribution and consequently the κ parameter and β coefficient in the context of Gamma and Inverse-Gamma PDSs (Eq. (4.24)). In the present work, the κ parameter is known as an input data and remains constant. Besides, it appears more convenient to consider an initial mean radius R 2 t=0 , instead of the initial β t=0 coefficient that may be difficult to apprehend.

The mean radius R 2 t=0 may be known either as an input data or with the help of a given relation depending on the situation to study. For instance, R 2 t=0 may be computed with the help of the single initial radius R t=0 2 considered by the monodisperse computation, in order to equate the initial specific interfacial area A t=0 I of the polydisperse and monodiperse computations and consequently study how the polydisperse solution departs from the simplified monodisperse one. This is the case for some numerical results provided in Section 5. In that context, after some algebraic manipulations detailed in Appendix C, the two initial mean radii are linked through the relation,

       R 2 t=0 = R t=0 2 κ (κ + 2) , Gamma, R 2 t=0 = R t=0 2 (κ -3) (κ -1) , Inverse Gamma, (4.27)
where the initial monodiperse radius R t=0 2 is given as an input data. As the mean radius consists of the first moment of the probability distribution (Eq. (4.17)), it reads in the present context (Eq. (4.24)):

R 2 t=0 =      κ β t=0
, Gamma,

β t=0 (κ -1)
, Inverse Gamma.

(4.28)

The initial β t=0 coefficient is then computed via Eq. (4.28), with the help of the constant κ parameter and the initial mean radius R 2 t=0 . The initial number N t=0 2 of particles is afterwards computed with the help of the fourth relation of (4.18), including the third moment of the probability distribution. In the present context, the initial number N t=0 2 reads,

N t=0 2 =        3α 2 4π β 3 t=0 (κ + 2)(κ + 1)κ , Gamma, 3α 2 4π (κ -3)(κ -2)(κ -1) β 3 t=0 , Inverse Gamma.

Illustration of the initial Gamma and Inverse Gamma PDFs

The Gamma and Inverse Gamma PDFs are displayed in Fig. 2 for various κ parameters and initial β t=0 coefficients. As described previously, the κ parameter and the initial mean radius R 2 t=0 are input data. The initial β t=0 coefficient is determined with the help of the initial mean radius R 2 t=0 via Eq. (4.28). The following Gamma and Inverse Gamma probability density functions consequently describe only the initial radius distribution of a cloud of droplets.

The initial mean radius is R 2 t=0 = 10 µm for all situations. Interesting behaviors appear. The bigger the κ parameter is, the more symmetrical and the sharper the two functions become, and consequently tend towards a monodisperse distribution. This monodisperse-like behavior will be recovered later when examining the specific interfacial area. However, when κ Gamma ≤ 1, the qualitative shape of the Gamma PDF changes drastically, and the first condition of Eq. (4. [START_REF] Olmos | Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up[END_REF]) is no longer satisfied. It then appears that in the present two-phase context, the previous general restriction: κ Gamma > 0 transforms to κ Gamma > 1, as a consequence of the droplet condition: lim R2→Rmin (f (R 2 )) = 0 demanded by Eq. (4. [START_REF] Olmos | Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up[END_REF]). Nevertheless, the present behavior observed with κ Gamma ≤ 1 may be interesting in other contexts. Modeling interfacial area in porous media may be a relevant example. In that context, the mechanistic model for shock initiations of solid explosions is presented for example in [START_REF] Massoni | A mechanistic model for shock initiation of solid explosives[END_REF] [START_REF] Massoni | A mechanistic model for shock initiation of solid explosives[END_REF].

Unlike the Gamma PDF, the Inverse Gamma PDF shows only minor changes regarding its shape when the κ Inverse Gamma parameter tends to its lowest admissible limit, i.e. κ Inverse Gamma → 3. Moreover, in this same limit, the probability density function tends to spread out along the mean radius value in the desired asymmetric way. The second condition of Eq. (4.19) (lim R2→Rmax (f (R 2 )) = 0) is satisfied for both Gamma and Inverse Gamma PDFs. Indeed, beyond a certain radius R max , both probability density functions tend to zero. In the present illustration, as the β t=0 coefficient corresponds to the initial mean radius R 2 t=0 of 10 µm through Eq. (4.28), the Gamma (4.20) and Inverse-Gamma (4.21) PDFs depend only the κ parameter. Consequently, the radius R max depends only on the κ parameter as well. As observed in Fig. 2, the maximum radius R max increases as the κ parameter decreases. Oppositely, the minimum radius R min decreases with the κ parameter. The interval [R min , R max ] gets wider as κ decreases. However, for κ Gamma ≤ 1 the minimum radius R min is zero as observed in Fig. 2, and is not admissible. A reasonable estimate of R min can be determined with the help of the critical Weber number. A study based on the concept of a critical Weber number is presented in [START_REF] Pilch | Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop[END_REF] [START_REF] Pilch | Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop[END_REF] which permits prediction of the maximum size of stable fragments.

As mentioned in Section 4.2, other PDFs are suitable for the present work as long as they satisfy the conditions introduced. In Appendix A, the Log-Normal and Rosin-Rammler PDFs are presented and depicted for various coefficients. Those display distributions similar to the ones observed in Fig. 2, provided that the various coefficients are well adjusted.

Impact of the probability density functions on the interfacial area of the two-phase flow and mean

radius of the polydisperse droplets Some examples of the Gamma (4.20) and Inverse-Gamma (4.21) PDFs are displayed in Fig. 2. The previous results represent the initial radius distribution of a cloud of droplets. In the following, the impact of the κ parameter on the evolution of the specific interfacial area A I of the two-phase flow, as well as the evolution of the mean radius R 2 of the liquid droplets, are examined. A single control volume where N 2 liquid droplets are present is considered. When the liquid phase of the two-phase flow is made of spherical monodisperse droplets, the specific interfacial area A I is provided by Eq. (3.8), that is recalled hereafter. As the droplets are supposed to be spherical with radius R 2 , the specific number of droplets N 2 is provided by Eq. (3.5). This last relation is here rewritten under the form,

A I = 4πR 2 2 N 2 , R 2 = 3 α 2 4 3 π N 2 , Monodisperse. (4.30)
For a given specific number of particles N 2 and liquid volume fraction α 2 , the specific interfacial area A I and radius R 2 in the context of monodisperse droplets, are determined with the help of those last two relations. In Sections 4.1 and 4.2, the specific interfacial area A I , and consequently the mean radius R 2 , have been reconsidered to account for polydisperse droplets through distribution moments, and are provided by the second and third relations of (4.18). As the Gamma (4.20) and Inverse-Gamma (4.21) PDFs are considered, the corresponding equations reduce, for the specific interfacial area A I , to Eq. (4.26) that is recalled hereafter. Moreover the combination of Eqs. (4.25) and (4.28) yields the following expression of the mean radius R 2 ,

A I =          4πN 2 (κ + 1)κ 4πN 2 3α 2 (κ + 2)(κ + 1)κ -2 3 , 4πN 2 (κ -2)(κ -1) 3α 2 4πN 2 (κ -3)(κ -2)(κ -1) 2 3
,

R 2 =            κ 3 4πN2 3α2 (κ + 2)(κ + 1)κ , Gamma, 3 3α2 4πN2 (κ -3)(κ -2)(κ -1) (κ -1) , Inverse Gamma.
(4.31) Figure 3 compares the specific interfacial area A I and mean radius R 2 , in the monodisperse (Eq. (4.30)) and polydisperse situations (Eq. (4.31)), for various liquid volume fractions α 2 . Recall that the specific interfacial A I provided by the polydisperse relation (4.31) results from the expression of the β coefficient (Eq. (4.25)). This last expression results from the combination of the third moments of probability (Eq. (4.24)) and the two-phase relation (4.18). Consequently, unlike the previous results, the initial mean radius R 2 t=0 is not used to compute the β coefficient.

The following results depend only on the κ parameter, the amount of liquid α 2 and the specific number of droplets N 2 in the control volume. In the present section, the specific number of droplets N 2 is considered constant and is set to N 2 = 10 12 droplets per unit volume. As predicted by the previous analysis, the specific interfacial area A I , as well as the mean radius R 2 , computed with the help of the Gamma (4.20) and Inverse-Gamma (4.21) PDFs tend to the ones provided by the simplified monodisperse relation (4.30) in the event of a large κ parameter. However, as expected, A I and R 2 depart significantly from the monodisperse limit when the κ parameter takes lower values. The evolution of A I and R 2 was tested with the previous κ Gamma and κ Inverse Gamma parameters but different specific number of droplets N 2 . No major changes were observed and the corresponding results are omitted for the sake of space restrictions. Note that in the present section, the specific number of droplets N 2 is arbitrary chosen as the previous results only intent to analyze the evolution of A I and R 2 for various volume fractions α 2 . In the next section, two-phase numerical results are presented and the specific number of droplets N 2 , as well as the volume fraction α 2 , are provided by the balance equations of (2.1).

Beforehand, another point of view is considered. In the following, both the specific number of droplets N 2 and liquid volume fraction α 2 remain constant. Those are set to 10 12 particles per unit volume and 0.2 respectively. Evolutions of the specific interfacial area A I and mean radius R 2 are examined according to various κ parameters. Figure 4 depicts the corresponding evolutions. It can be mathematically proven that the present A I expressions, based on the the Gamma (4.20) and Inverse-Gamma (4.21) PDFs, tend to monodisperse relation (4.30) in the limit κ → ∞. Demonstrations are provided in Appendix B. Nevertheless, Fig. 4 shows that for κ paramaters of the order of few decade units, the specific interfacial area A I gets quite close to the monodisperse limit. The same can be said for the mean radius R 2 . Consequently, in the present context, κ Gamma ∈]1, 50] and κ Inverse Gamma ∈]3, 50] appear to be a fair interval. The proposed method is quite simple. In order to account for the polydisperse character of the two-phase flow, the specific interfacial area A I is reconsidered with the help of a continuous probability distribution. Multiple distribution functions can be used, as long as their m 0 , m 1 , m 2 and m 3 distribution moments are available. The specific interfacial area A I is determined through the third relation of (4.18). The Gamma and Inverse Gamma PDFS have been considered previously. In that context, A I is computed by Eq. (4.26) yielding only few code modifications. In addition to the initial mean radius R 2 t=0 and initial volume fraction α t=0 2 , only the κ parameter, controlling the shape of the polydisperse distribution, is requested as an input data and remains constant. Those data are also used to determine the initial specific number of droplets N t=0 2 with the help of Eqs. (4.28) and (4.29).

Numerical results

Comparison between the conventional method, considering a single droplet radius in a control volume, and the present method accounting for the polydisperse effects is now addressed. The test case consists of a 1D simplification of the explosion situation depicted in Fig. 1. A layer of either monodiperse or polydisperse water droplets is initially present in the 1D domain and is surrounded by air on both sides. Material interfaces are then initially present. On the left of the liquid layer, air is initially at an elevated pressure, and represents initial explosion conditions. On the right, air is at atmospheric conditions. As time goes on, the present initial explosion conditions yield two-phase suspensions. Both material interfaces and two-phase suspensions are then present in this numerical test. The present test case is depicted in Fig. 5. A fractional step method is used. Computation of the source terms of (2.1) is decoupled from transport and wave propagation. The hyperbolic step, i.e. the resolution of the two-phase flow equation system without pressure and velocity relaxation source terms, is addressed with the first-order Godunov (1959) scheme [START_REF] Godunov | A finite difference scheme for numerical computation of the discontinuous wave solutions of equations of fluid dynamics[END_REF] including nonconservative terms (Saurel and Abgrall, 1999 [38]). Obviously, higher-order methods may be used but add unnecessary complexity, see for example the second-order MUSCL-type method including non-conservative terms presented in [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] in a similar two-phase context. The first-order method is then preferred as the analysis is free of extra ingredients such as flux limiters and gradient computation. The numerical scheme involves resolution of the Riemann problem and is stable under the conventional CFL condition. In the following, CFL = 0.8 for all test cases. The Riemann problem is solved with the help of the HLLC-type solver of [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF] [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. Stiff pressure relaxation is considered according to the method provided in Le [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF] [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF]. Coefficient µ, present in the right-hand side of (2.1), is then considered very large (µ → ∞).

For both methods, the viscous drag is treated as a velocity source term and is computed with the help of Eqs. (3.7) and (3.10) reminded hereafter,

N 2 F 1→2 = λ (u 1 -u 2 ) with λ = 1 8 C d ρ 1 |u 2 -u 1 | × A I . ( 5
Re 2 = 2R 2 ρ 1 |u 2 -u 1 | µ 1 .
(5.2)

A particle radius R 2 is then needed. The very purpose of the present computations is to highlight the effects of the polydisperse character of the liquid cloud on the specific interfacial area A I only. Consequently, the R 2 radius involved in the computation of the drag coefficient C d is determined with the help of Eq. (3.5) for both methods, Those relations are recalled hereafter,

R 2 = 3 α 2 4 3 π N 2 . ( 5
A I =                4πR 2 2 N 2 , Monodisperse, 4πN 2 (κ + 1)κ 4πN 2 3α 2 (κ + 2)(κ + 1)κ -2 3 , Polydisperse Gamma, 4πN 2 (κ -2)(κ -1) 3α 2 4πN 2 (κ -3)(κ -2)(κ -1) 2 3
, Polydisperse Inverse Gamma.

(5.4)

Recall that α 2 and N 2 are provided by the balance equations of (2.1). For the sake of clarity, only the Inverse Gamma PDF is used in the following. For both phases, the Stiffened-Gas equation of state (Le Métayer et al., 2004 [50]) is considered,

p k (ρ k , e k ) = (γ k -1) ρ k e k -γ k p ∞,k , (5.5)
with γ k the ratio of heat capacities (γ k = C p,k /C v,k ) and p ∞,k a constant related to the attractive effects of phase k. The dynamic viscosity of the gas phase is µ 1 = 1.8 × 10 -5 Pa. s. The various equation-of-state parameters are provided in Table 2. Table 2: Stiffened-Gas coefficients for air and liquid water [START_REF] Chiapolino | Numerical investigations of two-phase finger-like instabilities[END_REF].

Monodisperse

First, results provided by the monodisperse simplification are presented. In the following, three initial particle radii are considered, namely R t=0 . A single size of particle is considered in a control volume. However, the radius R 2 evolves with time through Eq. ( 5.3) as the liquid volume fraction α 2 and the specific number of particle N 2 are provided by the balance equations (2.1). The corresponding results are given in Figs. 6, 7 and8. Computation in the absence of viscous drag force is also considered for the sake of comparison. x (m) x (m) = 300 µm. Computation in the absence of viscous drag force is also considered. All results are given at time t = 1.2 ms. The air volume fraction α1, the pressure p1 = p2 = p, the air and liquid densities, respectively ρ1 and ρ2, as well as the air and liquid speeds, respectively u1 and u2, are presented. A close-up view of the liquid density ρ2 is also shown. x (m) x (m) Results provided by Fig. 6 indicate that the larger the liquid particles are, the faster the incident shock wave goes. Indeed, the initial amount of liquid is the same for all computations: α t=0 2 = 0.9999 in the liquid layer. The larger the droplets are in a control volume, the smaller their number N 2 is as verified in Fig. 7. Consequently, the interfacial A I in this very same control volume is lesser when the droplets are large as seen in Fig. 8.

α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u (m/s) ρ 2 (kg/m 3 )
N 2 (m -3 ) R 2 (m) N 2 (m -3 ) R 2 (m) N 2 (m -3 ) R 2 (m)
R 2 =30 m Initial R 2 =300 m Initial R 2 =30 m R 2 =300 m x (m) x (m) A I (m 2 /m 3 ) A I (m 2 /m 3 )
The specific interfacial area A I describes the available interaction surface between the gas and the liquid phases in a control volume and affects the viscous drag force (Eq. (5.1)). When the interfacial area is small, the viscous interactions are weak between the two phases and the shock wave initiated by the high-pressure air travels faster. In that event, the speed of the air u 1 and of speed of the liquid droplets u 2 are quite different as seen in Fig. 6. The largest difference is found in the absence of viscous drag where consequently no interaction (other than the one described by the non-conservative terms) is present. The incident shock wave is also the fastest in these circumstances. Oppositely, when the specific interfacial area is large as a result of the smallest particle radius (R t=0 2 = 3 µm), the viscous interactions between the two phases are strong and the incident shock wave is the slowest. The liquid droplets being small, those are easily dragged by the carrier air phase. Indeed, Fig. 6 shows that the speed of the droplets u 2 is quasi identical to the speed u 1 of the carrier phase.

The size R 2 of the liquid droplets clearly impacts the two-phase flow solution through the specific interfacial area A I and consequently through the viscous interactions between the gas and the liquid phases. The previous computations are simplified in the sense that a single size of particles is considered in a control volume. In the following, the polydisperse aspect of the liquid droplets is accounted for according to the method presented earlier. All corresponding results are plotted against those provided by the monodisperse simplification with an initial radius of R t=0 2 = 30 µm.

Polydisperse

The initialization of the two-phase flow conditions is done according to the method presented in Section 4.2.1. As previously with the monodisperse computations, the initial volume fraction α t=0 2 of the liquid phase is considered as an input data, as it is usually the case in two-phase flow computations. The initial mean radius, i.e. R 2 t=0 for the polydisperse computation or R t=0 2 for the monodisperse computation, is also known either as an input data or with the help of Eq. (4.27) ensuring a common initial specific interfacial area A t=0 I,mono. = A t=0 I,poly. as detailed in Appendix C. Indeed, for a proper comparison of the two computations, the following test cases consider either a common initial specific interfacial area A t=0 I,mono. = A t=0 I,poly. or a common mean radius R 2 t=0 = R t=0 2 between the two computations. For the polydisperse computations the Inverse Gamma probability density function is considered with various κ parameters. The resulting initial droplet distributions are depicted in Fig. 9 for the two events, i.e. common initial interfacial area or common initial mean radius between the monodisperse and polydisperse computations.

Common initial interfacial area A t=0

I between the polydisperse and monodisperse computations The analysis begins with the first situation, i.e. the initial specific interfacial area is common between the two computations A t=0 I,mono. = A t=0 I,poly. . Examining the evolution of the polydisperse solution when initialized with the same initial specific interfacial area as the one of the monodisperse situation is important. Indeed, as seen previously with the monodisperse computations, A I plays a major role in the viscous interactions between the phases. The following results show how the polydisperse solution departs from the monodisperse simplification and show the impact of the polydisperse character of the liquid droplets on their mean radius R 2 and volume fraction α 2 .

The initial mean radius of the monodisperse computation is R t=0 2 = 30 µm. Two κ parameters are considered, κ = 5 and κ = 3.1. As seen in Fig. 9 those two parameters lead to very different initial size distributions. Yet, those two distributions yield the same initial specific interfacial area A t=0 I that corresponds to one of the monodisperse computation. When κ = 3.1, the PDF indicates a high probability to get droplets of radii R 2 ≃ 1 µm. However, such a low value appears quite close to the minimum physical radius R min estimated with the help of the critical Weber number, and may consequently be problematic as regards to the physical representation of liquid droplets. The present size distribution mathematically provides A t=0 I,mono. = A t=0 I,poly. but appears fictitious. We will come back to this point a bit further. Beforehand, the monodiperse results and the polydisperse results obtained with κ = 5 are compared. The corresponding results are presented in Figs. 10 and11. x (m) x (m) The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 5. All results are given at time t = 1.2 ms and are presented in terms of α1, p, ρ1, ρ2, u1 and u2 for both computations. The two computed solutions are very close, only slight differences appear in the various profiles. The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 5. All results are given at time t = 1.2 ms and are presented in terms of AI , u, N2 and R2 for both computations. The two computed solutions are very close, only slight differences appear in the various profiles.

α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u 1 (m/s) u 2 (m/s)
In the present conditions, when κ Inverse Gamma ≥ 5, the different flow variables of the monodisperse and polydisperse computations are quasi merged, except for the mean radius R 2 and the specific number N 2 of the liquid droplets. However the resulting specific interfacial areas A I are very close. It then appears that with the present conditions, the polydisperse effects are correctly reproduced by the monodisperse simplification. The test is now repeated with a lower value of the κ parameter. This last one is now set to κ = 3.1 and is close to its admissible lower limit. Figures 12 and13 present the corresponding results. x (m) x (m) The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 3.1. All results are given at time t = 1.2 ms and are presented in terms of α1, p, ρ1, ρ2, u1 and u2 for both computations. As the κ parameter tends towards its lower limit, differences between the two solutions are visible. x (m) x (m)

α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u 1 (m/s) u 2 (m/s)
A I (m 2 /m 3 ) u (m/s) N 2 Mono. (m -3 ) N 2 Poly. (m -3 ) R 2 Mono. (m) R 2 Poly. (m)
Figure 13: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 3.1. All results are given at time t = 1.2 ms and are presented in terms of AI , u, N2 and R2 for both computations. As the κ parameter tends towards its lower limit, differences between the two solutions are visible.

The κ parameter being close to its lower limit, the impact of the PDF, taking into account the polydisperse effects, is now clearly seen on the various flow variables. Indeed, the present results show the impact of the particle size on the specific interfacial area A I and consequently on the flow variables (α, p, ρ, u) through the viscous interactions between the carrier gas phase and the dispersed liquid phase. Nevertheless, as indicated earlier by Fig. 9, the mean radius R 2 of the polydisperse droplets is physically questionable as it tends to the minimum physical radius R min ≃ 1 µm estimated with the help of the critical Weber number. Those very small particles then appear fictitious. Such an issue is controlled with the second initialization option that is now addressed. x (m) x (m) The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 50. All results are given at time t = 1.2 ms and are presented in terms of α1, p, ρ1, ρ2, u1 and u2 for both computations. The two solutions are very close. One more time, the κ parameter being quite large, the interfacial area A I of the polydisperse case is very close to the interfacial area provided by the monodisperse computation, as seen in Fig. 15. Consequently, the two solutions are in close agreement as well. The test is now repeated with a lower value: κ = 7. The corresponding results are provided in Figs. 16 and17. x (m) x (m) The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 7. All results are given at time t = 1.2 ms and are presented in terms of α1, p, ρ1, ρ2, u1 and u2 for both computations. Differences between the two solutions are visible. The κ parameter being lower, the specific interfacial areas A I of the two computations are clearly different, as seen in Fig. 17. The interfacial area provided by the polydisperse computation is less than the one provided by the monodisperse computation. The two-phase solutions are consequently different. Indeed, the available interaction surface between the gas and liquid phases is lesser and yields consequently a faster shock wave as suggested by the monodiperse analysis provided in Section 5.1. In the following, the κ parameter is lowered further, to κ = 3.5. The corresponding results are provided in Figs. 18 and19. x (m) x (m) The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 3.5. All results are given at time t = 1.2 ms and are presented in terms of α1, p, ρ1, ρ2, u1 and u2 for both computations. Differences between the two solutions are visible. The κ parameter is closer to its lower limit (κ = 3). The specific interfacial area A I provided by the polydisperse computation is significantly less than the one provided by the monodisperse computation as seen in Fig. 19. The polydisperse effects consequently affect significantly the shock wave that is much faster due to the lesser available interaction surface between the two phases. The flow variables (α, ρ, p, u) are consequently affected as well. The present results are similar to the one provided by the monodisperse simplification in the event of large particles yielding weak viscous interactions (Figs. 6, 7 and 8 of Section 5.1).

Common initial mean radius

α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u 1 (m/s) u 2 (m/s)
x (m) x (m) R 2 (m) u (m/s) A I (m 2 /m 3 ) N 2 (m -3 )
α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u 1 (m/s) u 2 (m/s)
x (m) x (m) R 2 (m) u (m/s) A I (m 2 /m 3 ) N 2 (m -3 )
α 1 p (Pa) ρ 1 (kg/m 3 ) ρ 2 (kg/m 3 ) u 1 (m/s) u 2 (m/s)
The previous results highlight the contribution of the polydisperse effects on the two phase flow, especially on the viscous interactions appearing in the relaxation zone between the gas and liquid phases. Those viscous interactions, and consequently the size of the relaxation zone, are controlled by the specific interfacial area, a key point in combustion and two-phase flow modeling. As seen in Fig. 19, the mean radius R 2 of the polydisperse solution is quite close to the single radius R 2 delivered by the simplified monodisperse computation. However, polydisperse droplets have a major impact on the specific interfacial area and consequently on the flow variables.

Two initialization options have been tested. The first one involves an identical initial specific interfacial area A t=0 I,mono. = A t=0 I,poly. for the monodisperse and polydisperse computations. The second involves an identical mean radius of the droplets R 2 t=0 = R t=0 2 . For both cases, the monodisperse and polydisperse solutions are quasi merged when the κ parameter is large enough as suggested by the analysis carried out in Section 4.3. For lower κ values, clear differences appear. Nevertheless, with the first set of results (identical initial A t=0 I ) the resulting mean radius R 2 t=0 of the polydisperse case may be very low and physically questionable. Yet, the second initialization option (identical initial R 2 t=0 = R t=0 2 ) gives control on the present issue and clearly shows the impact of the polydisperse effects on the relaxation zone between the gas and liquid phases.

Conclusion

Determination of the specific interfacial area A I is a key problem in two-phase flow modeling. In many situations the liquid phase is said to be polydisperse as it contains a substantial number of droplets of different sizes, making major effects on the two-phase flow. In the present paper, explosion situations are of particular interest. In such circumstances, material interfaces are present at the early times, as well as two-phase suspensions occurring at later timescales. Both interfacial and disperse flow conditions are then present and the two-phase flow model must be able to deal with both situations.

The specific interfacial area A I has been reconsidered to account for the polydisperse character of the liquid phase in a simplified way. Computation of A I relies on a continuous probability distribution. Gamma-like probability density functions have been used in the present work. However, the method may be used with various functions as long as their m 0 , m 1 , m 2 and m 3 moments are available. In the context of Gamma and Inverse Gamma distributions, only the κ constant is required as an input data, in addition to the initial mean radius R 2 t=0 and initial volume fraction

α t=0
2 that are necessary both for the conventional method and for the proposed method. The κ constant controls the shape of the polydisperse distribution and can be adapted to multiple situations. The overall method yields only few code modifications while taking into account the polydisperse aspect of the two-phase flow.

The two-phase equation system of [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF] [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shockbubble interaction[END_REF], variant of Baer and Nunziato's (BN) model (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], has been considered in the present work, as it is able to deal with both material interfaces and two-phase suspensions. The impact of the polydisperse character of the liquid phase has been highlighted with the help of a simplified 1D two-phase explosion test. This work can be continued in many directions. Among them are the introduction of fragmentation effects of the liquid droplets, as well as thermal effects and combustion processes. Comparison with experimental results is also part of future works. A.3). The βt=0, νt=0 and ηt=0 coefficients are computed via Eqs. (4.28), (A.3), (A.4), such that the initial mean radius R2 t=0 of the polydisperse particles is

R2 t=0 = 10 µm.
As observed in Section 4.2.2, the bigger the κ parameter is, the more symmetrical and the sharper the Gamma and Inverse-Gamma functions become, and consequently tend towards a monodisperse distribution, as seen in plot a) of As the droplets are supposed to be spherical, their mean radius, in the monodisperse case, is given by Eq. (4.30) that is recalled hereafter:

Inverse Gamma

The same reasoning is repeated for the Inverse-Gamma expression,

A Inverse I Gamma = 4πN 2 (κ -2)(κ -1) 3α 2 4πN 2 (κ -3)(κ -2)(κ -1) 2 3
(B.6) In Section 5, the initial solutions of the simplified 1D two-phase explosion test are determined according to either the first initialization option or the second one. The initial liquid volume fraction α t=0 2 is known. The constant κ parameter is also known as it is an input data. With the second initialization option, the initial mean radius R 2 t=0 of the polydisperse droplets is also considered as an input data, and is identical for both situations, i.e. the monodisperse and the polydisperse computations, R 2 t=0 = R t=0 2 . The present appendix is related to the first initialization option where the initial specific interfacial area A t=0 I is identical for both situations. In order to find a common initial A t=0 I for both computations, the initial mean radius R 2 t=0 of the polydisperse case is adapted and is provided by Eq. (4.27). The corresponding demonstration is addressed hereafter. When dealing with a monodisperse situation the initial specific interfacial area A t=0 I is computed via Eq. (3.8). However, when dealing with a polydisperse situation, the initial specific interfacial areas, combined with the Gamma and Inverse-Gamma PDFs, are provided by Eq. (4.26). With the help of Eq. (4.25), the initial interfacial areas A t=0 I become:

A Inverse I Gamma = 4πN 2 1 κ 2 1 -3 κ + 2 κ 2
A t=0 I =          4πR 2 t=0

(4. 29 )

 29 Finally, the initial interfacial area A t=0 I is computed via the third relation of(4.18), reducing to Eq. (4.26) in the context of Gamma and Inverse-Gamma PDFs. The initialization of the two-phase flow conditions consequently requires the initial volume fraction α t=0 2 and the initial mean radius of the droplets R 2 t=0 . When dealing with a monodisperse situation, the initial mean radius consists of the initial single and common radius R t=0 2 of the droplets. When dealing with polydisperse droplets, only the constant κ parameter is required as well, in addition to α t=0 2 and R 2 t=0 .

1 Figure 2 :

 12 Figure 2: Initial Gamma and Inverse Gamma probability density functions (Eqs. (4.20) and (4.21)) versus the particle radius, for various κ parameters. The βt=0 coefficient is computed via Eq. (4.28) such that the initial mean radius R2 t=0 of the polydisperse particles is R2 t=0 = 10 µm.

1 Figure 3 :

 13 Figure 3: Comparison of the interfacial area AI and mean radius R2, in the monodisperse and polydisperse situations, for various liquid volume fractions α2. The specific number of droplets N2 is constant and is set to N2 = 10 12 droplets per unit volume. For the monodisperse situation, AI and R2 are computed via Eq. (4.30). For the polydisperse situation, AI and R2 are computed via Eq. (4.31).

Figure 4 :

 4 Figure 4: Comparison of the interfacial area AI and mean radius R2, in the monodisperse and polydisperse situations, for various κ parameters. The specific number of droplets N2 and liquid volume fraction α2 are constant and are set to N2 = 10 12 droplets per unit volume and α2 = 0.2. For the monodisperse situation, AI and R2 are computed via Eq. (4.30). For the polydisperse situation, AI and R2 are computed via Eq. (4.31).

7 3 ρ 2 = 1050 kg/m 3 3 ρ 2 = 1050 kg/m 3 3 ρ 2 =

 732332332 Pa ρ 1 = 12 kg/m Atmospheric air α 1 = 0.9999 p = 10 5 Pa ρ 1 = 1.2 kg/m Liquid layer α 1 = 0.0001 p = 10 5 Pa ρ 1 = 1.2 kg/m 1050 kg/m 3

Figure 5 :

 5 Figure5: Simplified two-phase explosion test. A layer of liquid water is initially present in a 1D domain. The liquid volume fraction is α2 = 0.9999 in this zone, and atmospheric conditions are considered p = 10 5 Pa, ρ2 = 1050 kg/m -3 . The liquid layer is surrounded by air on both sides. Material interfaces are then initially present. On the left side, the air is initially dense and at elevated pressure: p = 10 7 Pa, ρ1 = 12 kg/m -3 , and represents initial explosion conditions. On the right side, the air is at atmospheric conditions: p = 10 5 Pa, ρ1 = 1.2 kg/m -3 . The air volume fraction is initially α1 = 0.9999 on both sides. The 1D domain is 2.5 m long. The liquid layer is placed at abscissa x = 1.4 m and is 0.025 m wide. The mesh consists on 2, 500 regular elements, yielding a space step of ∆x = 1 mm. Boundary conditions are non-reflective.

. 3 )

 3 Finally, the specific interfacial area A I is computed both through its monodisperse simplified formulation (3.8) and the proposed relation (4.18) accounting for the polydisperse character of the liquid cloud. As the Gamma (4.20) and Inverse-Gamma (4.21) PDFs are considered in the present work, the corresponding equations reduce to Eq. (4.26).

2 = 3 µm, R t=0 2 = 30 µm and R t=0 2 = 2 = (3α 2 )

 232222 300 µm. Those three radii lead to three initial specific numbers of particle N t=0 2 through Eq. (3.5): N t=0 / 4π R 3,t=0 2

Figure 6 :R t=0 2 = 3 µm, R t=0 2 = 30

 623230 Figure 6: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse results with various initial radii: R t=0 2 = 3 µm, R t=0 2 = 30 µm and R t=0 2

Figure 7 :R t=0 2 = 3 2 = 2 =

 72322 Figure 7: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse results with various initial radii: R t=0 2 = 3 µm (top plots), R t=0 2 = 30 µm (middle plots) and R t=0 2 = 300 µm (bottom plots). All results are given at time t = 1.2 ms. The specific number N2 and the radius R2 of the liquid particles are presented.

Figure 8 : 2 = 3 µm, R t=0 2 = 30 µm and R t=0 2 =

 82322 Figure 8: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse results with various initial radii: R t=0 2 = 3 µm, R t=0 2 = 30 µm and R t=0 2 = 300 µm. All results are given at time t = 1.2 ms. The specific interfacial area AI is presented for all radii on the left. On the right, a close-up view is shown.

2 Figure 9 : 2 =R t=0 2 =

 2922 Figure 9: Initial droplet distributions provided by the Inverse Gamma density function with various κ parameters. On the left, a common initial interfacial area A t=0 I,mono. = A t=0 I,poly. is considered between the monodisperse and polydisperse computations through an appropriate determination of the initial mean radius R2 t=0 of the polydisperse computation. The single radius R t=0 2 = 30 µm of the monodisperse computation is considered as an input data. The initial mean radius R2 t=0 of the polydisperse computation is determined with the help of Eq. (4.27) (see Appendix C). The initial mean radius R2 t=0 is used to find the βt=0 coefficient (Eq. (4.28) Section 4.2.1). On the right, the initial mean radius R2 t=0 is directly used as an input data. The initial monodisperse and polydisperse mean radii are considered equal: R2 t=0 = R t=0 2 = 30 µm.

Figure 10 :

 10 Figure 10: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results.The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0

AN 2 R 2 Figure 11 :

 2211 Figure 11: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results.The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0

Figure 12 :

 12 Figure 12: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results.The initial specific interfacial area is the same for both computations A t=0 I,mono. = A t=0 I,poly. and is computed via an initial monodisperse radius of R t=0

R 2 t=0 2 =

 22 between the polydisperse and monodisperse computations Unlike the previous results, the initial mean radius R 2 t=0 of the polydisperse computation and the initial single radius R t=0 2 of the monodisperse computation are now identical. The initial specific interfacial areas are consequently different A t=0 I,mono. = A t=0 I,poly. , as one is based on a polydisperse distribution of the droplets and the other supposes a single size (Eq. (5.4)). For the following test case, the initial mean radius is set to R 2 t=0 = R t=0 30 µm. Figures 14 and 15 compare the results provided by the monodisperse formulation of the specific interfacial area and those provided by the polydisperse relation (Eq. (5.4)). In this last situation, the Inverse Gamma PDF is used with κ = 50.

Figure 14 :

 14 Figure 14: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0

Figure 15 :R t=0 2 =

 152 Figure 15: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 50. All results are given at time t = 1.2 ms and are presented in terms of AI , u, N2 and R2 for both computations. The two solutions are very close.

Figure 16 :

 16 Figure 16: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0

Figure 17 :R t=0 2 =

 172 Figure 17: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 7. All results are given at time t = 1.2 ms and are presented in terms of AI , u, N2 and R2 for both computations. Differences between the two solutions are visible.

Figure 18 :

 18 Figure 18: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0

A I (m 2 /m 3 ) N 2 (m - 3 )Figure 19 :R t=0 2 =

 2323192 Figure 19: Simplified two-phase explosion test depicted in Fig. 5. Comparison of the monodisperse and polydisperse results. The initial mean radius is the same for both computations R2 t=0 = R t=0 2 = 30 µm. For the polydisperse computation, the Inverse Gamma distribution is used with κ = 3.5. All results are given at time t = 1.2 ms and are presented in terms of AI , u, N2 and R2 for both computations. Differences between the two solutions are visible.

Table A. 3 :Figure A. 20 :

 320 Figure A.20: Initial Gamma, Inverse Gamma, Log-Normal and Rosin-Rammler probability density functions (Eqs. (4.20), (4.21), (A.1), (A.2)) versus the particle radius, for various κ, σ and δ parameters (TableA.3). The βt=0, νt=0 and ηt=0 coefficients are computed via Eqs. (4.28), (A.3), (A.4), such that the initial mean radius R2 t=0 of the polydisperse particles is

Fig. A. 20 . 1 2 )

 2012 This behavior is recovered with the Rosin-Rammler PDF when the δ parameter is large enough as well. It is also recovered with the Log-Normal PDF but in the event of a low σ parameter. When the various parameters depart from their respective upper and lower limits, all PDFs present a similar behavior, adapted to the present two-phase context, as seen in plot b) of Fig.A.20. and can be rewritten as follows,A I Gamma = 4πN 2 κ 2The limit of Eq. (B.2) in the event κ → ∞ is,

R 2 = 3α 2 4πN 2 1 3 . (B. 4 )

 234 MonodisperseThe combination of Eqs. (B.3) and (B.4) yields the monodisperse expression (4.30) of the specific interfacial area,lim κ→∞ A I Gamma = 4πN 2 R 2 2 Monodisperse = A I Monodisperse . (B.5)

3 (B. 8 )

 38 The combination of Eqs. (B.4) and (B.8) yields the monodisperse expression (4.30) of the specific interfacial area, lim κ→∞ A Inverse I Gamma = 4πN 2 R 2 2 Monodisperse = A I Monodisperse . (B.9) Appendix C. Common initial specific interfacial area

2 (

 2 κ -2)(κ -1) β 2 t=0 , Inverse Gamma. (C.1)

Table 1

 1

	κ Gamma	R min	R max	κ Inverse Gamma	R min	R max
	30	4.5 µm	18 µm	30	5 µm	20 µm
	5	0.7 µm	33 µm	5	2 µm	50 µm
	1	0 µm	70 µm	3.1	1.2 µm	70 µm
	0.5	0 µm	70 µm	-	-	-

reports the radius R max necessary to satisfy f (R 2 ) = 0.001 × max (f (R 2 )).

Table 1 :

 1 Minimum and maximum radii for which the initial Gamma (4.20) and Inverse-Gamma (4.21) PDFs are numerically zero. Outside from the radius interval [Rmin, Rmax], the PDFs are less than ǫ = 0.001 × max (f (R2)).

Appendix A. Log-Normal and Rosin-Rammler probability density functions

In the present work, the Gamma and Inverse-Gamma distributions are considered. However, the method is not restricted to these functions and can be extended to other probability density functions such as the Log-Normal and Rosin-Rammler distributions. Those also consist of two-parameter families of continuous probability distributions. Their probability density functions read,

As for the Gamma and Inverse-Gamma PDFs, one parameter is considered as an input data and the second is computed with the help of the initial mean radius R 2 t=0 . Indeed, recall that for the Gamma and Inverse-Gamma PDFs, the κ parameter is given as an initial data and the initial β t=0 coefficient is computed from Eq. (4.28). For the Log-Normal and Rosin-Rammler PDFs, the σ and δ parameters are considered as input data respectively. The ν t=0 and η t=0 coefficients are computed from the initial mean radius R 2 t=0 as is consists of the first moment (m 1 ) of the corresponding probability function. The moments of the Log-Normal and Rosin-Rammler PDFs read:

The ν t=0 and η t=0 coefficients are then computed from Eqs. (A.3) and (A.4), while the σ and δ are free parameters. A.3. The initial mean radius is set to R 2 t=0 = 10 µm for all functions.

Appendix B. Monodisperse limit of the interfacial area

As mentioned in Section 4.3, it can be mathematically proven that the A I expressions (4.31), based on the the Gamma (4.20) and Inverse-Gamma (4.21) PDFs, tend to the monodisperse relation (4.30) in the limit κ → ∞. The demonstration is provided hereafter.

Gamma

The specific interfacial area reads, when combined to the Gamma PDF,

The initial number of droplets N t=0 2 is provided by Eq. (3.5) in the monodisperse case and Eq. (4.29) for the polydisperse case. Those are recalled hereafter,

The introduction of the initial number of droplets N t=0

The initial man radius R 2 t=0 of the polydisperse droplets is provided by Eq. (4.28) that is recalled hereafter,

The combination of Eqs. (C.3) and (C.4) yields,

The initial liquid volume fraction α t=0 2 is supposed to be the same for the three radius distribution cases. After some algebraic manipulations, equaling the initial specific interfacial area A t=0 I provided by the monodisperse distribution to the initial specific interfacial area provided by either the Gamma or Inverse-Gamma PDF leads to the following relation,

The κ parameter is an input data. When the first initialization option is considered, the initial single radius of the monodisperse situation R t=0 2 is also known as an input data. The initial mean radius R 2 t=0 of the droplets, in the polydisperse situation, is then computed with the help of Eq. (C.6), and ensures a common initial specific interfacial area A t=0 I between the monodisperse and polydisperse computations.