1. Introduction
1.1. Background
The modelling of compressible multiphase fluid mixtures is fundamental to a wide range of scientific and engineering applications, including geophysical flows, aerospace engineering and industrial processes. A key challenge in describing such systems lies in formulating a mathematical framework that consistently captures the interactions between phases while accounting for compressibility effects, diffusion and interfacial dynamics. (We utilise ‘phase’ to denote different fluid constituents (e.g. air and water).) A widely used class of partial differential equation (PDE)-based models for multiphase flows consists of sharp-interface models, which either explicitly track the phase boundaries and prescribe conditions at the interface (Ghias, Mittal & Dong Reference Ghias, Mittal and Dong2007; Shen, Ren & Ding Reference Shen, Ren and Ding2020) or are approximated by smooth interface (level-set) approximations (Sussman, Smereka & Osher Reference Sussman, Smereka and Osher1994; Sethian Reference Sethian1999; ten Eikelder & Akkerman Reference ten Eikelder and Akkerman2021). These models have been successfully applied to problems where the interface between phases remains well defined. Complex mixture problems, however, where phases are dispersed throughout the domain and interpenetrate one another, or emerge from phase separation, do not belong to the class of sharp-interface problems. As a consequence, these problems demand an alternative modelling approach. Below we discuss two classes of models that describe these complex mixtures: compressible two-phase models and phase-field models.
Compressible two-phase flow models play a fundamental role in describing mixture flows, where compressibility and interfacial dynamics are essential. Among these models, the Baer-Nunziato model is often regarded as one of the most comprehensive frameworks for describing two-phase flows (Baer & Nunziato Reference Baer and Nunziato1986). This diffuse-interface, non-equilibrium model consists of separate mass, momentum and energy balance equations for each phase, along with an additional equation governing the topology of the two-fluid interface. Over time, significant extensions have been proposed, including generalisations to
$N$
-component mixtures (see, e.g. Müller et al. Reference Müller, Hantke and Richter2016), broadening the applicability to more complex multiphase systems. We refer for details on the analysis and derivation of the Baer-Nunziato model to Hillairet (Reference Hillairet2019), Perrier & Gutiérrez (Reference Perrier and Gutiérrez2021) and, for numerical methods, to Andrianov & Warnecke (Reference Andrianov and Warnecke2004), Tokareva & Toro (Reference Tokareva and Toro2010), Crouzet et al. (Reference Crouzet, Daude, Galon, Helluy, Hérard, Hurisse and Liu2013), Coquel, Hérard & Saleh (Reference Coquel, Hérard and Saleh2017), Sirianni, Re & Abgrall (Reference Sirianni, Re and Abgrall2025).
Despite their generality, Baer-Nunziato-type models introduce considerable computational complexity due to their large number of equations and wave families. To address this, reduced models with fewer equations have been developed (Kapila et al. Reference Kapila, Menikoff, Bdzil, Son and Stewart2001; Murrone & Guillard Reference Murrone and Guillard2005). These reduced formulations are often derived through relaxation procedures that assume velocity, pressure or temperature equilibrium between phases (Saurel & Abgrall Reference Saurel and Abgrall1999). However, such simplifications come at a cost: reduced models may fail to strictly conserve energy, suffer from an ill-posed mixture speed of sound (Saurel, Petitpas & Berry Reference Saurel, Petitpas and Berry2009) or lose asymptotic consistency with the full nonlinear system (Saleh & Seguin Reference Saleh and Seguin2020). These limitations motivate the exploration of alternative modelling approaches for compressible mixtures, particularly those that incorporate interfacial dynamics and phase transition effects in a thermodynamically consistent manner.
Over the past decades, phase-field modelling has become a valuable tool for moving boundary problems in computational fluid mechanics (Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler1998; Gomez & van der Zee Reference Gomez and van der Zee2018). It addresses two key challenges: geometrical representation and physical modelling. The phase-field variable naturally captures topological changes, providing advantages over traditional interface-tracking methods for complex flows. Additionally, phase-field models offer a thermodynamically consistent framework for capillary effects, phase transitions and material composition. They also provide an effective description of wetting, making them well suited for complex wetting phenomena (Aland & Mokbel Reference Aland and Mokbel2021; Bhopalam, Bueno & Gomez Reference Bhopalam, Bueno and Gomez2022; Demont et al. Reference Demont, Stoter, Diddens and van Brummelen2024).
The state-of-the art phase-field models for fluids may roughly be classified into two categories. The first category describes single-fluid flow with phase change (Gomez et al. Reference Gomez, Hughes, Nogueira and Calo2010; Bresch, Gisclon & Lacroix-Violet Reference Bresch, Gisclon and Lacroix-Violet2019) through the Navier–Stokes Korteweg (NSK) model. The two phases correspond to the liquid and gaseous states of the same component. In the NSK model the density field itself acts as a phase-field variable. Areas with low density identify as vapour and high-density regions as liquid. Liquid–vapour transitions are often incorporated in the NSK model by means of the van der Waals capillarity model, named after the 1910 Nobel Laureate in physics (Van der Waals Reference Van der Waals1894).
The second category consists of models that describes viscous (incompressible) mixture flow through the Navier–Stokes Cahn–Hilliard (Allen–Cahn) (NSCHAC) model. The origin of NSCHAC models may be traced back to 1977 when Hohenberg and Halperin proposed a NSCHAC model via the coupling between the Navier–Stokes equations, describing viscous fluid flow, and the Cahn–Hilliard equation, describing spinoidal decomposition (Hohenberg & Halperin Reference Hohenberg and Halperin1977). The most important limitation of this model is that the density of the fluids are assumed equal. This precludes the applicability of the model to a large number of practical problems such as ink-droplet dynamics. Since the end of the last century, a number of efforts of extending the model to the case of non-matching densities have been made, e.g. Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998), Boyer (Reference Boyer2002), Ding, Spelt & Shu (Reference Ding, Spelt and Shu2007), Abels, Garcke & Grün (Reference Abels, Garcke and Grün2012), Shen, Yang & Wang (Reference Shen, Yang and Wang2013), Aki et al. (Reference Aki, Dreyer, Giesselmann and Kraus2014). Each of these works aims to describe the same physics, yet the proposed models are different. The occurrence of a many NSCH(AC) models shows that consensus in the realm of single-phase multicomponent flows has been missing. Note that this is in contrast to single-fluid flow with phase change, where the NSK model is the widely adopted model. In two recent articles we have proposed a unified mixture-theory framework of NSCHAC models as a resolution to this open problem (see ten Eikelder et al. (2023), ten Eikelder (Reference Brunk and ten Eikelder2025) for the theoretical framework and ten Eikelder & Schillinger (Reference ten Eikelder and Schillinger2024), Brunk & ten Eikelder (Reference Brunk and ten Eikelder2025) for supporting benchmark computations). The framework indicates a single NSCHAC model that is invariant to the choice of fundamental variables.
1.2. Objective and main results
Although significant progress has been made in modelling multiphase flows, a general phase-field theory for compressible mixtures remains largely undeveloped, despite some important contributions (see, e.g. Málek & Rajagopal (Reference Málek and Rajagopal2008); Freistühler & Kotschote (Reference Freistühler and Kotschote2017); Huang & Johnsen (Reference Huang and Johnsen2023, Reference Huang and Johnsen2024); Mukherjee & Gomez (Reference Mukherjee and Gomez2024) and the references therein). While compressible two-phase flow models exist and provide a well-established framework for describing compressible mixtures, they differ fundamentally from phase-field models in both structure and formulation. Phase-field models, on the other hand, have been extensively used for capturing interfacial dynamics in incompressible fluid mixtures but have not been systematically extended to fully compressible regimes. As a result, the connection between compressible two-phase flow models and phase-field approaches is largely unresolved. Although both modelling strategies aim to describe similar physical processes, their underlying formulations remain distinct. This disconnect has hindered the development of a unified theoretical framework that accommodates both compressibility effects and the diffuse-interface nature of phase-field models.
In this paper we take an initial step toward expanding the theoretical foundation of compressible mixture modelling by developing a thermodynamically consistent theory for compressible, isothermal
$N$
-phase mixtures. As part of this development, we draw connections to both phase-field models and compressible two-phase flow models, outlining key points of overlap and distinction. Our approach introduces a set of
$N$
mass balance equations coupled with a single momentum balance law. In particular, we derive the following multiphase-field model for constituents
${\alpha } = 1, \ldots , N$
:
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \sum _{{\beta }} \tilde {\rho }_{\beta }\boldsymbol{\nabla }\mu _{\beta } - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) - \textrm {div} \left (\sum _{{\beta }} \boldsymbol{M}_{{\alpha }{\beta }}\boldsymbol{\nabla }\mu _{\beta }\right ) + \sum _{{\beta }} m_{{\alpha }{\beta }} \mu _{\beta } &=0, \end{align}
\begin{align} \mu _{\alpha } - \dfrac {\partial \varPsi }{\partial \tilde {\rho }_{\alpha }} + \textrm {div} \left ( \dfrac {\partial \varPsi }{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right )&=0. \end{align}
Here
$\boldsymbol{v}$
is the fluid velocity,
$\tilde {\rho }_{\alpha }$
are the partial densities,
$\rho = \sum _{\beta }\tilde {\rho }_{\beta }$
is the mixture density,
$\boldsymbol{g}$
is the force vector,
$\nu$
is the dynamical viscosity,
$\nu \lambda$
is the second viscosity coefficient and
$\boldsymbol{D}$
is the symmetric velocity gradient. Furthermore,
$\varPsi$
is the free energy,
$\mu _\alpha$
are constituent chemical potentials, and
$\boldsymbol{M}_{{\alpha }{\beta }}$
and
$m_{{\alpha }{\beta }}$
are mobility parameters. Equation (1.1a
) describes the mixture momentum equation, (1.1b
) are the constituent mass balance laws and (1.1c
) define the chemical potentials.
A distinguishing feature of the model, in contrast to compressible two-phase flow models, is its preservation of thermodynamic consistency despite its reduced complexity (i.e. the combination of
$N$
mass balance laws with just a single momentum balance law). Furthermore, the model connects to phase-field models and distinguishes itself from compressible two-phase flow models as follows. For specific closure models, it extends the NSK framework to multiple fluids, where the last two terms in (1.1b
) vanish in the single-fluid case. Additionally, under particular closure choices, these terms take the form of Cahn–Hilliard and Allen–Cahn-type contributions, establishing a direct link to NSCH(AC) models. Key distinctions from existing compressible two-phase flow models lie in (i) the absence of an explicit evolution equation for the volume fraction, and (ii) the formulation of the equation of state, as given in (1.1c
).
We note that the proposed model offers the description of compressible spinodal-decomposition phenomena, which may occur in phase-transforming mixtures with significant density changes. A representative case is a water–air system undergoing vapour–liquid transformation, covering both condensation and evaporation (e.g. fog formation, evaporating droplets). A single mixture Helmholtz free energy encodes a van der Waals-type equation of state for the vapour–liquid transition and a demixing/mixing effect between the water and air components. The resulting thermodynamic quantities (pressure, chemical potentials, Gibbs free energy) determine whether phase change and separation/mixing occur. This contrasts with compressible two-phase flow models (e.g. Baer–Nunziato (BN)), where phase change is introduced via mass-transfer source terms proportional to chemical-potential differences. In the proposed model, mass-transfer source terms occur via Allen–Cahn contributions representing a chemical reaction.
Beyond the above example, a central advantage of the proposed framework in comparison with compressible two-phase flow models is that all thermodynamic driving forces are derived from a single Helmholtz free energy. As a consequence, pressure(s), chemical potentials, capillarity/Korteweg-type stresses, Cahn–Hilliard diffusive fluxes and Allen–Cahn-type mass-transfer contributions are mutually consistent and satisfy the same equilibrium conditions, which reduces ambiguity in closure choices and avoids internal inconsistencies that can arise when these ingredients are specified independently. This unified structure is particularly beneficial in the compressible multiphase setting, where constituent interactions influence both phase separation and mixture-acoustic response. A practical implication is that it facilitates the design of robust discretisations that respect the underlying energy structure, thereby reducing spurious interfacial forces. At the same time, the present framework is isothermal and adopts a single-velocity mixture description; therefore, it is not intended for regimes with significant slip or detailed thermal/latent-heat effects. In addition, the diffuse-interface approach implies an interfacial thickness set by the Helmholtz free energy. As such, reliable simulations must resolve this length scale to obtain accurate capillary forces and grid-independent results.
1.3. Plan of the paper
The remainder of the paper is structured as follows. In § 2 we present the continuum theory of rational mechanics for compressible isothermal fluid mixtures. Next, in § 3 we perform constitutive modelling via the Coleman–Noll procedure. Then, in § 4 we discuss some properties of the model. In § 5 we study the associated first-order system. Next, in § 6 we discuss the case of binary mixtures. Subsequently, in § 7 we discuss connections to existing models. Finally, in § 8 we conclude and provide further research directions.
2. Continuum mixture theory
This section outlines the continuum theory for mixtures of compressible, isothermal constituents, serving as a foundational framework. While the theory is well known, its inclusion here is essential, particularly as it differs from the starting point for the design of compressible two-phase flow models. The outlined mixture-theory description is closely aligned with ten Eikelder, Van Der Zee & Schillinger (Reference ten Eikelder, Van Der Zee and Schillinger2024), ten Eikelder (Reference ten Eikelder2025), underscoring its suitability as a first-principles basis for the development for a wide range of models.
Continuum mixture theory is based on the following three general principles outlined in the groundbreaking work of Truesdell & Toupin (Reference Truesdell and Toupin1960), Truesdell (Reference Truesdell1984).
-
(i) All properties of the mixture must be mathematical consequences of properties of the constituents.
-
(ii) So as to describe the motion of a constituent, we may in imagination isolate it from the rest of the mixture, provided we allow properly for the actions of the other constituents upon it.
-
(iii) The motion of the mixture is governed by the same equations as is a single body.
The first principle explains that the mixture is formed by its individual constituents. The second principle describes that the components of the physical model are interconnected through interaction terms. Lastly, the third principle states that the movement of the mixture cannot be distinguished from that of a single fluid.
In § 2.1 we provide the groundwork for the continuum theory of mixtures, covering essential kinematics. Following this, § 2.2 presents the evolution laws for individual constituents and mixtures.
2.1. Preliminaries and kinematics
The central concept of the continuum theory of mixtures posits that the material body consists of
$N$
constituent bodies
$\mathscr{B}_{\alpha }$
, where
$\alpha$
ranges from 1 to
$N$
. These constituent bodies
$\mathscr{B}_{\alpha }$
can simultaneously occupy the same region in space. Let
$\boldsymbol{X}_{\alpha }$
represent the spatial position of a particle of
$\mathscr{B}_{\alpha }$
in the Lagrangian (reference) configuration. The spatial position of a particle is described by the invertible deformation map:
\begin{align} \boldsymbol{x} := \boldsymbol{\chi }_{{\alpha }}(\boldsymbol{X}_{{\alpha }},t). \end{align}
The constituent partial mass density
$\tilde {\rho }_{{\alpha }}$
and specific mass density
$\tilde {\rho }_{\alpha }\gt 0$
are respectively defined as
\begin{align} \tilde {\rho }_{{\alpha }}({\boldsymbol{x}},t) :&= \displaystyle \lim _{ \vert V \vert \rightarrow 0} \dfrac {M_{{\alpha }}(V)}{\vert V \vert }, \end{align}
\begin{align} \rho _{\alpha }({\boldsymbol{x}},t) :&= \displaystyle \lim _{\vert V_{{\alpha }}\vert \rightarrow 0} \dfrac {M_{{\alpha }}(V)}{\vert V_{{\alpha }}\vert }, \end{align}
where
$V \subset \varOmega$
(measure
$\vert V \vert$
) is an arbitrary control volume around
$\boldsymbol{x}$
,
$V_{{\alpha }} \subset V$
(measure
$\vert V_{{\alpha }}\vert$
) is the volume of constituent
$\alpha$
so that
$V =\cup _{{\alpha }}V_{{\alpha }}$
. Furthermore, the constituents masses are
$M_{{\alpha }}=M_{{\alpha }}(V)$
, and the total mass in
$V$
is
$M=M(V)=\sum _{{\alpha }}M_{{\alpha }}(V)$
. The mixture density is the superposition of the partial mass densities:
\begin{align} \rho ({\boldsymbol{x}},t):&= \displaystyle \lim _{ \vert V \vert \rightarrow 0} \dfrac {M(V)}{\vert V \vert }=\displaystyle \sum _{{\alpha }}\tilde {\rho }_{\alpha }({\boldsymbol{x}},t). \end{align}
We introduce the mass fractions and volume fractions respectively as
\begin{align} Y_{\alpha }({\boldsymbol{x}},t):&= \displaystyle \lim _{ \vert V \vert \rightarrow 0} \dfrac {M_{{\alpha }}(V)}{M(V)}=\dfrac {\tilde {\rho }_{\alpha }}{\rho }, \end{align}
\begin{align} \phi _{\alpha }({\boldsymbol{x}},t):&= \displaystyle \lim _{ \vert V \vert \rightarrow 0} \dfrac {|V_{{\alpha }}|}{|V|}=\dfrac {\tilde {\rho }_{\alpha }}{\rho _{\alpha }}, \end{align}
which satisfy
\begin{align} \displaystyle \sum _{{\alpha }} Y_{\alpha }({\boldsymbol{x}},t)&= 1, \end{align}
\begin{align} \displaystyle \sum _{{\alpha }} \phi _{{\alpha }}({\boldsymbol{x}},t)&= 1. \end{align}
Next, we introduce the material time derivative
$\grave {\unicode {x03C8}}_{\alpha }$
of the differentiable constituent function
$\unicode {x03C8}_{\alpha }$
:
\begin{align} \grave {\unicode {x03C8}}_{\alpha }=\partial _t\unicode {x03C8}_{{\alpha }}(\boldsymbol{X}_{{\alpha }},t) \vert _{\boldsymbol{X}_{\alpha }}. \end{align}
Here the notation
$\vert _{\boldsymbol{X}_{\alpha }}$
is used to emphasise that
$\boldsymbol{X}_{\alpha }$
is held fixed. The constituent velocity is a constituent material derivative of the deformation map:
\begin{align} \boldsymbol{v}_{{\alpha }}(\boldsymbol{x},t)=\partial _t\boldsymbol{\chi }_{{\alpha }}(\boldsymbol{X}_{{\alpha }},t) \vert _{\boldsymbol{X}_{\alpha }} = \grave {\boldsymbol{\chi }}_{\alpha }. \end{align}
The mass-averaged mixture velocity
$\boldsymbol{v}$
(also called barycentric velocity) is defined via the following identification:
\begin{align} \rho ({\boldsymbol{x}},t) \boldsymbol{v}({\boldsymbol{x}},t) = \displaystyle \sum _{{\alpha }} \tilde {\rho }_{\alpha }({\boldsymbol{x}},t) \boldsymbol{v}_{\alpha }({\boldsymbol{x}},t). \end{align}
Furthermore, the (constituent) peculiar velocity is defined as
\begin{align} \boldsymbol{w}_{{\alpha }}({\boldsymbol{x}},t):=\boldsymbol{v}_{{\alpha }}({\boldsymbol{x}},t)-\boldsymbol{v}({\boldsymbol{x}},t) \end{align}
and specifies the constituent velocity relative to the mixture’s general motion. Introducing the (scaled) peculiar velocity
\begin{align} \boldsymbol{J}_{\alpha }({\boldsymbol{x}},t) :&= \tilde {\rho }_{\alpha }({\boldsymbol{x}},t) \boldsymbol{w}_{\alpha }({\boldsymbol{x}},t), \end{align}
we have the identity
\begin{align} \displaystyle \sum _{{\alpha }} \boldsymbol{J}_{\alpha }({\boldsymbol{x}},t) &= 0. \end{align}
Finally, we introduce the material derivative of the mixture as
\begin{align} \dot {\unicode {x03C8}}({\boldsymbol{x}},t) &= \partial _t \unicode {x03C8}({\boldsymbol{x}},t) + \boldsymbol{v}({\boldsymbol{x}},t)\boldsymbol{\cdot }\boldsymbol{\nabla }\unicode {x03C8}({\boldsymbol{x}},t). \end{align}
2.2. Constituent balance laws
The motion of each constituent in the continuum theory of mixtures is expressed by an individualised set of balance laws, as outlined in the second general principle. These laws incorporate interaction terms that depict how the different constituents interact. The motion of each constituent
${\alpha } = 1, \ldots , N$
is governed by the following set of local balance laws for all
$\boldsymbol{x}\in \varOmega$
and
$t \gt 0$
:
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}_{\alpha }) &= \gamma _{\alpha }, \end{align}
\begin{align} \partial _t (\tilde {\rho }_{\alpha }\boldsymbol{v}_{\alpha }) + \textrm {div} \left ( \tilde {\rho }_{\alpha }\boldsymbol{v}_{\alpha }\otimes \boldsymbol{v}_{\alpha } \right ) - \textrm {div} \boldsymbol{T}_{\alpha } - \tilde {\rho }_{\alpha } \boldsymbol{b}_{\alpha } &= \boldsymbol{\pi }_{\alpha }, \end{align}
\begin{align} \boldsymbol{T}_{\alpha }-\boldsymbol{T}_{\alpha }^T &=\boldsymbol{N}_{\alpha }. \end{align}
Here
$\varOmega$
is the spatial domain. Equation (2.13a
) describes the local constituent mass balance law, (2.13b
) the local constituent linear momentum balance law and (2.13c
) the local constituent angular momentum balance. The interaction terms are
$\gamma _{{\alpha }}, {\boldsymbol{\pi }}_{\alpha }$
and
$\boldsymbol{N}_{\alpha }$
. These denote respectively the mass supply of constituent
$\alpha$
due to chemical reactions with the other constituents, the momentum exchange rate of constituent
$\alpha$
and the intrinsic moment of momentum of constituent
$\alpha$
. Furthermore,
$\boldsymbol{T}_{\alpha }$
is the Cauchy stress tensor of constituent
$\alpha$
,
$\boldsymbol{b}_{\alpha }$
the constituent external body force. In this paper we assume equal body forces (
$\boldsymbol{b}_{\alpha }= \boldsymbol{b}$
for
${\alpha } = 1, \ldots , N$
) and restrict to body forces of gravitational type
$\boldsymbol{b} = -b \boldsymbol{\jmath } = -b \boldsymbol{\nabla }y$
, with
$y$
the vertical coordinate,
$\boldsymbol{\jmath }$
the vertical unit vector and
$b$
a constant.
In addition, one can infer that the constituent mass balance law can be written as
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) + \textrm {div} \boldsymbol{J}_{\alpha } &= \gamma _{\alpha }. \end{align}
We denote the kinetic and gravitational energies of the constituents respectively as
\begin{align} \mathscr{K}_{\alpha } &=\tilde {\rho }_{\alpha } \|\boldsymbol{v}_{\alpha }\|^2/2, \end{align}
\begin{align} \mathscr{G}_{\alpha } &=\tilde {\rho }_{\alpha } b y, \end{align}
where
$\|\boldsymbol{v}_{\alpha }\|=(\boldsymbol{v}_{\alpha } \boldsymbol{\cdot }\boldsymbol{v}_{\alpha })^{1/2}$
is the Euclidean norm of the velocity
$\boldsymbol{v}_{\alpha }$
.
Remark 2.1 (Volume-averaged mixture velocity). Alternative to the mass-averaged mixture velocity (2.8), another commonly adopted mixture velocity is the volume-averaged velocity
$\boldsymbol{u}$
defined as
${\boldsymbol{u}} = \sum _{\alpha } \phi _{\alpha }\boldsymbol{v}_{\alpha }$
. In the case of incompressible flow (
$\rho _{\alpha }(\boldsymbol{x},t)= \rho _{\alpha }$
) without mass transfer (
$\gamma _{\alpha }=0$
) one may deduce from (2.5b
) and (2.13a
) that it is divergence free, i.e.
$\textrm {div} {\boldsymbol{u}} = 0$
.
2.3. Mixture balance laws
Next, we proceed to the continuum balance laws of the mixtures. In agreement with the first general principle, the motion of the mixture results from that of the individual constituents. Summing the individual balance laws (2.13) across all constituents yields
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v} \otimes \boldsymbol{v} \right ) - \textrm {div} \boldsymbol{T} - \rho \boldsymbol{b} &=0, \end{align}
\begin{align} \boldsymbol{T}-\boldsymbol{T}^T &=0, \\[9pt] \nonumber \end{align}
where
\begin{align} \boldsymbol{T} :&= \sum _{\alpha } \boldsymbol{T}_{\alpha }-\tilde {\rho }_{\alpha }\boldsymbol{w}_{\alpha }\otimes \boldsymbol{w}_{\alpha }, \end{align}
\begin{align} \boldsymbol{b} :&=\frac {1}{\rho }\sum _{\alpha } \tilde {\rho }_{\alpha }\boldsymbol{b}_{\alpha }. \end{align}
In agreement with the third general principle, we have postulated the following balance conditions:
\begin{align} \displaystyle \sum _{\alpha } \gamma _{\alpha } &= 0, \end{align}
\begin{align} \displaystyle \sum _{\alpha } \boldsymbol{\pi }_{\alpha } &= 0, \end{align}
\begin{align} \displaystyle \sum _{\alpha } \boldsymbol{N}_{\alpha } &= 0. \end{align}
The first general principle of mixture theory additionally implies that the kinetic and gravitational energy of the mixture are the superposition of the constituent energies:
\begin{align} \mathscr{K} &= \displaystyle \sum _{{\alpha }} \mathscr{K}_{\alpha }, \end{align}
\begin{align} \mathscr{G} &= \displaystyle \sum _{{\alpha }} \mathscr{G}_{\alpha }. \end{align}
Remark 2.2 (Kinetic energy). The kinetic energy of the mixture can be decomposed as
\begin{align} \mathscr{K} &= \kern4pt \overline {\kern-4pt \mathscr{K}} + \displaystyle \sum _{{\alpha }} \frac {1}{2} \tilde {\rho }_{\alpha } \|\boldsymbol{w}_{\alpha }\|^2, \end{align}
\begin{align} \kern5pt \overline {\kern-5pt \mathscr{K}} &= \frac {1}{2} \rho \|\boldsymbol{v}\|^2, \end{align}
where
$\kern5pt \overline {\kern-5pt \mathscr{K}}$
is a kinetic energy quantity expressed in pure mixture variables, and where the second member is comprised of kinetic energies based on peculiar velocities.
3. Constitutive modelling
This section is devoted to the construction of constitutive models, employing a structured approach akin to ten Eikelder et al. (Reference ten Eikelder, Van Der Zee and Schillinger2024), ten Eikelder (Reference ten Eikelder2025). First, § 3.1 provides assumptions and details on the Coleman–Noll modelling procedure (Coleman & Noll Reference Coleman and Noll1974). Following this, in § 3.2 we determine the constitutive modelling restriction following from § 3.1. Finally, in § 3.3 we identify specific constitutive models that align with these restrictions.
3.1. Assumptions and modelling choices
Instead of working with the full set of balance laws (2.13a )–(2.13c ), we restrict to the reduced set
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) + \textrm {div} \boldsymbol{H}_{\alpha } &= \zeta _{\alpha }, \end{align}
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) - \textrm {div} \boldsymbol{T} - \rho \boldsymbol{b} &=0, \end{align}
\begin{align} \boldsymbol{T}-\boldsymbol{T}^T &=0, \end{align}
and where
${\alpha }=1,\ldots ,N$
in (3.1a
). We have decomposed the mass-transfer terms into conservative and possible non-conservative parts via
$\gamma _{\alpha } = \zeta _{\alpha } - \textrm {div} \boldsymbol{j}_{\alpha }$
with
$\boldsymbol{H}_{\alpha } := \boldsymbol{J}_{\alpha } + \boldsymbol{j}_{\alpha }$
. The state variables in the system (3.1) are
$\tilde {\rho }_{\alpha }$
(
${\alpha } = 1,\ldots ,N$
) and
$\boldsymbol{v}$
. We seek constitutive models for
$\boldsymbol{T}$
,
$\boldsymbol{H}_{\alpha }$
and
$\zeta _{\alpha }$
(
${\alpha } = 1,\ldots ,N$
). Hereby we discard the definition (2.10), but enforce the balance condition (2.11).
Remark 3.1 (Decomposition mass transfer). Rather than modelling
$\boldsymbol{H}_{\alpha }$
and
$\zeta _{\alpha }$
, one may work directly with the original terms
$\boldsymbol{J}_{\alpha }$
and
$\gamma _{\alpha }$
, as both approaches yield equivalent closure models. The decomposition
$\gamma _{\alpha } = \zeta _{\alpha } - \operatorname {div}\boldsymbol{j}_{\alpha }$
is introduced solely for interpretative purposes. For example, assuming identical constituent velocities,
$\boldsymbol{v}_{\alpha } = \boldsymbol{v}$
, does not conflict with the constitutive model for
$\boldsymbol{H}_{\alpha }$
.
We postulate the energy-dissipation law
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t} \mathscr{E} = \mathscr{W} - \mathscr{D} \end{align}
as a design principle for constructing closure models for
$\boldsymbol{T}, \boldsymbol{H}_{\alpha }$
and
$\zeta _{\alpha }$
(
${\alpha }=1,\ldots ,N$
). The total energy is composed of the Helmholtz free energy, the kinetic energy and the gravitational energy:
\begin{align} \mathscr{E} = \displaystyle \int _{\mathcal{R}(t)}(\varPsi + \kern5pt \overline {\kern-5pt \mathscr{K}} + \mathscr{G}\kern2pt)\,\textrm {d}v. \end{align}
In this formulation,
$\mathcal{R}(t) \subset \varOmega$
is an arbitrary, time-dependent control volume with volume element
$\textrm {d}v$
and unit outward normal
$\boldsymbol{\nu }$
, transported by the velocity field
$\boldsymbol{v}$
. Moreover,
$\mathscr{W}$
denotes the work rate applied on the boundary
$\partial \mathcal{R}(t)$
, while
$\mathscr{D}$
represents the internal dissipation, which is required to satisfy
$\mathscr{D}\geqslant 0$
.
Remark 3.2 (Energy-dissipation statement). In the energy-dissipation law (
3.2
), the kinetic energy term is taken as
$ \kern5pt \overline {\kern-5pt \mathscr{K}}$
, rather than the full kinetic energy of the mixture
$\mathscr{K}$
(see
Remark 2.2
). Consequently, (
3.2
) can be seen as a simplified formulation of the second law of thermodynamics for mixtures (see, e.g. ten Eikelder et al. Reference ten Eikelder, van der Zee, Akkerman and Schillinger2023, Reference ten Eikelder, Van Der Zee and Schillinger2024). Moreover, the use of control volumes is not essential; an equivalent modelling constraint, as discussed in §
3.2
, may be derived through analogous steps at the local PDE level.
In this paper we postulate the free energy to belong to the constitutive class
\begin{align} \varPsi = \hat {\varPsi }\left (\left \{\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N},\left \{\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N}\right )\!, \end{align}
where we note that both
$\left \{\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N}$
and
$\left \{\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N}$
consist of independent variables. In addition, we introduce the constituent chemical potentials
\begin{align} \hat {\mu }_{\alpha } &= \dfrac { \partial \hat {\varPsi }}{\partial \tilde {\rho }_{\alpha }} - \textrm {div}\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}. \end{align}
Remark 3.3 (Relation to compressible two-phase flow models). The modelling choices presented in this section distinguish our approach from conventional compressible two-phase flow models in several important ways. First, our formulation is based on the reduced set of balance laws (
3.1
). In contrast, compressible two-phase flow models typically include additional evolution equations for the volume fractions, which do not naturally arise from the continuum mixture theory described in §
2
. Consequently, while the state variables in compressible two-phase flow models consist of specific densities
$\{\rho _{\alpha }\}$
and volume fractions
$\{\phi _{\alpha }\}$
(subject to the saturation constraint (
2.5b
)), our approach employs partial densities
$\{\tilde {\rho }_{\alpha }\}$
instead. Finally, the constitutive class (
3.4
) depends on all partial densities, thereby enabling the modelling of interaction (van der Waals) forces – an aspect that is absent in standard compressible two-phase flow models.
Remark 3.4 (Constitutive class free energy). In many conventional phase-field models, including those of the Cahn–Hilliard type, the free energy is expressed in terms of order parameters such as volume fractions or mass fractions (often referred to as concentrations). In contrast, the constitutive class ( 3.4 ) is formulated in terms of partial densities. We explore the relationship between these approaches in the binary case in § 6 and provide the corresponding derivation in Appendix C .
3.2. Modelling restriction
We proceed with the evaluation of the evolution of the energy (3.3). By applying Reynolds transport theorem to the free energy
$\hat {\varPsi }$
we have
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {d}v = \displaystyle \int _{\mathcal{R}(t)} \partial _t \hat {\varPsi } \,\textrm {d}v + \displaystyle \int _{\partial \mathcal{R}(t)} \hat {\varPsi } \boldsymbol{v} \boldsymbol{\cdot }{\boldsymbol{\nu }} \,\textrm {d}a. \end{align}
We apply the divergence theorem and expand the derivatives:
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {d}v = \displaystyle \int _{\mathcal{R}(t)} &\, \hat {\varPsi } \,\textrm {div} \boldsymbol{v} + \sum _{\alpha }\dfrac {\partial \hat {\varPsi }}{\partial \tilde {\rho }_{\alpha }} \dot {\tilde {\rho }}_{\alpha }+ \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\boldsymbol{\cdot }\left (\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right )\dot \,\textrm {d}v. \end{align}
Substituting the identity for the material derivative
\begin{align} (\boldsymbol{\nabla }\omega )\dot = \boldsymbol{\nabla }(\dot {\omega }) - (\boldsymbol{\nabla }\omega )^T\boldsymbol{\nabla }\!\boldsymbol{v} \end{align}
with
$\omega = \tilde {\rho }_{\alpha }$
into (3.7) yields
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {d}v = \displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {div} \boldsymbol{v} + \sum _{\alpha }\dfrac {\partial \hat {\varPsi }}{\partial \tilde {\rho }_{\alpha }} \dot {\tilde {\rho }}_{\alpha }+ \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\boldsymbol{\cdot }\boldsymbol{\nabla }\bigl(\dot {\tilde {\rho }}_{\alpha }\bigr) - \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\boldsymbol{\cdot }(\boldsymbol{\nabla }\tilde {\rho }_{\alpha })^T\boldsymbol{\nabla }\!\boldsymbol{v} \,\textrm {d}v. \end{align}
By subsequently integrating by parts, we arrive at
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi }\,\textrm {div} \boldsymbol{v} + \sum _{\alpha } \hat {\mu }_{\alpha } \dot {\tilde {\rho }}_{\alpha } - \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}: \boldsymbol{\nabla }\!\boldsymbol{v} \,\textrm {d}v \nonumber \\ &\,+ \displaystyle \int _{\partial \mathcal{R}(t)}\sum _{\alpha } \dot {\tilde {\rho }}_{\alpha } \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a. \end{align}
Substituting the constituent mass balance laws (2.14) and applying integration by parts leads to
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi } \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \hat {\varPsi }\,\textrm {div} \boldsymbol{v} -\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } \textrm {div}\boldsymbol{v} +\sum _{\alpha } \boldsymbol{\nabla }\hat {\mu }_{\alpha } \boldsymbol{\cdot }\boldsymbol{H}_{\alpha } \nonumber \\ &\quad - \sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}: \boldsymbol{\nabla }\!\boldsymbol{v} +\sum _{\alpha } \hat {\mu }_{\alpha } \zeta _{\alpha }\,\textrm {d}v\nonumber \\ &\quad + \displaystyle \int _{\partial \mathcal{R}(t)}\left (\sum _{\alpha } \dot {\tilde {\rho }}_{\alpha } \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}-\hat {\mu }_{\alpha } \boldsymbol{H}_{\alpha }\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a. \end{align}
Subsequently, the evolution equations for the kinetic and gravitational energies follow straightforward from (3.1a ) and (3.1b ) (see ten Eikelder et al. Reference ten Eikelder, van der Zee, Akkerman and Schillinger2023 for details):
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \mathscr{K} \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} - \boldsymbol{\nabla }\!\boldsymbol{v} : \boldsymbol{T}+\rho \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{g} \,\textrm {d}v+ \displaystyle \int _{\partial \mathcal{R}(t)} \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{T} {\boldsymbol{\nu }} \,\textrm {d}a, \end{align}
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \mathscr{G} \,\textrm {d}v &=- \displaystyle \int _{\mathcal{R}(t)} \rho \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{g}\,\textrm {d}v. \end{align}
Summing (3.11) and (3.12) yields
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t} \mathscr{E} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T\boldsymbol{T}-\sum _{\alpha } \left (\hat {\mu }_{\alpha }\boldsymbol{H}_{\alpha } -\dot {\tilde {\rho }}_{\alpha } \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right )\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a \nonumber \\ &\quad - \displaystyle \int _{\mathcal{R}(t)} \left (\boldsymbol{T} +\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\right )\boldsymbol{I}\right ):\boldsymbol{\nabla }\!\boldsymbol{v}\nonumber \\ &\quad + \sum _{\alpha }\left (-\boldsymbol{\nabla }\hat {\mu }_{\alpha }\boldsymbol{\cdot }\boldsymbol{H}_{\alpha } - \hat {\mu }_{\alpha }\zeta _{\alpha } \right )\,\textrm {d}v, \end{align}
where we identify the rate of work and the dissipation as
\begin{align} \mathscr{W} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T\boldsymbol{T}-\sum _{\alpha } \left (\hat {\mu }_{\alpha }\boldsymbol{H}_{\alpha } -\dot {\tilde {\rho }}_{\alpha } \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right )\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a, \end{align}
\begin{align} \mathscr{D} &= \displaystyle \int _{\mathcal{R}(t)} \left (\boldsymbol{T} +\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\right )\boldsymbol{I}\right ):\boldsymbol{\nabla }\!\boldsymbol{v}\nonumber \\ &\quad + \sum _{\alpha }\left (-\boldsymbol{\nabla }\hat {\mu }_{\alpha }\boldsymbol{\cdot }\boldsymbol{H}_{\alpha } - \hat {\mu }_{\alpha }\zeta _{\alpha } \right )\,\textrm {d}v. \end{align}
Since the control volume
$\mathcal{R}=\mathcal{R}(t)$
is arbitrary, the energy-dissipation law is satisfied provided that the following local inequality holds:
\begin{align} \left (\boldsymbol{T} +\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\right )\boldsymbol{I}\right ):\boldsymbol{\nabla }\!\boldsymbol{v} - \sum _{\alpha } \boldsymbol{H}_{\alpha }\boldsymbol{\cdot }\boldsymbol{\nabla }\hat {\mu }_{\alpha } - \sum _{\alpha } \zeta _{\alpha }\hat {\mu }_{\alpha } &\geqslant 0. \end{align}
3.3. Selection of constitutive models
By means of the Colemann–Noll procedure, the inequality (3.15) permits us to restrict to mixture stress tensors
$\boldsymbol{T}$
, constituent diffusive fluxes
$\boldsymbol{H}_{\alpha }$
and constituent mass fluxes
$\zeta _{\alpha }$
that belong to the constitutive classes:
\begin{align} \boldsymbol{T} &= \hat {\boldsymbol{T}}\left (\boldsymbol{\nabla }\!\boldsymbol{v}, \left \{\tilde {\rho }_{\beta }\right \}\!, \left \{\boldsymbol{\nabla }\tilde {\rho }_{\beta }\right \}\!, \left \{\hat {\mu }_{\beta }\right \}\!,\left \{\boldsymbol{\nabla }\hat {\mu }_{\beta }\right \}\right )\!, \end{align}
\begin{align} \boldsymbol{H}_{\alpha } &= \hat {\boldsymbol{H}}_{\alpha }\left (\left \{\tilde {\rho }_{\beta }\right \}\!, \left \{\boldsymbol{\nabla }\tilde {\rho }_{\beta }\right \}\!, \left \{\hat {\mu }_{\beta }\right \}\!,\left \{\boldsymbol{\nabla }\hat {\mu }_{\beta }\right \}\right )\!, \end{align}
\begin{align} \zeta _{\alpha } &= \hat {\zeta }_{\alpha }\left (\left \{\tilde {\rho }_{\beta }\right \}\!, \left \{\hat {\mu }_{\beta }\right \}\right )\!. \\[9pt] \nonumber \end{align}
We seek for constitutive models (3.16) that render each of the three terms in (3.15) non-negative:
\begin{align} \left (\hat {\boldsymbol{T}} +\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\right )\boldsymbol{I}\right ):\boldsymbol{\nabla }\!\boldsymbol{v} &\geqslant 0, \end{align}
\begin{align} -\sum _{\alpha } \hat {\boldsymbol{H}}_{\alpha }\boldsymbol{\cdot }\boldsymbol{\nabla }\hat {\mu }_{\alpha } &\geqslant 0, \end{align}
\begin{align} -\sum _{\alpha } \hat {\zeta }_{\alpha }\hat {\mu }_{\alpha } &\geqslant 0. \end{align}
Constitutive choices: we make the following constitutive choices:
\begin{align} \hat {\boldsymbol{T}} &= -\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} - \left (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\right )\boldsymbol{I} + \nu \left (2\boldsymbol{D}+\lambda (\textrm {div}\,\boldsymbol{v})\boldsymbol{I}\right )\!, \end{align}
\begin{align} \hat {\boldsymbol{H}}_{\alpha } &= -\sum _{\beta } \boldsymbol{M}_{{\alpha }{\beta }} \boldsymbol{\nabla }\hat {\mu }_{\beta }, \end{align}
\begin{align} \hat {\zeta }_{\alpha } &= -\sum _{\beta } m_{{\alpha }{\beta }}\hat {\mu }_{\beta }. \end{align}
Here
$\nu \geqslant 0$
is the dynamic viscosity,
$\lambda \geqslant -2/d$
(with
$d$
denoting the number of dimensions),
$\boldsymbol{M}$
is a symmetric positive definite tensor, with the same dependencies as (3.16b
), satisfying
\begin{align} \sum _{\alpha } \boldsymbol{M}_{{\alpha }{\beta }} = \sum _{\beta } \boldsymbol{M}_{{\alpha }{\beta }} = 0,\quad \boldsymbol{M}_{{\alpha }{\beta }}|_{Y_{\gamma }=1}=0\quad \text{for all } {\gamma }=1,\ldots ,N, \end{align}
and
$m$
is a symmetric positive definite scalar mobility, with the same dependencies as (3.16c
), satisfying
\begin{align} \sum _{\alpha } m_{{\alpha }{\beta }} = \sum _{\beta } m_{{\alpha }{\beta }} = 0,\quad m_{{\alpha }{\beta }}|_{Y_{\gamma }=1}=0\quad \text{for all } {\gamma }=1,\ldots ,N. \end{align}
Moreover, compatibility with the angular momentum condition (3.1c ) requires that
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} = \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \otimes \boldsymbol{\nabla }\tilde {\rho }_{\alpha }. \end{align}
Lemma 3.5 (Compatibility of constitutive choices). The choices (3.18a )–(3.18c ) are compatible with the restrictions (3.17a )–(3.17c ), as well as with the balance of relative gross motion (2.11) and the balance of mass supply (2.18a ).
Proof. For the stress tensor, a direct calculation yields
\begin{align} \Biggl (\hat {\boldsymbol{T}} &+\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} + \Bigl (\sum _{\alpha } \hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } - \hat {\varPsi }\Bigr )\boldsymbol{I}\Biggr ):\boldsymbol{\nabla }\!\boldsymbol{v} \nonumber \\ &= 2 \nu \Biggl (\boldsymbol{D} - \frac {1}{d} (\textrm {div}\,\boldsymbol{v})\boldsymbol{I}\Biggr ):\Biggl (\boldsymbol{D} - \frac {1}{d} (\textrm {div}\,\boldsymbol{v})\boldsymbol{I}\Biggr ) + \nu \Biggl (\lambda + \frac {2}{d}\Biggr ) (\textrm {div}\,\boldsymbol{v})^2 \ge 0, \end{align}
which confirms compatibility with (3.17a
). For the diffusive flux, the condition
$\sum _{\alpha } \boldsymbol{M}_{{\alpha }{\beta }}=0$
guarantees compatibility with the balance (2.11), while the positive semi-definiteness of
$\boldsymbol{M}_{{\alpha }{\beta }}$
ensures the thermodynamical restriction (3.17b
). Similarly, for the mass-transfer term, the properties
$\sum _{\alpha } m_{{\alpha }{\beta }}=0$
and the positive definiteness of
$m_{{\alpha }{\beta }}$
ensure compatibility with the mass supply balance (2.18a
) and the restriction (3.17c
).
Remark 3.6 (Mass transfer in compressible two-phase flow models). We consider the special case:
\begin{align} m_{{\alpha }{\beta }} = \begin{cases} -\hat {m}_{{\alpha }{\beta }} & \text{if } {\alpha } \neq {\beta }, \\ \sum _{{\gamma }\neq {\alpha }} \hat {m}_{{\alpha }{\gamma }} & \text{if } {\alpha } = {\beta }. \end{cases} \end{align}
Under this assumption, we obtain
\begin{align} \hat {\zeta }_{\alpha } &= - \sum _{{\beta }\neq {\alpha }} m_{{\alpha }{\beta }} \boldsymbol{\nabla }\hat {\mu }_{\beta } - m_{{\alpha }{\alpha }} \boldsymbol{\nabla }\hat {\mu }_{\alpha } \nonumber \\ &= \sum _{{\beta }\neq {\alpha }} \hat {m}_{{\alpha }{\beta }} \boldsymbol{\nabla }\hat {\mu }_{\beta } - \sum _{{\gamma }\neq {\alpha }} \hat {m}_{{\alpha }{\gamma }} \boldsymbol{\nabla }\hat {\mu }_{\alpha } \nonumber \\ &= -\sum _{{\beta }} \hat {m}_{{\alpha }{\beta }} \boldsymbol{\nabla }(\hat {\mu }_{\alpha } - \hat {\mu }_{\beta }). \end{align}
This formulation is commonly used for binary flows in compressible two-phase models; see, e.g. Coquel et al. (Reference Coquel, Gallouët, Helluy, Hérard, Hurisse and Seguin2013), Saurel & Pantano (Reference Saurel and Pantano2018), Pelanti (Reference Pelanti2022).
This concludes the fundamental exploration of constitutive models compatible with the imposed energy-dissipative modelling restriction. Substitution provides the model
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p + \textrm {div} \left (\sum _{\beta } \boldsymbol{\nabla }\tilde {\rho }_{\beta } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\beta }} \right )& \nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) - \sum _{{\beta }} \textrm {div} \left (\boldsymbol{M}_{{\alpha }{\beta }}\boldsymbol{\nabla }\hat {\mu }_{\beta }\right ) + \sum _{{\beta }} m_{{\alpha }{\beta }} \hat {\mu }_{\beta } &=0, \end{align}
\begin{align} \hat {\mu }_{\alpha } - \dfrac {\partial \varPsi }{\partial \tilde {\rho }_{\alpha }} + \textrm {div} \left ( \dfrac {\partial \varPsi }{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right )&=0, \end{align}
where the pressure
$p$
takes the form
\begin{align} p=\sum _{\beta } (\hat {\mu }_{\beta }\tilde {\rho }_{\beta }) -\hat {\varPsi }. \end{align}
Remark 3.7 (Relation to Cahn–Hilliard/Allen–Cahn/Korteweg models). The framework connects to Cahn–Hilliard/Allen–Cahn/Korteweg-type models by selecting
$\hat {\varPsi } = \hat {\varPsi }_0+(1/2)\sum _{{\beta }{\gamma }} \lambda _{{\beta }{\gamma }} \boldsymbol{\nabla }\tilde {\rho }_{\beta }\boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {\rho }_{\gamma }$
, where
$\hat {\varPsi }_0=\hat {\varPsi }_0(\left \{\tilde {\rho }_{\beta }\right \})$
and where
$\lambda _{{\beta }{\gamma }}$
are model parameters.
Remark 3.8 (Alternative pressure variable). A common alternative approach is to decompose the stress tensor into its volumetric (isotropic) and deviatoric (traceless) components, which allows one to define the volumetric pressure
$p^{\textrm {vol}}$
:
\begin{align} \hat {\boldsymbol{T}} &= -p^{\textrm {vol}} \boldsymbol{I} + \hat {\boldsymbol{T}}^{\textrm {dev}}, \end{align}
\begin{align} p^{\textrm {vol}} &= \sum _{\beta } \left (\hat {\mu }_{\beta }\tilde {\rho }_{\beta } - \frac {1}{d} \boldsymbol{\nabla }\tilde {\rho }_{\beta } \boldsymbol{\cdot }\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\beta }}\right ) -\hat {\varPsi } - \nu \left (\lambda + \frac {2}{d}\right )\textrm {div} \boldsymbol{v}, \end{align}
\begin{align} \hat {\boldsymbol{T}}^{\textrm {dev}} &=\sum _{\beta }\left ( \boldsymbol{\nabla }\tilde {\rho }_{\beta } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\beta }}- \frac {1}{d}\boldsymbol{\nabla }\tilde {\rho }_{\beta } \boldsymbol{\cdot }\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\beta }} \boldsymbol{I}\right ) + 2 \nu \left ( \boldsymbol{D} - \frac {1}{d} (\textrm {div} \boldsymbol{v}) \boldsymbol{I}\right )\!. \end{align}
We explicitly state the compatibility with the energy-dissipation condition.
Theorem 3.9 (Compatibility energy dissipation). The model (3.25) is compatible with the energy-dissipation condition (3.2).
Proof. This follows from Lemma3.5, in particular we have
\begin{align} \mathscr{D} &= \displaystyle \int _{\mathcal{R}(t)} 2 \nu \left ( \boldsymbol{D} - \frac {1}{d} (\textrm {div} \boldsymbol{v}) \boldsymbol{I}\right ):\left (\boldsymbol{D} - \frac {1}{d} (\textrm {div} \boldsymbol{v}) \boldsymbol{I}\right )+ \nu \left (\lambda + \frac {2}{d}\right )\left (\textrm {div} \boldsymbol{v}\right )^2\nonumber \\ &\,\quad \quad + \displaystyle \sum _{{\alpha },{\beta }} M_{{\alpha }{\beta }}\boldsymbol{\nabla }\hat {\mu }_{\alpha }\boldsymbol{\cdot }\boldsymbol{\nabla }\hat {\mu }_{\beta } +\sum _{{\alpha },{\beta }} m_{{\alpha }{\beta }}\hat {\mu }_{\alpha }\hat {\mu }_{\beta }\,\textrm {d}v \geqslant 0. \end{align}
4. Properties
We first discuss compact formulations in § 4.1. Next, we discuss the decomposition of the free energy in § 4.2. Finally, in § 4.3 we present the equilibrium conditions of the model.
4.1. Compact formulations
Simplification of the model (3.25) can be achieved by employing the following identity.
Lemma 4.1 (Identity Korteweg stresses). We have the following identity for the pressure and Korteweg stresses:
\begin{align} \boldsymbol{\nabla }\!p + \textrm {div} \left ( \sum _{\alpha }\boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right ) = \sum _{\alpha }\tilde {\rho }_{\alpha }\boldsymbol{\nabla }\hat {\mu }_{\alpha }. \end{align}
Proof. See Lemma A.1.
Remark 4.2. The identity in Lemma 4.1 represents the multicomponent extension of the well-known identity for the Korteweg stress.
Applying Lemma 4.1 we find the more compact form
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \sum _{{\beta }} \tilde {\rho }_{\beta }\boldsymbol{\nabla }\hat {\mu }_{\beta } - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) - \sum _{{\beta }} \textrm {div} \left (M_{{\alpha }{\beta }}\boldsymbol{\nabla }\hat {\mu }_{\beta }\right ) + \sum _{{\beta }} m_{{\alpha }{\beta }} \hat {\mu }_{\beta } &=0, \end{align}
for
${\alpha } = 1,\ldots , N$
. Furthermore, we note that the model may alternatively be written in a form that more closely links to existing phase-field models, i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \rho \sum _{{\beta }} Y_{\beta }\boldsymbol{\nabla }\hat {\mu }_{\beta } - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &=0, \end{align}
\begin{align} \partial _t (\rho Y_{\alpha }) + \textrm {div}(\rho Y_{\alpha } \boldsymbol{v}) - \sum _{{\beta }} \textrm {div} \left (M_{{\alpha }{\beta }}\boldsymbol{\nabla }\hat {\mu }_{\beta }\right ) + \sum _{{\beta }} m_{{\alpha }{\beta }} \hat {\mu }_{\beta } &=0, \end{align}
for
${\alpha } = 1,\ldots ,N-1$
. We note that the combination of the mixture mass balance (4.3b
) and the
$N-1$
constituent balance laws (4.3c
) are equivalent to the
$N$
balance laws (4.2b
) or
$N$
balance laws (4.3c
).
We observe that, in the single-fluid regime, the model simplifies to the compressible Navier–Stokes equations.
Proposition 4.3 (Reduction to compressible Navier–Stokes). If the chemical potentials
$\hat {\mu }_{\alpha }$
are well defined for
$\tilde {\rho }_{\alpha } = 0$
(
${\alpha } = 1,\ldots ,N$
), the multiconstituent system (4.3) reduces to the compressible Navier–Stokes equations in the single constituent regime (
$Y_{\beta }= 1$
), i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p - \textrm {div} \left (\nu \left (2\boldsymbol{D} + \lambda \textrm {div}\boldsymbol{v}\right ) \right ) -\rho \boldsymbol{b}&= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
with
$\rho _{\beta } = \tilde {\rho }_{\beta } = \rho$
and
$p=p(\rho )$
.
4.2. Decomposition of the free energy
We proceed with decomposing the free energy into its components; in agreement with the first metaphysical principle of continuum mixture theory we have
\begin{align} \hat {\varPsi }\left (\left \{\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N},\left \{\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N}\right ) = \sum _{\beta } \hat {\varPsi }_{\beta }\left (\left \{\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N},\left \{\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,\ldots ,N}\right )\!, \end{align}
where
$\hat {\varPsi }_{\beta }$
is the volume-measure constituent free energy. The split of the free energy, (4.5), highlights that the constituent free energies may depend on quantities of other constituents (
$\tilde {\rho }_{\alpha }$
and
$\boldsymbol{\nabla }\tilde {\rho }_{\alpha }$
for all
${\alpha } = 1,\ldots ,N$
), allowing for the incorporation of, for example, attraction/repelling forces.
Next, we associate constituent chemical potentials and constituent partial pressures with each
$\hat {\varPsi }_{\beta }$
in the sense:
\begin{align} \hat {\mu }_{{\beta }{\alpha }} :&= \dfrac { \partial \hat {\varPsi }_{\beta }}{\partial \tilde {\rho }_{\alpha }} - \textrm {div}\dfrac {\partial \hat {\varPsi }_{\beta }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}, \end{align}
\begin{align} p_{\beta } &= \left (\sum _{\alpha }\hat {\mu }_{{\beta }{\alpha }} \tilde {\rho }_{\alpha }\right ) - \hat {\varPsi }_{\beta }. \end{align}
The constituent partial pressures
$p_{\alpha }$
are an extension of the classical partial pressure that includes gradient contributions essential to modelling surface tension effects. These constituent quantities satisfy the properties
\begin{align} \hat {\mu }_{\alpha } &= \sum _{{\beta }} \hat {\mu }_{{\beta }{\alpha }}, \end{align}
\begin{align} p &= \sum _{\beta } p_{\beta }. \end{align}
The latter may be referred to as Dalton’s law. As such, the split (4.5) reveals that the system may be written as
\begin{align} \sum _{\beta } \left ( \partial _t (\tilde {\rho }_{\beta } \boldsymbol{v}) + \textrm {div} \left ( \tilde {\rho }_{\beta } \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p_{\beta } + \textrm {div} \left (\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }_{\beta }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right )-\tilde {\rho }_{\beta } \boldsymbol{g} \right ) & \nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right ) &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) - \sum _{{\beta }} \textrm {div} \left (M_{{\alpha }{\beta }}\boldsymbol{\nabla }\hat {\mu }_{\beta }\right ) + \sum _{{\beta }} m_{{\alpha }{\beta }} \hat {\mu }_{\beta } &=0. \end{align}
Remark 4.4 (Absence of interaction forces). When interaction forces are absent, i.e.
$\hat {\varPsi }_{\beta } = \hat {\varPsi }_{\beta }(\tilde {\rho }_{\beta }, \boldsymbol{\nabla }\tilde {\rho }_{\beta })$
, we have
\begin{align} \hat {\mu }_{{\beta }{\alpha }}&=\hat {\mu }_{\alpha } \delta _{{\alpha }{\beta }}, \end{align}
\begin{align} p_{{\beta }}&= \hat {\mu }_{\beta }\tilde {\rho }_{\beta } - \hat {\varPsi }_{\beta }, \end{align}
so that
\begin{align} \sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }_{\beta }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} = \boldsymbol{\nabla }\tilde {\rho }_{\beta } \otimes \dfrac {\partial \hat {\varPsi }_{\beta }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\beta }}. \end{align}
Remark 4.5 (Pressure/temperature). We emphasise that the model is isothermal: a single (constant) temperature is assumed and the thermodynamic state of the mixture is determined by one Helmholtz free energy
$\hat {\varPsi }(\{\tilde {\rho }_\beta \}\!,\{\boldsymbol{\nabla }\tilde {\rho }_\beta \})$
. Consequently, despite the absence of explicit species energy equations, the model satisfies a standard mixture energy-dissipation balance derived from the second law. A single mixture pressure
$p$
is defined from
$\hat {\varPsi }$
; when
$\hat {\varPsi }=\sum _{\beta }\hat {\varPsi }_{\beta }$
, constituent pressures
$p_{\beta }$
can also be introduced for interpretation, with
$p=\sum _{\beta }p_{\beta }$
in general and
$p_{\beta }$
not necessarily equal.
4.3. Equilibrium conditions
We characterise the equilibrium conditions of the model (4.8) by
\begin{align} \mathscr{D} = 0. \end{align}
The identity (3.28) provides
\begin{align} 2 \nu ^{\textit{eq}} \left ( \boldsymbol{D}^{\textit{eq}} - \frac {1}{d} (\textrm {div} \boldsymbol{v}^{\textit{eq}}) \boldsymbol{I}\right ):\left (\boldsymbol{D}^{\textit{eq}} - \frac {1}{d} (\textrm {div} \boldsymbol{v}^{\textit{eq}}) \boldsymbol{I}\right )&=0, \end{align}
\begin{align} \nu ^{\textit{eq}} \left (\lambda + \frac {2}{d}\right )\left (\textrm {div} \boldsymbol{v}^{\textit{eq}}\right )^2 &=0, \end{align}
\begin{align} \displaystyle \sum _{{\alpha },{\beta }} \boldsymbol{M}_{{\alpha }{\beta }}^{\textit{eq}}\boldsymbol{\nabla }\hat {\mu }_{\alpha }^{\textit{eq}}\boldsymbol{\cdot }\boldsymbol{\nabla }\hat {\mu }_{\beta }^{\textit{eq}} &=0, \end{align}
\begin{align} \sum _{{\alpha },{\beta }} m_{{\alpha }{\beta }}^{\textit{eq}}\hat {\mu }_{\alpha }^{\textit{eq}}\hat {\mu }_{\beta }^{\textit{eq}} &=0, \end{align}
where the superscripts ‘
$\textrm {eq}$
’ denote the equilibrium of the quantity. Identities (4.12a
) and (4.12b
) dictate that
$\boldsymbol{v}^{\textit{eq}}$
are rigid motions. Since
$\boldsymbol{M}_{{\alpha }{\beta }}$
is symmetric positive definite, we decompose it as
$\boldsymbol{M}_{{\alpha }{\beta },ij}=\sum _q \boldsymbol{B}_{q,{\alpha },i}\boldsymbol{B}_{q,{\beta },j}$
for some
$\boldsymbol{B}_{q,{\alpha },i}$
with the same dependencies as
$\boldsymbol{M}_{{\alpha }{\beta },ij}$
. Hence, (4.12c
) implies that
$\sum _{\alpha } \boldsymbol{B}_{\alpha }^{\textit{eq}} \boldsymbol{\nabla }\hat {\mu }_{\alpha }^{\textit{eq}} =0$
. In a similar fashion we deduce that
$\sum _{\alpha } B_{{\alpha }, q}^{\textit{eq}} \hat {\mu }_{\alpha }^{\textit{eq}} =0$
for all
$q=1,\ldots ,N$
, where
$m_{{\alpha }{\beta }} = \sum _q B_{{\alpha }, q}B_{{\beta }, q}$
. These conditions imply that
$\hat {\boldsymbol{H}}^{\textit{eq}}_{\alpha } = 0, \hat {\zeta }^{\textit{eq}}_{\alpha }=0$
for
${\alpha }=1,\ldots ,N$
. Finally, restricting to zero velocities,
$\boldsymbol{v}^{\textit{eq}}=\boldsymbol{0}$
, the mass balance laws imply constant in time partial densities,
$\tilde {\rho }_{\alpha }^{\textit{eq}}(\boldsymbol{x},t) = \tilde {\rho }_{\alpha }^{\textit{eq}}(\boldsymbol{x})$
, while the momentum balance dictates that
$\sum _{\alpha } Y_{\alpha }^{\textit{eq}}\boldsymbol{\nabla }\hat {\mu }_{\alpha }^{\textit{eq}} = \boldsymbol{g}$
. Summarising, we have the equilibrium conditions
\begin{align} \boldsymbol{v}^{\textit{eq}}&=\boldsymbol{0}, \end{align}
\begin{align} \tilde {\rho }_{\alpha }^{\textit{eq}}&=\textrm {const }\, (\text{in time}), \end{align}
\begin{align} \sum _{\alpha } b_{\alpha }^{\textit{eq}} \hat {\mu }_{\alpha }^{\textit{eq}} &=0, \end{align}
\begin{align} \sum _{\alpha } \boldsymbol{B}_{\alpha }^{\textit{eq}} \boldsymbol{\nabla }\hat {\mu }_{\alpha }^{\textit{eq}} &=0, \end{align}
\begin{align} \sum _{\alpha } Y_{\alpha }^{\textit{eq}} \boldsymbol{\nabla }\hat {\mu }_{\alpha }^{\textit{eq}} &=\boldsymbol{g}. \end{align}
Remark 4.6 (No equilibrium conditions pressure). If the free energy
$\varPsi$
does not depend on
$\left \{\boldsymbol{\nabla }\tilde {\rho }_{\alpha }\right \}$
, the condition (4.13e
) converts into
$\boldsymbol{\nabla }\!p = \rho \boldsymbol{g}$
(see
Lemma 4.1
). Thus, there is no equilibrium condition on the partial pressures
$p_{\alpha }, {\alpha }=1,\ldots ,N$
. This is in contrast to some compressible two-phase models in which pressure equilibrium is assumed; see, e.g. Kapila et al. (Reference Kapila, Menikoff, Bdzil, Son and Stewart2001), Murrone & Guillard (Reference Murrone and Guillard2005), Daude et al. (Reference Daude, Galon, Potapov, Beccantini and Mianné2023) (and §
7.1
).
5. Analysis of the first-order system
In this section we examine the first-order system associated with (3.25). Beyond presenting the formulation, our aim is to clarify its structural properties and to check its consistency with established principles of fluid dynamics. We first introduce the system in § 5.1, then analyse its potential hyperbolic structure in § 5.2 to determine whether the dynamics admit finite-speed propagation and, finally, derive the Riemann invariants and jump conditions in § 5.3 to verify compatibility with classical interfacial balance laws.
5.1. First-order system
We study the first-order system associated with (3.25) in the absence of mass transfer (
$\hat {\zeta }_{\alpha } = 0$
) and body forces (
$\boldsymbol{b}=0$
):
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) &=0. \end{align}
The system satisfies the energy law (3.2) with
$\mathscr{D}=0$
. This may be written in the local form
\begin{align} \partial _t \left(\varPsi + \kern5pt \overline {\kern-5pt \mathscr{K}} \kern3.5pt \right) + \textrm {div}\left (\left(\varPsi + \kern5pt \overline {\kern-5pt \mathscr{K}} \kern3.5pt \right)\boldsymbol{v}\right ) = 0. \end{align}
Remark 5.1 (Derivation energy law). The energy law (5.2) follows directly from (5.1) by multiplying with the appropriate weights. In this regard, we note the following identities:
\begin{align} \sum _{{\alpha }} \dfrac {\partial \varPsi }{\partial \tilde {\rho }_{\alpha }} \left ( \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) \right ) &= \partial _t \varPsi + \textrm {div} (\varPsi \boldsymbol{v}) + p \textrm {div} \boldsymbol{v}, \end{align}
\begin{align} \sum _{{\alpha }} \left (-\frac {1}{2}|\boldsymbol{v}|^2\left ( \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) \right )\right ) &= -\frac {1}{2}|\boldsymbol{v}|^2\left ( \partial _t \rho + \textrm {div} (\rho \boldsymbol{v})\right )\!, \end{align}
\begin{align} \boldsymbol{v}\boldsymbol{\cdot }\left ( \partial _t (\rho \boldsymbol{v}) + \textrm {div}(\rho \boldsymbol{v} \otimes \boldsymbol{v}) + \boldsymbol{\nabla }\!p \right ) &= \partial _t \kern5pt \overline {\kern-5pt \mathscr{K}} + \textrm {div}(\kern5pt \overline {\kern-5pt \mathscr{K}} \boldsymbol{v}) + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\!p \nonumber \\ &\quad + \frac {1}{2}|\boldsymbol{v}|^2 \left (\partial _t \rho + \textrm {div} (\rho \boldsymbol{v})\right )\!. \end{align}
Remark 5.2 (Energy law for non-smooth solutions). In case of shocks or other regularities the stability condition takes the form
\begin{align} \partial _t \left(\varPsi + \kern5pt \overline {\kern-5pt \mathscr{K}} \kern3.5pt \right) + \textrm {div}\left (\left(\varPsi + \kern5pt \overline {\kern-5pt \mathscr{K}} + p \right)\boldsymbol{v}\right ) \leqslant 0. \end{align}
which may be established using vanishing-viscosity solutions.
In order to analyse the wave speeds of the model, we deduce the evolution equations of the partial pressures from (5.1b ), i.e.
\begin{align} \partial _t p_{\alpha } + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\!p_{\alpha } + \sum _{{\beta }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 \textrm {div} \boldsymbol{v} =0, \end{align}
where the speed of sound quantities
$a_{{\alpha }{\beta }}$
account for the interaction between constituents:
\begin{align} a_{{\alpha }{\beta }}^2 := \dfrac {\partial p_{\alpha }}{\partial \tilde {\rho }_{\beta }}. \end{align}
The evolution of the classical mixture pressure results from the superposition of the constituent pressure equations
\begin{align} \partial _t p + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\!p + \rho a^2 \textrm {div} \boldsymbol{v}= 0, \end{align}
where the mixture speed of sound
$a$
satisfies
\begin{align} \rho a^2 &= \sum _{{\beta }} \tilde {\rho }_{\beta } a_{{\beta }}^2, \end{align}
and the constituent speed of sound quantities are given by
\begin{align} a_{\beta }^2 &= \dfrac {\partial p}{\partial \tilde {\rho }_{\beta }}= \sum _{{\alpha }} a_{{\alpha }{\beta }}^2. \end{align}
Remark 5.3 (Speed of sound compressible two-phase models). The explicit dependence of the free energy on each constituent (see (3.4)) is reflected in the definitions of the sound speeds
$a_{\beta }$
and
$a_{{\alpha }{\beta }}$
. This approach contrasts with the sound speed definitions typically employed in compressible two-phase flow models (e.g. Baer & Nunziato Reference Baer and Nunziato1986). In the absence of interaction forces (cf.
Remark 4.4
), the relation reduces to
$a_{\beta }^2 = \mathrm{d}p_{\beta }/\mathrm{d}\tilde {\rho }_{\beta }$
, where the standard definition in those models is
$a_{\beta }^2 = \mathrm{d}p_{\beta }/\mathrm{d}\rho _{\beta }$
.
5.2. Hyperbolic structure
In the following we analyse the hyperbolic structure of the first-order model
\begin{align} \partial _t \boldsymbol{v} + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\!\boldsymbol{v} + \rho ^{-1}\boldsymbol{\nabla }\!p &=0, \end{align}
\begin{align} \partial _t p_{\alpha } + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\!p_{\alpha } + \sum _{{\beta }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 \textrm {div} \boldsymbol{v} &=0, \\[9pt] \nonumber \end{align}
where the momentum equation is a direct consequence from (5.1). In this section we consider the one-dimensional case
$d=1$
for conciseness and refer to Appendix D for dimension
$d=3$
. We consider the case
$a_{{\alpha }{\beta }},a_{{\beta }},a\gt 0$
. We write the system (5.10) in the matrix-vector form
\begin{align} \partial _t \boldsymbol{W} + \boldsymbol{A} \partial _x \boldsymbol{W} = 0, \end{align}
where
\begin{align} \boldsymbol{W} = \begin{bmatrix} p_1 \\[-3pt] \vdots \\ p_N\\ v \end{bmatrix}\!, \quad \boldsymbol{A} =\begin{bmatrix} v & 0 & \ldots & & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[5pt] 0 & v & 0 & \ldots & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2\\[-3pt] \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \ldots & 0 & v & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{(N-1){\beta }}^2\\[5pt] 0 & \ldots & & 0 & v & \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ \rho ^{-1} & & \ldots & & \rho ^{-1} & v \end{bmatrix}\!. \end{align}
The system admits
$N+1$
eigenvalues:
\begin{align} \lambda _1 = v - a \lt \lambda _2 = \ldots = \lambda _{N} = v \lt \lambda _{N+1} = v+a. \end{align}
The corresponding right eigenvectors are
\begin{align} \boldsymbol{r}_1 &= \begin{bmatrix} - \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ - \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ a \end{bmatrix}\!, \, \boldsymbol{r}_2 = \begin{bmatrix} - 1\\ 1 \\ 0\\[-6pt] \vdots \\ \\ 0 \end{bmatrix}\!, \, \boldsymbol{r}_3 = \begin{bmatrix} - 1\\ 0 \\ 1\\ 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\, \boldsymbol{r}_N = \begin{bmatrix} - 1\\ 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0 \end{bmatrix}\!, \nonumber\\ \boldsymbol{r}_{N+1} &=\begin{bmatrix} \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ a \end{bmatrix}\!, \end{align}
and the left eigenvectors are
\begin{align}& {\boldsymbol{l}}_1 = \begin{bmatrix} - \dfrac {1}{2 \rho a^2}\\[-3pt] \vdots \\ - \dfrac {1}{2 \rho a^2}\\[3pt] \dfrac {1}{2 a} \end{bmatrix}\!, \, {\boldsymbol{l}}_2 = \begin{bmatrix} - \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\ \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\ 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_3 = \begin{bmatrix} - \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{3{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\ 0\\ \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\ 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \ldots , \end{align}
\begin{align} &\qquad\qquad\quad {\boldsymbol{l}}_N = \begin{bmatrix} - \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\ 0\\[-3pt] \vdots \\ 0\\ \dfrac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{N+1} =\begin{bmatrix} \dfrac {1}{2 \rho a^2}\\[-3pt] \vdots \\ \dfrac {1}{2 \rho a^2}\\[3pt] \dfrac {1}{2 a} \end{bmatrix}\!. \end{align}
The eigenvectors are scaled so that
$\boldsymbol{r}_j\boldsymbol{\cdot }{\boldsymbol{l}}_j=1$
for all
$j=1,\ldots ,N+1$
. Both right and left eigenvectors are linearly independent and span
$\mathbb{R}^{N+1}$
. This implies that the system is hyperbolic if the eigenvalues are real valued, which holds when
$a_{{\alpha }{\beta }}^2 \gt 0$
for all
${\alpha },{\beta } = 1,\ldots ,N$
.
Remark 5.4 (Van der Waals equation of state). The van der Waals equation of state does not, in general, guarantee that
$a_{{\alpha }{\beta }}^2 \gt 0$
, and consequently, the resulting system is not hyperbolic.
Proposition 5.5 (Structure of the waves). The characteristic fields associated with eigenvalues
$\lambda _j, j=2,\ldots ,N$
are linearly degenerate (i.e.
$\boldsymbol{r}_j\boldsymbol{\cdot }\partial \lambda _j/\partial \boldsymbol{W}=0$
), and those associated with eigenvalues
$\lambda _j, j=1,N+1$
are genuinely nonlinear (i.e.
$\boldsymbol{r}_j\boldsymbol{\cdot }\partial \lambda _j/{}\partial \boldsymbol{W}\neq 0$
).
This means that fields associated with eigenvalues
$\lambda _j, j=2,\ldots ,N$
correspond to contact discontinuities, whereas fields associated with eigenvalues
$\lambda _j, j=1,N+1$
correspond to shock or rarefaction waves, respectively.
Proof.
A direct computation shows that the fields associated with eigenvalues
$\lambda _j, j=2,\ldots ,N$
are linearly degenerate (i.e.
$\boldsymbol{r}_j\boldsymbol{\cdot }\partial \lambda _j/\partial \boldsymbol{W}=0$
).
Next we focus on the field associated with eigenvalues
$\lambda _j, j=1,N+1$
. We have
\begin{align} \boldsymbol{r}_1\boldsymbol{\cdot }\partial \lambda _1/\partial \boldsymbol{W}=\boldsymbol{r}_{N+1}\boldsymbol{\cdot }\partial \lambda _{N+1}/\partial \boldsymbol{W}= a+\sum _{{\alpha }{\beta }} \dfrac {\partial a}{\partial p_{\alpha }}\tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2. \end{align}
Next we note that, for an arbitrary variable
$\omega$
, we have
\begin{align} \dfrac {\partial a}{\partial \omega } = \frac {1}{2a \rho }\left (\dfrac {\partial (\rho a^2)}{\partial \omega }-a^2 \frac {\partial \rho }{\partial \omega }\right )\!. \end{align}
Applying this identity for
$\omega = p_{\alpha }$
and substituting (5.8) and (5.9) yields
\begin{align} \dfrac {\partial a}{\partial p_{\alpha }} &= \frac {1}{2a \rho }\left (\sum _{{\gamma }{\tau }} \dfrac {\partial (\tilde {\rho }_{\tau } a_{{\gamma }{\tau }}^2)}{\partial p_{\alpha }}-a^2 \sum _{\tau } \frac {\partial \tilde {\rho }_{\tau }}{\partial p_{\alpha }}\right )\nonumber \\ &= \frac {1}{2a \rho }\left (\sum _{{\tau }{\gamma }} \tilde {\rho }_{\tau } \dfrac {\partial a_{{\gamma }{\tau }}^2}{\partial p_{\alpha }}+\sum _{{\tau }{\gamma }} a_{{\gamma }{\tau }}^2\bigl(a^2\bigr)^{-1}_{{\tau }{\alpha }}-a^2 \sum _{\tau } \bigl(a^2\bigr)^{-1}_{{\tau }{\alpha }}\right )\nonumber \\ &= \frac {1}{2a \rho }\left (\sum _{{\tau }} \tilde {\rho }_{\tau } \dfrac {\partial a_{{\tau }}^2}{\partial p_{\alpha }}+1-a^2 \sum _{{\tau }}\bigl(a^2\bigr)^{-1}_{{\tau }{\alpha }}\right )\!. \end{align}
Inserting into (5.17) gives
\begin{align} \boldsymbol{r}_1\boldsymbol{\cdot }\partial \lambda _1/\partial \boldsymbol{W}&=\boldsymbol{r}_{N+1}\boldsymbol{\cdot }\partial \lambda _{N+1}/\partial \boldsymbol{W}\nonumber \\&= a+\sum _{{\alpha }{\beta }} \dfrac {\partial a}{\partial p_{\alpha }}\tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2\nonumber \\ &= a+\frac {1}{2a \rho }\sum _{{\alpha }{\beta }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 \left (\sum _{{\tau }} \tilde {\rho }_{\tau } \dfrac {\partial a_{{\tau }}^2}{\partial p_{\alpha }}+1-a^2 \sum _{{\tau }}\bigl(a^2\bigr)^{-1}_{{\tau }{\alpha }}\right ) \nonumber \\ &= a+\frac {1}{2a \rho }\left (\sum _{{\alpha }{\beta }{\tau }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 \tilde {\rho }_{\tau } \dfrac {\partial a_{{\tau }}^2}{\partial p_{\alpha }}+\sum _{{\alpha }{\beta }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 - a^2 \sum _{{\tau }}\tilde {\rho }_{\tau }\right ) \nonumber \\ &= a +\frac {1}{2a \rho }\sum _{{\alpha }{\beta }{\tau }} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2 \tilde {\rho }_{\tau } \dfrac {\partial a_{{\tau }}^2}{\partial p_{\alpha }}, \end{align}
which is strictly positive by the aforementioned assumption
$a_{{\alpha }{\beta }},a\gt 0$
in combination with noting that the phase speed of sound increases with the pressure
$\partial a_{{\tau }}/\partial p_{\alpha }\gt 0$
(see also Murrone & Guillard (Reference Murrone and Guillard2005) for a similar assumption).
5.3. Riemann invariants and jump conditions
Next, we derive the Riemann invariants of the first-order system (5.1).
Proposition 5.6. The Riemann invariants associated with the various waves are
\begin{align} \lambda _1: & \quad \left \{\left \{Y_{\alpha }\right \}_{{\alpha }\neq {\beta }}, v+\int _p \frac {1}{\rho a} \textrm {d}p\right \}\!, \end{align}
\begin{align} \lambda _2,\ldots ,\lambda _{N}: & \quad \left \{v,p\right \} ,\end{align}
\begin{align} \lambda _{N+1}: & \quad \left \{\left \{Y_{\alpha }\right \}_{{\alpha }\neq {\beta }}, v-\int _p \frac {1}{\rho a} \textrm {d}p \right \}\!, \end{align}
for some fixed
${\beta } \in \{1,\ldots ,N\}$
.
Proof.
To find the Riemann invariants, we seek for
$\omega _i$
such that
$\partial _{\boldsymbol{W}} \omega _i \boldsymbol{\cdot }\boldsymbol{r}_j = 0$
. Focusing on the Riemann invariants associated with
$\lambda _1$
, we have
\begin{align} 0= \dfrac {\partial \omega _i}{\partial \boldsymbol{W}}\boldsymbol{\cdot }\boldsymbol{r}_1 = a \dfrac {\partial \omega _i}{\partial v} - \sum _{\beta } \tilde {\rho }_{\beta } a_{\beta }^2 \dfrac {\partial \omega _i}{\partial p_{\beta }}. \end{align}
A direct computation reveals that
$Y_{\alpha } = \tilde {\rho }_{\alpha }/\rho$
are Riemann invariants:
\begin{align} \dfrac {\partial Y_{\alpha }}{\partial \boldsymbol{W}}\boldsymbol{\cdot }\boldsymbol{r}_1 &= - \sum _{\beta } \tilde {\rho }_{\beta } a_{\beta }^2 \dfrac {\partial Y_{\alpha }}{\partial p_{\beta }}\nonumber \\ &= -\frac {1}{\rho } \sum _{\beta } \tilde {\rho }_{\beta } a_{\beta }^2\dfrac {\partial \tilde {\rho }_{\alpha }}{\partial p_{\beta }} + \frac {\tilde {\rho }_{\alpha }}{\rho ^2}\sum _{\beta } \tilde {\rho }_{\beta } a_{\beta }^2 \dfrac {\partial \rho }{\partial p_{\beta }}\nonumber \\ &= -\frac {1}{\rho } \tilde {\rho }_{\alpha } + \frac {\tilde {\rho }_{\alpha }}{\rho ^2}\sum _{\beta } \tilde {\rho }_{\beta } \nonumber \\ &= 0. \end{align}
Next, we verify that
$v+\int _p ( {1}/{\rho a}) \textrm {d}p$
is the final Riemann invariant:
\begin{align} \dfrac {\partial v+\int _p \frac {1}{\rho a} \textrm {d}p}{\partial \boldsymbol{W}}\boldsymbol{\cdot }\boldsymbol{r}_1 &= a \dfrac {\partial v+\int _p \frac {1}{\rho a} \textrm {d}p}{\partial v} - \sum _{\alpha } \tilde {\rho }_{\alpha } a_{\alpha }^2 \dfrac {\partial v+\int _p \frac {1}{\rho a} \textrm {d}p}{\partial p_{\alpha }}\nonumber \\ &=a -\sum _{\alpha } \tilde {\rho }_{\alpha } a_{\alpha }^2\dfrac {\partial \int _p \frac {1}{\rho a} \textrm {d}p}{\partial p_{\alpha }} \nonumber \\ &= a-\sum _{\alpha } \tilde {\rho }_{\alpha } a_{\alpha }^2\frac {1}{\rho a} \nonumber \\ &=0. \end{align}
In a similar fashion we find that the Riemann invariants associated with
$\lambda _{N+1}$
are
$Y_{\alpha }$
,
${\alpha } = 1,\ldots ,N$
and
$v-\int _p ( {1}/{\rho a}) \textrm {d}p$
. Finally, we compute the Riemann invariants associated with
$\lambda _2,\ldots ,\lambda _N$
:
\begin{align} 0= \dfrac {\partial \omega _i}{\partial \boldsymbol{W}}\boldsymbol{\cdot }\boldsymbol{r}_j = - \dfrac {\partial \omega _i}{\partial p_1} + \dfrac {\partial \omega _i}{\partial p_j}\quad \text{for }j = 2,\ldots ,N. \end{align}
This is equivalent to
\begin{align} \dfrac {\partial \omega _i}{\partial p_1} = \ldots = \dfrac {\partial \omega _i}{\partial p_N}. \end{align}
It is now easy to verify that
$\omega _i = v,p$
solve system (5.26).
Finally, we note that the conservative nature of the system permits us to uniquely define shock waves due to the Rankine–Hugoniot conditions, i.e.
\begin{align} \left [\tilde {\rho }_{\alpha } (v-\sigma )\right ] &= 0, \quad {\alpha } = 1,\ldots ,N, \end{align}
\begin{align} \left [\rho v(v-\sigma ) + p\right ] &= 0, \end{align}
with
$ [ \phi ] = \phi _R - \phi _L$
the jump between the left and right states of the shock wave and
$\sigma$
the speed of the shock wave.
6. Binary mixtures
In this section we restrict to binary mixtures (
${\alpha } = 1,2)$
for which we discuss alternative – equivalent – forms. The model (4.2) takes for binary mixtures the following form:
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \tilde {\rho }_1\boldsymbol{\nabla }\hat {\mu }_1 + \tilde {\rho }_2\boldsymbol{\nabla }\hat {\mu }_2 - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_1 + \textrm {div}(\tilde {\rho }_1 \boldsymbol{v}) - \textrm {div} \left (M\boldsymbol{\nabla }\left (\hat {\mu }_1-\hat {\mu }_2)\right )\right ) + m \left (\hat {\mu }_1-\hat {\mu }_2)\right ) &=0, \end{align}
\begin{align} \partial _t \tilde {\rho }_2 + \textrm {div}(\tilde {\rho }_2 \boldsymbol{v}) - \textrm {div} \left (M\boldsymbol{\nabla }\left (\hat {\mu }_2-\hat {\mu }_1)\right )\right ) + m \left (\hat {\mu }_2-\hat {\mu }_1)\right ) &=0. \end{align}
Here
$M =- M_{12}=-M_{21}=M_{11}=M_{22}$
and
$m=-m_{12}=-m_{21}=m_{11}=m_{22}$
. By means of variable transformation, we may express the model in terms of the standard mixture quantities. The diffusive fluxes and mass-transfer terms constitute a single quantity by means of the balance conditions (2.11) and (2.18a
):
\begin{align} \hat {\zeta } :&= \hat {\zeta }_1-\hat {\zeta }_2, \quad \quad \hat {\zeta }_1 = \frac {1}{2}\hat {\zeta }, \quad \quad \hat {\zeta }_2= -\frac {1}{2}\hat {\zeta }, \end{align}
\begin{align} \hat {\boldsymbol{H}} :&= \hat {\boldsymbol{H}}_1-\hat {\boldsymbol{H}}_2, \quad \quad \hat {\boldsymbol{H}}_1 = \frac {1}{2}\hat {\boldsymbol{H}}, \quad \quad \hat {\boldsymbol{H}}_2= -\frac {1}{2}\hat {\boldsymbol{H}}. \\[9pt] \nonumber \end{align}
Additionally, we introduce a mass-fraction order parameter
\begin{align} Y :&= Y_1-Y_2, \quad \quad Y_1=\frac {1+Y}{2}, \quad \quad Y_2=\frac {1-Y}{2}, \end{align}
and deduce the relation for the density of the mixture
$\rho$
:
\begin{align}&\qquad\qquad\quad \frac {2}{\rho } =\frac {Y_1}{\tilde {\rho }_1} + \frac {Y_2}{\tilde {\rho }_2}=\frac {1}{\tilde {\rho }_1}\frac {1+Y}{2} + \frac {1}{\tilde {\rho }_2}\frac {1-Y}{2}, \end{align}
\begin{align}& \tilde {\rho }_1 = \rho \frac {1+Y}{2}, \quad \tilde {\rho }_2 =\rho \frac {1-Y}{2}, \quad \rho = \tilde {\rho }_1+\tilde {\rho }_2, \quad Y = \dfrac {\tilde {\rho }_1-\tilde {\rho }_2}{\tilde {\rho }_1+\tilde {\rho }_2}. \end{align}
Furthermore, expressing the model solely in terms of mixture quantities requires the conversion of the component chemical potential quantities. To this purpose, we define
\begin{align} \hat {\mu }_\rho &= \dfrac { \partial \hat {\varPsi }}{\partial \rho } - \textrm {div}\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }, \end{align}
\begin{align} \hat {\mu }_Y &= \dfrac { \partial \hat {\varPsi }}{\partial Y} - \textrm {div}\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }Y}. \end{align}
Additionally, we introduce the following identities.
Lemma 6.1 (Transformation identities). The component chemical potential quantities may be expressed in mixture quantities by means of the identities:
\begin{align} \tilde {\rho }_1 \boldsymbol{\nabla }\hat {\mu }_1+\tilde {\rho }_2\boldsymbol{\nabla }\hat {\mu }_2&=\rho \boldsymbol{\nabla }\hat {\mu }_\rho - \hat {\mu }_Y \boldsymbol{\nabla }Y, \end{align}
\begin{align} \hat {\mu }_1-\hat {\mu }_2 &= 2\frac {\hat {\mu }_Y}{\rho }. \end{align}
Proof. See LemmaB.2.
Applying the variable transformation (6.4) and the identities of Lemma6.1 provides
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \rho \boldsymbol{\nabla }\hat {\mu }_\rho - \hat {\mu }_Y \boldsymbol{\nabla }Y - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\hat {M}\boldsymbol{\nabla }\frac {\hat {\mu }_Y}{\rho }\right ) + \hat {m} \frac {\hat {\mu }_Y}{\rho } &=0, \end{align}
where
$\hat {M}=4M$
and
$\hat {m}=4m$
. Identifying
$ \hat {\varPsi }(\tilde {\rho }_1,\boldsymbol{\nabla }\tilde {\rho }_1,\tilde {\rho }_2,\boldsymbol{\nabla }\tilde {\rho }_2) \equiv \check {\varPsi }(\rho ,\boldsymbol{\nabla }\!\rho , Y,\boldsymbol{\nabla }Y)$
so that
$ \hat {\mu }\rho = \check {\mu }\rho$
and
$ \hat {\mu }_Y = \check {\mu }_Y$
, and setting
$ \hat {M} = \check {M}$
, reveals that (6.7) matches with the models derived in Appendix C.
Furthermore, the model may be expressed using a mass-measured free energy
$\hat {\psi } = \rho ^{-1}\hat {\varPsi }$
. For this purpose, we introduce the chemical potentials
\begin{align} \hat {\nu }_\rho &= \dfrac { \partial \hat {\psi }}{\partial \rho } - \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \hat {\psi }}{\partial \boldsymbol{\nabla }\!\rho }\right )\!, \end{align}
\begin{align} \hat {\nu }_Y &= \dfrac { \partial \hat {\psi }}{\partial Y} - \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }Y}\right )\!, \\[9pt] \nonumber \end{align}
which are related to (6.5) via
$\hat {\mu }_\rho = \rho \hat {\nu }_\rho + \hat {\psi }$
and
$\hat {\mu }_Y = \rho \hat {\nu }_Y$
(see also the proof of PropositionC.2). Inserting these provides the alternative formulation:
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \rho \boldsymbol{\nabla }(\rho \hat {\nu }_\rho + \hat {\psi }) - \rho \hat {\nu }_Y \boldsymbol{\nabla }Y&\nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\hat {M}\boldsymbol{\nabla }\hat {\nu }_Y\right ) + \hat {m} \hat {\nu }_Y &=0. \end{align}
This formulation aligns with the model derived in Appendix C.2 through the identification
$ \hat {\psi }(\tilde {\rho }_1,\boldsymbol{\nabla }\tilde {\rho }_1,\tilde {\rho }_2,\boldsymbol{\nabla }\tilde {\rho }_2) \equiv \check {\check {\psi }}(\rho ,\boldsymbol{\nabla }\!\rho , Y,\boldsymbol{\nabla }Y)$
, which implies that
$ \hat {\nu }_\rho = \check {\check {\mu }}_\rho$
and
$ \hat {\nu }_Y = \check {\check {\mu }}_Y$
, and additionally identifying
$ \hat {M} = {\check {\check {M}}}$
.
7. Connections to existing models
This section aims to establish connections with existing models. Driven by showing parallels between models, we primarily focus on models with the same number of equations. In § 7.1 we compare with established compressible two-phase models, in § 7.2 we examine links to phase-field models for binary fluids and in § 7.3 we highlight similarities with an
$N$
-phase incompressible model.
7.1. Compressible multiphase models
In this subsection we compare with compressible multiphase three-equation models in (7.1.1) and with the Baer-Nunziato model in (7.1.2). For a more complete overview of compressible multiphase models, we refer to the review paper of Saurel & Pantano (Reference Saurel and Pantano2018).
7.1.1. Three equations models
Motivated by revealing similarities between models, here we describe the connection of the first-order system (5.1) to existing three-equation compressible two-phase models.
The model of Dumbser (Reference Dumbser2011) describing a two-phase flow scenario with phase 1 a liquid phase and phase 2 a gas phase, is given by (the model of Dumbser (Reference Dumbser2011) contains an additional term that represents gravitational forces; in this comparison we consider the version without this term):
\begin{align} \partial _t (\rho _1\phi _1 \boldsymbol{v}) + \textrm {div} \left ( \rho _1\phi _1 \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\breve {p}_1 &= 0, \end{align}
\begin{align} \partial _t (\rho _1\phi _1) + \textrm {div}( \rho _1\phi _1 \boldsymbol{v}) &=0, \end{align}
\begin{align} \partial _t \phi _1 + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi _1 &=0. \end{align}
The unknowns are the specific density of phase 1,
$\rho _1$
, the volume fraction of phase 1,
$\phi _1$
, the velocity
$\boldsymbol{v}$
and the pressure of phase 1,
$p_1$
. This pressure is given by the equation of state of the form
$\breve {p}_1 = \breve {p}_1(\rho _1)$
for which a Tait equation of state is employed. We compare this model to the binary version of the first-order system (5.1), i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_1 + \textrm {div}(\tilde {\rho }_1 \boldsymbol{v}) &=0. \end{align}
\begin{align} \partial _t \tilde {\rho }_2 + \textrm {div}(\tilde {\rho }_2 \boldsymbol{v}) &=0, \end{align}
with unknowns
$\tilde {\rho }_1, \tilde {\rho }_2, \boldsymbol{v}$
and
$p = p(\tilde {\rho }_1,\tilde {\rho }_2)$
. We observe the following similarities and differences between the models.
-
(i) Structure equations: the mass balance equation of phase 1 is present in both models: (7.1b ) and (7.2b ). In contrast, the momentum equations differ in the sense that (7.1a ) represents the momentum balance of phase 1, whereas (7.2a ) describes the momentum balance of the mixture. We remark that (7.1) contains a single (specific) density (
$\rho _1$
), for which one might argue that this quantity represents the mixture density (and may be denoted as
$\rho$
). Finally, we note that (7.2b
) may be written as(7.3)where the right-hand side is non-zero in general, showing that it is incompatible with (7.1c ).
\begin{align} \partial _t \phi _1 + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi _1 = - \frac {\phi _1}{\rho _1}(\partial _t \rho _1 + \textrm {div}(\rho _1 \boldsymbol{v})), \end{align}
-
(ii) Symmetry: model (7.1) is not symmetric with respect to the phases in the sense that it contains mass and momentum balance equations with respect to phase 1, but not of phase 2. Of course, the evolution equation (7.1c ) is symmetric due to the saturation constraint
$\phi _1+\phi _2=1$
. -
(iii) Conservative structure: in contrast to (7.2), the model (7.1) is of a non-conservative type due to the interface evolution equation (7.1c ).
-
(iv) Energy law: model (7.2) satisfies the energy law (3.2) with
$\mathscr{D}=0$
. No energy law is provided for (7.1). -
(v) Equation of state: the equations of state of the models,
$\breve {p}_1 = \breve {p}_1(\rho _1)$
for model (7.1) and
$p = p(\tilde {\rho }_1,\tilde {\rho }_2)$
for model (7.2), differ in two key aspects:-
(a) their dependency on a single constituent density quantity versus dependency on each constituent density quantity,
-
(b) the particular form of this/these density quantity/quantities, i.e. specific constituent density versus partial constituent densities.
-
-
(vi) Speed of sound: similarly as the equations of state, the constituent speed of sound of the models,
$a_1^2 = \textrm {d}_{\rho _1}\breve {p}_1$
for model (7.1) and
$a_{{\alpha }{\beta }}^2 = \partial _{\tilde {\rho }_{\beta }}p_{{\alpha }}$
for model (7.2), differ in two key aspects:-
(a) their dependency on a single constituent density quantity versus dependency on each constituent density quantity,
-
(b) the particular form of this/these density quantity/quantities, i.e. specific constituent density versus partial constituent densities.
-
-
(vii) Eigenvalues: the eigenvalues of (7.1) coincide with those of the first-order system (5.1).
-
(viii) Hyperbolicity: due to the adoption of the Tait equation of state in (7.1), this model is hyperbolic when
$\rho _1\gt 0$
and
$\phi _1\gt 0$
, whereas the system (5.1) is hyperbolic when
$a^2_{{\alpha }{\beta }}\gt 0$
for all
${\alpha },{\beta }=1,2$
. -
(ix) Riemann invariants: for model (5.1), these are given in (5.21), while for the model (7.1), these are not provided.
-
(x) Rankine–Hugoniot conditions: for model (5.1), these are given in (5.27), due to the non-conservative nature of model (7.1) these are not well defined.
Next, we study the similarity of model (5.1) with that of Daude et al. (Reference Daude, Galon, Potapov, Beccantini and Mianné2023), which is given by
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\breve {p} &=0, \end{align}
\begin{align} \partial _t (\rho _{\alpha } \phi _{\alpha } ) + \textrm {div}(\rho _{\alpha } \phi _{\alpha } \boldsymbol{v}) &=0\quad \text{ for } {\alpha }= 1,2,3, \\[9pt] \nonumber \end{align}
where the unknowns are
$\boldsymbol{v}, \phi _{\alpha }, \rho _{\alpha }, {\alpha }=1,2,3$
and
$p$
. These quantities are connected through the equations of state
$\breve {p}_{\alpha } = \breve {p}_{\alpha }(\rho _{\alpha })$
in which the constituent pressures
$\breve {p}_{\alpha }$
are assumed individually equal and equal to the mixture pressure
$\breve {p}_1=\breve {p}_2=\breve {p}_3=\breve {p}$
. We compare this model to the ternary version of the first-order system (5.1), i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) &=0 \quad \text{for} \ {\alpha } = 1,2,3, \end{align}
with unknowns
$\tilde {\rho }_{\alpha }, {\alpha } = 1,2,3, \boldsymbol{v}$
and
$p$
and the equation of state
$p = p(\left \{\tilde {\rho }_{\alpha }\right \}_{{\alpha }=1,2,3})$
. We observe the following similarities and differences between the models.
-
(i) Structure equations: through the identity
$\tilde {\rho }_{\alpha }=\rho _{\alpha }\phi _{\alpha }$
we observe that the form of the equations of (7.4) and (7.5) is identical. -
(ii) Symmetry: both model (7.4) and model (7.5) are symmetrical with respect to the phases.
-
(iii) Conservative structure: both model (7.4) and model (7.5) are conservative.
-
(iv) Energy law: for smooth solutions, model (7.4) satisfies the energy law
(7.6)where
\begin{align} \partial _t (\rho u) + \textrm {div}(\rho u \boldsymbol{v} + \breve {p}\boldsymbol{v}) = 0, \end{align}
$\rho u = \sum \tilde {\rho }_{\alpha } u_{\alpha }$
,
$u_{\alpha } = \epsilon _{\alpha } + |\boldsymbol{v}|^2/2$
and
$\textrm {d}\epsilon _{\alpha }/\textrm {d}\rho _{\alpha } = p_{\alpha }/\rho _{\alpha }^2= \breve {p}/\rho _{\alpha }^2$
for some constituent energy quantity
$\epsilon _{\alpha }$
. On the other hand, (7.2) satisfies the energy law (3.2) with
$\mathscr{D}=0$
. This may be written in the local form(7.7)
\begin{align} \partial _t (\rho \psi ) + \textrm {div}(\rho \psi \boldsymbol{v} + p\boldsymbol{v}) = 0. \end{align}
-
(v) Equation of state: we observe the following similarities and differences regarding pressure quantities.
-
(a) Both models contain constituent pressure quantities and a single mixture pressure quantity.
-
(b) The constituent pressure quantities of the models,
$\breve {p}_{\alpha }=\breve {p}_{\alpha }(\rho _{\alpha })$
for model (7.4) and
$p_{\alpha }(\left \{\tilde {\rho }_{\beta }\right \}_{{\beta }=1,2,3})$
for model (7.5), differ in two aspects:-
(bi) their dependency on a single constituent density quantity versus dependency on each constituent density quantity,
-
(bii) the particular form of this/these density quantity/quantities, i.e. specific constituent density versus partial constituent densities.
-
-
(c) The mixture pressure
$\breve {p}$
of (7.4) assumes a mechanical equilibrium, i.e.(7.8)whereas the mixture pressure
\begin{align} \breve {p} = \breve {p}_1=\breve {p}_2=\breve {p}_3, \end{align}
$p$
of model (7.5) is given by Dalton’s law:(7.9)
\begin{align} p = p_1+p_2+p_3. \end{align}
-
-
(vi) Speed of sound: We observe the following similarities and differences regarding speed of sound quantities.
-
(a) Both models contains both constituent speed of sound quantities and a single mixture speed of sound quantity.
-
(b) The constituent speed of sound quantities of the models,
$\hat {a}_{\alpha }^2 = \textrm {d}_{\rho _{\alpha }}\breve {p}_{\alpha }$
for model (7.4) and
$a_{{\alpha }{\beta }}^2 = \partial _{\tilde {\rho }_{\beta }}p_{{\alpha }}$
for model (7.5), differ in two aspects:-
(bi) their dependency on a single constituent density quantity versus dependency on each constituent density quantity,
-
(bii) the particular form of this/these density quantity/quantities, i.e. specific constituent density versus partial constituent densities.
-
-
(c) The mixture speed of sound
$\hat {a}$
of (7.4) satisfies the Wood formula(7.10)whereas the mixture speed of sound
\begin{align} \frac {1}{\rho \hat {a}^2} = \sum _{{\alpha }=1,2,3} \frac {\phi _{\alpha }}{\rho _{\alpha } \hat {a}_{\alpha }^2}, \end{align}
$a$
of model (7.5) is given by(7.11)
\begin{align} \rho a^2 = \sum _{{\alpha },{\beta }=1,2,3} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2. \end{align}
-
-
(vii) Eigenvalues: the eigenvalues of (7.4) and of (7.5) and of the same form, i.e.
(7.12)with
\begin{align} \lambda _1 = v-\bar {a}, \quad \lambda _{2,3} = v,\quad \lambda _4 = v + \bar {a}, \end{align}
$\bar {a} = \hat {a}$
for (7.4) and
$\bar {a}=a$
for (7.5).
-
(viii) Hyperbolicity: the models (7.4) and (7.5) are hyperbolic when the speeds of sound are positive.
-
(ix) Riemann invariants: the Riemann invariants associated with the various waves of the models (7.4) and (7.5) are of the same form, i.e.
(7.13a)
\begin{align} \lambda _1: & \quad \left \{Y_1,Y_2, v+\int _p \frac {1}{\rho \bar {a}} \textrm {d}p \right \} \!,\end{align}
(7.13b)
\begin{align} \lambda _2,\lambda _3: & \quad \left \{v,p\right \}\!, \end{align}
(7.13c)
\begin{align} \lambda _{4}: & \quad \left \{Y_1,Y_2, v-\int _p \frac {1}{\rho \bar {a}} \textrm {d}p \right \}\!, \end{align}
with
$\bar {a} = \hat {a}$
for (7.4) and
$\bar {a}=a$
for (7.5). -
(x) Rankine–Hugoniot conditions: the Rankine–Hugoniot conditions of both model (7.4) and model (7.5) are well defined and match, i.e.
(7.14a)
\begin{align} \left [\tilde {\rho }_{\alpha } (v-\sigma )\right ] &= 0, \quad {\alpha } = 1,2,3 ,\end{align}
(7.14b)
\begin{align} \left [\rho v(v-\sigma ) + \bar {p}\right ] &= 0, \end{align}
with
$\bar {p}=\breve {p}$
for model (7.4) and
$\bar {p}=p$
for model (7.5), where we recall the jump
$ [ \phi ] = \phi _R - \phi _L$
and the speed of the shock wave
$\sigma$
.
Remark 7.1 (Numerical difficulties Wood formula). The Wood formula (7.10) has a non-monotonic variation with volume fraction
$\phi _{\alpha }$
. In Saurel et al. (Reference Saurel, Petitpas and Berry2009) it is argued that this causes numerical difficulties due to inaccuracies in wave transmission across interfaces.
We provide an overview of the comparison in table 1.
Comparison of the first-order versions of the three compressible models with three equations (for binary fluids).
7.1.2. The BN model
Due to its importance, we consider here the seven-equation two-fluid model of Baer & Nunziato (Reference Baer and Nunziato1986) with the source terms
\begin{align} \partial _t \phi _{\alpha } +\boldsymbol{v}_I \!\boldsymbol{\cdot }\boldsymbol{\nabla }\phi _{\alpha } &= \varPhi _{\alpha }, \end{align}
\begin{align} \partial _t(\phi _{\alpha } \rho _{\alpha }) + \textrm {div}(\phi _{\alpha } \rho _{\alpha } \boldsymbol{v}_{\alpha }) &= \varGamma _{\alpha }, \end{align}
\begin{align} \partial _t(\phi _{\alpha } \rho _{\alpha } \boldsymbol{v}_{\alpha }) + \textrm {div}\big (\phi _{\alpha }(\rho _{\alpha } \boldsymbol{v}_{\alpha } \otimes \boldsymbol{v}_{\alpha } + \breve {p}_{\alpha } \boldsymbol I)\big ) - p_I \boldsymbol{\nabla }\phi _{\alpha } &= \boldsymbol{m}_{\alpha } + \varGamma _{\alpha } \boldsymbol{v}_I , \end{align}
\begin{align} \partial _t(\phi _{\alpha } \rho _{\alpha } e_{\alpha }) + \textrm {div} \big (\phi _{\alpha }(\rho _{\alpha } e_{\alpha } + \breve {p}_{\alpha })\boldsymbol{v}_{\alpha } \big ) - p_I \,\boldsymbol{v}_I \!\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha _{\alpha } &= Q_{\alpha } + \boldsymbol{m}_{\alpha }\!\boldsymbol{\cdot }\boldsymbol{v}_I +\varGamma _{\alpha } H_I - p_I \varPhi _{\alpha }, \end{align}
for
${\alpha }=1,2$
. The right-hand side are the source terms with
$\varPhi _{\alpha }$
the source term for the volume-fraction equation,
$\varGamma _{\alpha }$
the mass-transfer term,
$\boldsymbol{m}_{\alpha }$
the drag effect,
$Q_{\alpha }$
the interfacial heat transfer and
$H_I$
an interfacial enthalpy variable. These source terms add up to zero:
$\varPhi _1+\varPhi _2=0$
,
$\varGamma _1+\varGamma _2=0$
,
$\boldsymbol{M}_1+\boldsymbol{M}_2=0$
and
$Q_1+Q_2=0$
. We refer for alternative formulations of the source terms to Guillard & Labois (Reference Guillard and Labois2006) and Zein, Hantke & Warnecke (Reference Zein, Hantke and Warnecke2010). Furthermore,
$\boldsymbol{v}_I=\boldsymbol{v}_I(\boldsymbol{v}_1,\boldsymbol{v}_2)$
and
$p_I=p_I(\breve {p}_1,\breve {p}_2)$
are interfacial velocity and pressure quantities. The energies
$e_{\alpha }$
are given by
$e_{\alpha } = \epsilon _{\alpha } + (1/2)\|\boldsymbol{v}_{\alpha }\|^2$
, where
$\epsilon _{\alpha }=\epsilon _{\alpha }(\rho _{\alpha }, \breve {p}_{\alpha })$
are the internal energies. The volume fractions are assumed to add up to unity, i.e.
$\phi _1+\phi _2 = 1$
, and as such (7.15a
) represents a single equation, making the total number of equations seven. The seven unknowns of this system are
$\phi _1, \rho _1, \rho _2, \breve {p}_1, \breve {p}_2, \boldsymbol{v}_1, \boldsymbol{v}_2$
and the internal energies then follow from the two equations of state,
$\epsilon _1=\epsilon _1(\rho _1, \breve {p}_1)$
and
$\epsilon _2=\epsilon _2(\rho _2, \breve {p}_2)$
.
Remark 7.2 (Extensions of BN models). The BN framework admits extensions that incorporate viscous stresses and heat conduction, generalisations to
$N$
-fluid mixtures and surface tension effects; see ,e.g. Saurel & Pantano (
Reference Saurel and Pantano2018
). Since most formulations in the literature concern the classical inviscid two-fluid (seven-equation) system (7.15), we confine the discussion here to that setting.
We compare to the binary, first-order (inviscid) formulation of our model without gravity, but with mass transfer, i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_1 + \textrm {div}(\tilde {\rho }_1 \boldsymbol{v}) &= \zeta _1 = - m_{11} \hat {\mu }_1- m_{12} \hat {\mu }_2, \end{align}
\begin{align} \partial _t \tilde {\rho }_2 + \textrm {div}(\tilde {\rho }_2 \boldsymbol{v}) &= \zeta _2 = - m_{21} \hat {\mu }_1- m_{22} \hat {\mu }_2, \end{align}
with chemical potential
$\hat {\mu }_{\alpha }= \partial \varPsi / \partial \tilde {\rho }_{\alpha }$
and mixture pressure
$p= p(\tilde {\rho }_1,\tilde {\rho }_2)=\sum _{{\beta }=1,2} (\hat {\mu }_{\beta }\tilde {\rho }_{\beta }) -\hat {\varPsi }$
which reduces to
$p = \sum _{{\alpha }=1,2,{\beta }=1,2} \tilde {\rho }_{\alpha } \tilde {\rho }_{\beta } (\partial \hat {\psi }_{\beta }/\partial \tilde {\rho }_{\alpha })$
with
$\hat {\varPsi }_{\beta } = \tilde {\rho }_{\beta } \hat {\psi }_{\beta }$
.
We observe the following differences and similarities:
-
(i) Structure equations: as is common in compressible two-phase flow models, the model of Baer & Nunziato (Reference Baer and Nunziato1986) contains an evolution equation of the volume fraction. Equation (7.15a ) does not appear here in the proposed model (7.16), as detailed in (7.3). Next, the mass balance equations of the BN model may be written as
(7.17a)
\begin{align} \partial _t \tilde {\rho }_1 + \textrm {div}(\tilde {\rho }_1 \boldsymbol{v}) + \textrm {div}(\tilde {\rho }_1 \boldsymbol{w}_1) &= \varGamma _1, \end{align}
(7.17b)
\begin{align} \partial _t \tilde {\rho }_2 + \textrm {div}(\tilde {\rho }_2 \boldsymbol{v}) + \textrm {div}(\tilde {\rho }_2 \boldsymbol{w}_2) &= \varGamma _2, \end{align}
which reveal that the difference lies in (i) the presence of the peculiar velocity
$\boldsymbol{w}_{\alpha }$
, defined in (2.9); and (ii) possibly the mass-transfer terms. In order to compare the momentum equations, we add (7.15c
) for
${\alpha } = 1,2$
and find that(7.18)The difference due of the models regarding the terms with the peculiar velocities in the mass and momentum equations is standard when comparing models with multiple momentum equations with models with a single momentum equation; see, e.g. ten Eikelder et al. (Reference ten Eikelder, van der Zee, Akkerman and Schillinger2023). Apart from this term, the mass balance equations and momentum equations coincide between the two models when
\begin{align} \partial _t(\rho \boldsymbol{v}) + \textrm {div}\left (\rho \boldsymbol{v} \otimes \boldsymbol{v}\right) + \textrm {div}\left (\sum _{\alpha }\phi _{\alpha }\rho _{\alpha } \boldsymbol{w}_{\alpha } \otimes \boldsymbol{w}_{\alpha } \right ) + \boldsymbol{\nabla }\left ( \sum _{\alpha } \phi _{\alpha }\breve {p}_{\alpha } \right ) &= 0. \end{align}
$p$
coincides with
$\sum _{\alpha } \phi _{\alpha } \breve {p}_{\alpha }$
. This holds for the perfect gas equation of state. In this regard we recall (4.7b
) and refer to the discussion on the equation of state below. The absence of energy equations in the proposed model (7.16) is due to the isothermal assumption.
-
(ii) Symmetry: the symmetry of model (7.15) depends on the form of the source terms and the interfacial quantities
$\boldsymbol{v}_I$
and
$p_I$
. Usually the source terms
$\varPhi _{\alpha }, \varGamma _{\alpha }, \boldsymbol{m}_{\alpha }$
and
$Q_{\alpha }$
are symmetric with respect to the phases (i.e. a renumbering does not change the model). However, this does not always hold for the interfacial quantities
$\boldsymbol{v}_I$
and
$p_I$
. Namely, typical choices for these quantities include single phase values, i.e.
$\boldsymbol{v}_I = \boldsymbol{v}_1$
or
$ \boldsymbol{v}_I =\boldsymbol{v}_2$
and
$p_I = \breve {p}_1$
or
$p_I = \breve {p}_2$
; see Zein et al. (Reference Zein, Hantke and Warnecke2010). -
(iii) Conservative structure: we consider the systems without source terms. In contrast to the conservative model (7.2), the BN model (7.15) is of a non-conservative term in (i) to the interface evolution (7.15a ), (ii) the momentum equation (7.15c ), and (iii) the energy equation (7.15d ). As mentioned in the introduction, these non-conservative terms vanish across genuinely nonlinear waves.
-
(iv) Energy/entropy law: the BN model (7.15) satisfies a mixture entropy law; see, e.g. Zein et al. (Reference Zein, Hantke and Warnecke2010). The energy-dissipation law (3.2) (with
$\mathscr{D}=0$
) of model (7.16) also follows from the second law of thermodynamics for isothermal mixtures with a single momentum equation. -
(v) Equation of state: a key distinction between the models lies in their equations of state. Clearly, the equation of state of the (7.15) model is of the form
$\epsilon _{\alpha }=\epsilon _{\alpha }(\breve {p}_{\alpha },\rho _{\alpha })$
, whereas in the (7.16) model the barotropic equation of state
$p=p(\tilde {\rho }_1,\tilde {\rho }_2)$
is adopted. Other key differences are:-
(a) their dependency on single constituent quantities (in the BN model (7.15)) versus dependency on quantities of each constituent (in the proposed model (7.16)),
-
(b) the particular form of density quantity/quantities, i.e. specific constituent density (in the BN model (7.15)) versus partial constituent densities (in the proposed model (7.16)).
To emphasise the structural parallels between the BN model (7.15) and the proposed model (7.16), we restrict attention to the case
(7.19)which encompasses the perfect gas equation of state but excludes, for example, the van der Waals equation of state. In this setting, the constituent pressure becomes
\begin{align} \hat {\varPsi }_{\alpha } = \hat {\varPsi }_{\alpha }(\tilde {\rho }_{\alpha }), \end{align}
(7.20)Hence, this quantity is of the form
\begin{align} p_{{\alpha }} = \tilde {\rho }_{\alpha }^2 \,\dfrac {\partial \hat {\psi }_{\alpha }}{\partial \tilde {\rho }_{\alpha }}. \end{align}
$p_{\alpha } = \phi _{\alpha } \bar {p}_{\alpha }$
, which closely resembles the pressure quantity
$\phi _{\alpha } \breve {p}_{\alpha }$
in the BN model (7.15). Moreover, when the perfect gas equation of state is adopted the pressures match. For the BN model (7.15), the constituent perfect gas equation of state is represented by
$\breve {p}_{\alpha } = R_{{\alpha }} \theta \rho _{\alpha }$
, where the constituent perfect gas equation of state in the proposed model (7.16) is
$p_{\alpha } = R_{{\alpha }} \theta \tilde {\rho }_{\alpha }$
, with
$R_{\alpha }$
the specific gas constant of constituent
$\alpha$
. In this case we find that(7.21)
\begin{align} p_{\alpha } = \phi _{\alpha } \breve {p}_{\alpha }. \end{align}
-
-
(vi) Speed of sound: we observe the following similarities and differences regarding speed of sound quantities.
-
(a) Both models contain constituent speed of sound quantities, for model (7.15),
(7.22)and
\begin{align} \breve {a}_{\alpha }^2 = \rho _{\alpha }^{-1} (\partial _{\breve {p}_{\alpha }} \breve {\epsilon }_{\alpha })^{-1}\left ( \dfrac {p_{\alpha }}{\rho _{\alpha }} - \rho _{\alpha } \partial _{\rho _{\alpha }}\breve {\epsilon }_{\alpha } \right ) \end{align}
$a_{{\beta }}^2 = \partial _{\tilde {\rho }_{\beta }}p=\sum _{\alpha } a_{{\alpha }{\beta }}^2, a_{{\alpha }{\beta }}^2 = \partial _{\tilde {\rho }_{\beta }}p_{{\alpha }}$
for model (7.16), which differ in three aspects:-
(ai) their dependency on a single constituent density quantity versus dependency on each constituent density quantity,
-
(aii) the particular form of this/these density quantity/quantities, i.e. specific constituent density versus partial constituent densities,
-
(aiii) the form of the equation due to the isothermal assumption in the proposed model (7.16).
-
-
(b) The mixture speed of sound of model (7.15) is not defined/applicable, whereas the mixture speed of sound
$a$
of model (7.16) is given by(7.23)
\begin{align} \rho a^2 = \sum _{{\beta }=1,2} \tilde {\rho }_{\beta } a_{{\beta }}^2 = \sum _{{\alpha },{\beta }=1,2} \tilde {\rho }_{\beta } a_{{\alpha }{\beta }}^2. \end{align}
-
-
(vii) Eigenvalues: the number of eigenvalues in one spatial dimension of the two systems is different due to the number of equations; for model (7.15), we have
$\lambda _1=v_I, \lambda _{2,4}=v_1\pm c_1, \lambda _3=v_1, \lambda _{5,7}=v_2\pm c_2, \lambda _6=v_2$
, whereas those of the proposed model (7.16) are
$\lambda _{1,3} = v\pm c$
and
$\lambda _2 = v$
. -
(viii) Hyperbolicity: when the speed of sounds are positive, the proposed model (7.16) is hyperbolic and the BN model (7.15) is hyperbolic if additionally
$v_I \neq v_1,v_2$
. This property is sometimes called ‘weakly hyperbolic’; see, e.g. Coquel et al. (Reference Coquel, Hérard and Saleh2017). -
(ix) Mass transfer: the symmetry and vanishing-sum conditions on
$m_{{\alpha }{\beta }}$
yield the following mass-transfer terms: (7.24)
\begin{align} \hat {\zeta }_{1} &= -m_{12}\,(\hat {\mu }_1-\hat {\mu }_2), \end{align}
(7.25)
\begin{align} \hat {\zeta }_{2} &= -m_{12}\,(\hat {\mu }_2-\hat {\mu }_1). \end{align}
In contrast, the BN model (7.15) typically expresses mass transfer in the form
(7.26)
\begin{align} \varGamma _1 &= -c_{12}(\breve {g}_1-\breve {g}_2), \end{align}
(7.27)
\begin{align} \varGamma _2 &= -c_{12}(\breve {g}_2-\breve {g}_1), \end{align}
where the quantity
$c$
has the dependency
$c=c(\tilde {\rho }_1,\tilde {\rho }_2)$
and where(7.28)with
\begin{align} \breve {g}_{\alpha }= \epsilon _{\alpha } - \theta _{\alpha } \eta _{\alpha } + \frac {\breve {p}_{\alpha }}{\rho _{\alpha }}, \end{align}
$\theta _{\alpha }$
the temperature and
$\eta _{\alpha }$
the entropy of fluid
$\alpha$
(see, e.g. Zein et al. Reference Zein, Hantke and Warnecke2010). These two formulations are, in general, not equivalent. Again, to highlight structural similarities between the BN model (7.15) and the proposed model (7.16), we restrict attention to the case (7.19). In this scenario we have the identities (7.29)
\begin{align} \hat {\mu }_{{\alpha }} &= \hat {\psi }_{\alpha } + \tilde {\rho }_{\alpha } \dfrac {\partial \hat {\psi }_{\alpha }}{\partial \tilde {\rho }_{\alpha }}, \end{align}
(7.30)
\begin{align} \hat {p}_{{\alpha }} &= \tilde {\rho }_{\alpha }^2 \dfrac {\partial \hat {\psi }_{\alpha }}{\partial \tilde {\rho }_{\alpha }}, \end{align}
so that
(7.31)where we have inserted the Gibbs relation
\begin{align} \hat {\mu }_{{\alpha }} &= \epsilon _{\alpha } - \theta _{\alpha } \eta _{\alpha } + \frac {p_{\alpha }}{\tilde {\rho }_{\alpha }}, \end{align}
$\hat {\psi }_{\alpha } = \epsilon _{\alpha } - \theta _{\alpha } \eta _{\alpha }$
. We can now identify
$\hat {\mu }_{\alpha }$
with
$g_{\alpha }$
when the constituent pressure satisfies
$p_{\alpha } = \breve {p}_{\alpha } \phi _{\alpha }$
; see (7.21). Finally, by identifying
$m_{12}$
with
$c_{12}$
one can obtain identical mass transfer in this case.
-
(x) Rankine–Hugoniot conditions: for model (7.16), these are given in (7.14), while for the BN model (7.15), these are
(7.32a)
\begin{align} \left [ \phi _{\alpha }\right ] &=0 ,\end{align}
(7.32b)
\begin{align} \left [ \rho _{\alpha } (v_{\alpha }-\sigma )\right ] &= 0, \end{align}
(7.32c)
\begin{align} \left [ \rho _{\alpha } v_{\alpha }(v_{\alpha }-\sigma ) + \breve {p}_{\alpha } \right ] &= 0, \end{align}
(7.32d)
\begin{align} \left [ \rho _{\alpha } e_{\alpha } (v_{\alpha }-\sigma ) + \breve {p}_{\alpha } v_{\alpha } \right ] &= 0, \end{align}
for
${\alpha } = 1,2$
.
7.2. Phase-field models for binary fluids
We show resemblance of model (6.1) with two existing binary fluid models in the literature. By using LemmaA.2, the model (6.7) converts into
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\check {p} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \right )&\nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\hat {M}\boldsymbol{\nabla }\frac {\hat {\mu }_Y}{\rho }\right ) + \hat {m} \frac {\hat {\mu }_Y}{\rho } &=0, \end{align}
with
$p = \rho \hat {\mu }_\rho - \hat {\varPsi }$
. By setting
$\hat {m}=0$
and inserting arbitrary source terms on the right-hand side of (7.33b
) and (7.33c
), the model coincides with that of Mukherjee & Gomez (Reference Mukherjee and Gomez2024).
Next, through the identification
$\hat {\varPsi }=\rho \hat {\psi }$
, the model (6.9) converts into
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\hat {p} + \textrm {div} \left (\rho \boldsymbol{\nabla }Y\otimes \dfrac {\partial \hat { \psi }}{\partial \boldsymbol{\nabla }Y} +\rho \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \hat { \psi }}{\partial \boldsymbol{\nabla }\!\rho } \right )&\nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\hat {M}\boldsymbol{\nabla }\hat {\nu }_Y\right ) + \hat {m} \hat {\nu }_Y &=0, \end{align}
with
$\hat {p} = \rho \hat {\mu }_\rho - \hat {\varPsi }$
. By setting
$\hat {m}=0$
and taking
$\psi = \psi _0(\rho ,Y) + \epsilon \|\boldsymbol{\nabla }Y\|^2/2$
, the model converts into the compressible Navier–Stokes Cahn–Hilliard model proposed in Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998), i.e.
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\breve {p} + \textrm {div} \left (\rho \boldsymbol{\nabla }Y\otimes \boldsymbol{\nabla }Y \right )&\nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\hat {M}\boldsymbol{\nabla }\left ( \dfrac {\partial \psi _0}{\partial Y} - \frac {\epsilon }{\rho } \textrm {div}(\rho \boldsymbol{\nabla }Y)\right )\right ) &=0, \end{align}
with
$\breve {p} = \rho ^2 \partial \psi _0/\partial \rho$
. We emphasise here that the dependence of
$ \check {\check {\psi }}$
on
$\boldsymbol{\nabla }\!\rho$
is not taken into account in Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998). In other words, in terms of gradient dependency, in Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998) depends on
$\boldsymbol{\nabla }Y = \boldsymbol{\nabla }((\tilde {\rho }_1-\tilde {\rho }_2)/(\tilde {\rho }_1+\tilde {\rho }_2))$
, whereas in (6.9) the free energy depends on
$\boldsymbol{\nabla }\tilde {\rho }_1$
and
$\boldsymbol{\nabla }\tilde {\rho }_2$
.
Remark 7.3 (Matching first-order systems). The corresponding first-order systems of Mukherjee & Gomez ( Reference Mukherjee and Gomez2024 ) and Lowengrub & Truskinovsky ( Reference Lowengrub and Truskinovsky1998 ) have not been analysed in those papers. However, differences between model (6.1) and those models vanish for first-order systems, meaning that the properties in § 5 directly carry over.
We sketch an overview of the comparison of the binary models in table 2. The present model shares many structural similarities with the formulations of Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998) and Mukherjee & Gomez (Reference Mukherjee and Gomez2024), particularly at the first-order level, where distinctions largely arise from higher-order terms. Depending on the specific modelling choices (e.g. mobility), in some scenarios the differences may diminish as the interface thickness tends to zero; however, they remain relevant in problems with finite interface thickness or non-trivial fluid mixtures, where the models can yield significantly different predictions.
Comparison binary models.
$^*$
The model of Mukherjee & Gomez (Reference Mukherjee and Gomez2024) is equipped with an energy law when the source terms in the model are set to zero.
$^{**}$
The first-order systems match; see Remark7.3.
7.3. Incompressible
$N$
-phase flows
We compare model (4.2) with the incompressible
$N$
-phase model given in ten Eikelder (Reference ten Eikelder2025), which reads, for phases (constituents)
${\alpha } = 1, \ldots , N$
,
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \sum _{\beta } \tilde {\rho }_{\beta } \boldsymbol{\nabla }g_{\beta } - \textrm {div} \left ( \nu (2 \boldsymbol{\nabla} ^s \boldsymbol{v}+\bar {\lambda }(\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) +\textrm {div} \breve {\boldsymbol{H}}_{\alpha } - \breve {\zeta }_{\alpha }&=0, \end{align}
\begin{align} \textrm {div}\boldsymbol{v} + \displaystyle \sum _{{\beta }} \rho _{\beta }^{-1} \boldsymbol{\nabla }\boldsymbol{\cdot }\breve {\boldsymbol{H}}_{\beta } - \displaystyle \sum _{{\beta }}\rho _{\beta }^{-1}\breve {\zeta }_{\beta } &= 0, \end{align}
\begin{align} \breve {\boldsymbol{H}}_{\alpha } + \sum _{\beta } \breve {\boldsymbol{M}}_{{\alpha }{\beta }}\boldsymbol{\nabla }g_{\beta }&=0, \end{align}
\begin{align} \breve {\zeta }_{\alpha } + \sum _{\beta } \breve {m}_{{\alpha }{\beta }}g_{\beta }&=0, \end{align}
with the chemical potential quantities
\begin{align} g_{\alpha } &= \dfrac {\hat {\mu }_{\alpha }^\phi + \lambda }{\rho _{\alpha }}, \end{align}
\begin{align} \hat {\mu }_{\alpha }^\phi &= \dfrac { \partial \hat {\varPsi }}{\partial \phi _{\alpha }} - \textrm {div}\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\phi _{\alpha }}. \end{align}
Model (7.36) is incompressible in the sense that
$\rho _{\alpha }$
is constant is time and space, i.e.
$\rho _{\alpha }(\boldsymbol{x},t) = \rho _{\alpha }$
. As such, directly enforcing the saturation condition (2.5b
) reveals that
$\left \{\tilde {\rho }_{\alpha }(\boldsymbol{x},t)=\rho _{\alpha }\phi _{\alpha }(\boldsymbol{x},t)\right \}_{{\alpha }=1,\ldots ,N}$
constitute
$N-1$
independent variables. Instead, formulation (7.36) enforces the saturation constraint through (7.36c
), therefore avoiding ambiguities. The set of variables is completed, analogously to a single-phase system, with Lagrange multiplier
$\lambda$
that enforces (7.36c
).
Obviously, since models (3.25) and (4.2) have largely been derived in a similar fashion, they both share key features:
$N$
mass balance laws,
$1$
momentum balance law and similar closure models for diffusive fluxes and mass-transfer terms. The key difference lies in the presence of (7.36c
). This manifests itself in the form of the chemical potentials, namely
$g_{\alpha } = \hat {\mu }_{\alpha } + \rho _{\alpha }^{-1}\lambda$
. The model may be written into the form
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\!p + \textrm {div} \left (\sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right )& \nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \tilde {\rho }_{\alpha } + \textrm {div}(\tilde {\rho }_{\alpha } \boldsymbol{v}) +\textrm {div} \breve {\boldsymbol{H}}_{\alpha } - \breve {\zeta }_{\alpha } &=0, \end{align}
\begin{align} \textrm {div}\boldsymbol{v} + \displaystyle \sum _{{\beta }} \rho _{\beta }^{-1} \boldsymbol{\nabla }\boldsymbol{\cdot }\breve {\boldsymbol{H}}_{\beta } - \displaystyle \sum _{{\beta }}\rho _{\beta }^{-1}\breve {\zeta }_{\beta } &= 0, \end{align}
where
$p = \lambda + \sum _{\alpha } \tilde {\rho }_{\alpha } \hat {\mu }_{\alpha } - \hat {\varPsi }$
, which closely resembles (3.25). This shows that the pressure in (7.38) is the superposition of the thermodynamical pressure
$\sum _{\alpha } \tilde {\rho }_{\alpha } \hat {\mu }_{\alpha } - \hat {\varPsi }$
and the Lagrange multiplier
$\lambda$
.
The observed strong connection is in line with expectations (just as for single-phase flow), reflecting the natural link between compressible and incompressible models. However, this does not hold for existing compressible two-phase flow models, where establishing their incompressible counterparts is not straightforward.
8. Conclusion
In this work we have developed a thermodynamical framework for compressible, isothermal
$N$
-phase mixtures. The framework is formulated within the principles of continuum mixture theory, which serves as a fundamental basis for modelling multi-physics systems. We have established some connections with existing phase-field and compressible two-phase flow models in an effort to help bridge the gap between these fields of research. Specifically, the proposed formulation extends the NSK model to multiple fluids and adapts the Navier–Stokes-Cahn–Hilliard/Allen–Cahn (NSCH/AC) models to general
$N$
-phase compressible mixtures. A key distinction from existing compressible two-phase flow models is that reducing the complexity of those models (i.e. using fewer equations) often leads to the loss or degradation of essential thermodynamic properties. In contrast, the framework proposed in this work maintains thermodynamic consistency regardless of its level of complexity.
Several open challenges remain for future research. We highlight some key directions. First, the framework offers a foundation for studying spinodal decomposition in compressible mixtures, with potential applications in phase-change problems such as evaporation and condensation. Second, by incorporating free energies with explicit interface width parameters, one can investigate sharp-interface limits, both formally and rigorously. Other important extensions include the incorporation of chemical reactions and non-isothermal effects. Furthermore, an essential avenue for future work, beyond the scope of this paper, is the development of efficient numerical schemes that accurately capture interfacial dynamics while preserving key thermodynamic properties.
Acknowledgements
The authors thank the reviewers for their constructive comments and suggestions.
Funding
MtE acknowledges support from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) via project EI 1210/5-1, number 566600860. DS gratefully acknowledges support from the DFG via the Emmy Noether Award SCH 1249/2-1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Identities Korteweg tensors
Lemma A.1 (Identity Korteweg stresses). We have the following identity for the pressure and Korteweg stresses:
\begin{align} \boldsymbol{\nabla }\!p + \textrm {div} \left ( \sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right ) = \sum _{\alpha } \tilde {\rho }_{\alpha }\boldsymbol{\nabla }\hat {\mu }_{\alpha }. \end{align}
Here we recall that
$p=\sum _{\alpha }\hat {\mu }_{\alpha }\tilde {\rho }_{\alpha }-\hat {\varPsi }$
.
Proof. This follows from the sequence of identities
\begin{align} &\boldsymbol{\nabla }\!p + \textrm {div} \left ( \sum _{\alpha } \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }} \right )\nonumber \\ &=\boldsymbol{\nabla }\left (\sum _{\alpha }\hat {\mu }_{\alpha }\tilde {\rho }_{\alpha }-\hat {\varPsi }\right ) + \sum _{\alpha } \left (\boldsymbol{\nabla }\tilde {\rho }_{\alpha } \textrm {div} \left (\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right ) + (\boldsymbol{H} \tilde {\rho }_{\alpha })\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right )\nonumber \\ &=\sum _{\alpha } \left (\boldsymbol{\nabla }\left (\hat {\mu }_{\alpha }\tilde {\rho }_{\alpha } \right ) - \dfrac {\partial \hat {\varPsi }}{\partial \tilde {\rho }_{\alpha }}\boldsymbol{\nabla }\tilde {\rho }_{\alpha } - (\boldsymbol{H} \tilde {\rho }_{\alpha }) \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}+ \boldsymbol{\nabla }\tilde {\rho }_{\alpha } \textrm {div} \left (\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right ) + (\boldsymbol{H} \tilde {\rho }_{\alpha })\dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_{\alpha }}\right )\nonumber \\ &=\sum _{\alpha } \tilde {\rho }_{\alpha }\boldsymbol{\nabla }\hat {\mu }_{\alpha }, \end{align}
where
$\boldsymbol{H}\tilde {\rho }_{\alpha }$
denotes the Hessian of the partial density
$\tilde {\rho }_{\alpha }$
.
Lemma A.2 (Identity Korteweg stresses – binary fluids – volume measure). The following identity holds for the pressure and Korteweg stresses in binary fluids:
\begin{align} \boldsymbol{\nabla }\check {p} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \right )= \rho \boldsymbol{\nabla }\check {\mu }_\rho -\check {\mu }_Y\boldsymbol{\nabla }Y. \end{align}
Here we recall that
$\check {p} = \check {\mu }_\rho \rho - \check {\varPsi }$
.
Proof. This follows from the sequence of identities
\begin{align} &\boldsymbol{\nabla }\check {p} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \right ) \nonumber \\ &=\boldsymbol{\nabla }(\check {\mu }_\rho \rho ) - \boldsymbol{\nabla }\check {\varPsi } + \boldsymbol{\nabla }Y \textrm {div}\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} + (\boldsymbol{H}\!Y)\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} + \boldsymbol{\nabla }\!\rho \textrm {div}\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + (\boldsymbol{H}\!\rho )\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\nonumber \\ &=\boldsymbol{\nabla }(\check {\mu }_\rho \rho ) -\boldsymbol{\nabla }Y\dfrac {\partial \check { \varPsi }}{\partial Y}-\boldsymbol{\nabla }\!\rho \dfrac {\partial \check { \varPsi }}{\partial \rho } + \boldsymbol{\nabla }Y \textrm {div}\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} + \boldsymbol{\nabla }\!\rho \textrm {div}\dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \nonumber \\ &=\rho \boldsymbol{\nabla }\check {\mu }_\rho -\check {\mu }_Y\boldsymbol{\nabla }Y, \end{align}
with
$\boldsymbol{H}\!\rho$
and
$\boldsymbol{H}Y$
the Hessians
$\rho$
and
$Y$
.
Lemma A.3 (Identity Korteweg stresses – binary fluids – mass measure). The following identity holds for the pressure and Korteweg stresses in binary fluids:
\begin{align} \boldsymbol{\nabla }{\check {\check {p}}} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \rho \dfrac {\partial \check {\check {\psi }}} {\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \rho \dfrac {\partial \check {\psi }}{\partial \boldsymbol{\nabla }\!\rho } \right )= \rho \boldsymbol{\nabla }\left(\rho \check {\check {\mu }}_\rho + \check {\check {\psi }}\right)-\rho \check {\check {\mu }}_Y\boldsymbol{\nabla }Y. \end{align}
Here we recall that
$ {\check {\check {p}}} = \rho ^2 \check {\check {\mu }}_\rho$
.
Proof. This follows from the sequence of identities
\begin{align} &\boldsymbol{\nabla }{\check {\check {p}}} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \rho \dfrac {\partial \check {\psi }}{\partial \boldsymbol{\nabla }\!\rho } \right ) \nonumber \\ &=\rho ^2 \boldsymbol{\nabla }\check {\check {\mu }}_\rho + 2\rho \check {\check {\mu }}_\rho \boldsymbol{\nabla }\!\rho + \boldsymbol{\nabla }Y \textrm {div}\left (\rho \dfrac {\partial \check { \psi }}{\partial \boldsymbol{\nabla }Y} \right )+ (\boldsymbol{H}\!Y)\rho \dfrac {\partial \check { \psi }}{\partial \boldsymbol{\nabla }Y} + \boldsymbol{\nabla }\!\rho \textrm {div}\left (\rho \dfrac {\partial \check {\check { \psi }}}{\partial \boldsymbol{\nabla }\!\rho } \right ) \nonumber \\ &\quad + (\boldsymbol{H}\!\rho )\rho \dfrac {\partial \check {\check { \psi }}}{\partial \boldsymbol{\nabla }\!\rho }\nonumber \\ &=\rho ^2 \boldsymbol{\nabla }\check {\check {\mu }}_\rho + 2\rho \check {\check {\mu }}_\rho \boldsymbol{\nabla }\!\rho - \rho \check {\check {\mu }}_\rho \boldsymbol{\nabla }\!\rho + \rho \dfrac {\partial \check {\check { \psi }}}{\partial \rho } + \rho \boldsymbol{\nabla }\check {\check {\psi }} - \rho \dfrac {\partial \check {\check { \psi }}}{\partial \rho } \boldsymbol{\nabla }\!\rho - \rho \dfrac {\partial \check {\check { \psi }}}{\partial Y} \boldsymbol{\nabla }Y - \rho \check {\check {\mu }}_Y \boldsymbol{\nabla }Y \nonumber \\ &\quad + \rho \dfrac {\partial \check {\check { \psi }}}{\partial Y} \boldsymbol{\nabla }Y \nonumber \\ &=\rho \boldsymbol{\nabla }\left(\rho \check {\check {\mu }}_\rho + \check {\check {\psi }}\right)-\rho \check {\check {\mu }}_Y\boldsymbol{\nabla }Y. \end{align}
Appendix B. Variable transformation – binary mixture
In this section we provide details on the transformation of constituent quantities to mixture quantities for a binary mixture. The variable transformations are given by
\begin{align} \tilde {\rho }_1 = \rho \frac {1+Y}{2}, \quad \tilde {\rho }_2 =\rho \frac {1-Y}{2}, \quad \rho = \tilde {\rho }_1+\tilde {\rho }_2, \quad Y = \dfrac {\tilde {\rho }_1-\tilde {\rho }_2}{\tilde {\rho }_1+\tilde {\rho }_2}. \end{align}
The derivatives of the scalar quantities are
\begin{align} \dfrac {\partial \tilde {\rho }_1}{\partial \rho } = \dfrac {1+Y}{2}, \quad & \dfrac {\partial \tilde {\rho }_2}{\partial \rho } = \dfrac {1-Y}{2}, \end{align}
\begin{align} \dfrac {\partial \tilde {\rho }_1}{\partial Y} = \dfrac {\rho }{2}, \quad & \dfrac {\partial \tilde {\rho }_2}{\partial Y} =-\dfrac {\rho }{2}, \end{align}
\begin{align} \dfrac {\partial \rho }{\partial \tilde {\rho }_1} = 1, \quad &\dfrac {\partial \rho }{\partial \tilde {\rho }_2} = 1, \end{align}
\begin{align} \dfrac {\partial Y}{\partial \tilde {\rho }_1} = \dfrac {2\tilde {\rho }_2}{(\tilde {\rho }_1+\tilde {\rho }_2)^2}, \quad &\dfrac {\partial Y}{\partial \tilde {\rho }_2} = -\dfrac {2\tilde {\rho }_1}{(\tilde {\rho }_1+\tilde {\rho }_2)^2}. \end{align}
The derivatives of the gradient quantities are
\begin{align}&\qquad\qquad\qquad \dfrac {\partial \boldsymbol{\nabla }\!\rho }{\partial \tilde {\rho }_1} = \dfrac {\partial \boldsymbol{\nabla }\!\rho }{\partial \tilde {\rho }_2} = 0, \end{align}
\begin{align}& \dfrac {\partial \boldsymbol{\nabla }Y}{\partial \tilde {\rho }_1} = -4 \frac {\tilde {\rho }_2}{(\tilde {\rho }_1+\tilde {\rho }_2)^3}\boldsymbol{\nabla }\tilde {\rho }_1 + 2 \frac {\tilde {\rho }_1-\tilde {\rho }_2}{(\tilde {\rho }_1+\tilde {\rho }_2)^3}\boldsymbol{\nabla }\tilde {\rho }_2, \end{align}
\begin{align}& \dfrac {\partial \boldsymbol{\nabla }Y}{\partial \tilde {\rho }_2} = 2 \frac {\tilde {\rho }_1-\tilde {\rho }_2}{(\tilde {\rho }_1+\tilde {\rho }_2)^3}\boldsymbol{\nabla }\tilde {\rho }_1 + 4 \frac {\tilde {\rho }_1}{(\tilde {\rho }_1+\tilde {\rho }_2)^3}\boldsymbol{\nabla }\tilde {\rho }_2. \end{align}
Lemma B.1 (Relation Korteweg stresses). The Korteweg stresses are related via
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_1 \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_1}+\boldsymbol{\nabla }\tilde {\rho }_2 \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_2} = \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}. \end{align}
Proof. The chain rule and the variable transformations (B2) provide
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_1 &= \frac {1+Y}{2}\boldsymbol{\nabla }\!\rho + \frac {\rho }{2} \boldsymbol{\nabla }Y, \end{align}
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_2 &= \frac {1-Y}{2}\boldsymbol{\nabla }\!\rho - \frac {\rho }{2} \boldsymbol{\nabla }Y. \end{align}
Next, the chain rule and the variable transformations (B3) provide
\begin{align} \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_1} &= \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} \dfrac {1-Y}{\rho }, \end{align}
\begin{align} \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_2} &= \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } - \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} \dfrac {1+Y}{\rho }. \end{align}
Combining the above expressions gives
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_1 \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_1} &= \frac {1+Y}{2} \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \frac {1-Y}{2}\boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}\nonumber \\ &\quad + \frac {1-Y^2}{2\rho } \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} + \frac {\rho }{2} \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }, \end{align}
\begin{align} \boldsymbol{\nabla }\tilde {\rho }_2 \otimes \dfrac {\partial \hat {\varPsi }}{\partial \boldsymbol{\nabla }\tilde {\rho }_2} &= \frac {1-Y}{2} \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \frac {1+Y}{2}\boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}\nonumber \\ &\quad - \frac {1-Y^2}{2\rho } \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} - \frac {\rho }{2} \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }. \end{align}
The addition of these two identities provides the result.
Lemma B.2 (Relations chemical potentials). The chemical potentials are related via
\begin{align} \hat {\mu }_1\tilde {\rho }_1+\hat {\mu }_2\tilde {\rho }_2 &= \hat {\mu }_\rho \rho , \end{align}
\begin{align} \tilde {\rho }_1 \boldsymbol{\nabla }\hat {\mu }_1+\tilde {\rho }_2\boldsymbol{\nabla }\hat {\mu }_2&=\rho \boldsymbol{\nabla }\hat {\mu }_\rho - \hat {\mu }_Y \boldsymbol{\nabla }Y, \end{align}
\begin{align} \hat {\mu }_1-\hat {\mu }_2 &= 2\frac {\hat {\mu }_Y}{\rho }. \end{align}
Proof. The first and third identity follow from the chain rule:
\begin{align} \hat {\mu }_\rho &= \hat {\mu }_1 \frac {\partial \tilde {\rho }_1}{\partial \rho }+\hat {\mu }_2 \frac {\partial \tilde {\rho }_2}{\partial \rho } = \frac {1}{\rho }(\hat {\mu }_1 \tilde {\rho }_1+\hat {\mu }_2 \tilde {\rho }_2), \end{align}
\begin{align} \hat {\mu }_Y &= \hat {\mu }_1 \frac {\partial \tilde {\rho }_1}{\partial Y}+\hat {\mu }_2 \frac {\partial \tilde {\rho }_2}{\partial Y} = \frac {\rho }{2} (\hat {\mu }_1 -\hat {\mu }_2). \end{align}
Taking the gradient of the first identity yields
\begin{align} \hat {\mu }_1 \boldsymbol{\nabla }\tilde {\rho }_1+\hat {\mu }_2 \boldsymbol{\nabla }\tilde {\rho }_2 + \tilde {\rho }_1\boldsymbol{\nabla }\hat {\mu }_1 +\tilde {\rho }_2\boldsymbol{\nabla }\hat {\mu }_2&= \hat {\mu }_\rho \boldsymbol{\nabla }\!\rho +\rho \boldsymbol{\nabla }\hat {\mu }_\rho . \end{align}
Noting the identity
\begin{align} \boldsymbol{\nabla }\hat {\varPsi } = \hat {\mu }_1 \boldsymbol{\nabla }\tilde {\rho }_1+\hat {\mu }_2 \boldsymbol{\nabla }\tilde {\rho }_2 = \hat {\mu }_\rho \boldsymbol{\nabla }\!\rho + \hat {\mu }_Y \boldsymbol{\nabla }Y \end{align}
completes the proof.
Appendix C. Constitutive modelling binary mixture
This section presents alternative derivations of the modelling framework for binary fluids. In Appendix C.1 we detail the derivation using a volume-measure-based free energy (
$\varPsi$
), while Appendix C.2 covers the derivation based on a mass-measure-based free energy (
$\psi$
), which are related via
$\varPsi = \rho \psi$
. We utilise the system of balance laws (3.1) in mixture variables, i.e.
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) + \textrm {div} \boldsymbol{H} &= \zeta , \end{align}
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) - \textrm {div} \boldsymbol{T} - \rho \boldsymbol{b} &=0, \end{align}
\begin{align} \boldsymbol{T}-\boldsymbol{T}^T &=0, \end{align}
with
$Y=Y_1-Y_2$
,
$\boldsymbol{H} = \boldsymbol{H}_1-\boldsymbol{H}_2$
and
$\zeta = \zeta _1-\zeta _2$
. In both cases we use the energy-dissipative design condition (3.2) where the energy is given in (3.3).
C.1. Volume-measure free energy
We postulate the free energy to pertain to the constitutive class
\begin{align} \varPsi = \check {\varPsi }(\rho , \boldsymbol{\nabla }\!\rho , Y,\boldsymbol{\nabla }Y), \end{align}
and introduce the chemical potentials
\begin{align} \check {\mu }_Y &= \dfrac { \partial \check {\varPsi }}{\partial Y} - \textrm {div}\dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}, \end{align}
\begin{align} \check {\mu }_\rho &= \dfrac { \partial \check {\varPsi }}{\partial \rho } - \textrm {div}\dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }. \end{align}
We proceed with the evaluation of the evolution of the energy (C2). Following the same steps as in § 3.2, we apply the Reynolds transport theorem, the divergence theorem and subsequently expand the derivatives to find that
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \check {\varPsi } \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \dfrac {\partial \check {\varPsi }}{\partial \rho } \dot {\rho } + \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\boldsymbol{\cdot }\dot {\overline {\boldsymbol{\nabla }\!\rho }} + \dfrac {\partial \check {\varPsi }}{\partial Y} \dot {Y} + \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}\boldsymbol{\cdot }\dot {\overline {\boldsymbol{\nabla }Y}} + \check {\varPsi } \,\textrm {div} \boldsymbol{v}\,\textrm {d}v, \end{align}
where we recall the notation
\begin{align} \dot {\overline {\boldsymbol{\nabla }\omega }}:=\partial _t \boldsymbol{\nabla }\omega + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\left (\boldsymbol{\nabla }\omega \right ) \end{align}
for the material derivative of the quantity
$\boldsymbol{\nabla }\omega$
. By substituting the identity (3.8) for
$\omega = Y, \rho$
, and subsequently integrating by parts, we arrive at
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \check {\varPsi } \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \check {\mu }_\rho \dot {\rho } + \check {\mu }_Y \dot {Y} - \left (\boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} + \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\right ): \boldsymbol{\nabla }\!\boldsymbol{v} + \check {\varPsi }\,\textrm {div} \boldsymbol{v}\,\textrm {d}v \nonumber \\ &\quad+ \displaystyle \int _{\partial \mathcal{R}(t)}\left (\dot {Y} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}+\dot {\rho } \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a. \end{align}
Substituting the mass balance (C1a ) and (C1b ) and applying integration by parts leads to
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \check {\varPsi } \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} -\check {\mu }_\rho \rho \textrm {div}\boldsymbol{v}+\boldsymbol{\nabla }\left (\!\dfrac {\check {\mu }_Y}{\rho }\!\right ) \boldsymbol{\cdot }\boldsymbol{H} - \left (\!\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }+\boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}\!\right ): \boldsymbol{\nabla }\!\boldsymbol{v} \nonumber \\ &\quad + \check {\varPsi }\textrm {div} \boldsymbol{v} + \left (\!\dfrac {\check {\mu }_Y}{\rho }\!\right )\zeta \textrm {d}v+ \!\displaystyle \int _{\partial \mathcal{R}(t)}\!\left (\!-\rho ^{-1}\check {\mu }_Y \boldsymbol{H}+\dot {\rho } \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }+\dot {Y} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}\!\right )\boldsymbol{\cdot }\boldsymbol{\nu } \textrm {d}a. \end{align}
Taking the sum of (C7) and (3.12), we obtain
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t} \mathscr{E} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T\boldsymbol{T}-\rho ^{-1}\check {\mu }_Y \boldsymbol{H}+\dot {Y} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}+\dot {\rho } \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a \nonumber \\ &\quad - \displaystyle \int _{\mathcal{R}(t)} \left ( \boldsymbol{T} +\check {\mu }_\rho \rho \textrm {div}\boldsymbol{v}-\boldsymbol{\nabla }\left (\dfrac {\check {\mu }_Y}{\rho }\right )\boldsymbol{\cdot }\boldsymbol{H} + \left (\boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}+\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\right ): \boldsymbol{\nabla }\!\boldsymbol{v} \right . \nonumber \\ &\quad \left .- \check {\varPsi }\,\textrm {div} \boldsymbol{v} - \left (\dfrac {\check {\mu }_Y}{\rho }\right )\zeta \right )\,\textrm {d}v, \end{align}
in which we identify the rate of work and the dissipation as
\begin{align} \mathscr{W} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T\boldsymbol{T}-\left (\dfrac {\check {\mu }_Y}{\rho }\right ) \boldsymbol{H}+\dot {Y} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}+\dot {\rho } \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho }\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a, \end{align}
\begin{align} \mathscr{D} &= \displaystyle \int _{\mathcal{R}(t)} \left ( \boldsymbol{T} + \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y}+\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \left (\check {\mu }_\rho \rho -\check {\varPsi }\right )\boldsymbol{I} \right ): \boldsymbol{\nabla }\!\boldsymbol{v} \nonumber \\ &\quad -\boldsymbol{\nabla }\left (\dfrac {\check {\mu }_Y}{\rho }\right ) \boldsymbol{\cdot }\boldsymbol{H} - \left (\dfrac {\check {\mu }_Y}{\rho }\right ) \zeta \,\textrm {d}v. \end{align}
Due to the arbitrariness of the control volume, the energy-dissipation law holds when the local inequality is satisfied:
\begin{align} \left ( \boldsymbol{T} + \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} + \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } + \left (\check {\mu }_\rho \rho -\check {\varPsi }\right )\boldsymbol{I} \right ): \boldsymbol{\nabla }\!\boldsymbol{v} -\boldsymbol{\nabla }\left (\!\dfrac {\check {\mu }_Y}{\rho }\!\right ) \boldsymbol{\cdot }\boldsymbol{H} - \left (\!\dfrac {\check {\mu }_Y}{\rho }\!\right ) \zeta \geqslant 0. \end{align}
Based on the inequality (C10), we now restrict ourselves to stress tensors
$\boldsymbol{T}$
, diffusive fluxes
$\boldsymbol{h}$
and mass fluxes
$\zeta$
belonging to the constitutive classes:
\begin{align} \boldsymbol{T} &= \check {\boldsymbol{T}}(\boldsymbol{\nabla }\!\boldsymbol{v}, \rho ,\boldsymbol{\nabla }\!\rho ,Y, \boldsymbol{\nabla }Y, \check {\mu }_\rho , \boldsymbol{\nabla }\check {\mu }_\rho , \check {\mu }_Y, \boldsymbol{\nabla }\check {\mu }_Y), \end{align}
\begin{align} \boldsymbol{H} &= \check {\boldsymbol{H}}(\rho ,\boldsymbol{\nabla }\!\rho ,Y, \boldsymbol{\nabla }Y, \check {\mu }_\rho , \boldsymbol{\nabla }\check {\mu }_\rho , \check {\mu }_Y, \boldsymbol{\nabla }\check {\mu }_Y), \end{align}
\begin{align} \zeta &= \check {\zeta }(\rho ,Y, \check {\mu }_\rho , \check {\mu }_Y). \\[9pt] \nonumber \end{align}
We select closure models that render each of the terms in (C10) separately non-negative, i.e.
\begin{align} \check {\boldsymbol{T}} &= - \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} - \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } - \left (\check {\mu }_\rho \rho -\check {\varPsi }\right )\boldsymbol{I} + \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}), \end{align}
\begin{align} \check {\boldsymbol{H}} &= -\check {\boldsymbol{M}} \boldsymbol{\nabla }\left (\dfrac {\check {\mu }_Y}{\rho }\right )\!, \end{align}
\begin{align} \check {\zeta } &= -\check {m}\left (\dfrac {\check {\mu }_Y}{\rho }\right )\!, \end{align}
a positive definite mobility tensor
$ \check {\boldsymbol{M}}$
with the same dependencies as (C11b
), satisfying
$ \check {\boldsymbol{M}}|_{Y=0,1} = 0$
, and a non-negative scalar mobility
$ \check {m}$
with the same dependencies as (C11c
), also satisfying
$ \check {m}|_{Y=0,1} = 0$
.
This completes the constitutive modelling, and the system is given by
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }\check {p} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \right ) & \nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\check {\boldsymbol{M}}\boldsymbol{\nabla }\left (\dfrac {\check {\mu }_Y}{\rho }\right )\right ) + \check {m} \left (\dfrac {\check {\mu }_Y}{\rho }\right ) &=0, \end{align}
\begin{align} \check {\mu }_Y - \dfrac {\partial \check {\varPsi }}{\partial Y}+\textrm {div} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} &=0, \end{align}
\begin{align} \check {\mu }_\rho - \dfrac { \partial \check {\varPsi }}{\partial \rho } + \textrm {div}\dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } &=0, \end{align}
with the pressure
$\check {p} =\check {\mu }_\rho \rho -\check {\varPsi }$
. Analogously to Lemma4.1, the Korteweg stress and pressure contributions in the momentum equation may be expressed using the chemical potential quantities.
Lemma C.1 (Identity Korteweg stresses – binary fluids – volume measure). The following identity holds for the pressure and Korteweg stresses in binary fluids:
\begin{align} \boldsymbol{\nabla }\check {p} + \textrm {div} \left ( \boldsymbol{\nabla }Y\otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }Y} +\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check { \varPsi }}{\partial \boldsymbol{\nabla }\!\rho } \right )= \rho \boldsymbol{\nabla }\check {\mu }_\rho -\check {\mu }_Y\boldsymbol{\nabla }Y. \end{align}
Proof. See LemmaA.2.
Applying the lemma converts the model into the compact form:
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \rho \boldsymbol{\nabla }\check {\mu }_\rho -\check {\mu }_Y\boldsymbol{\nabla }Y - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left (\check {\boldsymbol{M}}\boldsymbol{\nabla }\left (\dfrac {\check {\mu }_Y}{\rho }\right )\right ) + \check {m} \left (\dfrac {\check {\mu }_Y}{\rho }\right ) &=0, \end{align}
\begin{align} \check {\mu }_Y - \dfrac {\partial \check {\varPsi }}{\partial Y}+\textrm {div} \dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }Y} &=0, \end{align}
\begin{align} \check {\mu }_\rho - \dfrac { \partial \check {\varPsi }}{\partial \rho } + \textrm {div}\dfrac {\partial \check {\varPsi }}{\partial \boldsymbol{\nabla }\!\rho } &=0. \end{align}
C.2. Mass-measure free energy
The mass-measure free energy is postulated to belong to the class
\begin{align} \psi = \check {\check {\psi }}(\rho ,\boldsymbol{\nabla }\!\rho , Y,\boldsymbol{\nabla }Y), \end{align}
and the chemical potential quantities defined as
\begin{align} \check {\check {\mu }}_Y &= \dfrac { \partial \check {\check {\psi }}}{\partial Y} - \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right )\!, \end{align}
\begin{align} \check {\check {\mu }}_\rho &= \dfrac { \partial \check {\check {\psi }}}{\partial \rho } - \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }\right )\!. \end{align}
We proceed with the evaluation of the evolution of the energy (C16) and find, via the Reynolds transport theorem and the mass balance (C1a ), that
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \rho \check {\check {\psi }} \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \rho \dot { \check {\check {\psi }}} \,\textrm {d}v\nonumber \\ &= \displaystyle \int _{\mathcal{R}(t)} \rho \left (\dfrac {\partial \check {\check {\psi }}}{\partial \rho } \dot {\rho } + \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }\boldsymbol{\cdot }\dot {\overline {\boldsymbol{\nabla }\!\rho }} +\dfrac {\partial \check {\check {\psi }}}{\partial Y} \dot {Y} + \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\boldsymbol{\cdot }\dot {\overline {\boldsymbol{\nabla }Y}}\right ) \,\textrm {d}v. \end{align}
On the account of the relation
\begin{align} \dot {\overline {\boldsymbol{\nabla }\omega }} = \boldsymbol{\nabla }(\dot {\omega }) - (\boldsymbol{\nabla }\omega )^T\boldsymbol{\nabla }\!\boldsymbol{v}, \end{align}
for
$\omega = \rho , Y$
, we arrive at
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \rho \check {\check {\psi }} \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \rho \check {\check {\mu }}_Y \dot {Y} +\rho \check {\check {\mu }}_\rho \dot {\rho } - \left (\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \rho \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\boldsymbol{\nabla }Y \otimes \dfrac {\partial \rho \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right ): \boldsymbol{\nabla }\!\boldsymbol{v} \,\textrm {d}v \nonumber \\ &\quad+ \displaystyle \int _{\partial \mathcal{R}(t)} \left (\dot {\rho } \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\dot {Y} \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a. \end{align}
By substituting the mass balance laws (C1a ) and (C1b ) we find that
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t}\displaystyle \int _{\mathcal{R}(t)} \rho \check {\check {\psi }} \,\textrm {d}v &= \displaystyle \int _{\mathcal{R}(t)} \boldsymbol{\nabla }\check {\check {\mu }}_Y \boldsymbol{\cdot }\boldsymbol{H} -\rho ^2 \check {\check {\mu }}_\rho \textrm {div}\boldsymbol{v} + \check {\check {\mu }}_Y\zeta \nonumber \\ &\quad - \left (\boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \rho \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\boldsymbol{\nabla }Y \otimes \dfrac {\partial \rho \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right ): \boldsymbol{\nabla }\!\boldsymbol{v}\,\textrm {d}v \nonumber \\ &\quad+ \displaystyle \int _{\partial \mathcal{R}(t)}\left (- \check {\check {\mu }}_Y \boldsymbol{H}+\dot {\rho } \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\dot { c}\rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }c}\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a. \end{align}
Following the same procedure as in § 3 we find that
\begin{align} \dfrac {\textrm {d}}{\textrm {d}t} \mathscr{E} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T {\check {\check {\boldsymbol{T}}}}- \check {\check {\mu }}_Y \boldsymbol{H}+\dot {\rho } \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\dot {Y} \rho \dfrac {\partial \check {\check {\psi }}|}{\partial \boldsymbol{\nabla }Y}\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a \nonumber \\ &\quad -\! \displaystyle \int _{\mathcal{R}(t)} \!\left (\! {\check {\check {\boldsymbol{T}}}} + \rho ^2 \check {\check {\mu }}_\rho \boldsymbol{I} + \rho \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho } + \rho \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} \!\right )\!: \boldsymbol{\nabla }\!\boldsymbol{v} -\boldsymbol{\nabla }\check {\check {\mu }}_Y \boldsymbol{\cdot }\boldsymbol{H} - \check {\check {\mu }}_Y\zeta \textrm {d}v, \end{align}
where we identify the rate of work and dissipation term:
\begin{align} \mathscr{W} &= \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^T \check {\check {\boldsymbol{T}}}- \check {\check {\mu }}_Y\boldsymbol{H}+\dot {\rho } \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+\dot {Y} \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right )\boldsymbol{\cdot }\boldsymbol{\nu } \,\textrm {d}a, \end{align}
\begin{align} \mathscr{D} &= \displaystyle \int _{\mathcal{R}(t)} \left ( {\check {\check {\boldsymbol{T}}}} +\rho ^2 \check {\check {\mu }}_\rho \boldsymbol{I}+ \rho \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+ \rho \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} \right ): \boldsymbol{\nabla }\!\boldsymbol{v} -\boldsymbol{\nabla }\check {\check {\mu }}_Y \boldsymbol{\cdot }\boldsymbol{H} - \check {\check {\mu }}_Y\zeta \,\textrm {d}v. \end{align}
From the design condition
$\mathscr{D}\geqslant 0$
we obtain the local modelling restriction:
\begin{align} \left ( {\check {\check {\boldsymbol{T}}}} +\rho ^2 \check {\check {\mu }}_\rho \boldsymbol{I} + \rho \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }+ \boldsymbol{\nabla }Y \otimes \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} \right ): \boldsymbol{\nabla }\!\boldsymbol{v} -\boldsymbol{\nabla }\check {\check {\mu }}_Y \boldsymbol{\cdot }\boldsymbol{H} - \check {\check {\mu }}_Y\zeta \geqslant 0. \end{align}
Based on the modelling restriction (C24) we introduce with the constitutive classes:
\begin{align} \boldsymbol{T} &= {\check {\check {\boldsymbol{T}}}}(\boldsymbol{\nabla }\!\boldsymbol{v}, \rho ,\boldsymbol{\nabla }\!\rho ,Y, \boldsymbol{\nabla }Y, \check {\check {\mu }}_\rho , \boldsymbol{\nabla }\check {\check {\mu }}_\rho , \check {\check {\mu }}_Y, \boldsymbol{\nabla }\check {\check {\mu }}_Y), \end{align}
\begin{align} \boldsymbol{H} &= {\check {\check {\boldsymbol{H}}}}(\rho ,\boldsymbol{\nabla }\!\rho ,Y, \boldsymbol{\nabla }Y, \check {\check {\mu }}_\rho , \boldsymbol{\nabla }\check {\check {\mu }}_\rho , \check {\check {\mu }}_Y, \boldsymbol{\nabla }\check {\check {\mu }}_Y), \end{align}
\begin{align} \zeta &= {\check {\check {\zeta }}}(\rho , Y, \check {\check {\mu }}_\rho , \check {\check {\mu }}_Y). \end{align}
We choose closure models that ensure each term in (C10) remains separately non-negative:
\begin{align} \check {\boldsymbol{T}} &= - \rho \boldsymbol{\nabla }Y \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} - \rho \boldsymbol{\nabla }\!\rho \otimes \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho } -\rho ^2 \check {\check {\mu }}_\rho \boldsymbol{I} + \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}), \end{align}
\begin{align} {\check {\check {\boldsymbol{H}}}} &= - {\check {\check {\boldsymbol{M}}}} \boldsymbol{\nabla }\check {\check {\mu }}_Y, \end{align}
\begin{align} {\check {\check {\zeta }}} &= - {\check {\check {m}}} \check {\check {\mu }}_Y. \end{align}
Here the mobility tensor
$ {\check {\check {\boldsymbol{M}}}}$
is positive definite, depends on the same variables as (C25b
) and satisfies
$ {\check {\check {\boldsymbol{M}}}}|_{Y=0,1}=0$
. The scalar mobility
$ {\check {\check {m}}} \geqslant 0$
shares the same dependencies as (C25c
) and satisfies
$ {\check {\check {m}}}|_{Y=0,1}=0$
.
This concludes the constitutive modelling, with the system given by
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \boldsymbol{\nabla }{\check {\check {p}}} + \textrm {div} \left ( \boldsymbol{\nabla }\!\rho \otimes \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho } +\boldsymbol{\nabla }Y\otimes \rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y} \right ) & \nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left ( {\check {\check {\boldsymbol{M}}}}\boldsymbol{\nabla }\check {\check {\mu }}_Y\right ) + {\check {\check {m}}} \check {\check {\mu }}_Y &=0, \end{align}
\begin{align} \check {\check {\mu }}_Y -\dfrac { \partial \check {\check {\psi }}}{\partial Y} + \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }Y}\right )&=0, \end{align}
\begin{align} \check {\check {\mu }}_\rho - \dfrac { \partial \check {\check {\psi }}}{\partial \rho } + \frac {1}{\rho }\textrm {div}\left (\rho \dfrac {\partial \check {\check {\psi }}}{\partial \boldsymbol{\nabla }\!\rho }\right ) &=0, \end{align}
with the pressure
$ {\check {\check {p}}} = \rho ^2 \check {\check {\mu }}_\rho$
.
Proposition C.2 (Equivalence binary mixture models). The mixture models (C15) and (C27) are equivalent.
Proof.
This follows from identifying
$ \check {\varPsi } = \rho \check {\check {\psi }}$
,
$ \check {\boldsymbol{M}} = {\check {\check {\boldsymbol{M}}}}$
and
$ \check {m} = {\check {\check {m}}}$
. In particular, the first identification provides
$ \check {\mu }_Y = \rho \check {\check {\mu }}_Y$
and the identification of the pressure follows from
\begin{align} \check {p} = \rho \check {\mu }_\rho - \check {\varPsi } = \rho (\rho \check {\check {\mu }}_\rho + \check {\check {\psi }}) - \rho \check {\check {\psi }} = \rho ^2 \check {\check {\mu }}_\rho = {\check {\check {p}}}, \end{align}
with
$\check {\mu }_\rho =\rho \check {\check {\mu }}_\rho + \check {\check {\psi }}$
.
Finally, utilising LemmaA.3 the model may be written in the compact form:
\begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm {div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \rho \boldsymbol{\nabla }(\rho \check {\check {\mu }}_\rho + \check {\check {\psi }}) - \rho \check {\check {\mu }}_Y \boldsymbol{\nabla }Y&\nonumber \\ - \textrm {div} \left ( \nu (2\boldsymbol{D}+\lambda (\textrm {div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{g} &= 0, \end{align}
\begin{align} \partial _t \rho + \textrm {div}(\rho \boldsymbol{v}) &= 0, \end{align}
\begin{align} \partial _t (\rho Y) + \textrm {div}(\rho Y \boldsymbol{v}) - \textrm {div} \left ( {\check {\check {M}}}\boldsymbol{\nabla }\check {\check {\mu }}_Y\right ) + {\check {\check {m}}} \check {\check {\mu }}_Y &=0. \end{align}
Appendix D. Hyperbolicity in
$d$
dimensions
We write the system in the matrix-vector form
\begin{align} \partial _t \boldsymbol{W} + \boldsymbol{A}_{x_1} \partial _{x_1} \boldsymbol{W} + \boldsymbol{A}_{x_2} \partial _{x_2} \boldsymbol{W} + \boldsymbol{A}_{x_3} \partial _{x_3} \boldsymbol{W} = 0, \end{align}
where
\begin{align} \boldsymbol{W} = \begin{bmatrix} p_1 \\[3pt] \vdots \\[3pt] p_N\\[3pt] v_{x_1} \\[3pt] v_{x_2} \\[3pt] v_{x_3} \end{bmatrix}\!, \quad &\boldsymbol{A}_{x_1} =\begin{bmatrix} v_{x_1} & 0 & \ldots & & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2 & 0 & 0\\[3pt] 0 & v_{x_1} & 0 & \ldots & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2& 0 & 0\\[3pt] \vdots & \ddots & \ddots & \ddots & \vdots & \vdots & 0 & 0\\[3pt] 0 & \ldots & 0 & v_{x_1} & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{(N-1){\beta }}^2& 0 & 0\\[3pt] 0 & \ldots & & 0 & v_{x_1} & \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2& 0 & 0\\[3pt] \rho ^{-1} & & \ldots & & \rho ^{-1} & v_{x_1} & 0 & 0 \\[3pt] 0 & & & &\ldots & 0& v_{x_1} & 0 \\[3pt] 0 & & & &\ldots & & 0& v_{x_1} \end{bmatrix}\!,\nonumber \\[3pt] &\boldsymbol{A}_{x_2} =\begin{bmatrix} v_{x_2} & 0 & \ldots & & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2 & 0\\[3pt] 0 & v_{x_2} & 0 & \ldots & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2 & 0\\[3pt] \vdots & \ddots & \ddots & \ddots & \vdots & 0 & \vdots & 0\\[3pt] 0 & \ldots & 0 & v_{x_2} & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{(N-1){\beta }}^2 & 0 \\[3pt] 0 & \ldots & & 0 & v_{x_2} & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2 & 0\\[3pt] 0 & & \ldots & & 0 &v_{x_2} &0 & 0 \\[3pt] \rho ^{-1} & & \ldots & & \rho ^{-1} & 0 & v_{x_2} & 0 \\[3pt] 0 & & & &\ldots & & 0& v_{x_2} \end{bmatrix}\!,\nonumber \\[3pt] &\boldsymbol{A}_{x_3} =\begin{bmatrix} v_{x_3} & 0 & \ldots & & 0 & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[3pt] 0 & v_{x_3} & 0 & \ldots & 0 & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2\\[3pt] \vdots & \ddots & \ddots & \ddots & \vdots & 0 & 0 & \vdots \\[3pt] 0 & \ldots & 0 & v_{x_3} & 0 & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{(N-1){\beta }}^2\\[3pt] 0 & \ldots & & 0 & v_{x_3} & 0 & 0 & \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\[3pt] 0 & & \ldots & &0 &v_{x_3} &0 & 0 \\[3pt] 0 & & &\ldots & &0 &v_{x_3} & 0 \\[3pt] \rho ^{-1} & & \ldots & & \rho ^{-1} & 0 & 0 & v_{x_3} \end{bmatrix}\!. \end{align}
The matrices
$\boldsymbol{A}_{x_i}, i=1,2,3$
have
$N+d$
eigenvalues:
\begin{align} \lambda _{i,1} = v_i - a \lt \lambda _{i,2} = \ldots = \lambda _{i,N+d-1} = v_i \lt \lambda _{i,N+d} = v_i+a. \end{align}
The corresponding right eigenvectors are
\begin{align} \boldsymbol{r}_{x_1,1} &= \begin{bmatrix} - \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ - \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ a\\0\\0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_1,2} = \begin{bmatrix} - 1\\ 1 \\ 0\\[-6pt] \vdots \\ \\ 0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_1,3} = \begin{bmatrix} - 1\\ 0 \\ 1\\ 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\, \boldsymbol{r}_{x_1,N} = \begin{bmatrix} - 1\\ 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0\\0\\0 \end{bmatrix}\!, \nonumber \\&\qquad \boldsymbol{r}_{x_1,N+1} = \begin{bmatrix} 0\\[-6pt] \vdots \\ 0\\ 1\\ 0 \end{bmatrix}\!,\, \boldsymbol{r}_{x_1,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 0\\ 1 \end{bmatrix}\!, \, \boldsymbol{r}_{x_1,N+3} =\begin{bmatrix} \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ a\\0\\0 \end{bmatrix}\!, \end{align}
\begin{align} \boldsymbol{r}_{x_2,1} &= \begin{bmatrix} - \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ - \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ 0\\a\\0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_2,2} = \begin{bmatrix} - 1\\ 1 \\ 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_2,3} = \begin{bmatrix} - 1\\ 0 \\ 1\\ 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\, \boldsymbol{r}_{x_2,N} = \begin{bmatrix} - 1\\ 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0\\0\\0 \end{bmatrix}\!, \nonumber \\&\qquad \boldsymbol{r}_{x_2,N+1} = \begin{bmatrix} 0\\[-6pt] \vdots \\ 1\\ 0\\ 0 \end{bmatrix}\!,\, \boldsymbol{r}_{x_2,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 0\\ 1 \end{bmatrix}\! , \, \boldsymbol{r}_{x_2,N+3} =\begin{bmatrix} \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ 0\\a\\0 \end{bmatrix}\!,\\ \nonumber \end{align}
\begin{align} \boldsymbol{r}_{x_3,1} &= \begin{bmatrix} - \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ - \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ 0\\0\\a \end{bmatrix}\!, \, \boldsymbol{r}_{x_3,2} = \begin{bmatrix} - 1\\ 1 \\ 0\\[-6pt] \vdots \\ \\ 0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_3,3} = \begin{bmatrix} - 1\\ 0 \\ 1\\ 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\, \boldsymbol{r}_{x_3,N} = \begin{bmatrix} - 1\\ 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0\\0\\0 \end{bmatrix}\!, \nonumber \\&\qquad \boldsymbol{r}_{x_3,N+1} = \begin{bmatrix} 0\\[-6pt] \vdots \\ 1\\ 0\\ 0 \end{bmatrix}\!,\, \boldsymbol{r}_{x_3,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0 \end{bmatrix}\!, \, \boldsymbol{r}_{x_3,N+3} =\begin{bmatrix} \sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2\\[-3pt] \vdots \\ \sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2\\ 0\\0\\a \end{bmatrix}\!, \end{align}
and the left eigenvectors are
\begin{align} {\boldsymbol{l}}_{x_1,1} = \begin{bmatrix} - \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ - \frac {1}{2 \rho a^2}\\[6pt] \frac {1}{2 a}\\[3pt] 0\\0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_1,2} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_1,3} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{3{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\nonumber \\ {\boldsymbol{l}}_{x_1,N} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\ 0\\0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_1,N+1} = \begin{bmatrix} 0\\[-3pt] \vdots \\ 0\\ 1\\ 0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_1,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 0\\ 1 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_1,N+3} =\begin{bmatrix} \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ \frac {1}{2 \rho a^2}\\[6pt] \frac {1}{2 a}\\[3pt] 0\\0 \end{bmatrix}\!,\\ \nonumber \end{align}
\begin{align} {\boldsymbol{l}}_{x_2,1} = \begin{bmatrix} - \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ - \frac {1}{2 \rho a^2}\\[3pt] 0\\ \frac {1}{2 a}\\[3pt] 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_2,2} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_2,3} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{3{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\nonumber \\ {\boldsymbol{l}}_{x_2,N} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\[-3pt] \vdots \\ 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\ 0\\0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_2,N+1} = \begin{bmatrix} 0\\[-3pt] \vdots \\ 1\\ 0\\ 0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_2,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 0\\ 1 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_2,N+3} =\begin{bmatrix} \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ \frac {1}{2 \rho a^2}\\[3pt] 0\\ \frac {1}{2 a}\\[3pt] 0 \end{bmatrix}\!,\\ \nonumber \end{align}
\begin{align} {\boldsymbol{l}}_{x_3,1} = \begin{bmatrix} - \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ - \frac {1}{2 \rho a^2}\\[3pt] 0\\0\\\frac {1}{2 a} \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_3,2} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{2{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{2{\beta }}^2\right)}\\[12pt] 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_3,3} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{3{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{3{\beta }}^2\right)}\\[12pt] 0\\[-6pt] \vdots \\ 0 \end{bmatrix}\!, \ldots ,\nonumber \\ {\boldsymbol{l}}_{x_3,N} = \begin{bmatrix} - \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{N{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\[-6pt] \vdots \\ 0\\[3pt] \frac {\sum _{\beta } \tilde {\rho }_{\beta } a_{1{\beta }}^2}{\sum _{\beta } \tilde {\rho }_{\beta } \left(a_{1{\beta }}^2+a_{N{\beta }}^2\right)}\\[12pt] 0\\ 0\\0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_3,N+1} = \begin{bmatrix} 0\\[-6pt] \vdots \\ 1\\ 0\\ 0 \end{bmatrix}\!,\, {\boldsymbol{l}}_{x_3,N+2} = \begin{bmatrix} 0 \\[-3pt] \vdots \\ 0\\ 1\\ 0 \end{bmatrix}\!, \, {\boldsymbol{l}}_{x_3,N+3} =\begin{bmatrix} \frac {1}{2 \rho a^2}\\[3pt] \vdots \\ \frac {1}{2 \rho a^2}\\[3pt] 0\\0\\\frac {1}{2 a} \end{bmatrix}\!. \end{align}
The eigenvectors are scaled so that
$\boldsymbol{r}_j\boldsymbol{\cdot }{\boldsymbol{l}}_j=1$
for all
$j=1,\ldots ,N+1$
. Both right and left eigenvectors are linearly independent and span
$\mathbb{R}^{N+d}$
. This implies that the system is hyperbolic if the eigenvalues are real valued, which holds when
$a_{{\alpha }{\beta }}^2=\textrm {d}p_{\alpha }/\textrm {d}\tilde {\rho }_{\beta } \gt 0$
. Similar to the one-dimensional case, fields associated with
$\lambda _2,\ldots ,\lambda _{N+{2}}$
correspond to contact discontinuities and those associated with
$\lambda _1,\lambda _{N+{3}}$
correspond to shock or rarefaction waves, respectively.
Open access
