1 Introduction – context
Two-phase flows involve a large range of physical scales. Thus, one generally identifies different regimes (Ishii Reference Ishii1975) according to the topology and the geometry of the interfaces separating both fluids. In flows called separated-phase flows, the variation of the material interface is described at the bulk fluid scale, while for dispersed flows, like bubbly flows or sprays, one phase is disperse within a carrier fluid and the interface involves scales that are much smaller than the macroscopic dynamics of the flow. In the latter case, for the purpose of a statistical modelling, both the topology and the geometry of the interface are often fixed by the statistical configuration space such as the polydisperse spherical droplets of a spray described by a Williams–Boltzmann equation (Williams Reference Williams1958, Reference Williams1985). However, industrial applications may often involve flows where different regimes are at play successively or even simultaneously. For example in combustion chambers, as described in Kah (Reference Kah2010), Le Chenadec (Reference Le Chenadec2012), a liquid fuel is injected at high velocity and pressure. In the vicinity of the injector, the fluid forms a jet, which falls into the category of separated-phase flows. Yet, further away from the injection point, the liquid fuel turns into a cloud of liquid drops. As there is a strong connection between the distribution of the droplets in size and velocity and the quality of the combustion in such turbulent flows (Reveillon & Vervisch Reference Reveillon and Vervisch2005; Sabat et al. Reference Sabat, Vié, Larat and Massot2019), the accurate prediction of the birth of this cloud in the chamber is a key challenge for the industrial community (Le Touze Reference Le Touze2015; Fiorina et al. Reference Fiorina, Vié, Franzelli, Darabiha, Massot, Dayma, Dagaut, Moureau, Vervisch and Berlemont2016). Another example is given by the boiling crisis issue in nuclear primary circuits. There, the opposite process occurs: the liquid phase vaporizes and generates a disperse flow of small bubbles. These may coalesce and hinder an efficient cooling of the nuclear core, once a critical volume fraction of vapour is reached (see a description of the physical context and the safety stakes in Fleau et al. (Reference Fleau, Mimouni, Mérigoux and Vincent2015)).
The large range of scales and the difference in the nature of the regimes involved make two-phase flows hard to simulate. Direct numerical simulation, i.e. trying to capture the scales of the smallest interfacial variations, is restricted today to academic configurations (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Ling et al. Reference Ling, Fuster, Zaleski and Tryggvason2017; Vaudor et al. Reference Vaudor, Ménard, Aniszewski, Doring and Berlemont2017) and is too expensive in an industrial context. Consequently, an important need in model development has arisen. New modelling approaches should enable the description of a range of scales to describe features of the interface, below a cutoff threshold that cannot be resolved. These unresolved scales will be referred to as small scales (with respect to the bulk scale) in the following. Accounting for small scales requires us to provide a model for them directly at the bulk scale. We shall call this sub-scale modelling.
For disperse-phase flows, where all interfacial structures are sub-scale, various approaches have been proposed. On the one hand, two-fluid systems of equations are obtained by averaging flow models involving fluids separated by sharp interfaces (Ishii Reference Ishii1975). This process yields equations for a two-phase mixture. However, standard averaging leads to limited information about the interfacial structures: one only gets access to quantities like the volume fraction or the mean interfacial area. Other higher-order quantities such as mean curvatures may be introduced (Drew Reference Drew1990), but the closure of the associated equations necessitates additional assumptions on the flow topology. The mathematical properties of such systems are not straightforward: dissipative terms may not be clearly identified, neither may the equilibrium and asymptotic regimes. On the other hand, one can consider a statistical approach inspired from kinetic theory, where the modelling relies on a probability density function defined over configuration space that spans features of the flow like: the size of the particles, temperature, additional shape criteria. Such an approach has led to many contributions, see for example Marchisio & Fox (Reference Marchisio and Fox2005), Masi, Simonin & Bédat (Reference Masi, Simonin and Bédat2011), Yuan, Laurent & Fox (Reference Yuan, Laurent and Fox2012), Vié, Laurent & Massot (Reference Vié, Laurent and Massot2013), Masi et al. (Reference Masi, Simonin, Riber, Sierra and Gicquel2014), Sabat et al. (Reference Sabat, Vié, Larat and Massot2019) and references therein. Extensions to other flow regimes are not straightforward and are a current research topic (Essadki et al. Reference Essadki, de Chaisemartin, Laurent and Massot2018, Reference Essadki, Drui, de Chaisemartin, Larat, Ménard and Massot2019).
Separated-phase regimes intend to describe the interface at the same scale as the scale of the bulk flow resolution. Some methods are based on a sharp description of the interface like the front tracking method (Chern et al. Reference Chern, Glimm, Mcbryan, Plohr and Yaniv1986; Tryggvason et al. Reference Tryggvason, Bunner, Esmaeeli, Juric, Al-Rawahi, Tauber, Han, Nas and Jan2001; Terashima & Tryggvason Reference Terashima and Tryggvason2010), the level set method (Mulder, Osher & Sethian Reference Mulder, Osher and Sethian1992; Fedkiw et al. Reference Fedkiw, Aslam, Merriman and Osher1999; Osher & Fedkiw Reference Osher and Fedkiw2001; Liu, Khoo & Wang Reference Liu, Khoo and Wang2005) or the volume of fluid (VOF) method (Hirt & Nichols Reference Hirt and Nichols1981; Ménard, Tanguy & Berlemont Reference Ménard, Tanguy and Berlemont2007). Other methods may also use mixture models for capturing the interface in a transition zone (Saurel & Abgrall Reference Saurel and Abgrall1999; Allaire, Clerc & Kokh Reference Allaire, Clerc and Kokh2002; Chanteperdrix, Villedieu & Vila Reference Chanteperdrix, Villedieu and Vila2002; Kokh & Lagoutiére Reference Kokh and Lagoutiére2010; Grenier, Vila & Villedieu Reference Grenier, Vila and Villedieu2013; Saurel & Pantano Reference Saurel and Pantano2018). Such techniques usually do not involve sub-scale modelling. Recent efforts have been conducted in the community (Vincent et al. Reference Vincent, Larocque, Lacanette, Toutant, Lubin and Sagaut2008; Anez et al. Reference Anez, Ahmed, Hecht, Duret, Reveillon and Demoulin2019; Kedelty, Uglietta & Herrmann Reference Kedelty, Uglietta and Herrmann2018) in the spirit of large eddy simulations (LES) to model the interaction between interfaces and small-scale phenomena.
In two-phase flows, many phenomena can be expressed at different scales like pressure effects, capillarity or heat transfer. In the present work we wish to explore the coupling between these scales within the restricted focus on kinematic phenomena by considering a classic bulk motion, supplemented with small-scale effects. At first glance, it would seem reasonable to aim for models involving a bulk velocity for each component, like the well-known models of Ishii (Reference Ishii1975), Drew (Reference Drew1983), Drew & Passman (Reference Drew and Passman1999), Morel (Reference Morel2015) for example. However, such models are still a very active research area that covers a wide range of matters from physical modelling to mathematical ill-posedness issues. Indeed, guaranteeing unconditional hyperbolicity for such models is not straightforward (Stuhmiller Reference Stuhmiller1977; Sha & Soo Reference Sha and Soo1979; Sursock Reference Sursock1982; Drew & Passman Reference Drew and Passman1999; Ndjinga Reference Ndjinga2007; Lhuillier, Chang & Theofanous Reference Lhuillier, Chang and Theofanous2013; Panicker, Passalacqua & Fox Reference Panicker, Passalacqua and Fox2018), nor is the definition of weak solutions due to the presence of non-conservative terms in the equations (Andrianov & Warnecke Reference Andrianov and Warnecke2004). Although these problems are crucial in the general case, we prefer to trade generality for a much simpler framework by assuming that each component has the same bulk velocity, in order to strictly focus on the two-scale analysis. Obviously, these assumptions do not match realistic flows pertaining to complex applications. However, they are consistent with our analysis, that aims at singling out the kinematics effects of each scale, and to study how they couple the evolution of these scales.
In order to cautiously introduce this two-scale kinematics hypothesis, we choose to rely on Hamilton’s principle of stationary action. Indeed, after choosing a set of parameters that will describe the system and its physical constraints, like mass conservation, the stationary action principle (see e.g. Finlayson (Reference Finlayson1972), Gavrilyuk & Gouin (Reference Gavrilyuk and Gouin1999), Berdichevsky (Reference Berdichevsky2009) and references therein) only requires us to define scalar-valued functions describing energies associated with the system evolution. In our case, this boils down to considering the classic kinetic and potential energies involved in a simple barotropic two-phase model, and to add an extra term for the small-scale kinetic energy, as in Gavrilyuk & Saurel (Reference Gavrilyuk and Saurel2002), Gavrilyuk (Reference Gavrilyuk, dell’Isola and Gavrilyuk2012). Then, we obtain a system of five conservation laws that is checked to be hyperbolic. In a second step, we complete the system with a dissipative structure, by adding new terms that will be compatible with an energy (mathematical entropy) inequality.
The resulting complete 5-equation model does not a priori rely on hypotheses regarding the topology of the interface. It involves a convective part and some stiff sources terms. These source terms can be interpreted as relaxation terms. By supposing infinitely fast relaxation processes, two equilibrium systems can be derived from our 5-equation model. These models have been used in the literature for the modelling of separated-phase flows (Chanteperdrix et al. Reference Chanteperdrix, Villedieu and Vila2002; Caro et al. Reference Caro, Coquel, Jamet and Kokh2006; Grenier et al. Reference Grenier, Vila and Villedieu2013; Bernard-Champmartin & De Vuyst Reference Bernard-Champmartin and De Vuyst2014). In this sense, this 5-equation system can be considered as a parent model for simple separated-phase flows models. Nevertheless, this system is also compatible with the description of a dispersed flow model. More specifically, under adequate hypotheses, it is possible to retrieve from this model an evolution equation for the bubble dynamics that is similar to the Rayleigh–Plesset equation (Lord Rayleigh Reference Lord Rayleigh1917; Plesset & Prosperetti Reference Plesset and Prosperetti1977). Thus, we have obtained a hierarchy of two-phase models that includes models suitable for both separated- and disperse-phase flows. This inter-model connection enables two different estimates for two relaxation parameters related to mechanical equilibrium between materials in the separated-phase model. This result is notable, as these estimates are often replaced by infinitely fast relaxation processes or heuristic values (Chanteperdrix et al. Reference Chanteperdrix, Villedieu and Vila2002; Gallouet, Hérard & Seguin Reference Gallouet, Hérard and Seguin2004; Hérard & Hurisse Reference Hérard and Hurisse2005; Saurel, Petitpas & Berry Reference Saurel, Petitpas and Berry2009). Few works have already addressed this matter and proposed similar estimates. The choice of this parameter has a significant influence on the damping effects in the system, as will be shown in the ending sections. The different choices for these fluid parameters are then tested by studying the acoustic regime of the models and comparing their dispersion relations with reference data issued from both experiments and the bubbly flows model studied in Cheng, Drew & Lahey (Reference Cheng, Drew and Lahey1985). This study will also shed some light on the similarities and differences between the models of the hierarchy, by considering their acoustic behaviour as a benchmark tool.
Before going any further, let us mention that this work aims at contributing at a wider effort in the numerical simulation of two-phase flows involving different scales and different regimes. Indeed, as both disperse- and separated-flow models may be connected through a parent model via limits of relaxation processes, this suggests we investigate numerical methods that comply with such a property at the discrete level. First, the community has provided many important contributions with the development of asymptotic preserving (AP) techniques (see e.g. Caflisch, Jin & Russo Reference Caflisch, Jin and Russo1997; Jin Reference Jin1999; Gosse & Toscani Reference Gosse and Toscani2004; Buet & Cordier Reference Buet and Cordier2007; Dumbser, Enaux & Toro Reference Dumbser, Enaux and Toro2008; Berthon & Turpault Reference Berthon and Turpault2011; Chalons, Girardin & Kokh Reference Chalons, Girardin and Kokh2013). Such discretization methods for systems of equations are endowed with stability and accuracy properties that are preserved even when considering the limits to asymptotic regimes. Then the availability of an entropy inequality that accounts for source terms related to the small-scale kinematic effects suggests that it can be used as a stability criterion including multi-scale phenomena (Gallice Reference Gallice2003; Bouchut Reference Bouchut2004). Moreover, due to the key role of the stiff source terms in the system, is important to consider stiff ordinary differential equation (ODE) integrators like those presented in Hairer & Wanner (Reference Hairer and Wanner1996). Finally, the acoustic study provided in this work may also be used as a benchmark for numerical development. These important matters do not fit the scope of the present work and will not be discussed further in the following.
The paper is structured as follows. First, we apply classic modelling guidelines: the conservative part of the 5-equation model is derived thanks to the stationary action principle. Next, the system is equipped with a dissipative structure by enabling an entropy budget. After completing the definition of the parent 5-equation model, two sub-models are derived by considering the different cases of vanishing parameters. In a second part of the paper, after having led the first modelling part in a quite agnostic manner with respect to the two-phase flow regime, we shall consider the specific case of a bubbly fluid. In this context, we see that it is possible to identify the two parameters of the parent 5-equation model as a micro-viscosity and a micro-inertia. Furthermore, estimates for these parameters will be proposed. Finally, we study each model of the hierarchy in the acoustic regime, by expressing their respective dispersion relation. These acoustic behaviours are then compared practically with experimental and reference data.
2 A 5-equation hyperbolic two-fluid system
Before going into exposing the paradigm in the following subsection, we will present the four main simplifying assumptions used in this work:
(H1) there is no mass transfer between both phases;
(H2) there is no slip velocity between both phases;
(H3) we use a barotropic equation of state (EOS) for each pure phase and will not consider the evolution of energy and energy exchanges of the phases;
(H4) interfacial forces are neglected.
These assumptions are restrictive but they enable a simple framework that allows us to focus on our two-scale study. Extensions will be discussed in the conclusion. Both fluids are governed by a barotropic EOS characterized by

where
$\unicode[STIX]{x1D70C}_{k}$
,
$f_{k}$
and
$p_{k}$
are respectively the densities, the specific barotropic energies and the partial pressures of each fluid,
$k=1,2$
. Then,
$c_{k}^{2}=(\text{d}p_{k}/\text{d}\unicode[STIX]{x1D70C}_{k})$
denote the sound velocities within each pure material
$k$
. If we note by
$Y_{k}$
the mass fraction of each fluid, we have
$Y_{1}+Y_{2}=1$
. We moreover postulate that fluids
$1$
and
$2$
are immiscible, thus we can also define the volume fractions
$\unicode[STIX]{x1D6FC}_{k}$
, such that
$\unicode[STIX]{x1D6FC}_{1}+\unicode[STIX]{x1D6FC}_{2}=1$
. For convenience, in the rest of this paper we set
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{1}$
,
$Y=Y_{1}$
, which allows us to define the density
$\unicode[STIX]{x1D70C}$
of the medium by

Finally, we suppose given a function

that is the barotropic energy of the medium. The choice of
$f$
will be specified later.
Thanks to hypotheses (H1) and (H2), the total mass and the mass fractions of each fluid are conserved under the same velocity field
$\boldsymbol{u}$
, what can be written as


or equivalently as
$\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k})+\text{div}(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\boldsymbol{u})=0$
,
$k=1,2$
.
2.1 Variational principle for two-phase models
The first step in our modelling work aims at deriving a system of conservation laws equipped with a mathematical entropy evolution equation and checking that the convective structure of the system is hyperbolic. Following the lines of Berdichevsky (Reference Berdichevsky2009) for barotropic fluids, of Gavrilyuk & Saurel (Reference Gavrilyuk and Saurel2002) for a system of equations for two compressible fluids and two temperatures and of Caro et al. (Reference Caro, Coquel, Jamet and Kokh2006) for a homogeneous isothermal two-fluid model, we propose to use Hamilton’s principle of stationary action to derive the conservative structure of our two-phase model. Before going any further, we introduce a few notations: if
$\boldsymbol{a}$
and
$\boldsymbol{b}$
are two column vectors whose components are
$a_{i}$
and
$b_{i}$
,
$1\leqslant i\leqslant d$
, then
$\boldsymbol{a}^{T}$
is a row vector,
$\boldsymbol{a}^{\text{T}}\boldsymbol{b}=\sum _{i}a_{i}b_{i}$
is a real number and
$\boldsymbol{a}\boldsymbol{b}^{\text{T}}$
is a square matrix of size
$d$
with
$(\boldsymbol{a}\boldsymbol{b}^{\text{T}})_{i,j}=a_{i}b_{j}$
. If
$\boldsymbol{x}\in \mathbb{R}^{d}\mapsto A$
is a field of square matrices of size
$d$
then
$\text{div}(\unicode[STIX]{x1D608})$
is a column vector of size
$d$
, where
$\text{div}(\unicode[STIX]{x1D63C})_{i}=\sum _{j}\unicode[STIX]{x2202}_{x_{j}}\unicode[STIX]{x1D608}_{ij}$
. Finally, we note by
$D_{t}(\cdot )=\unicode[STIX]{x2202}_{t}(\cdot )+\boldsymbol{u}^{\text{T}}\unicode[STIX]{x1D735}(\cdot )$
the material derivative.
Now, we need to define the kinetic energy
$E_{kin}$
and the potential energy
$E_{pot}$
that are involved in the conservative transformations of our system. We postulate that the specific potential energy of our system is the free energy
$f$
. Also, following Gavrilyuk & Saurel (Reference Gavrilyuk and Saurel2002), a key element of our study consists in considering a two-scale kinetic energy composed of a bulk kinetic energy
$\unicode[STIX]{x1D70C}|\boldsymbol{u}|^{2}/2$
and a small-scale kinetic energy
$\unicode[STIX]{x1D708}(D_{t}\unicode[STIX]{x1D6FC})^{2}/2$
, where
$\unicode[STIX]{x1D708}>0$
is a function of the volume fraction
$\unicode[STIX]{x1D6FC}$
only and whose physical interpretation and characterization is given later on. This allows a finer representation of general two-phase flows and the consequences of such a hypothesis will be further discussed in §§ 3 and 4.
Now, the Lagrangian
${\mathcal{L}}=E_{kin}-E_{pot}$
of our system is fully defined, namely,

Let
$\mathscr{V}(t)$
be the volume occupied by a portion of the fluid at time
$t\in [t_{0},t_{1}]$
and
$\unicode[STIX]{x1D6FA}=\{(\boldsymbol{x},t)\in \mathbb{R}^{3}\times [t_{0},t_{1}]~|~\boldsymbol{x}\in \mathscr{V}(t)\}$
be the subset of all space–time points covered by
$\mathscr{V}(t)$
during
$t_{0}\leqslant t\leqslant t_{1}$
. The space variable
$\boldsymbol{X}\in \mathscr{V}(t_{0})$
denotes the Lagrangian coordinates associated with the reference frame at instant
$t=t_{0}$
. If
$(\boldsymbol{X},t)\mapsto \unicode[STIX]{x1D753}$
is the mapping that gives the position
$\unicode[STIX]{x1D753}(\boldsymbol{X},t)$
at instant
$t$
of a fluid element that was located at
$\boldsymbol{X}$
at time
$t=t_{0}$
, then obviously
$\unicode[STIX]{x1D6FA}=\{(\unicode[STIX]{x1D753}(\boldsymbol{X},t),t)~|~\boldsymbol{X}\in \mathscr{V}(t_{0}),~t_{0}\leqslant t\leqslant t_{1}\}$
.
From a pure Eulerian point of view, a transformation of the medium is fully characterized by the fields
$(\boldsymbol{x},t)\mapsto (\unicode[STIX]{x1D70C},\boldsymbol{u},Y,\unicode[STIX]{x1D6FC})$
. Equivalently, we can say that a transformation of the medium is fully characterized by the Eulerian fields
$(\boldsymbol{x},t)\mapsto (Y,\unicode[STIX]{x1D6FC})$
and the Lagrangian mapping
$(\boldsymbol{X},t)\mapsto \unicode[STIX]{x1D753}$
under the hypothesis that
$\unicode[STIX]{x1D753}$
is compliant with the mass conservation.
If
$(\boldsymbol{x},t)\mapsto (Y,\unicode[STIX]{x1D6FC})$
and
$(\boldsymbol{X},t)\mapsto \unicode[STIX]{x1D753}$
is a given transformation of the medium, we consider a family of transformations
$(\boldsymbol{x},t,\unicode[STIX]{x1D701})\mapsto (\widehat{Y},\widehat{\unicode[STIX]{x1D6FC}})$
and
$(\boldsymbol{X},t,\unicode[STIX]{x1D701})\mapsto \widehat{\unicode[STIX]{x1D753}}$
parametrized by
$\unicode[STIX]{x1D701}$
which lies in the vicinity of
$0$
, such that:
(i)
$(\widehat{Y},\widehat{\unicode[STIX]{x1D6FC}})(\boldsymbol{x},t,\unicode[STIX]{x1D701}=0)=(Y,\unicode[STIX]{x1D6FC})(\boldsymbol{x},t)$ and
$\widehat{\unicode[STIX]{x1D753}}(\boldsymbol{X},t,\unicode[STIX]{x1D701}=0)=\unicode[STIX]{x1D753}(\boldsymbol{X},t)$ ;
(ii)
$\widehat{Y}$ and
$\widehat{\unicode[STIX]{x1D753}}$ verify the mass and partial mass conservation (2.4) and (2.5), for
$\unicode[STIX]{x1D701}$ close to
$0$ .
We adopt the classic definition of the infinitesimal transformations that act on the medium by introducing the infinitesimal displacement
$(\boldsymbol{x},t)\mapsto \unicode[STIX]{x1D743}$
, where

and by setting for any Eulerian field
$(\boldsymbol{x},t,\unicode[STIX]{x1D701})\mapsto \widehat{b}$

Using the lines of Gavrilyuk & Saurel (Reference Gavrilyuk and Saurel2002), Berdichevsky (Reference Berdichevsky2009), Gavrilyuk (Reference Gavrilyuk, dell’Isola and Gavrilyuk2012), our hypotheses provide the infinitesimal variations of
$\unicode[STIX]{x1D70C}$
,
$Y$
,
$u$
and
$D_{t}\unicode[STIX]{x1D6FC}$
. We obtain indeed that

and

Now, we can define the Hamiltonian action
${\mathcal{A}}(\unicode[STIX]{x1D701})$
by setting

and compute the infinitesimal variations of the action,
$\unicode[STIX]{x1D6FF}{\mathcal{A}}=(\text{d}{\mathcal{A}}/\text{d}\unicode[STIX]{x1D701})(\unicode[STIX]{x1D701}=0)$
:

We make the classic assumption (Gavrilyuk & Gouin Reference Gavrilyuk and Gouin1999; Gavrilyuk & Saurel Reference Gavrilyuk and Saurel2002) that for our transformation family,
$\unicode[STIX]{x1D743}$
and
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}$
vanish on
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$
. Thanks to Green’s formula, we obtain after tedious calculations that

We follow the stationary action principle that boils down to the postulate that a physical transformation of the medium should extremize the Hamiltonian action
${\mathcal{A}}$
. In our case, this yields


For the choice of
${\mathcal{L}}$
expressed by (2.6), we have


Relations (2.13) and (2.14) will respectively provide the evolution equations for the momentum and the volume fraction. Indeed, reinjecting (2.15) into (2.13)–(2.14) provides


We define the medium pressure
$p$
and a new variable
$w$
by setting

Using the mass and partial mass conservation hypotheses, we obtain that the fluid transformations are governed by the following system of equations





Here, the stationary action principle only provides the conservative elements for our model: the momentum (2.18c
), the evolution equation for the volume fraction (2.18d
) and the small-scale pulsation evolutions (2.18e
). These equations supplement masses conservation, equations (2.18a
)–(2.18b
), which has been postulated. Let us note that the evolution of
$\unicode[STIX]{x1D6FC}$
is driven by (2.18d
) and (2.18e
) which involves more complex features than transport phenomena. Equation (2.18e
) can be interpreted as a small-scale momentum conservation equation. It is not surprising: the main role of the stationary action principle is to derive a momentum conservation equation from the data of kinetic energy and potential energy. In our case, we supplied kinetic energies for both the bulk and the small scales. Consequently, we recover momentum conservation equations for the bulk and the small scales that are coupled together. For the bulk-scale momentum, the small-scale kinematics accounts for small-scale inertia. In the following,
$\unicode[STIX]{x1D708}$
will be referred to as micro-inertia.
We shall examine the dissipative structures of system (2.18) in the next section.
2.2 Dissipation and second principle of thermodynamics
While we are not interested in bulk dissipation phenomena, we aim at describing small-scale dissipation associated with the small-scale kinetic energy
$\unicode[STIX]{x1D708}(D_{t}\unicode[STIX]{x1D6FC})^{2}/2$
, see equation (2.6). This way, the damping of small-scale pulsations of the interface due to various dissipative phenomena like the liquid viscosity can be taken into account.
Irreversible damping effects can be introduced in our two-phase system, by adding terms to (2.18) that are compatible with a mathematical entropy evolution equation. To do so, we first suppose that (2.18a
)–(2.18d
) are valid, but we discard (2.18e
) and consider that
$D_{t}w$
is now a quantity to be defined.
In the barotropic case, it is classic (Godlewski & Raviart Reference Godlewski and Raviart1996; Serre Reference Serre2007) to consider the total free energy of the system as a mathematical entropy

We now seek for an entropy flux function
$\boldsymbol{G}$
and a proper evolution principle for
$D_{t}w$
, such that

or equivalently

Relation (2.21) reads

Using (2.18a
)–(2.18d
) to express
$D_{t}\boldsymbol{u}$
,
$D_{t}\unicode[STIX]{x1D70C}$
,
$D_{t}Y$
,
$D_{t}\unicode[STIX]{x1D6FC}$
and

altogether with (2.17), we obtain that

A simple choice for ensuring (2.24) consists in setting

where
$\unicode[STIX]{x1D700}>0$
is a constant. Let us underline that this choice is guided by our knowledge of the non-dissipative structure of the system of equation resulting from the stationary action principle. This yields a definition of
$\boldsymbol{G}$
and a new evolution equation for
$w$
that reads

Consequently, the generic form of our two-phase flow system reads





In order to complete the definition of our model, we need to specify the free energy of the medium. We consider that
$f$
is the sum of a bulk mixture free energy and a compaction energy
$\unicode[STIX]{x1D6FC}\mapsto e(\unicode[STIX]{x1D6FC})$
, where
$e$
is a given function (Gavrilyuk & Saurel Reference Gavrilyuk and Saurel2002). We set

For this choice, granted that
$\unicode[STIX]{x1D70C}_{k}^{2}\unicode[STIX]{x2202}f_{k}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{k}=p_{k}$
,
$k=1,2$
, a straightforward calculation gives that

and system (2.27) reads here






Let us conclude this section by stating well-posedness properties of our two-phase model with micro-inertia
$\unicode[STIX]{x1D708}$
. We consider the sole convective part of system (2.27) for one-dimensional problems by discarding the source terms. The resulting system is hyperbolic and its characteristic velocities are

where

Details related to the eigenstructure of system (2.27) are presented in appendix B.
2.3 Submodels and traversing the hierarchy
Based on our assumptions, the richest model (2.30) involves both parameters
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
. We examine two limit flow regimes, obtained for vanishing values of the parameters
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
. First, we study the case of a negligible micro-inertia compared to the internal dissipation effects, i.e.
$\unicode[STIX]{x1D708}\rightarrow 0$
and
$\unicode[STIX]{x1D700}=O(1)$
. Second, we consider the case when both micro-inertia and internal dissipation tend to zero, with
$\unicode[STIX]{x1D708}\rightarrow 0$
and
$\unicode[STIX]{x1D700}\rightarrow 0$
. Both cases allow to recover the two-phase systems presented in (Chanteperdrix et al.
Reference Chanteperdrix, Villedieu and Vila2002). These systems are both composed of a conservative part and a (possibly null) stiff source term.
2.3.1 A 4-equation model for
$\unicode[STIX]{x1D708}\rightarrow 0$
and
$\unicode[STIX]{x1D700}=O(1)$
We consider system (2.30) and suppose here that
$\unicode[STIX]{x1D708}\rightarrow 0$
for a fixed value of
$\unicode[STIX]{x1D700}$
. We also assume here that
$\unicode[STIX]{x1D708}^{\prime }(\unicode[STIX]{x1D6FC})/\unicode[STIX]{x1D708}$
is proportional to
$\unicode[STIX]{x1D6FC}^{-1}$
. This assumption will be justified in § 3.1 where an identification of
$\unicode[STIX]{x1D708}$
is proposed. Then, we can see that (2.30e
) provides that

By using equation (2.30d ), we obtain that the limit regime is governed by







The stiff source term in (2.34d ) drives all the dissipation effects within system (2.34).
2.3.2 A 3-equation model for
$\unicode[STIX]{x1D708}\rightarrow 0$
and
$\unicode[STIX]{x1D700}\rightarrow 0$
The last model of our hierarchy can be obtained either by considering model (2.30) and the two vanishing coefficients
$\unicode[STIX]{x1D708}\rightarrow 0$
and
$\unicode[STIX]{x1D700}\rightarrow 0$
, or by taking
$\unicode[STIX]{x1D700}\rightarrow 0$
in (2.34). Formally, we obtain







In the specific case
$\text{d}e/\text{d}\unicode[STIX]{x1D6FC}=0$
, we recover the classic partial pressure equilibrium closure relation
$p_{1}=p_{2}$
that was studied in Chanteperdrix et al. (Reference Chanteperdrix, Villedieu and Vila2002): the resulting system is strictly hyperbolic as it possesses three distinct real-valued characteristic velocities
$\{u-c_{\mathit{Wood}},u,u+c_{\mathit{Wood}}\}$
, where
$c_{\mathit{Wood}}$
is defined by

Concerning the mathematical properties of (2.36), we refer the reader to Chanteperdrix et al. (Reference Chanteperdrix, Villedieu and Vila2002). One can note that for both relations (2.32) and (2.37), the subcharacteristic condition is verified since
$c_{\mathit{Wood}}\leqslant c_{\mathit{Frozen}}\leqslant c$
. Let us remark that these characteristic velocities do not match the propagation velocities of acoustic waves as the latter depend on the the wave frequency. This will be made clear in § 3.3.
3 Application to bubbly fluids in the small perturbation regime
The derivation of model (2.27) relies on general mechanics and thermodynamics principles that do not involve specific hypotheses on the flow topology. Thus, the hierarchy of two-fluid models described in § 2.3 is generic and the interpretation of some specific terms and variables such as the small-scale kinetic energy
$\frac{1}{2}\unicode[STIX]{x1D708}(D_{t}\unicode[STIX]{x1D6FC})^{2}$
, the variable
$w$
or the constant
$\unicode[STIX]{x1D700}$
, depend on the precise physical context within which the hierarchy is used.
Hereafter, we focus on the case of an incompressible liquid transporting a disperse phase of compressible spherical homogeneous bubbles. Even if this assumption is not realistic in most cases, we still assume that there is no slip velocity between the phases. Besides, since we will consider bubbles, which are only pulsating in volume while remaining spherical, the impact of assumption (H4) is not too strong an assumption. In Prosperetti (Reference Prosperetti1977), Prosperetti studies thermal effects in bubble oscillations and proves that, in the small perturbation limit, the bubble oscillation mechanism is driven by three dimensionless numbers,

where
$\bar{\unicode[STIX]{x1D706}}$
is the mean free path in the gas,
$\unicode[STIX]{x1D706}_{g}$
is the wavelength within the bubble,
$R_{0}$
is the bubble radius at equilibrium and
$\unicode[STIX]{x1D700}_{g}^{th}$
and
$\unicode[STIX]{x1D700}_{l}^{th}$
are the characteristic thicknesses of thermal conduction within the gas and the liquid;
$G_{1}$
is shown to measure the ratio of the characteristic penetration thicknesses of acoustic and thermal phenomena and
$G_{1}\ll 1$
for a broad range of acoustic frequencies. Next,
$G_{2}$
compares the thermal penetration characteristic thickness with the bubble radius and this allows us to discriminate two main regimes. When
$G_{2}\ll 1$
, thermal equilibrium is always reached within an acoustic period and the bubble oscillation can be considered as isothermal. When
$G_{2}\gg 1$
, thermal conduction between liquid and gas is negligible and the transformation can be considered as adiabatic. In both cases, internal energy is simply driven by the other state variables and the energy equation is redundant. Hypothesis (H3) thus extends this last statement to all of the regimes in
$G_{2}$
, meaning that, with
$G_{2}$
, bubble oscillation goes from an isothermal to an adiabatic regime and no energy equation is needed. We make three additional modelling assumptions:
(H5) the mass of each bubble remains constant during a medium transformation (no break-up nor collapse);
(H6) the surrounding liquid is incompressible and thus has a constant density
$\unicode[STIX]{x1D70C}_{2}=\overline{\unicode[STIX]{x1D70C}}_{2}$ ;
(H7) for a bubble located at
$(\boldsymbol{x},t)$ , the pressure is uniform within the bubble and equal to
$p_{1}(\boldsymbol{x},t)$ and the pressure of the surrounding liquid is equal to
$p_{2}(\boldsymbol{x},t)$ .
Under these simple hypotheses, we obtain a specific model for bubbly fluids, compatible with system (2.30). By identification and comparison with standard models for bubbly flows, this allows us to give a physical meaning to the specific new terms listed above, to provide analytical expressions of the various parameters introduced during the modelling process and eventually to perform both physical and mathematical analyses of the models in this context, by looking at their behaviour in the acoustic regime.
3.1 Connection with the Rayleigh–Plesset equation
Let us consider a bubbly fluid which is monodisperse and characterized as follows: at each position and instant
$(\boldsymbol{x},t)$
, the distribution of the number of bubbles is defined by the density function
$n(\boldsymbol{x},t)$
and all bubbles are spherical with a radius
$R(\boldsymbol{x},t)$
. Given the phasic gas density
$\unicode[STIX]{x1D70C}_{1}(\boldsymbol{x},t)$
, all the bubbles have the same mass
${\mathcal{M}}_{b}(\boldsymbol{x},t)=4/3\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{1}(\boldsymbol{x},t)R^{3}(\boldsymbol{x},t)$
. Then, the gas volume fraction
$\unicode[STIX]{x1D6FC}$
and partial mass
$\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{1}$
can be related to the flow structure parameters by

Hypothesis (H5) boils down to
$D_{t}{\mathcal{M}}_{b}=0$
. Using the conservation of the partial mass
$\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{1}=\unicode[STIX]{x1D70C}Y$
for the gas (2.30b
) and
$n=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{1}/{\mathcal{M}}_{b}$
, we obtain the following conservation law for
$n$

The conservation of the partial mass
$(1-\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D70C}_{2}=\unicode[STIX]{x1D70C}(1-Y)$
for the surrounding fluid and hypothesis (H6) imply that
$\text{div}(\boldsymbol{u})$
is constrained by
$D_{t}\unicode[STIX]{x1D6FC}$
through the relation

Now, let us express
$D_{t}\unicode[STIX]{x1D6FC}$
,
$w$
and
$D_{t}w$
in terms of
$R(t)$
and its material derivatives. We have
$\unicode[STIX]{x1D6FC}=(4\unicode[STIX]{x03C0}/3)R^{3}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{1}/{\mathcal{M}}_{b}$
, which yields

Using (3.4), we obtain

Expressing
$m_{1}$
thanks to (3.2) and using (2.30d
) leads to

Relations (3.7) and (3.5) provide

Finally, combining (3.8) and (2.30e
) gives the evolution equation for
$R$
,

We see that the two-phase model (2.30), supplemented by our bubbly flow structure assumptions (H5–H6), provide the evolution (3.9) for a spatial distribution of bubble radii.
In fact, we clearly see that the bubble radii evolve only along the trajectories of the flow. Therefore, by means of the following notation:
${\dot{R}}=D_{t}R$
and
$\ddot{R}=D_{tt}^{2}R$
, we obtain an ODE of the radii along the trajectory of the flow

This ODE is analogous to a nonlinear oscillator with damping and forcing terms. This analogy also shows that the micro-inertia
$\unicode[STIX]{x1D708}$
is indeed connected to inertial effects, while
$\unicode[STIX]{x1D700}$
is related to damping and will be referred to as ‘micro-viscosity’. Now, the evolution of
$R(t)$
along the trajectories can be compared with other existing models that account for bubble vibrations in specific flow regimes, in order to provide estimates for
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D700}$
in (2.30).
First, we propose to proceed with the Rayleigh–Plesset equation (A 4g
) of the Cheng, Drew and Lahey (CDL) system (A 4). To do so, we examine a flow regime involving a bubble radius and a volume fraction small enough such that
$\unicode[STIX]{x1D6FC}\ll 1$
in (2.30) and
$kR\ll 1$
in (A 4). We neglect the surface tension in (A 4g
) by setting
$\unicode[STIX]{x1D70E}=0$
and relate
$p_{1}-p_{2}-\unicode[STIX]{x1D70C}\,\text{d}e/\text{d}\unicode[STIX]{x1D6FC}$
in (3.9) to
$p_{1i}-p_{2}$
in (A 4g
), recalling the hypothesis of uniform pressures (H7). By identifying the terms in
${\dot{R}}$
,
$\ddot{R}$
and
${\dot{R}}^{2}$
, we respectively obtain

The first two relations of (3.11) allow us to identify respectively
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
as

The third relation of (3.11) is redundant but compatible with the definition of
$\unicode[STIX]{x1D708}$
expressed in (3.12. Thus we see that for this specific flow regime, it is possible to insert a monodisperse bubble flow model into our two-phase model, for which the dynamics of the bubble radii degenerates to the Rayleigh–Plesset equation. Also, the resulting values of
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D700}$
match the results of Gavrilyuk & Saurel (Reference Gavrilyuk and Saurel2002) derived in the context of a Baer–Nunziato-type two-fluid model.
Next, let us remark that in the context of bubbly fluids, similar models for micro-inertia are available. In Temkin (Reference Temkin2005), a pulsational energy is considered in the form:
$E_{puls}={\textstyle \frac{1}{2}}M_{puls}{\dot{R}}^{2}$
and
$M_{puls}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{2}R^{3}$
. An alternative approach that allows us to incorporate small-scale bubble velocity is also proposed in Lhuillier, Theofanous & Liou (Reference Lhuillier, Theofanous and Liou2010), by accounting for pseudo-turbulent kinetic energies generated by particle pulsations
$K_{c}={\textstyle \frac{1}{2}}Q(\unicode[STIX]{x1D6FC})(D_{t}R)^{2}$
, where
$Q(\unicode[STIX]{x1D6FC})$
is proportional to
$3\unicode[STIX]{x1D6FC}$
in the dilute limit of the dispersed phase. In Gavrilyuk (Reference Gavrilyuk, dell’Isola and Gavrilyuk2012), a micro-scale kinetic energy of the form
$3\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{2}(D_{t}R)^{2}$
is used for modelling vibrations of bubbles within a bubbly fluid.
3.2 Connection with the linearized Rayleigh equation
Further, we examine the bubbly fluid model of § 3.1 by considering the regime of small variations of the bubbles radius. Let us assume that

where
$\overline{R}$
,
$\overline{p_{2}}$
,
$\overline{{\mathcal{M}}_{b}}$
,
$\overline{n}$
are constant values and
$0<r\ll 1$
is a small parameter. If one notes
$\overline{\unicode[STIX]{x1D6FC}}=4\unicode[STIX]{x03C0}\overline{R}^{3}\overline{n}/3$
,
$\overline{\unicode[STIX]{x1D70C}_{1}}=\overline{n}\overline{{\mathcal{M}}_{b}}/\overline{\unicode[STIX]{x1D6FC}}$
,
$\overline{p_{1}}=p_{1}(\overline{\unicode[STIX]{x1D70C}_{1}})$
,
$\overline{c_{1}}=c_{1}(\overline{\unicode[STIX]{x1D70C}_{1}})$
,
$\unicode[STIX]{x1D708}(\overline{\unicode[STIX]{x1D6FC}})=\overline{\unicode[STIX]{x1D708}}$
and
$\overline{\unicode[STIX]{x1D70C}}=\overline{\unicode[STIX]{x1D6FC}}\overline{\unicode[STIX]{x1D70C}_{1}}+(1-\overline{\unicode[STIX]{x1D6FC}})\overline{\unicode[STIX]{x1D70C}_{2}}$
then (3.13) yield



Identifying same-order terms with respect to
$r$
yields

and

Equation (3.17) is a second-order linear ODE in
$z$
that is consistent with the evolution of a linear harmonic oscillator with damping and forcing terms. This type of equation is classic in the literature for describing the motion of vibrating bubbles in the linear regime. Indeed, following Prosperetti (Reference Prosperetti1977), Cheng, Drew & Lahey (Reference Cheng, Drew and Lahey1983), we have

where
$\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{vis}+\unicode[STIX]{x1D6FE}_{th}+\unicode[STIX]{x1D6FE}_{ac}$
drives the damping intensity of the system. The coefficient
$\unicode[STIX]{x1D6FE}_{vis}$
pertains to viscous effects due to the surrounding liquid and is defined by
$\unicode[STIX]{x1D6FE}_{vis}=2\unicode[STIX]{x1D707}_{2}/(\overline{\unicode[STIX]{x1D70C}_{2}}\,\overline{R}^{2})$
,
$\unicode[STIX]{x1D6FE}_{th}$
is related to thermal exchanges between the gas and the liquid and
$\unicode[STIX]{x1D6FE}_{ac}$
concerns acoustic scattering by the bubbles. Expressions for
$\unicode[STIX]{x1D6FE}_{th}$
and
$\unicode[STIX]{x1D6FE}_{ac}$
are available in Prosperetti (Reference Prosperetti1977), Cheng et al. (Reference Cheng, Drew and Lahey1983), and references from previous studies therein. They involve several intermediate parameters but also characteristic parameters of the forcing term, like the frequency of the perturbation
$\unicode[STIX]{x1D6FF}p_{2}$
. Identifying terms in (3.17) and (3.18) yields the following definitions for
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D700}$

Thus we see that the above analysis and the resulting relation (3.19) provide a definition for
$\unicode[STIX]{x1D700}$
that is different from (3.12). More specifically, as
$\unicode[STIX]{x1D700}_{\mathit{Lin}}>\unicode[STIX]{x1D700}_{\mathit{RP}}$
we can see that (3.19) yields greater damping than (3.12). The discrepancy between
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$
and
$\unicode[STIX]{x1D700}_{\mathit{RP}}$
can be explained by simplifying hypotheses at the core of model (2.30) and our simple bubbly fluid model. Indeed, (H3) does not allow us to describe thermal exchange in the fluids. Also, we have completely neglected pressure fluctuations within the bubbles although they make an important contribution to damping effects.
In the following, we will rely on the following relations:
(i) Pfriem’s expression that can be found in Cheng et al. (Reference Cheng, Drew and Lahey1983) for
$\unicode[STIX]{x1D6FE}_{th}$ ,
(3.20)where$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{th}^{P40}(\unicode[STIX]{x1D714})=\unicode[STIX]{x1D714}_{n}\,\frac{3(\unicode[STIX]{x1D6FE}_{1}-1)\sqrt{2a_{1}}}{2\sqrt{\unicode[STIX]{x1D714}}R}, & & \displaystyle\end{eqnarray}$$
$a_{1}$ is the thermal diffusivity of the gas and
$\unicode[STIX]{x1D6FE}_{1}$ its ratio of specific heats.
(ii) The natural frequency
$\unicode[STIX]{x1D714}_{n}$ is given by the relation from Cheng et al. (Reference Cheng, Drew and Lahey1983),
(3.21)with$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{n}^{2}={\displaystyle \frac{3\unicode[STIX]{x1D705}_{1}p_{1}}{\unicode[STIX]{x1D70C}_{1}R^{2}}}, & & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D705}_{1}$ being the polytropic gas exponent.
(iii) For
$\unicode[STIX]{x1D6FE}_{ac}$ , we consider the relation found in Prosperetti (Reference Prosperetti1977),
(3.22)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{ac}^{P77}(\unicode[STIX]{x1D714})=0.5\,{\displaystyle \frac{\unicode[STIX]{x1D714}^{2}Rc_{2}}{c_{2}^{2}+(\unicode[STIX]{x1D714}R)^{2}}}. & & \displaystyle\end{eqnarray}$$
Finally, for practical purposes, we want to get rid of the dependence of these damping parameters on the frequency of the considered acoustic perturbation. A way to get simple and constant values for these damping effects is to evaluate expressions (3.20) and (3.22) at the natural frequency (3.21),

3.3 Dispersion relations for a plane and monochromatic wave
In order to test the relevance of our models and of the identification made for
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
, we compare the behaviour of systems (2.30), (2.34) and (2.36) in the acoustic regime, to experimental measures of sound waves dispersion in bubbly fluids. The linearization of systems (2.30), (2.34) and (2.36) and their acoustic properties are presented hereafter. The practical comparison with experimental data and reference dispersion relation will be given in § 4.
Considering smooth solutions of one-dimensional problems, all these systems can be expressed using the generic quasilinear form

Following standard lines (Whitham Reference Whitham1974; Burman & Sainsaulieu Reference Burman and Sainsaulieu1995), we seek a monochromatic wave solution of (3.24) by writing
$\boldsymbol{W}$
in the form

where
$\unicode[STIX]{x1D714}$
is the angular frequency,
$k$
the wavelength and
$r$
is a small amplitude parameter. The states
$\widehat{\boldsymbol{W}}^{(1)}$
and
$\boldsymbol{W}^{(0)}$
are both constant. The fluid parameters involved with
$\boldsymbol{W}^{(0)}$
are noted with the superscript
$^{(0)}$
and for the sake of simplicity, we suppose that
$\boldsymbol{W}^{(0)}$
is always a rest state, i.e.
$u^{(0)}=0$
.
Injecting (3.25) into (3.24) and identifying terms with respect to the powers of
$r$
yields

Consequently
$\unicode[STIX]{x1D714}$
and
$k(\unicode[STIX]{x1D714})$
are bound by the so-called dispersion relation

which allows us to defined the phase velocity and the spatial attenuation of the acoustic wave respectively by
$\text{Re}[\unicode[STIX]{x1D714}/k(\unicode[STIX]{x1D714})]$
and
$\text{Im}[k(\unicode[STIX]{x1D714})]$
. Now, let us detail the results for each system of our hierarchy. Let us note

(i) For the two-phase model with micro-inertia (2.30), the dispersion relation, the associated phase velocity
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ and the spatial attenuation
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ read
(3.29a )$$\begin{eqnarray}\displaystyle \left({\displaystyle \frac{k^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})}{\unicode[STIX]{x1D714}}}\right)^{2}={\displaystyle \frac{\unicode[STIX]{x1D708}\unicode[STIX]{x1D714}^{2}-\text{i}\unicode[STIX]{x1D700}\unicode[STIX]{x1D714}-c_{\mathit{Wood}}^{-2}\mathsf{H}(\unicode[STIX]{x1D70C}^{(0)},Y^{(0)},\unicode[STIX]{x1D6FC}^{(0)})}{\unicode[STIX]{x1D708}c_{\mathit{Frozen}}^{2}\unicode[STIX]{x1D714}^{2}-\text{i}\unicode[STIX]{x1D700}\,c_{\mathit{Frozen}}^{2}\unicode[STIX]{x1D714}-\mathsf{H}(\unicode[STIX]{x1D70C}^{(0)},Y^{(0)},\unicode[STIX]{x1D6FC}^{(0)})}}, & & \displaystyle\end{eqnarray}$$
(3.29b,c )$$\begin{eqnarray}\displaystyle c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})=\text{Re}\left[{\displaystyle \frac{\unicode[STIX]{x1D714}}{k^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})}}\right],\quad \unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})=\text{Im}[k^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})]. & & \displaystyle\end{eqnarray}$$
(ii) For the micro-inertia free model (2.34) obtained by the limit
$\unicode[STIX]{x1D708}\rightarrow 0$ ,
$\unicode[STIX]{x1D700}=O(1)$ , we get the dispersion relation, phase velocity
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$ and attenuation
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700}}$ defined by
(3.30a,b )$$\begin{eqnarray}\displaystyle \left({\displaystyle \frac{k^{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D714})}{\unicode[STIX]{x1D714}}}\right)^{2}={\displaystyle \frac{\text{i}\unicode[STIX]{x1D700}\unicode[STIX]{x1D714}+c_{\mathit{Wood}}^{-2}\mathsf{H}(\unicode[STIX]{x1D70C}^{(0)},Y^{(0)},\unicode[STIX]{x1D6FC}^{(0)})}{\text{i}\unicode[STIX]{x1D700}c_{\mathit{Frozen}}^{2}\unicode[STIX]{x1D714}+\mathsf{H}(\unicode[STIX]{x1D70C}^{(0)},Y^{(0)},\unicode[STIX]{x1D6FC}^{(0)})}},\quad c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D714})=\text{Re}\left[{\displaystyle \frac{\unicode[STIX]{x1D714}}{k^{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D714})}}\right],\qquad \quad & & \displaystyle\end{eqnarray}$$
(3.30c )$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D714})=\text{Im}[k^{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D714})]. & & \displaystyle\end{eqnarray}$$
(iii) Finally, for the full equilibrium model (2.36) when
$\unicode[STIX]{x1D708}\rightarrow 0$ ,
$\unicode[STIX]{x1D700}\rightarrow 0$ , the dispersion relation reads
(3.31a-c )$$\begin{eqnarray}\displaystyle {\displaystyle \frac{k(\unicode[STIX]{x1D714})^{2}}{\unicode[STIX]{x1D714}^{2}}}={\displaystyle \frac{1}{c_{\mathit{Wood}}^{2}}},\quad c_{\mathit{Phase}}(\unicode[STIX]{x1D714})=c_{\mathit{Wood}},\quad \unicode[STIX]{x1D6FD}=0. & & \displaystyle\end{eqnarray}$$
We recall that
$c_{\mathit{Wood}}$
and
$c_{\mathit{Frozen}}$
are defined by (2.37) and (2.32).
Let us first remark a distinctive behaviour of the acoustic waves for the system (2.30) equipped with both damping and micro-inertia: in the case of low internal dissipation i.e. small values of
$\unicode[STIX]{x1D700}$
, the dispersion relation (3.29) leads to resonance in the vicinity of the frequency

We observe that when one accounts for internal damping with
$\unicode[STIX]{x1D700}>0$
and with micro-inertia (respectively without micro-inertia), the phase velocity of the acoustic wave
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
(respectively
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$
) is not equal to the sound velocity
$c$
(respectively
$c_{\mathit{Frozen}}$
) issued from the characteristic velocities when one discards the source terms in the system. This underlines the fact that the micro-inertia and the damping source terms have a substantial influence on the phase velocity of the acoustic waves.

Figure 1. Phase velocity and spatial attenuation for the 5-equation model and influence of the value of
$\unicode[STIX]{x1D708}$
.
Let us now examine the variations of the dispersion relations across the hierarchy. As
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D700}$
reach their asymptotic limit, the transition from one model to another materializes through the dispersion relations (3.29), (3.30) and (3.31). Indeed, we see that

Let us note that these limits are not uniform over all frequencies: this is illustrated in figure 2. Indeed, we can observe that, whatever the value of
$\unicode[STIX]{x1D700}$
, there is a frequency above which
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$
will be close to
$c_{\mathit{Frozen}}$
. However, this critical frequency increases when
$\unicode[STIX]{x1D700}$
decreases. On the contrary, the transition from
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
to
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$
is more uniform, as illustrated in figure 1. Indeed, when
$\unicode[STIX]{x1D708}$
decreases, the damping effects due to
$\unicode[STIX]{x1D700}$
prevail: although the resonant frequency only depends on
$\unicode[STIX]{x1D708}$
, the effects of resonance are completely attenuated because of
$\unicode[STIX]{x1D700}$
. One can also note that the asymptotic behaviour for spatial attenuation at low frequencies does not depend on the value of
$\unicode[STIX]{x1D708}$
, but depends on the value of
$\unicode[STIX]{x1D700}$
, see figure 2 for the 4-equation model.
The model hierarchy also shows through when one spans frequency values
$\unicode[STIX]{x1D714}$
(and can also be noticed in figures 1 and 2). Indeed, if one considers acoustic waves at low frequencies
$\unicode[STIX]{x1D714}\ll 1$
, then the dispersion relations (3.29), (3.30) and (3.31) yield

In terms of phase velocities and attenuation we obtain

In the limit
$\unicode[STIX]{x1D714}\rightarrow 0$
, the phase velocity of the acoustic waves for all models tends to
$c_{\mathit{Phase}}=c_{\mathit{Wood}}$
and the spatial attenuation vanishes. Let us now turn to high frequencies
$\unicode[STIX]{x1D714}\gg 1$
. From (3.29), (3.30) and (3.31) we have

Thus the acoustic waves of both systems equipped with internal damping are such that

Consequently, in the limit
$\unicode[STIX]{x1D714}\rightarrow +\infty$
, the phase velocities of the acoustic waves associated with (2.30) and (2.18) tend to
$c_{\mathit{Frozen}}$
, while it remains constant and equal to
$c_{\mathit{Wood}}$
for (2.36). In all cases the spatial attenuation tends to
$0$
.

Figure 2. Phase velocity and spatial attenuation for the 4-equation model and influence of the value of
$\unicode[STIX]{x1D700}$
.
4 Comparison with bubbly fluid reference data
From now on, we shall suppose that the compaction energy is null
$e(\unicode[STIX]{x1D6FC})=0$
, and that the barotropic EOS for each pure fluid has the form

where
$c_{k}$
,
$\unicode[STIX]{x1D70C}_{k}^{0}$
and
$p_{k}^{0}$
are real constants chosen as follows:


Thanks to the dispersion relations (3.29), (3.30) and (3.31), we can now study the responses of systems (2.30), (2.34) and (2.36) in the acoustic regime under a forced pressure oscillation and compare them with experimental results obtained for bubbly fluids. The data we shall use rely on two experimental works by Silberman (Reference Silberman1957) and Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008). Let us first briefly outline the framework of these studies. In the next paragraphs, elements related to the experimental data will be denoted with the superscript
$^{ref}$
.
The data of Silberman (Reference Silberman1957) are considered as reference experimental results in the domain of acoustic wave propagation for bubbly flows. They have been used in comparison with several models in the literature (van Wijngaarden Reference van Wijngaarden1972; Cheng et al.
Reference Cheng, Drew and Lahey1985; Commander & Prosperetti Reference Commander and Prosperetti1989; Drew & Passman Reference Drew and Passman1999). The method proposed in Silberman (Reference Silberman1957) consists in generating standing waves in various length steel pipes. The sound is generated at one end of the pipe while small hydrophones measure sound pressure at the other end. Measurements are performed between two nodes or antinodes that allow us to compute the phase velocity and the spatial attenuation. The size distribution of the bubbles is estimated using photographs. The resulting measurements were very accurate, except near the resonance frequency
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
. Indeed, for
$\unicode[STIX]{x1D714}$
close to
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
evaluation of the phase velocity was not possible due to the severe attenuation of the acoustic waves. In order to obtain data in this range of frequencies, we shall use the work of Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008) that involves acoustic wave propagating within a thin hair gel sample containing air bubbles. The sound waves are produced at one end of the system by a transducer and measurements are performed thanks to a hydrophone at the other end. The advantage of using the gel is that the distribution of bubble radii and volume fractions are accurately known. According to Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008), the difference in terms of acoustic behaviour between water and gel is negligible regarding the wave dispersion. Thanks to this set-up, the results of Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008) provide accurate data for both phase velocity and attenuation in the vicinity of
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
.
4.1 Influence of micro-viscosity and micro-inertia in the acoustic regime
We shall examine the behaviour of the models when
$\unicode[STIX]{x1D714}$
spans the possible frequencies and distinguish three main ranges of frequency for characterizing the phase velocity and the spatial attenuation. Then we will focus specifically on the near-resonance frequencies.
4.1.1 Comparison across the whole spectrum of frequencies
$\unicode[STIX]{x1D714}$
We consider a set of measures from Silberman (Reference Silberman1957) that involves a flow characterized by
$R=2.5\,\text{mm}$
and
$\unicode[STIX]{x1D6FC}=5.84\times 10^{-4}$
. Using relations (3.12) and (3.19), we obtain the following values for the parameters of the 5-equation model (2.30) and the 4-equation model (2.34)


Figure 3 displays both phase velocity and spatial attenuation for all models of the hierarchy, for the model of Cheng, Drew and Lahey (CDL model) (Cheng et al. Reference Cheng, Drew and Lahey1985), superimposed on the experimental results (Silberman Reference Silberman1957).

Figure 3. Dispersion relations for the CDL model and the different models of the hierarchy (numbered lines) and Silberman’s measures (symbols) for radii of bubbles around
$R=2.0\times 10^{-3}~\text{m}$
,
$\unicode[STIX]{x1D6FC}=5.84\times 10^{-4}$
. The model parameters are thus:
$\unicode[STIX]{x1D708}=3.6$
,
$\unicode[STIX]{x1D700}_{\mathit{RP}}=2.3$
,
$\unicode[STIX]{x1D700}_{\mathit{Lin}}=1.6\times 10^{3}$
.
(i) Range
$\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$ . For low frequencies, the pressure perturbation is very slow and we can thus expect the bubbles of the system to remain at an equilibrium state with respect to both mechanics and thermodynamics. Little internal dissipation is involved with this regime, which is visible through the measures that show a low spatial attenuation. The evaluation of
$c_{\mathit{Phase}}^{\mathit{ref}}$ in the experiment provides values that are close to
$c_{\mathit{Wood}}$ . All the models of the hierarchy show a good agreement with these data, as seen in figure 3. As
$\unicode[STIX]{x1D714}$ increases, the spatial attenuation
$\unicode[STIX]{x1D6FD}^{\mathit{ref}}$ also increases. This trend is correctly followed by the models of the hierarchy that account for internal dissipation. Nevertheless, Prosperetti (Reference Prosperetti1977) underlines that, in this regime and up to a certain frequency, the thermal dissipation is the dominant internal dissipation effect. For all models,
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ and
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700}}$ are lower than experimental data. However, if the 5-equation model (2.30) importantly underestimates
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ when the damping is identified with Rayleigh–Plesset
$(\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{RP}})$ , the match with the reference data is very good when
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{Lin}}$ .
(ii) Range
$\unicode[STIX]{x1D714}$ close to
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$ . Near resonance,
$\unicode[STIX]{x1D6FD}^{\mathit{ref}}$ increases with
$\unicode[STIX]{x1D714}$ and becomes very large. On the contrary,
$c_{\mathit{Phase}}^{\mathit{ref}}$ decreases as
$\unicode[STIX]{x1D714}$ increases. It reaches
$c_{\mathit{Phase}}^{\mathit{ref}}=0$ for some frequency
$\unicode[STIX]{x1D714}_{\mathit{ext}}^{\mathit{ref}}<\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$ . For
$\unicode[STIX]{x1D714}\in [\unicode[STIX]{x1D714}_{\mathit{ext}}^{\mathit{ref}},\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}]$ , there is a good agreement between
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ ,
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$ and
$c_{\mathit{Phase}}^{\mathit{ref}}$ . On the contrary, there is an important discrepancy between the phase velocity predicted by the CDL model and
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ ,
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700}}$ for this range of frequencies. It is worth noting that
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ is much closer to the reference value for
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{Lin}}$ than for
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{RP}}$ . Concerning the spatial attenuation, we can see that the models of the hierarchy that account for micro-inertia, namely (2.30) and (2.18) fit quite well the reference results. In contrast, models (2.34) and (2.36) yield a poor estimate of the spatial attenuation. This suggests that the source terms related to
$\unicode[STIX]{x1D708}$ in (2.27d ) and (2.27e ) play a key role in the system behaviour when
$\unicode[STIX]{x1D714}\approx \unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$ . Moreover it also hints that our estimate for
$\unicode[STIX]{x1D708}$ in (3.12) and (4.5) is coherent. Finally, the results suggest that the source terms related to
$\unicode[STIX]{x1D700}$ in (2.18) do not have a great influence on the values of
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$ in this range of frequencies.
(iii) Range
$\unicode[STIX]{x1D714}\gg \unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$ . For high frequencies, very few experimental data are available and thus the main comparison elements are given by the CDL model. In this regime, one can presume that acoustic radiation effects cannot be neglected. Indeed, the bubbles start emitting acoustic waves that are transmitted to the liquid. This process will remove energy from the bubbles and therefore will become the main damping effect of the system. For all systems of our hierarchy except (2.36), the phase velocity will tend to
$c_{\mathit{Frozen}}$ , which agrees with the behaviour of the CDL model. The spatial attenuation coefficients of the model hierarchy do not match well the reference data of the CDL model in this range of frequencies: (2.34) clearly overestimates damping (purple number 2 line), when model (2.30) with matching coefficients
$\unicode[STIX]{x1D708}_{\mathit{RP}}$ and
$\unicode[STIX]{x1D700}_{\mathit{RP}}$ provides too low a dissipation (red number 3 line). For the same model, the micro-viscosity choice
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$ clearly increases the damping effect but still at a much lower level than the dissipation of the CDL model.
4.1.2 Finer comparison near resonance
In the experiment of Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008), the set of measures is very dense for
$\unicode[STIX]{x1D714}$
close to
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
. In this paragraph we discard both the 4-equation model (2.34) and the 3-equation model (2.36) as they cannot produce resonant behaviour. The bubbles in Leroy et al. (Reference Leroy, Strybulevych, Page and Scanlon2008) are smaller than those of Silberman (Reference Silberman1957). Thus, we consider different values of
$(R,\unicode[STIX]{x1D6FC})$
by setting
$R\approx 8.1\times 10^{-5}~\text{m}$
and
$\unicode[STIX]{x1D6FC}=1.5\times 10^{-4}$
. Thanks to (3.12) and (3.19) we obtain


The results we obtain with this set of parameters is consistent with the previous comparison. Indeed, in figure 4 we can see that for
$\unicode[STIX]{x1D714}$
close to
$\unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
the phase velocity
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
is clearly overestimated for
$\unicode[STIX]{x1D700}_{\mathit{RP}}$
but the match with experimental data for
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$
seems more accurate. Regarding the attenuation
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
, the choice of
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$
gives clearly a better match with the reference data than
$\unicode[STIX]{x1D700}_{\mathit{RP}}$
.

Figure 4. Dispersion relations for the different models from the hierarchy (numbered lines) and the Leroy et al. measures (symbols) for radii of bubbles around
$R=8.1\times 10^{-5}~\text{m}$
,
$\unicode[STIX]{x1D6FC}=1.5\times 10^{-4}$
. The model parameters are thus:
$\unicode[STIX]{x1D708}=1.5\times 10^{-3}$
,
$\unicode[STIX]{x1D700}_{\mathit{RP}}=8.9\times 10^{-1}$
,
$\unicode[STIX]{x1D700}_{\mathit{Lin}}=8.7\times 10^{1}$
.
4.2 On the evaluation of the micro-viscosity and the micro-inertia
We succeeded in identifying the micro-inertia
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{\mathit{Lin}}=\unicode[STIX]{x1D708}_{\mathit{RP}}$
and the micro-viscosity
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{RP}}$
(respectively
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{Lin}}$
) in the 5-equation model (2.30) thanks to a comparison with the Rayleigh–Plesset equation (3.10) (respectively the linear radius evolution equation (3.17)). In this sense the 5-equation model (2.30) can be considered as a host for different bubbly fluid models that will characterize the micro-inertia and micro-viscosity.
In our case, the values of
$\unicode[STIX]{x1D700}_{\mathit{RP}}$
and
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$
are significantly different and for a given
$\unicode[STIX]{x1D714}\geqslant \unicode[STIX]{x1D714}_{\mathit{res}}^{\mathit{ref}}$
,
$\unicode[STIX]{x1D700}_{\mathit{RP}}$
and
$\unicode[STIX]{x1D700}_{\mathit{Lin}}$
yield different acoustic behaviours for system (2.30). Nevertheless, the agreement with reference data is better for
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D700}_{\mathit{Lin}}$
. We believe that this discrepancy comes from the over-simplifying hypotheses of § 3.1: (H6) and (H7) do not allow us to account for pressure fluctuations within the bubble, neither for pressure fluctuations in the liquid that act as inertial terms with respect to the bubble motion. These phenomena are the grounds for thermal and acoustic damping effects, that are respectively driven by
$\unicode[STIX]{x1D6FE}_{ac}$
and
$\unicode[STIX]{x1D6FE}_{th}$
(Prosperetti Reference Prosperetti1977; Cheng et al.
Reference Cheng, Drew and Lahey1983). Indeed, what we observed in § 4.1 is consistent with the analysis of Prosperetti (Reference Prosperetti1977): for the flows settings we chose, both thermal and acoustic damping are the main damping effects involved in the bubble vibrations. According to Prosperetti (Reference Prosperetti1977), when the radius is close to the value of Silberman’s experiment,
$R=1~\text{mm}$
, one can estimate that
$\unicode[STIX]{x1D6FE}_{th}/\unicode[STIX]{x1D6FE}_{vis}\in [5,5000]$
for
$\unicode[STIX]{x1D714}\leqslant 10^{5}~\text{s}^{-1}$
and
$\unicode[STIX]{x1D6FE}_{ac}/\unicode[STIX]{x1D6FE}_{vis}\in [5,5000]$
for
$\unicode[STIX]{x1D714}\geqslant 5\times 10^{3}~\text{s}^{-1}$
. When
$R=8.1\times 10^{-2}~\text{mm}$
, a value close to the bubbles radii in the Leroy et al. experiment, we have that
$\unicode[STIX]{x1D6FE}_{th}/\unicode[STIX]{x1D6FE}_{vis}\in [10,500]$
for
$\unicode[STIX]{x1D714}\leqslant 10^{6}~\text{s}^{-1}$
and
$\unicode[STIX]{x1D6FE}_{ac}/\unicode[STIX]{x1D6FE}_{vis}\in [3,500]$
for
$\unicode[STIX]{x1D714}\geqslant 10^{5}~\text{s}^{-1}$
.
Finally, we saw that, by increasing the value of
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
, it is possible to enforce an equivalent damping
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}\approx \unicode[STIX]{x1D6FD}_{th}$
that better fits the reference data in this range of frequencies, especially close to resonance. Thus, this enables an alternate means for determining
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
: one could tune
$(\unicode[STIX]{x1D700},\unicode[STIX]{x1D708})$
in a heuristic way to better fit
$c_{\mathit{Phase}}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
and
$\unicode[STIX]{x1D6FD}^{\unicode[STIX]{x1D700},\unicode[STIX]{x1D708}}$
with respect to reference data, and thus use the acoustic regime behaviour of model (2.30) as an evaluation tool.
5 Conclusion
In this work, we have derived a model for a compressible barotropic two-phase medium that accounts for both bulk and sub-scale vibrational kinematic phenomena. A notable feature of this model is that it is agnostic with respect to the considered regime, in the sense that there is no a priori assumption concerning the topology of the interfaces. The model may host the description of either separated-phase or disperse-phase flows. Following classic modelling guidelines, we have used Hamilton’s principle of stationary action and the elaboration of an entropy budget, to obtain respectively the conservative part and the dissipative structures of our model.
The resulting system is a 5-equation model whose convective part is hyperbolic. This system features two parameters:
$\unicode[STIX]{x1D700}$
that is related to internal dissipation effects and
$\unicode[STIX]{x1D708}$
that pertains to small-scale kinematic effects. Two reduced models can be obtained by considering the regimes
$\unicode[STIX]{x1D708}\rightarrow 0$
,
$\unicode[STIX]{x1D700}=O(1)$
and
$\unicode[STIX]{x1D708}\rightarrow 0$
,
$\unicode[STIX]{x1D700}\rightarrow 0$
. These limit regimes lead to models from the literature that have been used for describing compressible separated-phase flows (Chanteperdrix et al.
Reference Chanteperdrix, Villedieu and Vila2002; Bernard-Champmartin & De Vuyst Reference Bernard-Champmartin and De Vuyst2014). In this sense, we obtained a hierarchy of three compressible two-phase flow models.
As the 5-equation model is neutral with respect to the topology of the flow, it can be equipped with extra hypotheses to qualify the flow structure: the model has been analysed as a simple bubbly fluid model. In this context, we have shown that this 5-equation model is compatible with a sub-scale vibrational motion of bubbles, governed by the Rayleigh–Plesset equation in a linear or a nonlinear regime. These connections have allowed us to identify the parameters
$\unicode[STIX]{x1D700}$
and
$\unicode[STIX]{x1D708}$
with micro-viscosity and micro-inertial effects respectively. More precisely, it has enabled two possible evaluations for
$(\unicode[STIX]{x1D700},\unicode[STIX]{x1D708})$
: a first evaluation that is based on the nonlinear Rayleigh–Plesset equation and a second that uses the linearized Rayleigh–Plesset equation. The second evaluation suggests a way to compensate for neglected effects that may even be predominant in some flow configurations.
Finally, we have considered the acoustic regime, in order to compare the behaviour of monochromatic acoustic waves between each model of the hierarchy, a reference two-phase model proposed in Cheng et al. (Reference Cheng, Drew and Lahey1983) and some experimental data. This study has allowed us to further discriminate the domain of validity of the models of the hierarchy and their ability to obtain a resonance regime. Good matches with the reference data were obtained with the richest model. This validates the evaluations of
$(\unicode[STIX]{x1D700},\unicode[STIX]{x1D708})$
in this particular situation and shows that rich acoustic phenomena can be obtained with a single bulk velocity model.
The present study is a first building block to a larger effort to merge separated- and disperse-phase models. The present results leave now many open questions. First, many effects like gravity, thermal effects and capillary effects have been neglected in our study and obviously need to be accounted for. Next, we believe that supplementary information related to small scales must be added to the modelling process in order to enrich the small-scale description. Such an extension is under investigation for pulsating objects, that are homeomorphic to spheres and described by their surface and averaged curvatures (Cordesse et al. Reference Cordesse, Kokh, Di Battista, Drui and Massot2018). Finally, even though we have shown that the resulting models are compatible with different flow regimes, we did not provide any element that allows different regimes to interact within the same model. Addressing this question requires us to characterize the different regimes within the same system. This matter is the subject of ongoing work.
Acknowledgements
The help in the writing process of the paper and insightful comments of M. Lance and D. Legendre are gratefully acknowledged. The PhD of F.D. is funded by a CEA/DGA (Direction Générale de l’Armement – French Department of Defense) grant.
Appendix A. A bubbly flow model
We recall hereafter the two-phase bubbly system studied in Drew (Reference Drew1983), Cheng et al. (Reference Cheng, Drew and Lahey1985). For the sake of simplicity, we shall only consider one-dimensional problems. This model is composed of balance equations for mass, momentum and energy for each fluid. These equations are derived following the lines of Ishii (Reference Ishii1975), using averaging methods. Constitutive equations are also used for closure purposes, and finally, an additional interaction law relates the pressures of each component and closes definitely the system. The derivation of this interaction law is presented below. It is based on the same approach as the one used to derive the Rayleigh–Plesset equation (Lord Rayleigh Reference Lord Rayleigh1917; Plesset & Prosperetti Reference Plesset and Prosperetti1977).
First, consider the motion of a single spherical bubble of radius
$R$
: one makes the assumption that the bubble undergoes oscillations that are driven by a velocity potential
$\unicode[STIX]{x1D711}$
of the form

where
$r$
is the distance to the centre of the bubble (see Landau & Lifshitz (Reference Landau and Lifshitz1966)) and
$k$
is the wavenumber associated with the disturbance of the velocity field coming from the surrounding liquid. The evolution of
$R$
is derived by supposing that, in comparison with the gas inside the bubble, the liquid is almost incompressible and that the pressure far from the bubble interface is
$p_{2}$
. Then one supposes the dynamics of the bubble to verify the Bernoulli equation as follows:

where
$p_{2i}$
is the pressure of the liquid at the interface between the two fluids and
$\unicode[STIX]{x1D70C}_{2}$
is the density of the liquid. Moreover, one supposes that the motion of the bubble is constrained by the Laplace relation

where
$\unicode[STIX]{x1D70E}$
is the surface tension,
$\unicode[STIX]{x1D707}_{2}$
is the dynamic viscosity of the liquid and
$p_{1i}$
is the pressure of the gas at the bubble boundary. Supposing the amplitude of the bubble oscillations to be small, (A 3) and (A 2) are complemented by an additional relation that connects the gas interfacial pressure variations
$\unicode[STIX]{x1D6FF}p_{1i}$
to the radius variations
$\unicode[STIX]{x1D6FF}R$
. This relation accounts for thermal effects occurring at the interface and also for the thermodynamic properties of the gas. It involves complex expressions that will not be detailed here, we refer the reader to Prosperetti (Reference Prosperetti1977), Cheng et al. (Reference Cheng, Drew and Lahey1983) for a detailed view on this topic.
The other part of the model studied in Cheng et al. (Reference Cheng, Drew and Lahey1985) shows six bulk balance equations (A 4a
)–(A 4f
) detailed below. They characterize the evolution of the system parameters at the macroscopic scale and are derived by applying an ensemble averaging (Drew & Passman Reference Drew and Passman1999). Let us note by
$\unicode[STIX]{x1D70C}_{q}$
,
$u_{q}$
,
$p_{q}$
,
$h_{q}$
, respectively the averaged values of density, velocity, partial pressure and specific enthalpy of the fluid
$q=1,2$
. Each fluid is supposed to be a compressible material that is equipped with a pressure law of the form
$(\unicode[STIX]{x1D70C}_{k},h_{k})\mapsto p_{k}$
. In the case of a dispersed bubbly flow, the volume fraction of gas is defined by setting
$\unicode[STIX]{x1D6FC}=4\unicode[STIX]{x03C0}nR^{3}/3$
, where
$n$
is the bubble number density. The partial masses are given by
$m_{q}=\unicode[STIX]{x1D70C}_{q}\unicode[STIX]{x1D6FC}_{q}$
, where
$\unicode[STIX]{x1D6FC}_{1}=\unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D6FC}_{2}=1-\unicode[STIX]{x1D6FC}$
. Neglecting wall-shear effects and gravity, for one-dimensional problems the system reads as follows:













Finally, the system (A 4) has to be complemented with an evolution equation for the number density
$n$
. If one supposes that no coalescence or break-up occurs,
$n$
verifies the conservation equation

Appendix B. Eigenstructure of the two-phase model with micro-inertia
We consider the sole convective part of system (2.27) for one-dimensional problems by discarding the source terms. For smooth solutions, the obtained system may be expressed using the variable
$\boldsymbol{V}=(\unicode[STIX]{x1D70C},u,Y,\unicode[STIX]{x1D6FC},w)^{\text{T}}$
as follows:

The matrix
$A(\boldsymbol{V})$
possesses three distinct eigenvalues:
$u\pm c$
and
$u$
associated respectively with the eigenvectors

