2.1. Two-scale modelling assumptions

We propose successive sets of assumptions. First, regarding the fluid dynamics:

(H1a) \begin{align} {(a)}\ &\textrm{all phases have the same mean velocity;} \end{align}

(H1b) \begin{align} {(b)}\ &\textrm{the small-scale velocity fluctuations and capillarity forces } \nonumber\\ & \textrm{are modelled by a small-scale energy term}; \end{align}

(H1c) \begin{align} {(c)}\ & \textrm{the small-scale liquid phase is incompressible.} \end{align}

Following hypotheses (H1a) and (H1b), we account for velocities that naturally appear at small scale (e.g. droplets’ drag in the disperse regime) only through the advection of a local small-scale energy for the sake of clarity of the exposition of the new ideas of the paper. Potential extensions are discussed in the conclusion. The incompressibility of the small-scale is discussed in Appendix A to justify (H1c). It affects the physics at two levels. The first level pertains to sound propagation that is cancelled at the small scale. It is especially relevant in the disperse regime where no sound can propagate from one inclusion to another as underlined in Saurel et al. (Reference Saurel, Chinnayya and Carmouze2017). The second level is related to the volume occupation of the disperse phase and the evolution of its volume fraction.

Second, the following hypotheses are formulated solely for the sake of simplicity, but are not limiting for the models developed in the present work:

(H2a) \begin{align} {(a)}\ &\textrm{each phase is isothermal with barotropic equations of state;} \end{align}

(H2b) \begin{align} {(b)}\ &\textrm{large-scale capillarity effects are not modelled;} \end{align}

(H2c) \begin{align} {(c)}\ &\textrm{no mass exchanges between the phases are taken into account;} \end{align}

(H2d) \begin{align} {(d)}\ &\textrm{body forces are neglected.} \end{align}

Let us remark that (H1a) and (H2a) correspond to the thermal and kinematic equilibria as studied in Chanteperdrix, Villedieu & Vila (Reference Chanteperdrix, Villedieu and Vila2002), Caro et al. (Reference Caro, Coquel, Jamet and Kokh2005) and Drui et al. (Reference Drui, Larat, Kokh and Massot2019). The effects neglected by (H2a), (H2b) and (H2c) have been studied, for instance, in Perigaud & Saurel (Reference Perigaud and Saurel2005), Gaillard (Reference Gaillard2015), Schmidmayer et al. (Reference Schmidmayer, Petitpas, Daniel, Favrie and Gavrilyuk2017), Cordesse (Reference Cordesse2020), Loison et al. (Reference Loison, Kokh, Pichard and Massot2024), Caro et al. (Reference Caro, Coquel, Jamet and Kokh2005) and Pelanti (Reference Pelanti2022). The model can be extended to include these effects, but they are not necessary for the two-scale modelling framework at stake in this work. The particular modelling of mass transfer between the large scale and the small scale is not considered in this section but accounted in § 5.

We denote with the subscript $k=1$ the liquid and $k=2$ the gaseous phase. An additional superscript $d$ identifies the small-scale liquid phase, where the letter $d$ stands for ‘droplets’ or ‘disperse’ liquid phase. Following (H2a), each large-scale fluid $k=1,2$ is equipped with a barotropic equation of state of the form $\rho _k \mapsto e_k (\rho _k)$, where $\rho _k$ and $e_k$ are the density and the barotropic potential of the phase $k=1,2$. Denoting $({\cdot })\mapsto ({\cdot })'$ the derivative with respect to the density, the pressure $p_k$ of the large-scale phase $k=1,2$ is then defined by $p_k = \rho _k^2 e_k'(\rho _k)$. We assume that $p_k'(\rho _k)>0$ so that the sound speed $c_k$ associated with the phase $k=1,2$ is $c_k=\sqrt {p_k'(\rho _k)}$. With $\alpha _k$ the volume fraction of the phase $k=1,2,1^d$, we consider a two-scale mixture such that it obeys the volume constraint

(2.1)\begin{equation} \alpha_1 + \alpha_2 + \alpha_1^d = 1, \end{equation}

and its density $\rho$ is obtained from the effective phase densities $m_k := \alpha _k\rho _k$,

(2.2)\begin{equation} m_1 + m_2 + m_1^d = \rho. \end{equation}

Following Loison et al. (Reference Loison, Kokh, Pichard and Massot2024), we introduce the large-scale volume fractions defined by

(2.3)\begin{equation} \bar{\alpha}_k:=\frac{\alpha_k}{1-\alpha_1^d} \quad \textrm{for}\ k=1,2, \end{equation}

extending the usual definition at the large-scale perspective such that

(2.4)\begin{equation} \bar{\alpha}_1 + \bar{\alpha}_2=1. \end{equation}

We denote by $Y_k$ the mass fraction such that the specific barotropic potential $e$ of the two-phase material is defined by

(2.5)\begin{equation} e = Y_1 e_1(\rho_1) + Y_2 e_2(\rho_2) + Y_1^d e_1(\rho_1^d). \end{equation}

Under assumptions (H1a) and (H1b), there is a unique velocity $\boldsymbol {u}$ describing the mixture

(2.6)\begin{equation} \boldsymbol{u}:= \boldsymbol{u}_1 = \boldsymbol{u}_2 = \boldsymbol{u}_1^d. \end{equation}

Then, (H2c) provides the mass conservation of each phase

(2.7)\begin{equation} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u}) = 0 \quad \textrm{for} \ k=1,2,1^d. \end{equation}

Summing these equations provides the total mass conservation equation

(2.8)\begin{equation} \partial_t \rho + \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot} (\rho \boldsymbol{u})=0. \end{equation}

Denoting $D_t({\cdot })=\partial _t({\cdot }) + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } ({\cdot })$ the material time derivative, the incompressibility of the small-scale liquid phase (H1c) reads

(2.9)\begin{equation} D_t \rho_1^d = 0. \end{equation}

Then liquid volume fraction of the small scale is conserved

(2.10)\begin{equation} \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0. \end{equation}

The small-scale energy density accounting for both capillarity energy and kinetic energy of velocity fluctuations is advected at the mixture averaged velocity

(2.11)\begin{equation} D_t \zeta = 0. \end{equation}

2.2. Conservative motion with Hamilton's SAP

The equation of the two-scale fluid mixture's motion is derived similarly as Gavrilyuk & Gouin (Reference Gavrilyuk and Gouin1999) and Drui et al. (Reference Drui, Larat, Kokh and Massot2019) by using Hamilton's SAP. We set the kinetic and potential energies of the phase $k$ as

(2.12a,b)\begin{equation} E_k^{kin}:= \tfrac{1}{2}m_k \boldsymbol{u}^2, \quad E_k^{pot}:= m_k e_k\left(\frac{m_k}{\alpha_k}\right), \end{equation}

and the mixture kinetic and potential energies are $\sum _k E_k^{kin}$ and $\sum _k E_k^{pot}$. This enables the definition of the Lagrangian energy ${\mathcal {L}}_k=E_k^{kin}-E_k^{pot}$ of the phases $k=1,2$, ${\mathcal {L}}_1^d=E_1^{d,kin}-E_1^{d,pot}-m_1^d\zeta$ and a mixture Lagrangian energy ${\mathcal {L}}=\sum _k{\mathcal {L}}_k$. Capillarity energy at large scale is purposely not modelled for simplifying purposes (see (H2b)) and no particular flow regime is assumed for now at small scale. These effects are considered later in § 3, where we refine the small-scale model. Thereby, the Lagrangian associated with our system reads

(2.13)\begin{equation} {\mathcal{L}} = {\mathcal{L}}_1(\alpha_1, m_1, \boldsymbol{u}) +{\mathcal{L}}_2(\alpha_2, m_2, \boldsymbol{u}) +{\mathcal{L}}_1^d(m_1^d, \rho_1^d, \zeta, \boldsymbol{u}). \end{equation}

With mass conservation constraints, Hamilton's SAP leads to the following set of equations (see Appendix B for details)

(2.14)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0,\quad k=1,2,1^d, \\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u})=0,\\ \partial_t \zeta + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \zeta=0,\\ \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u} + p \boldsymbol{\mathsf{I}})=0, \end{array}\right\} \end{equation}

where $\boldsymbol{\mathsf{I}}$ is the identity matrix, and $p$ is the equilibrium pressure obtained for the mixture thermodynamic closure $p:=p_1=p_2$ given by the second line of (B16). This algebraic equation gives $\alpha _1$ and $p$, respectively, as the solution and the value of the equilibrium for given values of $m_1$, $m_2$ and $\alpha _1^d$. System (2.14) also admits a supplementary conservation equation

(2.15)\begin{equation} \partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot}(({\mathcal{H}}+p)\boldsymbol{u})=0, \end{equation}

where the total energy ${\mathcal {H}} = \rho \Vert \boldsymbol {u}\Vert ^2 - {\mathcal {L}}$ is a mathematical entropy (Godlewski & Raviart Reference Godlewski and Raviart1991) for (2.14). The mathematical entropy of the system is studied as the isothermal limit of the Euler–Fourier model in Serre (Reference Serre2010). It is shown to be convex and linked to the physical entropy of the mixture $s$ with ${\mathcal {H}}=\rho (\varepsilon -Ts) + \tfrac {1}{2}\rho \Vert \boldsymbol {u}\Vert ^2$ where $\varepsilon = e + Ts$ and $T$ are the internal energy and the temperature of the mixture, respectively.

Remark that $\alpha _1$ is not a variable of the set of conservation laws (2.14), but its effective dynamics for smooth solution can be written explicitly by applying the material time derivative to the algebraic equation $p_1(m_1/\alpha _1)=p_2(m_2/(1-\alpha _1-\alpha _1^d))$, and reads

(2.16) \begin{align} D_t \alpha_1 = \varLambda \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \quad \textrm{with} \ \varLambda = \alpha_1(1-\alpha_1)\tfrac{\rho_2 c_2^2\,-\,\rho_1c_1^2}{\alpha_1\rho_2c_2^2\,+\,(1-\alpha_1)\rho_1c_1^2}+\alpha_1^d \tfrac{\alpha_1\rho_2c_2^2}{\alpha_1\rho_2c_2^2\,+\,(1\,-\,\alpha_1)\rho_1c_1^2}. \end{align}

Moreover, system (2.14) is shown hyperbolic with a mixture's velocity $c_W^d:=c_W/(1-\alpha _1^d)$ where $c_W$ is the Wood sound velocity such that $(c_W^2)^{-1}=\rho (\alpha _1/(\rho _1c_1^2)+\alpha _2/(\rho _2c_2^2))$. This extends, in particular, the three-equation model proposed by Chanteperdrix et al. (Reference Chanteperdrix, Villedieu and Vila2002), which is recovered when $\alpha _1^d\rightarrow 0$.

Remark also that the small-scale specific energy $\zeta$ is chosen to be modelled here as a passively advected quantity. Therefore, it does not affect the dynamics of the two-scale mixture in absence of inter-scale mass transfer which is further modelled in § 5. The modelling of the effect of small-scale velocity fluctuations on large-scale dynamics, possibly via pressure contributions, is however out of the scope of this work. Because of this temporary passive role and for clarity purposes, this small-scale energy is removed from the modelling up to § 5.

2.3. Dissipative dynamics

Now that we have set the conservative structure of the model with Hamilton's SAP, we add dissipative processes to the model. We propose here to model both a pressure disequilibrium between the two large-scale phases as it is characteristic of two-phase flow models as in the model of Baer & Nunziato (Reference Baer and Nunziato1986) and fluid viscosity.

First, the pressure disequilibrium can encompass several physical phenomena with different timescales [e.g. bubbly flows (Drui et al. Reference Drui, Larat, Kokh and Massot2019) and stochastic thermodynamic relaxation (Perrier & Gutiérrez Reference Perrier and Gutiérrez2021)]. Let us mention that the limit of instantaneous relaxation can also be used to build a numerical scheme that solves the pressure equilibrium of (2.14) (Chanteperdrix et al. Reference Chanteperdrix, Villedieu and Vila2002).

In order to make this pressure relaxation a large-scale process that acts on each large-scale phase symmetrically, we define the large-scale volume fractions

(2.17)\begin{equation} \bar{\alpha}_k:=\frac{\alpha_k}{1-\alpha_1^d} \quad \textrm{for}\ k=1,2, \end{equation}

such that for $\alpha _1^d< 1$,

(2.18)\begin{equation} \alpha_1+\alpha_2+\alpha_1^d=1 \ \iff\ \bar{\alpha}_1 + \bar{\alpha}_2=1. \end{equation}

Then, we relax the pressure equilibrium of (2.14) into

(2.19)\begin{equation} D_t \bar{\alpha}_1 = \frac{p_1-p_2}{\epsilon}, \end{equation}

where $\epsilon >0$ has the dimension of a dynamic viscosity. Since the pressure is no longer unique with this pressure relaxation, the momentum equation can be expressed as follows:

(2.20)\begin{equation} \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+\bar{p} \boldsymbol{\mathsf{I}})=0, \end{equation}

where $\bar {p}$ is chosen to provide a signed dissipation of ${\mathcal {H}}$. Following Appendix B,

(2.21)\begin{align} \varsigma&:= \partial_t {\mathcal{H}} +\boldsymbol{\nabla}\boldsymbol{\cdot}(({\mathcal{H}}+\bar{p}) \boldsymbol{u}) \nonumber\\ &= (\kern0.7pt \bar{p}-\alpha_1p_1-(\alpha_2+\alpha_1^d)p_2)\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} - (p_1-p_2) D_t \alpha_1 \nonumber\\ &= (\kern0.7pt \bar{p}-\bar{\alpha}_1p_1-\bar{\alpha}_2p_2)\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} - \epsilon (1-\alpha_1^d) (D_t\bar{\alpha}_1)^2. \end{align}

Therefore, choosing $\bar {p}:=\bar {\alpha }_1p_1 + \bar {\alpha }_2p_2$ gives a signed mathematical entropy production $\varsigma \le 0$ and the relaxed model reads

(2.22)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k +\boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0\quad \textrm{for} \ k=1,2,1^d, \\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u})=0,\\ \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+ \bar{p}\boldsymbol{\mathsf{I}})=0,\\ D_t \bar{\alpha}_1= \epsilon^{{-}1}(p_1-p_2). \end{array}\right\} \end{equation}

The model (2.22) is obtained for the specific dissipative process (2.19) that gives a relaxed model for (2.14). It is also shown hyperbolic with a mixture's velocity $c_F^d:=c_F/(1-\alpha _1^d)$ where $c_F$ is the frozen sound velocity such that $c_F^2 = Y_1c_1^2+Y_2c_2^2$. The reader is referred to Loison et al. (Reference Loison, Kokh, Pichard and Massot2024) for a more thorough discussion on the pair of two-scale systems (2.14) and (2.22).

Second, the fluid viscosity in the context of a two-phase diffuse interface model with one velocity can be modelled with a Stokes source term on the momentum equation

(2.23)\begin{equation} \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+ \bar{p}\boldsymbol{\mathsf{I}})= \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{T}_{stokes}, \end{equation}

where $\boldsymbol {T}_{stokes}:=\mu (\boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u})^\textrm {T}) + \lambda (\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u})\boldsymbol {I}$ is the viscous stress, $\mu$ the dynamic viscosity and $\lambda$ the second viscosity with $\lambda + \frac {2}{3}\mu \ge 0$. After some calculations, this additional term modifies the mathematical entropy production (2.21) to yield

(2.24)\begin{align} &\partial_t {\mathcal{H}} +\boldsymbol{\nabla}\boldsymbol{\cdot}(({\mathcal{H}}+\bar{p}) \boldsymbol{u}-\boldsymbol{T}_{stokes}\boldsymbol{\cdot} \boldsymbol{u})\nonumber\\ &\quad ={-} \epsilon (1-\alpha_1^d) (D_t\bar{\alpha}_1)^2 - \!\left( 2\mu \left(\frac{\boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla}\boldsymbol{u})^T}{2}-\frac{1}{3}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\right)^2 + \left(\lambda + \frac{2}{3}\mu\right)(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})^2 \right)\!. \end{align}

The right-hand side is negatively signed and thus provides a dissipative term to the dynamics.

Including these two viscosity terms, a dimensionless analysis of the dissipative system is proposed in Appendix C. In the remainder of this paper, such viscous terms are neglected for readability purposes, but their absence could also correspond to regimes with high fluid velocities.

We remark that other dissipative processes could have been considered, but are not modelled here, such as drag for two-velocity models (Saurel et al. Reference Saurel, Chinnayya and Carmouze2017; Gavrilyuk Reference Gavrilyuk2020) or turbulent dissipation (Saurel, Gavrilyuk & Renaud Reference Saurel, Gavrilyuk and Renaud2003).

3.1. Kinetic-based model in the disperse regime

We start the modelling of the spray by assuming that:

(H3a) \begin{align} {(a)}\ &\textrm{the small-scale phase is composed of spherical droplets;} \end{align}

(H3b) \begin{align} {(b)}\ &\textrm{the droplets do not break nor coalesce.} \end{align}

Indeed, (H1b) shows that only small velocity differences are locally allowed at the small scale, and droplets are then very close to a spherical shape. The droplets are then described by a static minimal surface, which remains spherical. The kinetic modelling of the disperse small scale relies on the number density function (NDF) $f$ that counts the number of droplets within a small volume of the phase space around a point of the space–time domain. For generality purpose, we adopt a description for compressible inclusions whose sizes can change in time under either mass exchange or compressibility processes. Therefore, we consider the NDF $f(\boldsymbol {x}, t, \boldsymbol {v}, \hat {m})$ to be a distribution of velocities $\boldsymbol {v}$ and droplet masses $\hat {m}$ instead of a size parameter such as radius or surface. We also denote $f_{\boldsymbol {v},m}$ for compactness to underline the phase-space dependencies only. Then, the equation of evolution of the NDF expresses the conservation of the droplets in both the real and phase spaces (Williams Reference Williams1958; Marchisio & Fox Reference Marchisio and Fox2013) and reads

(3.1)\begin{equation} \partial_t\, f_{\boldsymbol{v},m} + \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot}(f_{\boldsymbol{v},m} \boldsymbol{v}) + \boldsymbol{\nabla}_{\boldsymbol{v}}\boldsymbol{\cdot}(\boldsymbol{F}(\boldsymbol{v}, \hat{m})f_{\boldsymbol{v},m}) + \partial_{\hat{m}}(R_{m}(\boldsymbol{v}, \hat{m})f_{\boldsymbol{v},m})=\varGamma, \end{equation}

where $\boldsymbol {F}$ and $R_{m}$ are unclosed rates of change of the velocity and the mass corresponding to body forces and mass transfer, while $\varGamma =\varGamma _{bu}+\varGamma _{coal}$ accounts for source terms such as break-up and coalescence phenomena. Hypotheses (H3b) and (H2c) imply $\varGamma =0$ and $R_{m}=0$, respectively. Moreover, the distribution of velocities is out of the scope of this work as assumed by (H1a) (see Vié et al. Reference Vié, Laurent and Massot2013 for an Eulerian model coupling the effects of size and velocity distributions). Therefore, we discard the velocity dependency by considering $n_m = \int _{\boldsymbol {v}} f_{\boldsymbol {v}, m} \,\textrm {d}\boldsymbol {v}$, the mass-based NDF in the limit of a vanishing Stokes number $St \rightarrow 0$, i.e. when the inertial timescale of inclusions is negligible compared to the timescale of the flow. Following the works of Jabin (Reference Jabin2002) and Massot (Reference Massot2007), the dynamics of $n_m$ reads

(3.2)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t n_{m} + \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot}(n_{m} \boldsymbol{u})=0,\\ \partial_t (n_{m}\boldsymbol{u}) + \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\cdot}(n_{m} \boldsymbol{u}\otimes\boldsymbol{u})=\boldsymbol{0}, \end{array}\right\} \end{equation}

where $\boldsymbol {u}$ is the unique velocity of every phase as defined by (2.6).

3.2. GeoMOM: definition and application to the incompressible spray

GeoMOM aims at enriching the models (2.14)–(2.22) with quantities describing the geometry of the interface at the small scale. This method introduced in Essadki et al. (Reference Essadki, de Chaisemartin, Massot, Laurent, Larat and Jay2016, Reference Essadki, de Chaisemartin, Laurent and Massot2018, Reference Essadki, Drui, de Chaisemartin, Larat, Ménard and Massot2019) is a reduced-order moment model for the polydisperse droplet distribution where the chosen moments are geometric quantities that can be defined for any flow regime and interface topology.

First, let us define those quantities. We consider a surface ${\mathcal {S}}$ defined by mapping a set ${\mathcal {U}}\subset \mathbb {R}^2$ onto ${\mathcal {S}}\subset \mathbb {R}^3$ such that we denote by $A(u,v)\,\textrm {d}u\,\textrm {d}v$ the infinitesimal surface element over ${\mathcal {S}}$. Then, the surface area and a surface-average operator are defined on ${\mathcal {S}}$ by

(3.3a,b)\begin{equation} S:= \int_{{\mathcal{U}}} A(u,v)\,{\rm d}u,{\rm d}v, \quad \langle {\cdot}\rangle:=\frac{1}{S}\int_{{\mathcal{U}}}({\cdot}) A(u,v)\,{\rm d}u\,{\rm d}v. \end{equation}

Considering the surface of the droplets in the two-scale mixture as described in § 2, the sum of their areas defines the small-scale interface area density $\varSigma$. Thanks to (3.3), one can also define surface-average Gauss and mean curvatures $\varSigma \langle G\rangle$ and $\varSigma \langle H\rangle$ as in Drew (Reference Drew1983) and Pope (Reference Pope1988). We use here these quantities to enhance the description of the disperse small scale, but their definitions are a priori independent of the surface geometry (see an application of these geometric quantities to the interface resulting from for the collision of two spherical droplets in Appendix D). This property is particularly convenient as it indicates that the kinetic-based model can be interpreted as a limit case for the mixed regime, and we discuss in conclusion how it allows coupling of the interface geometry between scales.

Second, we reduce the complexity of the spray dynamics, modelled by the distribution $n_m$ given in (3.2), by introducing a finite set of moments indexed by ${\mathcal {I}}$

(3.4)\begin{equation} M_i^m := \int_{\hat{m}} \hat{m}^i n_m \, {\rm d}\hat{m}, \quad i\in {\mathcal{I}} \textrm{ finite }\subset\mathbb{N}. \end{equation}

These scalars retain statistical information about the NDF and their dynamics can be obtained by integrating (3.2) against the corresponding monomials,

(3.5)\begin{equation} \partial_t M_i^m + \boldsymbol{\nabla}\boldsymbol{\cdot}(M_i^m\boldsymbol{u})=0, \quad i\in{\mathcal{I}}. \end{equation}

The method of moments then gives a reduced-order model for the small-scale spray in comparison with (3.2). In general, the equations are unclosed, and a reconstruction of the NDF based on the selected moments must be provided. The selection of these moments is not obvious and is usually motivated on the basis of mathematical properties of the resulting moment model.

Finally, the specificity of GeoMOM consists in the selection of moments related to geometric quantities defined with (3.3) to construct a model interpretable even out of the disperse regime. For the spherical droplets described by $n_m$, the local Gauss and mean curvatures on the sphere are constant and equal to $R^{-1}$ and $R^{-2}$, respectively, where $R=(3m/(4{\rm \pi} \rho _1^d))^{1/3}$ is the radius of the sphere of density $\rho _1^d$ and mass $m$. Therefore, considering a population of spherical droplets, we express the geometric quantities $\varSigma$, $\varSigma \langle G\rangle$ and $\varSigma \langle H\rangle$ as moments of the distribution $n_m$

(3.6a)$$\begin{gather} \varSigma\langle G\rangle = \int_{\hat{m}} 4{\rm \pi} n_m \, {\rm d}\hat{m}=4{\rm \pi} M^m_0, \end{gather}$$

(3.6b)$$\begin{gather}\varSigma\langle H\rangle = \int_{\hat{m}} 4{\rm \pi} \left(\frac{3}{4{\rm \pi} \rho_1}\right)^{1/3} \hat{m}^{1/3} n_m \, {\rm d}\hat{m}=4{\rm \pi} \left(\frac{3}{4{\rm \pi} \rho_1^d}\right)^{1/3} M^m_{1/3}, \end{gather}$$

(3.6c)$$\begin{gather}\varSigma = \int_{\hat{m}} 4{\rm \pi} \left(\frac{3}{4{\rm \pi} \rho_1^d}\right)^{2/3} \hat{m}^{2/3} n_m \, {\rm d}\hat{m}=4{\rm \pi} \left(\frac{3}{4{\rm \pi} \rho_1^d}\right)^{2/3} M^m_{2/3}. \end{gather}$$

Remark the special role of $\varSigma \langle G\rangle$ proportional to the zeroth-order moment of $n_m$. This results from the Gauss–Bonnet theorem (Kreyszig Reference Kreyszig1991), which indicates a geometric invariant $S\tilde {G}=4{\rm \pi}$ for continuous deformations of the sphere, where $\widetilde {({\cdot })}$ denotes the surface average as (3.3) applied to a unique droplet. It is also shown in Essadki et al. (Reference Essadki, Drui, de Chaisemartin, Larat, Ménard and Massot2019) that such an invariant enables to write the dynamics of geometric quantities $\varSigma$, $\varSigma \langle {\cdot }\rangle$ independently of the flow regime. Furthermore, the small-scale volume fraction is not surface-related in general, but in the specific case of spherical droplets, it is also linked to a moment of $n_m$ with

(3.7)\begin{equation} \alpha_1^d=\int_{\hat{m}} \frac{1}{\rho_1^d} \hat{m} n_{m} \, {\rm d}\hat{m} = \frac{1}{\rho_1^d} M^m_1. \end{equation}

From (3.6) and (3.7), we select ${\mathcal {I}}=\{0,1/3,2/3,1\}$ in (3.5) to obtain the dynamics of the geometric quantities

(3.8)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t (\varSigma\langle G\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u})=0,\\ \partial_t ((\rho_1^d)^{1/3}\varSigma\langle H\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot} ((\rho_1^d)^{1/3} \varSigma\langle H\rangle \boldsymbol{u})=0,\\ \partial_t ((\rho_1^d)^{2/3}\varSigma) + \boldsymbol{\nabla} \boldsymbol{\cdot} ((\rho_1^d)^{2/3}\varSigma \boldsymbol{u})=0,\\ \partial_t m_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_1^d \boldsymbol{u})=0. \end{array}\right\} \end{equation}

Up to here, the model has been derived for compressible spherical inclusions. We focus back on the description of a spray of droplets by assuming the incompressibility of the small scale thanks to the constraint $D_t \rho _1^d = 0$. Then, the geometric variables are governed by

(3.9)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t (\varSigma\langle G\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u}) = 0,\\ \partial_t (\varSigma\langle H\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle H\rangle \boldsymbol{u}) = 0,\\ \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u}) = 0,\\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0. \end{array}\right\} \end{equation}

System (3.9) corresponds to the system with surface-based moments of Essadki et al. (Reference Essadki, de Chaisemartin, Laurent and Massot2018) when no evaporation nor condensation is accounted for. Indeed, when the droplets are incompressible, the mass-based NDF $n_m$ relates to the surface-based NDF $n_{S}$ defined by

(3.10)\begin{equation} n_S(\boldsymbol{x}, t, \hat{S}):=n_m(\boldsymbol{x}, t, \hat{m})\delta(\hat{m}-m(\hat{S})), \quad m(\hat{S})=\rho_1^d\frac{\hat{S}}{3\sqrt{4{\rm \pi}}}. \end{equation}

As $\rho _1^d$ is a constant along the streamlines, there is no variation of surface area for the droplets and the dynamics is then driven by

(3.11)\begin{equation} \partial_t n_S + \boldsymbol{\nabla}_{\boldsymbol{x}}(n_S \boldsymbol{u})=0. \end{equation}

Again, when the droplets have the same density, the geometric quantities are also expressed through the moments $M_i^S$ of $n_S$,

(3.12a–d)\begin{equation} \varSigma\langle G\rangle = 4{\rm \pi} M^S_0, \quad \varSigma\langle H\rangle = \sqrt{4{\rm \pi}} M^S_{1/2}, \quad \varSigma = M^S_1, \quad \alpha_1^d = \frac{1}{3 \sqrt{4{\rm \pi}}} M^S_{3/2}. \end{equation}

The change of variables (3.10) only modifies the dimensions and orders of the moments. Consequently, the four initial half-integer moments can be recovered as in Essadki et al. (Reference Essadki, de Chaisemartin, Laurent and Massot2018). In the following, we stick with the surface-based moments for better comparison with the work of Essadki et al. In the end, the relations (3.12) between moments and geometric quantities provide a reduced-order model of the small-scale disperse phase (see Appendix D for an illustration of the reduced-order modelling applied to droplets resulting from the collision of two spherical droplets).

3.3. A two-scale model with a polydisperse small-scale model

The two-scale models derived in § 2 thanks to Hamilton's SAP can now be enhanced by incorporating the new geometric quantities in the Lagrangian and specifying whether these parameters are constrained by new conservation laws. We add then a small-scale capillarity energy $\sigma \varSigma$ along with the conservation constraints given by (3.9). Let us introduce a new variable $z=\varSigma /m_1^d$, which is transported, that is, $D_t z = 0$. As the large-scale capillarity is neglected (H2b), the energies associated with the large scales of the flow are left unchanged whereas the small-scale Lagrangian ${\mathcal {L}}_1^d$ becomes

(3.13)\begin{equation} {\mathcal{L}}_1^d(m_1^d, \rho_1^d, z:=\varSigma/m_1^d, \boldsymbol{u}):=\tfrac{1}{2}m_1^d \vert\boldsymbol{u}\vert^2 -m_1^d e_1(\rho_1^d) -\sigma m_1^d z. \end{equation}

Eventually, the system resulting from Hamilton's SAP (see details in Appendix E) together with the constraints (3.9) yields

(3.14)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0,\quad k=1,2,1^d, \\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) =0, \\ \partial_t \varSigma +\boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma\boldsymbol{u})=0,\\ \partial_t (\varSigma\langle G\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u}) =0, \\ \partial_t (\varSigma\langle H\rangle ) + \boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma\langle H\rangle \boldsymbol{u})=0,\\ \partial_t (\rho \boldsymbol{u})+ \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+p\boldsymbol{\mathsf{I}})=0, \end{array}\right\} \end{equation}

with $p:=p_1=p_2$. Let us remark that the new small-scale capillarity term in the Lagrangian does not alter momentum flux in this formulation compared to (2.14). Similarly to the previous models but with the updated definition of the Lagrangian ${\mathcal {L}}$, this model also admits a supplementary conservation equation

(3.15)\begin{equation} \partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot}(({\mathcal{H}}+p)\boldsymbol{u})=0, \quad {\mathcal{H}} = \sum_k m_k \Vert\boldsymbol{u}\Vert^2 - {\mathcal{L}}. \end{equation}

Both (3.14) and its relaxed version are hyperbolic with similar structure of eigenvalues and eigenvectors.

4.1. Kinetic-based model in the disperse regime

Let us model the spray of oscillating droplets by first discarding assumption (H3a) of spherical shapes, and rather assume:

(H4a) \begin{align} {(a)}\ &\textrm{the droplets' internal flow is irrotational; } \end{align}

(H4b) \begin{align} {(b)}\ &\textrm{the amplitude of the interface deformation is small; } \end{align}

(H4c) \begin{align} {(c)}\ &\textrm{the deformed droplets are either randomly orientated or they all } \nonumber \\ & \textrm{orientate to the same privileged direction;} \end{align}

(H4d) \begin{align} {(d)}\ &\textrm{the droplets' interface is deformed along the second } \nonumber \\ & \textrm{axisymmetric spherical harmonic.} \end{align}

The dynamics of a droplet satisfying (H4a) and (H4b) has been studied thoroughly (Rayleigh Reference Rayleigh1879; Prosperetti Reference Prosperetti1977; Plümacher et al. Reference Plümacher, Oberlack, Wang and Smuda2020). The choice of orientation (H4c) of the deformation depends on the physics, whether it is isotropic or not, and it is not discussed here. Denote by $\boldsymbol {r}$ the position of the droplet interface in spherical coordinates $(r,\theta,\phi )$ and $(\boldsymbol {e}_r$, $\boldsymbol {e}_{\theta }$, $\boldsymbol {e}_{\phi })$ the spherical orthonormalised basis. Following (H4d), we denote by $Y_2(\theta ) = \sqrt {5}/(4\sqrt {{\rm \pi} })(3\cos ^2\theta -1)$ the second axisymmetric spherical harmonic, and $R_0$ the radius of the non-deformed spherical droplet, such that the position can be written as

(4.1)\begin{equation} \boldsymbol{r} = (R_0 + x_2 Y_2) \boldsymbol{e}_r, \end{equation}

where $x_2$ denotes the amplitude of the deformation. The dynamics of this motion follows the harmonic oscillation

(4.2)\begin{equation} \ddot{x}_2+\omega^2 x_2 = 0,\quad \omega^2 = 8\frac{\sigma}{\rho_1^d R_0^3}=\tilde{\omega}^2S_0^{{-}3/2}, \end{equation}

where $\tilde {\omega }^2=(8(4{\rm \pi} )^{3/2}\sigma /\rho _1^d)$ is a constant along the streamlines. The dynamics of the droplets is then characterised by $(S_0,x_2,\dot {x}_2)$. For computational reasons, we introduce $\chi = (2/S_0)^{1/2}x_2$ which satisfies the same dynamics (4.2) as $x_2$. Assuming (H3b) and (H2c), the balance equation for the NDF $n_{\xi }$ in the phase space $\hat {\boldsymbol {\xi }}=(\hat {S}_0, \hat {\chi }, \hat {\dot {\chi }})$ yields

(4.3)\begin{equation} \partial_t n_{\xi}+ \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} (n_{\xi}\boldsymbol{u}) + \partial_{\hat{\chi}} (R_{\chi}n_{\xi}) + \partial_{\hat{\dot{\chi}}}(R_{\dot{\chi}}n_{\xi})=0, \end{equation}

where $R_{\chi }$ and $R_{\dot {\chi }}$ are rates of change that remain to be closed. This model can be seen as the particular case satisfying (H1a)–(H2c)–(H3b)–(H2a) of the kinetic model considered by O'Rourke & Amsden (Reference O'Rourke and Amsden1987), Amsden et al. (Reference Amsden, O'Rourke and Butler1989). With the oscillator model (4.2), we close (4.3) by fixing $R_{\chi } = \hat {\dot {\chi }}$ and $R_{\hat {\dot {\chi }}} = -\tilde {\omega }^2\hat {S}_0^{-3/2}\hat {\chi }$,

(4.4)\begin{equation} \partial_t n_{\xi}+ \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} (n_{\xi}\boldsymbol{u}) + \partial_{\hat{\chi}} (\hat{\dot{\chi}} n_{\xi}) + \partial_{\hat{\dot{\chi}}}(-\tilde{\omega}^2\hat{S}_0^{{-}3/2}\hat{\chi} n_{\xi})=0. \end{equation}

4.2. Modelling a spray of asynchronously oscillating droplets

With the new PBE (4.4), GeoMOM leads to a non-trivial closure problem and the geometric quantities chosen in (3.6) and (3.7) cannot model an arbitrary distribution. With such a choice, we propose in this section to approximate the distribution with

(4.5)\begin{equation} n_{\xi}(\boldsymbol{x}, t, \hat{S}_0, \hat{\chi}, \hat{\dot{\chi}}) = n_{S_0}(\boldsymbol{x}, t, \hat{S}_0) \frac{1}{\vert{\mathcal{E}}(\hat{S}_0)\vert} \mathbb{1}_{{\mathcal{E}}(\hat{S}_0)}(\hat{\chi}, \hat{\dot{\chi}}), \end{equation}

where the amplitudes and their rates of change $(\hat {\chi }, \hat {\dot {\chi }})$ are uniformly distributed on a compact space ${\mathcal {E}}\subset \mathbb {R}^2$ of area $\vert {\mathcal {E}}(\hat {S}_0)\vert$ that allows a maximal energy for a given droplet size $\hat {S}_0$. This approximation of the NDF discards situations where droplets oscillate synchronously, and consequently a macroscopic oscillation motion of the spray.

4.2.1. GeoMOM based on the classical surface-average operator

We apply GeoMOM with the same surface-average operators $\langle {\cdot }\rangle$ and $\widetilde {({\cdot })}$ defined in § 3.2, and geometric quantities $\varSigma \langle G\rangle$, $\varSigma \langle H\rangle$, $\varSigma$, $\alpha _1^d$. We establish then the relations between these geometric quantities and moments of the NDF in the context of a disperse regime,

(4.6)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \varSigma\langle G\rangle = \int_{\hat{\boldsymbol{\xi}}} S\tilde{G} n_{\xi}\,{\rm d}\hat{\boldsymbol{\xi}}, \quad \varSigma\langle H\rangle = \int_{\hat{\boldsymbol{\xi}}} S\tilde{H} n_{\xi}\, {\rm d}\hat{\boldsymbol{\xi}}, \\[14pt] \displaystyle \varSigma = \int_{\hat{\boldsymbol{\xi}}} S n_{\xi}\, {\rm d}\hat{\boldsymbol{\xi}}, \quad \alpha_1^d = \int_{\hat{\boldsymbol{\xi}}} V n_{\xi}\, {\rm d}\hat{\boldsymbol{\xi}}. \end{array}\right\} \end{equation}

When the oscillations are small, the integrands are functions of $\hat {\boldsymbol {\xi }}$ (see Appendix F) and write

(4.7a–d)\begin{equation} S\tilde{G}=4{\rm \pi}, \quad S\tilde{H} - S_0\tilde{H}_0 = \sqrt{4{\rm \pi}} S_0^{1/2}\chi^2, \quad V = \frac{1}{3\sqrt{4{\rm \pi}}} S_0^{3/2}, \quad S - S_0 = S_0 \chi^2, \end{equation}

where $\tilde {H}_0=\sqrt {4{\rm \pi} }S_0^{-1/2}$ is the surface-averaged mean curvature when the droplet is a sphere, i.e. $\chi =0$. We remark that $S\tilde {G}$ and $\alpha _1^d$ are constant despite the oscillation thanks to the Gauss–Bonnet theorem (Kreyszig Reference Kreyszig1991) and the incompressibility assumption. Finally, relations between geometric quantities and moments $M_{i,j,k}^{\xi }:=\int _{\hat {\boldsymbol {\xi }}} \hat {S}^i \hat {\chi }^j \hat {\dot {\chi }}^k n_{\xi }\, \textrm {d}\hat {\boldsymbol {\xi }}$ of $n_{\xi }$ are obtained

(4.8)\begin{equation} \left.\begin{array}{@{}c@{}} \varSigma\langle G\rangle = 4{\rm \pi} M_{0,0,0}^{\xi},\quad \varSigma\langle H\rangle = \sqrt{4{\rm \pi}} (M_{1/2, 0, 0}^{\xi} + M_{1/2,2,0}^{\xi}),\\[12pt] \varSigma = M_{1,0,0}^{\xi}+M_{1,2,0}^{\xi},\quad \alpha_1^d = \dfrac{1}{3\sqrt{4{\rm \pi}}}M_{3/2,0,0}^{\xi}. \end{array}\right\} \end{equation}

These relations extend those of (3.12) with moments dedicated to the oscillatory dynamics. Such decomposition leads us to define and choose the following geometric quantities and moments for our model

(4.9)\begin{equation} \left.\begin{array}{@{}c@{}} \varSigma\langle H\rangle_0 := \sqrt{4{\rm \pi}} M_{1/2, 0, 0}^{\xi},\quad \Delta\varSigma\langle H\rangle := \sqrt{4{\rm \pi}} M_{1/2,2,0}^{\xi},\\ \varSigma_0 := M_{1, 0, 0}^{\xi},\quad \Delta\varSigma := M_{1,2,0}^{\xi}, \end{array}\right\} \end{equation}

instead of just $\varSigma$ and $\varSigma \langle H\rangle$ which can be reconstructed with (4.8). Integrating (4.4) against $(1,\hat {S}_0^{1/2},\hat {S}_0,\hat {S}^{3/2}_0,\hat {S}^{1/2}_0\hat {\chi }^2,\hat {S}_0\hat {\chi }^2)$ provides

(4.10)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t (\varSigma\langle G\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u}) = 0,\\ \partial_t (\varSigma\langle H\rangle_0) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle H\rangle_0 \boldsymbol{u}) = 0,\\ \partial_t \varSigma_0 + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma_0 \boldsymbol{u}) = 0,\\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0,\\ \partial_t (\Delta\varSigma\langle H\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\Delta\varSigma\langle H\rangle \boldsymbol{u}) = 2 \sqrt{4{\rm \pi}} M_{1/2,1,1}^{\xi},\\ \partial_t (\Delta\varSigma) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\Delta\varSigma \boldsymbol{u}) = 2 M_{1,1,1}^{\xi}. \end{array}\right\} \end{equation}

We obtained that $\varSigma \langle G\rangle$, $\varSigma \langle H\rangle _0$, $\varSigma _0$ and $\alpha _1^d$ are conserved similarly as (3.9) with two additional equations for the oscillatory components $\Delta \varSigma \langle H\rangle$ and $\Delta \varSigma$.

4.2.2. Energies of the spray

In the context of two-scale modelling with Hamilton's SAP, we are specifically interested in defining the energies of the spray with the geometric quantities. For the oscillatory motion described by (4.2), the kinetic and potential energies of a droplet can be expressed as function of $\boldsymbol {\xi }$ (Rayleigh Reference Rayleigh1879, Appendix II)

(4.11a,b)\begin{equation} E^{kin,1d} = \frac{1}{2} \frac{\rho_1^d}{4(4{\rm \pi})^{3/2}}S_0^{5/2}\dot{\chi}^2, \quad E^{pot,1d} = \sigma S = \sigma S_0 + \sigma S_0 \chi^2. \end{equation}

It is then straightforward to obtain the energies of the spray $E^{kin,d}$ and $E^{pot,d}$ from moments of $n_{\xi }$ by integrating the expressions above,

(4.12a,b)\begin{equation} E^{kin,d} = \frac{1}{2} \frac{\rho_1^d}{4(4{\rm \pi})^{3/2}}M_{5/2,0,2}^{\xi}, \quad E^{pot,d} = \sigma \varSigma = \sigma M^{\xi}_{1,0,0} + \sigma M^{\xi}_{1,2,0}. \end{equation}

Remark that the moment $M_{5/2,0,2}^{\xi }$ is not linked to any of the selected geometric quantities of the model. We can find a closure for this moment using the approximation of $n^{\xi }$ given in (4.5). It requires to provide a definition for ${\mathcal {E}}$ which is chosen following the break-up criterion of O'Rourke & Amsden (Reference O'Rourke and Amsden1987). Therefore, we authorise the droplets to oscillate with an energy lower than a fraction $c\in [0,1]$ of the maximal energy $E_{max}$ before break-up,

(4.13)\begin{equation} {\mathcal{E}}(\hat{S}_0):=\{ (\hat{\chi},\hat{\dot{\chi}})\in\mathbb{R}^2 \textrm{ s.t. } E^{kin,1d}(\hat{S}_0, \hat{\chi}, \hat{\dot{\chi}}) + E^{pot,1d}(\hat{S}_0, \hat{\chi}, \hat{\dot{\chi}})\le c E_{max}(\hat{S}_0) \}. \end{equation}

The break-up energy $E_{max}(S_0)$ corresponds to a deformation up to an equatorial radius reaching half the value of the spherical radius

(4.14)\begin{equation} R_{eq} := \Vert\boldsymbol{r}({\rm \pi}/2)\Vert = \frac{R_0}{2}. \end{equation}

Then, the subset ${\mathcal {E}}$ boils down to ${\mathcal {E}}=\{\chi ^2 + (\dot {\chi }/\omega )^2 \le 2c/5\}$ whose area is $\vert {\mathcal {E}}\vert = \frac {2}{5} {\rm \pi}\omega c$. We can now close the expression of $E^{kin,d}$, $E^{pot,d}$ in (4.11) and find the dynamics of the parameter $c$ using (4.10) and the approximation of the NDF (4.5). It yields

(4.15a–c)\begin{equation} E^{kin,d} = \tfrac{1}{10} c \sigma \varSigma_0, \quad E^{pot,d} = \sigma\varSigma_0 + \tfrac{1}{10} c \sigma \varSigma_0,\quad D_t c = 0. \end{equation}

We remark that, as $(\chi,\dot {\chi })$ are uniformly distributed in ${\mathcal {E}}$, the mechanical energy is evenly distributed between kinetic and potential energies. Regarding the dynamics of $c$, the energy of the oscillation is advected along the streamline. Moreover, using the relation $\Delta \varSigma = M^{\xi }_{1,2,0}$, one can replace $c$ by geometric quantities with $c=10\Delta \varSigma /\varSigma _0$.

4.2.3. Two-scale model with the small-scale spray model of asynchronous droplets

Energies related to the small-scale oscillation are negatively signed in Hamilton's SAP as the kinetic energy is not a quadratic form of any kind of velocity. Indeed, the energy of the droplets is here considered as a whole energetic contribution without following their specific dynamics. Remark that the specific oscillatory energies for the small-scale are advected by the flow and then partly accounts for the velocity fluctuations introduced in § 2. While not being accounted here, the fluctuation energy term would be defined as the energy contribution of the velocity fluctuations which is not associated with droplets oscillations. We keep then the expression of the Lagrangian defined in (2.13), where we only modify the small-scale energies of ${\mathcal {L}}_1^d$ following

(4.16)\begin{equation} {\mathcal{L}}_1^d := \frac{1}{2}m_1^d \Vert\boldsymbol{u}\Vert^2 - m_1^d e_1(\rho_1^d) - \sigma \varSigma_0 \left(1 + \frac{c}{5}\right). \end{equation}

The constraints are the same as in § 3 with the additional advection constraint $D_t c=0$. The derivation of the two-scale mixture's dynamics with Hamilton's SAP is very similar to that of Appendix E and yields

(4.17)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0,\quad k=1,2,1^d, \\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) =0, \\ \partial_t \varSigma_0 +\boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma_0\boldsymbol{u})=0,\\ \partial_t (\Delta\varSigma) + \boldsymbol{\nabla} \boldsymbol{\cdot}(\Delta\varSigma\boldsymbol{u})=0,\\ \partial_t (\varSigma\langle G\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u}) =0, \\ \partial_t (\varSigma\langle H\rangle _0) + \boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma\langle H\rangle _0\boldsymbol{u})=0,\\ \partial_t (\rho \boldsymbol{u})+ \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+p\boldsymbol{\mathsf{I}})=0, \end{array}\right\} \end{equation}

with $p:=p_1=p_2$, and the additional equation on total energy as a mathematical entropy for (4.17),

(4.18)\begin{equation} \partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (({\mathcal{H}}+p)\boldsymbol{u})=0, \quad {\mathcal{H}} = \sum_k m_k \Vert\boldsymbol{u}\Vert^2-{\mathcal{L}}. \end{equation}

Now considering dissipation processes, we could also add the pressure relaxation as in the previous two-scale models. We propose here to focus on the dissipation associated to the small-scale oscillatory dynamics. Indeed, the oscillation motion eventually decreases and the associated energy dissipates into thermal energy or small-scale kinetic energy of the gas phase. Each of these last two energies are not modelled here, so the system loses this energy, and we only consider the following source term $R_c$ on the dynamics of $c$,

(4.19)\begin{equation} D_t c = R_c. \end{equation}

The mathematical entropy production may be written as

(4.20)\begin{equation} \varsigma:=\partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (({\mathcal{H}}+p)\boldsymbol{u}) ={-} (\partial_c {\mathcal{L}}) R_c. \end{equation}

As $\partial _c {\mathcal {L}} = - \sigma \varSigma _0/5$, we choose $R_c = - c/\tau$ to sign the mathematical entropy production $\varsigma \le 0$ and model the dissipation with an exponential decrease of characteristic time $\tau >0$. The final system reads

(4.21)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})= 0,\quad k=1,2,1^d, \\ \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0,\\ \partial_t \varSigma_0 +\boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma_0\boldsymbol{u})= 0,\\ \partial_t (\Delta\varSigma) + \boldsymbol{\nabla} \boldsymbol{\cdot}(\Delta\varSigma\boldsymbol{u})={-}\tau^{{-}1}\Delta\varSigma,\\ \partial_t (\varSigma\langle G\rangle) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma\langle G\rangle \boldsymbol{u}) = 0, \\ \partial_t (\varSigma\langle H\rangle _0) + \boldsymbol{\nabla} \boldsymbol{\cdot}(\varSigma\langle H\rangle _0\boldsymbol{u})= 0,\\ \partial_t (\rho \boldsymbol{u})+ \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+p \boldsymbol{\mathsf{I}})= 0, \end{array}\right\} \end{equation}

with $p:=p_1=p_2$. One can see that the component of the interface area density modelling the droplets’ dynamics have now its own dynamics with an additional source term that makes it dissipate over time.

4.3. Modelling a spray of synchronously oscillating droplets

With the choice of geometric quantities of the previous section, the moments are not suited to obtain a macroscopic oscillation of the spray. Therefore, we introduce now new geometric quantities with GeoMOM based on another surface-average operator, and the following NDF approximation

(4.22)\begin{equation} n_{\xi}(\hat{S}_0,\hat{\chi}, \hat{\dot{\chi}}) = \sum_{i} n_i \delta(\hat{S}_0 - (S_0)_i)\delta(\hat{\chi}-\chi_i)\delta(\hat{\dot{\chi}}-\dot{\chi}_i). \end{equation}

This corresponds to several populations of droplets which share the same size and oscillate synchronously.

4.3.1. GeoMOM based on an oriented surface-average operator

Similarly to § 3.2, we consider a surface ${\mathcal {S}}$ and its mapping ${\mathcal {U}}\subset \mathbb {R}^2$ onto ${\mathcal {S}}\subset \mathbb {R}^3$ such that $A(u,v)\, \textrm {d}u\,\textrm {d}v$ is the infinitesimal surface element over ${\mathcal {S}}$. We decompose the surface local area into two contributions, one related to a preferred direction $\boldsymbol {N}$ given by the large-scale dynamics as suggested by (H4c). In this sense, we decompose this infinitesimal surface element $A(u,v)\, \textrm {d}u\,\textrm {d}v$ using its definition using tangential vectors $\boldsymbol {e}_u$ and $\boldsymbol {e}_v$

(4.23a–c)\begin{equation} A = \Vert \boldsymbol{e}_u \times \boldsymbol{e}_v \Vert, \quad A_{{\parallel}} = | (\boldsymbol{e}_u \times \boldsymbol{e}_v) \boldsymbol{\cdot} \boldsymbol{N} |, \quad A_{{\perp}} = A - A_{{\parallel}}. \end{equation}

From them, we decompose the surface into parallel and perpendicular components $S_{\parallel }$ and $S_{\perp }$ such that

(4.24)\begin{equation} S = \int_{{\mathcal{U}}}A(u,v)\, {\rm d}u\,{\rm d}v = \int_{{\mathcal{U}}}A_{{\parallel}}(u,v)\, {\rm d}u\,{\rm d}v + \int_{{\mathcal{U}}}A_{{\perp}}(u,v)\, {\rm d}u\,{\rm d}v = S_{{\parallel}} + S_{{\perp}}. \end{equation}

We also define new surface-average operators similarly to (3.3)

(4.25a,b)\begin{equation} \langle {\cdot}\rangle_{{\parallel}}:=\frac{1}{S_{{\parallel}}}\int_{{\mathcal{U}}}({\cdot}) A_{{\parallel}}(u,v)\, {\rm d}u\,{\rm d}v, \quad \langle {\cdot}\rangle_{{\perp}}:=\frac{1}{S_{{\perp}}}\int_{{\mathcal{U}}}({\cdot}) A_{{\perp}}(u,v)\, {\rm d}u\,{\rm d}v, \end{equation}

and the surface-average operators $\widetilde {({\cdot })}^{\parallel }$, $\widetilde {({\cdot })}^{\perp }$ for one closed inclusion. Considering now the oscillatory motion (4.2), we have the following dynamics for the geometric quantities of one droplet (see details in Appendix F)

(4.26)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle S_{{\parallel}} = \dfrac{1}{2}S_0 - \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}}S_0\chi, \quad S_{{\parallel}}\tilde{H}^{{\parallel}} = \dfrac{1}{2}\sqrt{4{\rm \pi}}\sqrt{S_0} - \dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}}\sqrt{S_0}\chi,\\ \displaystyle S_{{\perp}} = \dfrac{1}{2}S_0 + \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}}S_0\chi, \quad S_{{\perp}}\tilde{H}^{{\perp}} = \dfrac{1}{2}\sqrt{4{\rm \pi}}\sqrt{S_0} + \dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}}\sqrt{S_0}\chi. \end{array}\right\} \end{equation}

Defining the oriented interface area densities $\varSigma _{\perp }$ and $\varSigma _{\parallel }$ within the mixture, we remark that $\varSigma = \varSigma _{\perp } + \varSigma _{\parallel }$ and $\varSigma \langle {\cdot }\rangle = \varSigma _{\perp }\langle {\cdot }\rangle _{\perp } + \varSigma _{\parallel }\langle {\cdot }\rangle _{\parallel }$. Integrating (4.26) against $n_{\xi }$ provides

(4.27)\begin{align} \left.\begin{array}{@{}c@{}} \displaystyle \varSigma_{{\parallel}} = \dfrac{1}{2}M^{\xi}_{1,0,0} - \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1,1,0}, \quad (\varSigma\langle H\rangle )_{{\parallel}} = \dfrac{1}{2}\sqrt{4{\rm \pi}}M^{\xi}_{1/2,0,0}-\dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1/2,1,0},\\ \displaystyle \varSigma_{{\perp}} = \dfrac{1}{2}M^{\xi}_{1,0,0} + \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1,1,0}, \quad (\varSigma\langle H\rangle )_{{\perp}} = \dfrac{1}{2}\sqrt{4{\rm \pi}}M^{\xi}_{1/2,0,0}+\dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1/2,1,0}. \end{array}\right\} \end{align}

Remark that both the parallel and perpendicular geometric variables represent the same information for the spray as they are related to the same moments of the NDF. Consequently, we here retain the perpendicular components only, and we extract the part dedicated to the dynamics by defining

(4.28)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \Sigma_{{\perp},0}:= \tfrac{1}{2}M^{\xi}_{1,0,0},\quad (\Delta\varSigma_{{\perp}}):= \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1,1,0}, \\ \displaystyle (\varSigma\langle H\rangle)_{{\perp},0}:=\tfrac{1}{2}\sqrt{4{\rm \pi}}M^{\xi}_{1/2,0,0},\quad (\Delta\varSigma\langle H\rangle)_{{\perp}}:=\dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}}M^{\xi}_{1/2,1,0}. \end{array}\right\} \end{equation}

We intentionally factored some indexes in the two last geometric quantities to lighten the notations. This splitting allows us to retain more information on the NDF in our model by including two moments rather than the sum of two moments. We know from (4.10) that $\varSigma _{\perp,0}$ and $(\varSigma \langle H\rangle )_{\perp,0}$ are conserved, while we integrate (4.4) to obtain the dynamics of the two other geometric quantities

(4.29)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \partial_t (\Delta\varSigma_{{\perp}}) + \boldsymbol{\nabla} \boldsymbol{\cdot} ( (\Delta\varSigma_{{\perp}}) \boldsymbol{u}) = \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}} M^{\xi}_{1,0,1},\\ \displaystyle \partial_t (\Delta\varSigma\langle H\rangle)_{{\perp}} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( (\Delta\varSigma\langle H\rangle)_{{\perp}} \boldsymbol{u}) = \dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}} M^{\xi}_{1/2,0,1}. \end{array}\right\} \end{equation}

With the moments considered up to now, we cannot model a macroscopic oscillatory spray as we lack some information regarding the distribution of $\hat {\dot {\chi }}$ in our model. We propose then to add two new geometric quantities that correspond to the unclosed moments of (4.29)

(4.30)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle(\varDelta_t\Sigma_{{\perp}}) := \varSigma_{{\perp}}\left\langle \dfrac{\partial_t A_{{\perp}}}{A_{{\perp}}}\right\rangle_{{\perp}} = \dfrac{1}{\sqrt{4{\rm \pi}}} \sqrt{\dfrac{5{\rm \pi}}{2}} M^{\xi}_{1,0,1},\\ \displaystyle (\varDelta_t\varSigma\langle H\rangle)_{{\perp}} := \varSigma_{{\perp}}\left\langle \partial_t H + H \dfrac{\partial_t A_{{\perp}}}{A_{{\perp}}}\right\rangle _{{\perp}} = \dfrac{1}{2}\sqrt{\dfrac{5{\rm \pi}}{2}} M^{\xi}_{1/2,0,1}. \end{array}\right\} \end{equation}

The above definitions with the oriented surface-average operators show that these quantities are also well-defined regardless of the flow regime. Finally, we consider the following eight geometric quantities to describe the spray of oscillating droplets: $\alpha _1^d$, $\Sigma _{\perp,0}$, $(\Delta \varSigma _{\perp })$, $(\varDelta _t\Sigma _{\perp })$, $(\varSigma \langle H\rangle )_{\perp,0}$, $(\Delta \varSigma \langle H\rangle )_{\perp }$, $(\varDelta _t\varSigma \langle H\rangle )_{\perp }$ and$\varSigma \langle G\rangle$.

4.3.2. Amplitude-based closure

With the eight moments given by the eight corresponding geometric quantities, we propose to look for a two-point quadrature which corresponds to two populations of droplets as defined in (4.22),

(4.31)\begin{equation} n_{\xi}(\hat{S}_0,\hat{\chi}, \hat{\dot{\chi}}) = \sum_{i=1,2} n_i \delta(\hat{S}_0 - (S_0)_i)\delta(\hat{\chi}-\chi_i)\delta(\hat{\dot{\chi}}-\dot{\chi}_i), \end{equation}

where $n_i$ are the weights or numbers of droplets that share the same abscissas $\chi _i$ and $\dot {\chi }_i$. The quadrature above admits a unique solution (excluding symmetry) under some realisability conditions that ensure that $\varSigma \langle G\rangle$, $(\varSigma \langle H\rangle )_{\perp,0}$, $\varSigma _{\perp,0}$ and $\alpha _1^d$ are positive with additional geometric constraints (see Appendix G). Remark that $n_i$ and $(S_0)_i$ depend only on $\varSigma \langle G\rangle$, $(\varSigma \langle H\rangle )_{\perp,0}$, $\Sigma _{\perp,0}$, $\alpha _1^d$ and, for $(i,j)\in \{(1,2),(2,1)\}$,

(4.32)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle\chi_i =\sqrt{\dfrac{2}{5{\rm \pi}}}\dfrac{\sqrt{4{\rm \pi}}(\Delta\varSigma_{{\perp}})-2\sqrt{(S_{0})_j}(\Delta\varSigma\langle H\rangle)_{{\perp}}}{n_i((S_{0})_i-\sqrt{(S_{0})_i(S_{0})_j})},\\ \displaystyle\dot{\chi}_i =\sqrt{\dfrac{2}{5{\rm \pi}}}\dfrac{\sqrt{4{\rm \pi}}(\varDelta_t\Sigma_{{\perp}})-2\sqrt{(S_{0})_j}(\varDelta_t\varSigma\langle H\rangle)_{{\perp}}}{n_i((S_{0})_i-\sqrt{(S_{0})_i(S_{0})_j})}. \end{array}\right\} \end{equation}

We easily obtain that $\alpha _1^d$, $\varSigma \langle G\rangle$, $(\varSigma \langle H\rangle )_{\perp,0}$ and $\Sigma _{\perp,0}$ are conserved such that, together with (4.29), it implies

(4.33)\begin{equation} \partial_t n_i + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\boldsymbol{u})= 0,\quad D_t (S_{0})_i = 0, \quad D_t\chi_i = \dot{\chi}_i, \quad i=1, 2. \end{equation}

This dynamics is expected as the two populations of oscillating droplets are advected by the flow at velocity $\boldsymbol {u}$.

4.3.3. Two-scale model with the small-scale spray model of synchronous droplets

Denote $\nu := \rho _1^d/(4(4{\rm \pi} )^{3/2})$ and $\gamma := 2\sigma$ such the kinetic and potential energies of the spray of oscillating droplets given in (4.11) can be written

(4.34a,b)\begin{equation} E^{kin,d} = \tfrac{1}{2} \nu M^{\xi}_{5/2,0,2}, \quad E^{pot,d} = \sigma M^{\xi}_{1,0,0} + \tfrac{1}{2}\gamma M^{\xi}_{1,2,0}, \end{equation}

and $\tilde {\omega }^2 = \gamma /\nu$. The closure (4.31) then yields

(4.35a,b)\begin{equation} E^{kin,d} = \sum_{i=1,2} \frac{1}{2}\nu n_i (S_0)_i^{5/2}\dot{\chi}_i^2,\quad E^{pot,d} = \sum_{i=1,2} \sigma n_i (S_0)_i + \frac{1}{2}\gamma n_i (S_0)_i \chi_i^2. \end{equation}

We extend the two-scale Lagrangian (2.13) by adding the energies above to account for the small-scale oscillation

(4.36)\begin{align} {\mathcal{L}} &= {\mathcal{L}}_1(\alpha_1, m_1, \boldsymbol{u}) +{\mathcal{L}}_2(\alpha_2, m_2, \boldsymbol{u}) +{\mathcal{L}}_{1}^d(m_1^d, \rho_1^d, \boldsymbol{u})\nonumber\\ &\quad +{\mathcal{L}}^{vib}_1(n_1, (S_0)_1, \chi_1, \dot{\chi}_1) +{\mathcal{L}}^{vib}_2(n_2, (S_0)_2, \chi_2, \dot{\chi}_1), \end{align}

where ${\mathcal {L}}^{vib}_i$ is defined using vibrating energies of (4.34)

(4.37)\begin{equation} {\mathcal{L}}^{vib}_i = \tfrac{1}{2}\nu n_i (S_0)_i^{5/2}\dot{\chi}_i^2 - \sigma n_i (S_0)_i - \tfrac{1}{2}\gamma n_i (S_0)_i \chi_i^2. \end{equation}

Recall that the kinetic energy is here positively signed as it is a quadratic form of $\dot {\chi }_1$ and $\dot {\chi }_2$ with their associated momentum equations after Hamilton's SAP. We also remark that similarly to the asynchronous case, a similar remark as that stated in § 4.2.3 holds regarding the modelling of small-scale energy fluctuations. We provided an expression using the quadrature's abscissas rather than the geometric quantities for computational convenience, but the dynamics of the geometric quantities is equivalently obtained using the quadrature expression of Appendix G. In addition to $\alpha _1$, the quantities $\chi _1$ and $\chi _2$ are also free variables which results in the following system with two additional momentum equations for each population of droplets of same size in the spray (see Appendix H):

(4.38)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0, \quad k=1,2,1^d,\\ \partial_t n_i + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i \boldsymbol{u})=0 ,\quad i=1,2,\\ \partial_t (n_i(S_0)_i) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i(S_0)_i \boldsymbol{u}) =0, \quad i=1,2,\\ \partial_t(n_i\dot{\chi}_i)+ \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\dot{\chi}_i \boldsymbol{u})={-}\omega_i^2 n_i\chi_i, \quad i=1,2, \\ \partial_t(n_i\chi_i)+ \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\chi_i \boldsymbol{u}) =n_i\dot{\chi}_i, \quad i=1,2,\\ \partial_t (\rho \boldsymbol{u})+ \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+p)=0, \end{array}\right\} \end{equation}

with $p:=p_1=p_2$ and $\omega _i^2 = \tilde {\omega }^2 (S_0)_i^{-3/2}$. Combining the equations at the fifth and sixth lines of the system above, one can recognise the equations of harmonic oscillators advected along the streamlines

(4.39)\begin{equation} D_t(D_t \chi_i) + \omega_i^2 \chi_i = 0,\quad i=1,2. \end{equation}

Similarly to the previous models, this system admits an additional conservation equation on the total energy ${\mathcal {H}}$ defined hereafter

(4.40)\begin{align} \partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (({\mathcal{H}}+p)\boldsymbol{u})=0,\quad {\mathcal{H}} = \sum_{k=1,2,1^d} \frac{1}{2}m_k\vert\boldsymbol{u}\vert^2 + \sum_{i=1,2} \frac{1}{2}\nu n_i (S_0)_i^{5/2}\dot{\chi}_i^2 -{\mathcal{L}}. \end{align}

Once again, we could consider a pressure relaxation model, but we would like to focus again on the dissipation associated to the oscillation process, and we introduce source terms $R_{\chi _i}$ in the new momentum equations

(4.41)\begin{equation} \partial_t(n_i\dot{\chi}_i) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\dot{\chi}_i \boldsymbol{u}) ={-}\omega_i^2 n_i\chi_i + R_{\chi_i},\quad i=1,2. \end{equation}

This source term provides the following mathematical entropy production:

(4.42)\begin{equation} \varsigma := \partial_t {\mathcal{H}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (({\mathcal{H}}+p)\boldsymbol{u}) =\sum_{i=1,2}\nu (S_0)_i^{5/2}R_{\chi_i}D_t\chi_i. \end{equation}

As the closure (4.31) groups droplets by size, we can now model the first-order size-dependent damping term of the viscous droplet in a light carrier phase (Prosperetti Reference Prosperetti1977; Plümacher et al. Reference Plümacher, Oberlack, Wang and Smuda2020). Then, we expect each population of oscillators to be damped following

(4.43)\begin{equation} D_{t}(D_t(\chi_i))+\omega_i^2 \chi_i={-}\beta_i D_t\chi_i \ \iff D_{t}\dot{\chi}_i+\omega_i^2 \chi_i={-}\beta_i \dot{\chi}_i, \end{equation}

with $\beta _i=4{\rm \pi} \nu _{vis}/(S_0)_i>0$ with $\nu _{vis}$ the liquid kinematic viscosity. One can then choose $R_{\chi _i}=-n_i\beta _i\dot {\chi }_i$ and recover both the above dissipation process for both populations of droplets and a signed production of mathematical entropy

(4.44)\begin{equation} \varsigma ={-}\beta_1n_1\nu(S_0)_1^{5/2}\dot{\chi}_1^2 -\beta_2n_2\nu(S_0)_2^{5/2}\dot{\chi}_2^2 \le 0. \end{equation}

This leads to the following dissipative model:

(4.45)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t m_k + \boldsymbol{\nabla} \boldsymbol{\cdot} (m_k \boldsymbol{u})=0, \quad k=1,2,1^d,\\ \partial_t n_i + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i \boldsymbol{u})=0, \quad i=1,2,\\ \partial_t (n_i(S_0)_i) + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i(S_0)_i \boldsymbol{u}) =0, \quad i=1,2,\\ \partial_t(n_i\dot{\chi}_i)+ \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\dot{\chi}_i \boldsymbol{u})={-}\omega_i^2 n_i\chi_i-\beta_in_i\dot{\chi}_i, \quad i=1,2,\\ \partial_t(n_i\chi_i)+ \boldsymbol{\nabla} \boldsymbol{\cdot} (n_i\chi_i \boldsymbol{u}) =n_i\dot{\chi}_i, \quad i=1,2,\\ \partial_t (\rho \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} \otimes \boldsymbol{u}+p \boldsymbol{\mathsf{I}})=0, \end{array}\right\} \end{equation}

with $p:=p_1=p_2$. This model shows a macroscopic synchronous oscillation through the two momentum equations on $\dot {\chi }_1$ and $\dot {\chi }_2$ and includes a physics-based dissipation rate. The system is here written using the weights and abscissas, but it can also be written using the geometric quantities of the model.

6.1. Geometry of the disperse regime through the classical averaging approach

First, the time evolution of the geometry for a generic two-phase flow is derived with an averaging process (Drew & Passman Reference Drew and Passman1999; Lhuillier Reference Lhuillier2004; Morel Reference Morel2015). This approach relies on the kinematics of the interface where the phases $k=1,2$ are located using a phase indicator function $\boldsymbol {X}_k$ and the interface is identified through its derivative $\boldsymbol {\nabla }\boldsymbol {X}_k$ in the sense of generalised function. Denoting by $\boldsymbol {u}_I$ the interface velocity, the kinematics of the interface reads

(6.1)\begin{equation} \partial_t \boldsymbol{X}_k + \boldsymbol{u}_I \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{X}_k = 0. \end{equation}

The ensemble-average operator $\langle {\cdot }\rangle _E$ is introduced such that the volume fraction and the interface area density are defined by $\alpha _k~=~\langle \boldsymbol {X}\rangle _E$ and $\varSigma =\langle \delta _I\rangle _E$ where $\delta _I=-\boldsymbol {n}\boldsymbol {\nabla }\boldsymbol {X}_k$ is the interface generalised function. Then, the evolution of the interface area density is derived in Lhuillier (Reference Lhuillier2004) as the trace of the interface tensor $\delta _I \boldsymbol {n}\otimes \boldsymbol {n}$. It eventually yields

(6.2)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t\alpha + \langle \boldsymbol{u}_I\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{X}\rangle _E=0,\\ \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} \langle \delta_I\boldsymbol{u}_I\rangle _E =\langle \boldsymbol{\mathsf{I}}-\boldsymbol{n}\otimes\boldsymbol{n}:\boldsymbol{\nabla}\boldsymbol{u}_I \delta_I\rangle _E. \end{array}\right\} \end{equation}

Such an averaging approach is compatible with the kinetic-based small-scale models of this work in the statistical sense. Similarly, the characteristics (size, oscillation amplitude, etc.) are following the probabilistic law given by the NDF. Focusing on the small-scale disperse regime while assuming (H1a), we have $\alpha =\alpha _1^d$, and we decompose $\boldsymbol {u}_I = \boldsymbol {u} + v_n \boldsymbol {n}$. Moreover, the symmetry of either spherical inclusions (H3a) or oscillation motion (H4d) gives $\langle v_n\boldsymbol {n}\delta _I\rangle _E=0$. The averaged equations then become

(6.3)\begin{align} \left.\begin{array}{@{}c@{}} D_t \alpha_1^d ={-}\langle v_n\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{X}\rangle _E,\\ \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u})= \dfrac{2}{3}\varSigma\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} -\langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}\boldsymbol{u}\rangle _E +\dfrac{2}{3}\langle \boldsymbol{\nabla}\boldsymbol{\cdot}(v_n\boldsymbol{n})\delta_I\rangle _E -\langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}(v_n\boldsymbol{n})\rangle _E, \end{array}\right\} \end{align}

where $\boldsymbol{\mathsf{q}} = \langle (\boldsymbol {n}\otimes \boldsymbol {n}-\boldsymbol{\mathsf{I}}/3)\delta _I\rangle _E$ is the anisotropic tensor. In this last equation, one can identify the different contributions to the interface area density evolution (from left to right of the right-hand side): the large-scale isotropic and anisotropic terms, the small-scale isotropic and anisotropic terms.

6.2. Comparison with the small-scale models

From the different small-scale models in this work, as well as proper choice of variables, we show that the interface area density equation of evolution can be thoroughly closed. We compare it with the unclosed kinematic set of equations (6.3) and identify the various contributions. The question of the choice of fundamental variables to describe the dynamics of the interface is eventually discussed.

6.2.1. Spray of compressible spherical inclusions

Before considering an incompressible small scale dedicated to the description of droplets, we focus on the case of compressible inclusions such as bubbles, which has been partially treated in § 3. We recast the unclosed dynamics for $\varSigma$ obtained in (3.8) into

(6.4)\begin{equation} \partial_t\varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u}) = \frac{2}{3}\varSigma \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} + \frac{2}{3} \varSigma \frac{D_t \alpha_1^d}{\alpha_1^d}. \end{equation}

Even though the previous equation is unclosed as the dynamics of $\alpha _1^d$ is not specified, one can still identify the isotropic contribution of the large scale of (6.3). As we consider spherical shapes, $\boldsymbol{\mathsf{q}}=\boldsymbol {0}$ and there are no anisotropic contributions. Therefore, the small-scale isotropic term is

(6.5)\begin{equation} \langle \boldsymbol{\nabla}\boldsymbol{\cdot}(v_n\boldsymbol{n})\delta_I\rangle _E = \varSigma \frac{D_t \alpha_1^d}{\alpha_1^d}. \end{equation}

A closed model is thus obtained with the dynamics of $\alpha _1^d$ as in the next case dealing with incompressible inclusions.

6.2.2. Spray of incompressible spherical droplets

The incompressible case has been treated in § 3 after assuming $D_t \rho _1^d = 0$. We have then the following conservation equations:

(6.6)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0,\\ \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u}) = 0. \end{array}\right\} \end{equation}

The anisotropic tensor $\boldsymbol{\mathsf{q}}$ is still nil and provides trivial anisotropic closures, while the isotropic closures read

(6.7a,b)\begin{equation} \langle v_n\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{X}\rangle _E = \alpha_1^d \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}, \quad \langle \boldsymbol{\nabla}\boldsymbol{\cdot}(v_n\boldsymbol{n})\delta_I\rangle _E ={-}\varSigma\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}. \end{equation}

This shows that the small-scale isotropic contributions balance the large-scale ones to maintain the incompressibility of the small scale.

6.2.3. Two-scale model with the small-scale spray model of asynchronous droplets

This case has been treated in §§ 4.2 and 5, where the oscillatory dynamics of the incompressible droplets is formulated in a decomposed form of the small-scale interface area density following $\varSigma = \varSigma _0 + \Delta \varSigma$. Considering the case with inter-scale mass-transfer, this results in the following set of equations:

(6.8)\begin{equation} \left.\begin{array}{@{}c@{}} \partial_t \alpha_1^d + \boldsymbol{\nabla} \boldsymbol{\cdot} (\alpha_1^d \boldsymbol{u}) = 0,\\ \displaystyle \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u}) ={-} \dfrac{\Delta\varSigma}{\tau} -\left(\dfrac{S_{0,avg}}{m_{avg}}+\dfrac{4}{5}q^2 \dfrac{S_{0,avg}}{m_{avg}}\right) R_{m_1}. \end{array}\right\} \end{equation}

Following the assumption (H4c) made on the orientation of the oscillating droplets, one can consider the right-hand side term as either isotropic or anisotropic

(6.9)\begin{align} \left.\!\begin{array}{@{}c@{}} \langle v_n\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{X}\rangle _E= \alpha_1^d \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u},\\ \langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}\boldsymbol{u}\rangle _E = 0,\\ \displaystyle \dfrac{2}{3}\langle \boldsymbol{\nabla}\boldsymbol{\cdot}(v_n\boldsymbol{n})\delta_I\rangle _E -\langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}(v_n\boldsymbol{n})\rangle _E ={-}\dfrac{2}{3}\varSigma\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} -\dfrac{\Delta\varSigma}{\tau} -\left(\dfrac{S_{0,avg}}{m_{avg}}+\dfrac{4}{5}q^2 \dfrac{S_{0,avg}}{m_{avg}}\right) R_{m_1}. \end{array}\!\right\} \end{align}

The excess of interface area density $\Delta \varSigma$ results from the balance between the inter-scale creation term and the exponentially dissipative decreasing. It proposes an additional perspective to the balance between coalescence and break-up as can be found in Vallet & Borghi (Reference Vallet and Borghi1999) and Anez et al. (Reference Anez, Ahmed, Hecht, Duret, Reveillon and Demoulin2019). This last phenomenon cannot be recovered here as we have not introduced any turbulence model for the small scale. This is the subject of a current work in progress.

6.2.4. Two-scale model with the small-scale spray model of synchronous droplets

This case has been studied in § 4.3, where incompressible droplets oscillate synchronously for a given size. The dynamics has been derived respectively for each size of droplets using the weights and abscissas rather than the geometric quantities. Without losing the properties of the polydisperse case, we propose here to focus on the monodisperse case detailed in Appendix I. It allows a more compact equation of evolution for the small-scale interfacial area density $\varSigma = M_{1,0,0}^{\xi } + M_{1,2,0}^{\xi }$

(6.10)\begin{equation} \partial_t \varSigma + \boldsymbol{\nabla} \boldsymbol{\cdot} (\varSigma \boldsymbol{u}) = \frac{16}{5} \frac{(\Delta\varSigma_{{\perp}})(\varDelta_t\Sigma_{{\perp}})}{(n_1^d)^{1/3}(6\sqrt{\rm \pi}\alpha_1^d)^{2/3}}. \end{equation}

Similarly to the previous cases, the incompressibility condition balances the isotropic terms to conserve $\alpha _1^d$. Moreover, the right-hand side of (6.10) only accounts for the small-scale anisotropic motion such that

(6.11)\begin{equation} \left.\begin{array}{@{}c@{}} \langle v_n\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{X}\rangle _E =\alpha_1^d \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u},\\ \langle \boldsymbol{\nabla}\boldsymbol{\cdot}(v_n\boldsymbol{n})\delta_I\rangle _E={-}\varSigma\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u},\\ \langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}\boldsymbol{u}\rangle _E=0,\\ \displaystyle \langle \boldsymbol{\mathsf{q}}:\boldsymbol{\nabla}(v_n\boldsymbol{n})\rangle _E ={-}\dfrac{2}{(3\sqrt{4{\rm \pi}})^{2/3}} \dfrac{(\Delta\varSigma_{{\perp}})(\varDelta_t\Sigma_{{\perp}})}{(n_1^d)^{1/3}(\alpha_1^d)^{2/3}}. \end{array}\right\} \end{equation}

This last identification of closures concludes our approach where the dynamics is first assessed, and the geometry kinematics is obtained secondly. Our approach shows that for this last case, the dynamics cannot be simply expressed relying on the usual quantities $\alpha _1^d$ and $\varSigma$. The SAP modelling strategy allows the proper set of variables and conservation equations to be identified in the dynamical case, and a closed conservation equation to be obtained for the interface area density in particular, an important building block for further modelling such as evaporation or heat transfer.