# A unified framework for $N$-phase Navier-Stokes Cahn-Hilliard Allen-Cahn
mixture models with non-matching densities

M.F.P. ten Eikelder e-mail<EMAIL_ADDRESS>
(†Institute for Mechanics, Computational Mechanics Group, Technical University
of Darmstadt )

###### Abstract

Over the past few decades, numerous $N$-phase incompressible diffuse-interface
flow models with non-matching densities have been proposed. Despite aiming to
describe the same physics, these models are generally distinct, and an
overarching modeling framework is absent. This paper provides a unified
framework for $N$-phase incompressible Navier-Stokes Cahn-Hilliard Allen-Cahn
mixture models with a single momentum equation. The framework naturally
emerges from continuum mixture theory, exhibits an energy-dissipative
structure, and is invariant to the choice of fundamental variables. This opens
the door to exploring connections between existing $N$-phase models and
facilitates the computation of $N$-phase flow models rooted in continuum
mixture theory.

Key words. $N$-phase flows, Navier-Stokes Cahn-Hilliard model; phase-field
models, incompressible flow, mixture theory, thermodynamic consistency.

AMS Subject Classification: Primary: 76T30, Secondary: 35Q35, 35R37, 76D05,
76D45, 80A99

##  1 Introduction

###  1.1 Background

Incompressible multi-phase flows are ubiquitous in nature, science and
engineering, with a wide range of application111In this work, the term “phase”
denotes the different fluid materials/constituents (e.g. air and water).. The
development of continuum models (and corresponding methods) that describe
these flows has been an active field of research for the last few decades.
This research can be (roughly) divided into (i) sharp interface models [1, 2,
3, 4], and (ii) diffuse-interface models. Within the diffuse-interface
category, phase-field models constitute a well-known class [5, 6, 7, 8]. While
we acknowledge the importance of each of the above approaches, the current
article focuses on phase-field models.

Phase-field models have gained popularity over the last decades, and have
become a versatile modeling technology with a wide range of applications in
science and engineering. They offer resolutions to challenging moving boundary
problems by simultaneously addressing geometrical representation and the
physical model, see e.g. [7, 9]. By representing interfaces implicitly through
continuous field variables, phase-field models eliminate the need for explicit
boundary tracking, enabling accurate and efficient simulations of phenomena
such as solidification [10], crack propagation in fracture mechanics [11], and
two-fluid flow dynamics [12].

The vast majority of incompressible, viscous, multi-phase flow models in the
literature is restricted to two fluids. In the realm of phase-field modeling,
a prototypical model is the Navier-Stokes Cahn-Hilliard Allen-Cahn (NSCHAC)
model. The first model of this kind, now known as model H, was proposed by
Hohenberg and Halperin in 1977 [13]. This model may be understood as a
simplification of the more complete two-component NSCHAC model in the sense
that (i) it is restricted to matching fluid densities, and (ii) it does not
permit mass transfer between phases (i.e. it does not contain an Allen-Cahn
type term). The foundation of this model is largely based on empirical
arguments; a derivation based on the concept of microforces (see [14]) was
established by Gurtin in 1996 [15]. In subsequent years, several efforts have
been made to relax the matching-density restriction, see e.g. [16, 17, 18],
and see e.g. [19, 20, 21, 22] for numerical simulations. Initially, these
models were classified into two distinct categories: (i) models with a mass-
averaged mixture velocity, and (ii) models with a volume-averaged mixture
velocity. In a recent article, we proposed a unified framework, rooted in
continuum mixture theory, which leads to a single Navier-Stokes Cahn-Hilliard
(NSCH) model that is invariant to the set of fundamental variables [23]; see
[22] for a divergence-conforming discretization with benchmarks. Contrary to
the above-mentioned classification, the framework indicates that
aforementioned classes models coincide, up to minor modifications.

Although most research in the field of multi-phase flows focuses on $N=2$
phases, there are various $N$-phase ($N>2$) incompressible flow models.
Similar to the two-phase case, the literature on $N$-phase models that
(partly) utilize continuum mixture theory is divided into two categories: (i)
models with a mass-averaged mixture velocity, and (ii) models with a volume-
averaged mixture velocity. Without attempting to be complete, we mention the
$N$-phase mass-averaged models [24] ($N=3$) and [25, 26] ($N\geq 2$), and the
$N$-phase volume-averaged models [27, 28, 29] ($N\geq 2$). Furthermore, there
are incompressible $N$-phase NSCH models that are not (partly) based on
continuum mixture theory, rather these models are established via coupling a
multi-phase Cahn-Hilliard (CH) model to the Navier-Stokes equations [30, 31,
32, 33, 34, 35]. We also refer to several theoretical considerations of Allen-
Cahn/Cahn-Hilliard (AC/CH) systems in isolation (ignoring inertial phenomena
present in fluid mechanic systems), see e.g. [36, 37, 38, 39, 40].

###  1.2 Objective and main results

A number of the existing $N$-phase phase-field models, mentioned above, and in
the references therein, provide different models (alongside with computational
methodologies) for the same physical situation: the dynamics of viscous,
incompressible (isothermal) $N$-phase mixture flows. Naturally, there is some
leeway in constitutive modeling, and not all models have the same complexity
level222To organize the various existing models one can adopt the
classification introduced in [41]. This classification is for example utilized
in [42, 43].. However, one can infer that models within the same complexity
class are already distinct before constitutive modeling. The above
observations raise questions regarding differences and connections between the
models. While the aforementioned unified framework of NSCHAC models [23] is
presented for the two-phase case, the adopted modeling principles therein are
at the core not restricted to two phases. There are however a number of non-
trivial considerations that come into play when examining the more general
case $N\geq 2$. Important elements to consider are (i) symmetry properties
with respect to the numbering of the phases, (ii) the reduction-consistency
property (an $N$-phase system reduces to an $(N-M)$-phase system in absence of
$M$ phases), and (iii) and the saturation constraint (volume-
fractions/concentrations add up to one).

The primary objective of this article is to lay down a unified framework of
$N$-phase NSCHAC mixture models. We limit our focus to isothermal phases, and
utilize continuum mixture theory [44] as point of departure. In particular, we
derive the following multi-phase-field model for phases (constituents)
$\alpha=1,...,N$:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\sum_{\beta}\phi_{\beta}\nabla(\mu_{\beta}+\lambda)$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{b}$ $\displaystyle=~{}0,$ (1.1a)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}(\hat{\mathbf{J}}_{\alpha}+\hat{\mathbf{j}}_{\alpha})-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0,$ (1.1b)
$\displaystyle\hat{\mathbf{J}}_{\alpha}-\sum_{\beta}\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle=~{}0,$ (1.1c)
$\displaystyle\hat{\mathbf{j}}_{\alpha}-\sum_{\beta}\mathbf{K}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle=~{}0,$ (1.1d)
$\displaystyle\hat{\zeta}_{\alpha}-\sum_{\beta}m_{\alpha{\beta}}\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle=~{}0,$ (1.1e)

for $\alpha=1,...,N$, where $\mathbf{v}$ is the fluid velocity,
$\rho_{\alpha}$ and $\tilde{\rho}_{\alpha}$ are constituent mass densities,
$\rho=\sum_{\beta}\tilde{\rho}_{\beta}$ the mixture density, $\mathbf{b}$ the
force vector, $\nu$ the dynamical viscosity, $\nu\bar{\lambda}$ the second
viscosity coefficient, $\nabla^{s}\mathbf{v}$ the symmetric velocity gradient,
and $\lambda$ the Lagrange multiplier pressure. Additionally, $\Psi$ is the
constituent free energy, $\mu_{\alpha}$ are constituent chemical potentials,
and $\mathbf{M}_{\alpha{\beta}}$ and $m_{\alpha{\beta}}$ are mobility
parameters. The model is composed of equation (1.1a) that details the mixture
momentum equation, $N$ constituent mass balance equations (1.1b), and models
for peculiar velocities (1.1c), and conservative and non-conservative mass
transfer models (1.1d)-(1.1e). Model (1.1) is expressed in terms of the mass-
averaged mixture velocity $\mathbf{v}$; an alternative – but equivalent –
formulation emerges when adopting the volume-averaged mixture velocity
$\mathbf{u}$:

$\displaystyle\partial_{t}\left(\rho\mathbf{v}\right)+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\sum_{\beta}\phi_{\beta}\nabla(\mu_{\beta}+\lambda)$
$\displaystyle-{\rm div}\left(\nu\left(2\nabla^{s}\mathbf{v}+\bar{\lambda}{\rm
div}\mathbf{v}\mathbf{I}\right)\right)-\rho\mathbf{b}$ $\displaystyle=~{}0,$
(1.2a) $\displaystyle{\rm
div}\mathbf{u}-\sum_{{\beta}}\rho_{\beta}^{-1}(\hat{\zeta}_{\beta}+{\rm
div}\hat{\mathbf{j}}_{\beta})$ $\displaystyle=~{}0,$ (1.2b)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}\left(\phi_{\alpha}\mathbf{v}\right)+\rho_{\alpha}^{-1}{\rm
div}(\hat{\mathbf{J}}_{\alpha}+\hat{\mathbf{j}}_{\alpha})-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0,$ (1.2c)

with
$\mathbf{v}=\mathbf{u}-\displaystyle\sum_{\beta}\rho_{\beta}^{-1}\hat{\mathbf{J}}_{\beta}$,
for $\alpha=1,...,N-1$, where
$\hat{\mathbf{J}}_{\alpha},\hat{\mathbf{j}}_{\alpha}$ and
$\hat{\zeta}_{\alpha}$ are defined in (1.1c), (1.1d) and (1.1e), respectively.
Analogously to the above formulation, the model is comprised of a mixture
momentum equation (1.1a), and $N$ constituent mass balance laws (1.1b). We
provide precise definitions of all quantities in the remainder of the article.

The classification as an NSCHAC model is evident in the combination of a
momentum equation with $N$ mass balance laws. The Cahn-Hilliard components
appear in the third members of the mass balance laws, whereas the Allen-Cahn
character materializes in the latter terms of the mass balance laws.
Furthermore, the model – in both formulations – displays a strong coupling
between the various equations; through the constituent densities
$\tilde{\rho}_{\alpha}$, the velocity $\mathbf{v}$ (or $\mathbf{u}$) and the
Lagrange multiplier pressure $\lambda$.

The secondary objective of this article is reveal connections between model
(1.1)-(1.2) and existing models in the literature. First, we compare model
(1.1)-(1.2) with the unified NSCHAC model [23] for the situation of two
constituents. Subsequently, we compare the framework to the model of Dong
[28]. Finally, we discuss the connections of the proposed framework with the
mixture-theory-compatible $N$-component model [45].

###  1.3 Plan of the paper

The remainder of the paper is organized as follows. In Section 2 we present
the continuum theory of rational mechanics for incompressible isothermal fluid
mixtures, highlighting the connections between different quantities and
formulations of evolution equations. Next, in Section 3, we conduct
constitutive modeling using the Coleman-Noll procedure. Following that,
Section 4 addresses the properties of the model. Subsequently, in Section 5 we
explore the connections of the novel model with existing models in the
literature. Finally, in Section 6, we provide a conclusion and outlook.

##  2 Continuum mixture theory

The purpose of this section is to outline the continuum theory of mixtures for
incompressible constituents, excluding thermal effects. This section aligns
with [45] at several points.

The continuum theory of mixtures is grounded in three general principles
introduced in the pioneering work of Truesdell and Toupin [44]:

  1. 1.

All properties of the mixture must be mathematical consequences of properties
of the constituents.

  2. 2.

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.

  3. 3.

The motion of the mixture is governed by the same equations as is a single
body.

The first principle communicates that the mixture is made up of its
constituent parts. The second principle asserts the connection of the
different components of the physical model through interaction terms. Lastly,
the latter principle states that one can not distinguish the motion of a
mixture from that of a single fluid.

In Section 2.1 we introduce the fundamentals of the continuum theory of
mixtures and the necessary kinematics. Then, in Sections 2.2 and 2.3, we
provide balance laws of individual constituents and associated mixtures.

###  2.1 Preliminaries

In the continuum theory of mixtures the material body is comprised of $N$
constituent bodies $\mathscr{B}_{\alpha}$, with $\alpha=1,\dots,N$. The bodies
$\mathscr{B}_{\alpha}$ are permitted to simultaneously occupy a shared region
in space. Denoting by $\mathbf{X}_{\alpha}$ the spatial position of a particle
of $\mathscr{B}_{\alpha}$ in the Lagrangian (reference) configuration, the
(invertible) deformation map defines the spatial position of a particle:

$\displaystyle\mathbf{x}:=\boldsymbol{\chi}_{\alpha}(\mathbf{X}_{\alpha},t),$
(2.1)

where $\mathbf{x}\in\Omega$, with $\Omega\in\mathbb{R}^{d}$ the domain
(dimension $d$). We introduce the constituent partial mass density
$\tilde{\rho}_{\alpha}$ and specific mass density $\rho_{\alpha}>0$
respectively as:

$\displaystyle\tilde{\rho}_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\displaystyle\lim_{|V|\rightarrow
0}\dfrac{M_{\alpha}(V)}{|V|},$ (2.2a)
$\displaystyle\rho_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\displaystyle\lim_{|V_{\alpha}|\rightarrow
0}\dfrac{M_{\alpha}(V)}{|V_{\alpha}|},$ (2.2b)

where $V\subset\Omega$ (measure $|V|$) is an arbitrary control volume around
$\mathbf{x}$, $V_{\alpha}\subset V$ (measure $|V_{\alpha}|$) is the volume of
constituent $\alpha$ so that $V=\cup_{\alpha}V_{\alpha}$. Here, 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 sum of the
partial mass densities:

$\displaystyle\rho(\mathbf{x},t):=$
$\displaystyle~{}\displaystyle\lim_{|V|\rightarrow
0}\dfrac{M(V)}{|V|}=\displaystyle\sum_{\alpha}\tilde{\rho}_{\alpha}(\mathbf{x},t).$
(2.3)

Additionally, we introduce the mass concentrations and volume fractions
respectively as:

$\displaystyle c_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\displaystyle\lim_{|V|\rightarrow
0}\dfrac{M_{\alpha}(V)}{M(V)}=\dfrac{\tilde{\rho}_{\alpha}}{\rho},$ (2.4a)
$\displaystyle\phi_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\displaystyle\lim_{|V|\rightarrow
0}\dfrac{|V_{\alpha}|}{|V|}=\dfrac{\tilde{\rho}_{\alpha}}{\rho_{\alpha}},$
(2.4b)

which sum up to one:

$\displaystyle\displaystyle\sum_{\alpha}c_{\alpha}(\mathbf{x},t)=$
$\displaystyle~{}1,$ (2.5a)
$\displaystyle\displaystyle\sum_{\alpha}\phi_{\alpha}(\mathbf{x},t)=$
$\displaystyle~{}1.$ (2.5b)

We assume that the constituents are incompressible, meaning that the specific
mass densities are (constituent-wise) constant:

$\displaystyle\rho_{\alpha}(\mathbf{x},t)=\rho_{\alpha}.$ (2.6)

By means of the incompressibility of the constituents, (2.6), and the
definitions (2.4), the volume fractions and concentrations are related by:

$\displaystyle\phi_{\alpha}=$
$\displaystyle~{}\frac{c_{\alpha}}{\rho_{\alpha}}\left(\sum_{\beta}\frac{c_{\beta}}{\rho_{\beta}}\right)^{-1},$
(2.7a) $\displaystyle c_{\alpha}=$
$\displaystyle~{}\rho_{\alpha}\phi_{\alpha}\left(\sum_{\beta}\rho_{\beta}\phi_{\beta}\right)^{-1},$
(2.7b)

for $\alpha=1,...,N$.

######  Remark 2.1 (Incompressibility $N$-phase model).

The assumption that the constituents are incompressible, definition (2.6),
poses challenges for $N$-phase systems. For example, the mappings (2.7) are
not invertible. We discuss these challenges throughout the article, and in
Appendix B.

######  Remark 2.2 (Alternative definitions incompressible mixtures).

Besides the current definition of incompressible constituents (2.6), which is
frequently adopted in the literature (see e.g. [26, 27, 28, 29]), there exist
other notions of incompressibility in mixture flows. We refer for an
alternative to [46] and the references therein.

We proceed with the introduction of the material time derivative
$\grave{\uppsi}_{\alpha}$ of the differentiable constituent function
$\uppsi_{\alpha}$:

$\displaystyle\grave{\uppsi}_{\alpha}=\partial_{t}\uppsi_{\alpha}(\mathbf{X}_{\alpha},t)|_{\mathbf{X}_{\alpha}}.$
(2.8)

Here we adopt the notation $|_{\mathbf{X}_{\alpha}}$ to indicate that
$\mathbf{X}_{\alpha}$ is held fixed. The constituent velocity now follows as
the constituent material derivative of the deformation map:

$\displaystyle\mathbf{v}_{\alpha}(\mathbf{x},t)=\partial_{t}\boldsymbol{\chi}_{\alpha}(\mathbf{X}_{\alpha},t)|_{\mathbf{X}_{\alpha}}=\grave{\boldsymbol{\chi}}_{\alpha}.$
(2.9)

In contrast to the mixture density, there appear various mixture velocities in
the literature. The two most popular are the mass-averaged velocity, denoted
$\mathbf{v}$, and the volume-averaged velocity, denoted $\mathbf{u}$, which
are respectively given by:

$\displaystyle\mathbf{v}(\mathbf{x},t)=$
$\displaystyle~{}\displaystyle\sum_{\alpha}c_{\alpha}(\mathbf{x},t)\mathbf{v}_{\alpha}(\mathbf{x},t),$
(2.10a) $\displaystyle\mathbf{u}(\mathbf{x},t)=$
$\displaystyle~{}\displaystyle\sum_{\alpha}\phi_{\alpha}(\mathbf{x},t)\mathbf{v}_{\alpha}(\mathbf{x},t).$
(2.10b)

We introduce peculiar velocities of the constituents relative to both mixture
velocities:

$\displaystyle\mathbf{w}_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\mathbf{v}_{\alpha}(\mathbf{x},t)-\mathbf{v}(\mathbf{x},t),$
(2.11a) $\displaystyle\boldsymbol{\omega}_{\alpha}(\mathbf{x},t):=$
$\displaystyle~{}\mathbf{v}_{\alpha}(\mathbf{x},t)-\mathbf{u}(\mathbf{x},t).$
(2.11b)

Additionally, we define the following (scaled) peculiar velocities (that
depend on $\mathbf{x}$ and $t$):

$\displaystyle\mathbf{J}_{\alpha}:=$
$\displaystyle~{}\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha},$ (2.12a)
$\displaystyle\mathbf{h}_{\alpha}=$
$\displaystyle~{}\phi_{\alpha}\mathbf{w}_{\alpha},$ (2.12b)
$\displaystyle\mathbf{J}_{\alpha}^{u}=$
$\displaystyle~{}\tilde{\rho}_{\alpha}\boldsymbol{\omega}_{\alpha}$ (2.12c)
$\displaystyle\mathbf{h}_{\alpha}^{u}=$
$\displaystyle~{}\phi_{\alpha}\boldsymbol{\omega}_{\alpha}.$ (2.12d)

######  Remark 2.3 (Terminology peculiar velocities).

The quantities (2.11) and (2.12) are in the literature often referred to as
“diffusion velocities” and “diffusive fluxes”, respectively. This terminology
is natural because the terms (2.12) appear in constituent mass balance laws
(see Section 2.2) as flux terms, and their constitutive models (see Section
3.4) have a diffusive character. However, utilizing constitutive models for
(2.12) is not essential (see [45]), and therefore we use the terminology
“(scaled) peculiar velocity” to reflect their original definitions (2.11) and
(2.12).

Direct consequences of (2.11), (2.12a), and (2.12d) are the properties:

$\displaystyle\displaystyle\sum_{\alpha}\mathbf{J}_{\alpha}=$
$\displaystyle~{}0,$ (2.13a)
$\displaystyle\displaystyle\sum_{\alpha}\mathbf{h}_{\alpha}^{u}=$
$\displaystyle~{}0.$ (2.13b)

The relation between the mass-averaged and volume-averaged velocities is
specified in the following lemma.

######  Lemma 2.4 (Relation mass-averaged and volume-averaged velocities).

The mass-averaged and volume-averaged velocity variables are related via:

$\displaystyle\mathbf{u}=$
$\displaystyle~{}\mathbf{v}+\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}=\mathbf{v}+\sum_{\alpha}\mathbf{h}_{\alpha},$
(2.14a) $\displaystyle\mathbf{v}=$
$\displaystyle~{}\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}.$
(2.14b)

###### Proof.

These relations result from the sequences of identities:

$\displaystyle\mathbf{u}=\sum_{\alpha}\phi_{\alpha}\mathbf{v}_{\alpha}=\sum_{\alpha}\phi_{\alpha}\mathbf{w}_{\alpha}+\sum_{\alpha}\phi_{\alpha}\mathbf{v}=\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}+\mathbf{v},$
(2.15a) $\displaystyle
0=\sum_{\alpha}\mathbf{J}_{\alpha}=\sum_{\alpha}\tilde{\rho}_{\alpha}(\mathbf{v}_{\alpha}-\mathbf{v})=\sum_{\alpha}\tilde{\rho}_{\alpha}(\mathbf{v}_{\alpha}-\mathbf{u}+\mathbf{u}-\mathbf{v})$
$\displaystyle=\sum_{\alpha}\mathbf{J}_{\alpha}^{u}+\rho(\mathbf{u}-\mathbf{v}).$
(2.15b)

∎

The relation between the scaled peculiar velocities is displayed in the next
lemma.

######  Lemma 2.5 (Relation scaled peculiar velocities).

The scaled peculiar velocities are related via:

$\displaystyle\mathbf{J}_{\alpha}=$
$\displaystyle~{}\mathbf{J}_{\alpha}^{u}-c_{\alpha}\sum_{\beta}\mathbf{J}_{\beta}^{u},$
(2.16a) $\displaystyle\mathbf{J}_{\alpha}^{u}=$
$\displaystyle~{}\mathbf{J}_{\alpha}-\tilde{\rho}_{\alpha}\sum_{\beta}\rho_{\beta}^{-1}\mathbf{J}_{\beta},$
(2.16b) $\displaystyle\mathbf{h}_{\alpha}=$
$\displaystyle~{}\mathbf{h}_{\alpha}^{u}-\phi_{\alpha}\rho^{-1}\sum_{\beta}\rho_{\beta}\mathbf{h}_{\beta}^{u},$
(2.16c) $\displaystyle\mathbf{h}_{\alpha}^{u}=$
$\displaystyle~{}\mathbf{h}_{\alpha}-\phi_{\alpha}\sum_{\beta}\mathbf{h}_{\beta}.$
(2.16d)

###### Proof.

These identities are a direct consequence of Lemma 2.5. ∎

Lastly, we define the material derivative of the mixture relative to the mass-
averaged velocity:

$\displaystyle\dot{\uppsi}(\mathbf{x},t)=$
$\displaystyle~{}\partial_{t}\uppsi(\mathbf{x},t)+\mathbf{v}(\mathbf{x},t)\cdot\nabla\uppsi(\mathbf{x},t).$
(2.17)

###  2.2 Constituent balance laws

In the continuum theory of mixtures, each constituent moves according to a
distinct set of balance laws, as specified by the second general principle.
These laws incorporate terms that model the interactions among the different
constituents. The following local balance laws apply to the motion of each
constituent $\alpha=1,\dots,N$ for all $\mathbf{x}\in\Omega$ and $t>0$:

$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha})$
$\displaystyle=~{}\gamma_{\alpha},$ (2.18a)
$\displaystyle\partial_{t}(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha})+{\rm
div}\left(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha}\otimes\mathbf{v}_{\alpha}\right)-{\rm
div}\mathbf{T}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{b}_{\alpha}$
$\displaystyle=~{}\boldsymbol{\pi}_{\alpha},$ (2.18b)
$\displaystyle\mathbf{T}_{\alpha}-\mathbf{T}_{\alpha}^{T}$
$\displaystyle=~{}\mathbf{N}_{\alpha},$ (2.18c)

Equation (2.18a) describes the local constituent mass balance law, where the
interaction term $\gamma_{\alpha}$ denotes the mass supply of constituent
$\alpha$ due to chemical reactions with the other constituents. Then, (2.18b)
represents the local constituent linear momentum balance law, where
$\mathbf{T}_{\alpha}$ is the Cauchy stress tensor of constituent $\alpha$,
$\mathbf{b}_{\alpha}$ is the constituent external body force, and
$\boldsymbol{\pi}_{\alpha}$ is the momentum exchange rate of constituent
$\alpha$ with the other constituents. We assume equal body forces
($\mathbf{b}_{\alpha}=\mathbf{b}$ for $\alpha=1,\dots,N$) throughout the
article. Additionally, we restrict to gravitational body forces:
$\mathbf{b}=-b\boldsymbol{\jmath}=-b\nabla y$, with $y$ being the vertical
coordinate, $\boldsymbol{\jmath}$ the vertical unit vector, and $b$ a
constant. Finally, (2.18c) describes the local constituent angular momentum
balance with $\mathbf{N}_{\alpha}$ the intrinsic moment of momentum.

We introduce a split of the mass transfer term into a conservative part and a
potentially non-conservative contribution via:

$\displaystyle\gamma_{\alpha}=\zeta_{\alpha}-{\rm div}\mathbf{j}_{\alpha}.$
(2.19)

The mass balance laws (3.1a) take the form:

$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha})+{\rm
div}\mathbf{j}_{\alpha}=\zeta_{\alpha}.$ (2.20)

By invoking the definitions in Section 2.1, one can deduce various alternative
– equivalent – formulations of the constituent mass balance laws (2.18a), such
as:

$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{v})+{\rm
div}(\mathbf{J}_{\alpha}+\mathbf{j}_{\alpha})$
$\displaystyle=~{}\zeta_{\alpha},$ (2.21a)
$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{u})+{\rm
div}(\mathbf{J}_{\alpha}^{u}+\mathbf{j}_{\alpha})$
$\displaystyle=~{}\zeta_{\alpha},$ (2.21b)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+{\rm
div}\mathbf{h}_{\alpha}+\rho_{\alpha}^{-1}{\rm div}\mathbf{j}_{\alpha}$
$\displaystyle=~{}\rho_{\alpha}^{-1}\zeta_{\alpha},$ (2.21c)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}\left(\phi_{\alpha}\mathbf{u}\right)+{\rm
div}\mathbf{h}_{\alpha}^{u}+\rho_{\alpha}^{-1}{\rm div}\mathbf{j}_{\alpha}$
$\displaystyle=~{}\rho_{\alpha}^{-1}\zeta_{\alpha},$ (2.21d)
$\displaystyle\rho\partial_{t}c_{\alpha}+\rho\mathbf{v}\cdot\nabla
c_{\alpha}+{\rm div}(\mathbf{J}_{\alpha}+\mathbf{j}_{\alpha})$
$\displaystyle=~{}\zeta_{\alpha},$ (2.21e)
$\displaystyle\rho\partial_{t}c_{\alpha}+\rho\mathbf{u}\cdot\nabla
c_{\alpha}+{\rm
div}(\mathbf{J}_{\alpha}^{u}+\mathbf{j}_{\alpha})-c_{\alpha}{\rm
div}\left(\sum_{\beta}\mathbf{J}_{\beta}^{u}\right)$
$\displaystyle=~{}\zeta_{\alpha}.$ (2.21f)

Additionally, by invoking the relation (2.7) we can deduce numerous
alternative – equivalent – formulations; for example, by inserting (2.7) into
(2.21c) we arrive at an uncommon formulation:

$\displaystyle\partial_{t}\left(\frac{c_{\alpha}}{\rho_{\alpha}}\left(\sum_{\beta}\frac{c_{\beta}}{\rho_{\beta}}\right)^{-1}\right)+{\rm
div}\left(\frac{c_{\alpha}}{\rho_{\alpha}}\left(\sum_{\beta}\frac{c_{\beta}}{\rho_{\beta}}\right)^{-1}\mathbf{v}\right)$
$\displaystyle+{\rm div}\mathbf{h}_{\alpha}+\rho_{\alpha}^{-1}{\rm
div}\mathbf{j}_{\alpha}$ $\displaystyle=\rho_{\alpha}^{-1}\zeta_{\alpha}.$
(2.22)

Similarly, one can write the constituent momentum balance laws (2.18b):

$\displaystyle\partial_{t}(\tilde{\rho}_{\alpha}\mathbf{v}+\mathbf{J}_{\alpha})+{\rm
div}\left(\tilde{\rho}_{\alpha}\mathbf{v}\otimes\mathbf{v}+\mathbf{J}_{\alpha}\otimes\mathbf{v}+\mathbf{v}\otimes\mathbf{J}_{\alpha}\right)$
$\displaystyle-~{}{\rm
div}\left(\mathbf{T}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)-\tilde{\rho}_{\alpha}\mathbf{b}_{\alpha}$
$\displaystyle=~{}\boldsymbol{\pi}_{\alpha},$ (2.23a)
$\displaystyle\partial_{t}(\tilde{\rho}_{\alpha}\mathbf{u}+\mathbf{J}_{\alpha}^{u})+{\rm
div}\left(\tilde{\rho}_{\alpha}\mathbf{u}\otimes\mathbf{u}+\mathbf{J}_{\alpha}^{u}\otimes\mathbf{u}+\mathbf{u}\otimes\mathbf{J}_{\alpha}^{u}\right)$
$\displaystyle+~{}{\rm
div}\left(\tilde{\rho}_{\alpha}\boldsymbol{\omega}_{\alpha}\otimes\boldsymbol{\omega}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)$
$\displaystyle-{\rm
div}\left(\mathbf{T}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)-\tilde{\rho}_{\alpha}\mathbf{b}_{\alpha}$
$\displaystyle=~{}\boldsymbol{\pi}_{\alpha}.$ (2.23b)

Finally, we introduce the constituent kinetic and gravitational energies,
respectively, as:

$\displaystyle\mathscr{K}_{\alpha}=$
$\displaystyle~{}\tilde{\rho}_{\alpha}\|\mathbf{v}_{\alpha}\|^{2}/2,$ (2.24a)
$\displaystyle\mathscr{G}_{\alpha}=$
$\displaystyle~{}\tilde{\rho}_{\alpha}by,$ (2.24b)

where
$\|\mathbf{v}_{\alpha}\|=(\mathbf{v}_{\alpha}\cdot\mathbf{v}_{\alpha})^{1/2}$
is the Euclidean norm of the velocity $\mathbf{v}_{\alpha}$.

###  2.3 Mixture balance laws

The standard formulation of mixture balance laws is well-known and follows
from summing the balance laws (2.18) over all constituents. To establish the
precise form, one can, for example, utilize the formulations (2.21a) and
(2.23a) and invoke the identity (2.13a) to obtain:

$\displaystyle\partial_{t}\rho+{\rm div}(\rho\mathbf{v})$
$\displaystyle=~{}0,$ (2.25a) $\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}$ $\displaystyle=~{}0,$ (2.25b)
$\displaystyle\mathbf{T}-\mathbf{T}^{T}$ $\displaystyle=~{}0,$ (2.25c)

where the mixture stress and mixture body force are given by, respectively:

$\displaystyle\mathbf{T}=$
$\displaystyle~{}\sum_{\alpha}\mathbf{T}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha},$
(2.26a) $\displaystyle\mathbf{b}=$
$\displaystyle~{}\frac{1}{\rho}\sum_{\alpha}\tilde{\rho}_{\alpha}\mathbf{b}_{\alpha},$
(2.26b)

and where we have postulated the following balance conditions to hold as
follows:

$\displaystyle\displaystyle\sum_{\alpha}\gamma_{\alpha}=$ $\displaystyle~{}0,$
(2.27a) $\displaystyle\displaystyle\sum_{\alpha}\boldsymbol{\pi}_{\alpha}=$
$\displaystyle~{}0,$ (2.27b)
$\displaystyle\displaystyle\sum_{\alpha}\mathbf{N}_{\alpha}=$
$\displaystyle~{}0,$ (2.27c)

and where we invoke (2.27a) via:

$\displaystyle\sum_{\alpha}\zeta_{\alpha}=$ $\displaystyle~{}0,$ (2.28a)
$\displaystyle\sum_{\alpha}\mathbf{j}_{\alpha}=$ $\displaystyle~{}0.$ (2.28b)

This formulation is compatible with the first general principle: the motion of
the mixture is derived from the motion of its individual constituents. In
addition, the postulate (2.27) is essential to ensure general principle three.
Even though the forms presented in (2.21) and (2.23) are equivalent, the
summation of these laws over the constituents does not provide a suitable
system of mixture balance laws for each of the formulations. Namely general
principle three communicates that the resulting equations of the mixture are
indistinguishable from that of a single body. Complying with this principle
restricts the forms of the mass balance law to (2.21a) and (2.21b), and
requires the identification of suitable mixture variables. These variables are
$\rho$, $\mathbf{v}$, $\mathbf{T}$ and $\mathbf{b}$, as defined above. In this
sense, the framework of continuum mixture theory serves as a guideline for
defining mixture variables. However, one can work with other variables as
well; and this is fully compatible with the framework.

We discuss other formulations that emerge from (2.21) and (2.23). Summation of
(2.21b)-(2.21f) over the constituents provides:

$\displaystyle\partial_{t}\rho+{\rm
div}\left(\rho\left(\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)$
$\displaystyle=~{}\sum_{\alpha}\gamma_{\alpha}=0,$ (2.29a) $\displaystyle{\rm
div}\left(\mathbf{v}+\sum_{\alpha}\mathbf{h}_{\alpha}\right)$
$\displaystyle=~{}\sum_{\alpha}\rho_{\alpha}^{-1}\gamma_{\alpha},$ (2.29b)
$\displaystyle{\rm div}\mathbf{u}$
$\displaystyle=~{}\sum_{\alpha}\rho_{\alpha}^{-1}\gamma_{\alpha}$ (2.29c)
$\displaystyle 0$ $\displaystyle=~{}\sum_{\alpha}\gamma_{\alpha},$ (2.29d)

where (2.29a)-(2.29c) follow from (2.21b)-(2.21d), respectively, and (2.29d)
results from both (2.21e) and (2.21f). We observe from (2.29a) that the term
in the inner brackets in the second term represents the mixture velocity.
Obviously, this matches the mass averaged velocity by invoking Lemma 2.4.
Next, note that (2.29b) also follows from the summation over the constituents
of (2.2). With the aid of Lemma 2.4, one can infer that (2.29b) and (2.29c)
are identical. Furthermore,
$\mathbf{v}+\sum_{\alpha}\mathbf{h}_{\alpha}=\mathbf{u}$ is a divergence-free
velocity whenever either (i) mass transfer is absent ($\gamma_{\alpha}=0$ for
all $\alpha$), or (ii) the constituent densities match
($\rho_{\alpha}=\rho_{\beta}$ for all $\alpha,{\beta}$). The equation (2.29d)
complies with the balance condition (2.27a). Finally, the summation of (2.23)
yields:

$\displaystyle\partial_{t}\left(\rho\left(\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)+{\rm
div}\left(\rho\left(\mathbf{u}\otimes\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\mathbf{u}+\rho^{-1}\mathbf{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)$
$\displaystyle+~{}{\rm
div}\left(\sum_{\alpha}\tilde{\rho}_{\alpha}\boldsymbol{\omega}_{\alpha}\otimes\boldsymbol{\omega}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}=~{}0.\hskip 56.9055pt$ (2.30)

Invoking Lemma 2.4 and Lemma 2.5 this may be written as:

$\displaystyle\partial_{t}\left(\rho\left(\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)+{\rm
div}\left(\rho\left(\mathbf{u}\otimes\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\mathbf{u}+\rho^{-1}\mathbf{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)$
$\displaystyle+~{}{\rm
div}\left(\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}=~{}0.\hskip 56.9055pt$ (2.31)

One can infer equivalence with the mass-averaged momentum equation by noting
the identities:

$\displaystyle{\rm
div}\left(\rho\left(\mathbf{u}\otimes\mathbf{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\mathbf{u}+\rho^{-1}\mathbf{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)\right)=\hskip
56.9055pt$ $\displaystyle{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}-\rho(\mathbf{v}-\mathbf{u})\otimes(\mathbf{v}-\mathbf{u})\right),\hskip
56.9055pt$ (2.32a) $\displaystyle{\rm
div}\left(\sum_{\alpha}\tilde{\rho}_{\alpha}\boldsymbol{\omega}_{\alpha}\otimes\boldsymbol{\omega}_{\alpha}-\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)={\rm
div}\left(\rho(\mathbf{v}-\mathbf{u})\otimes(\mathbf{v}-\mathbf{u})\right).\hskip
56.9055pt$ (2.32b)

In summary, an – equivalent – formulation of mixture balance laws (2.25) in
terms of the volume-averaged velocity is:

$\displaystyle\partial_{t}\rho+{\rm
div}\left(\rho\mathbf{u}+\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)$
$\displaystyle=~{}0,$ (2.33a)
$\displaystyle\partial_{t}\left(\rho\mathbf{u}+\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)+{\rm
div}\left(\rho\mathbf{u}\otimes\mathbf{u}+\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\mathbf{u}\right.$
$\displaystyle+\left.\mathbf{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}+\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\otimes\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}$ $\displaystyle=~{}0,$ (2.33b)
$\displaystyle\mathbf{T}-\mathbf{T}^{T}$ $\displaystyle=~{}0.$ (2.33c)

The various forms presented in this section show that the set of balance laws,
on both constituent level (Section 2.2) and mixture level (Section 2.3), is
invariant to the set of fundamental variables.

We close this section with a remark on the kinetic and gravitational energies.
According to the first metaphysical principle of mixture theory, the kinetic
and gravitational energies of the mixture equal the summation of the
constituent energies:

$\displaystyle\mathscr{K}=$
$\displaystyle~{}\displaystyle\sum_{\alpha}\mathscr{K}_{\alpha},$ (2.34a)
$\displaystyle\mathscr{G}=$
$\displaystyle~{}\displaystyle\sum_{\alpha}\mathscr{G}_{\alpha}.$ (2.34b)

The kinetic energy of the mixture can be decomposed as:

$\displaystyle\mathscr{K}=$
$\displaystyle~{}\bar{\mathscr{K}}+\displaystyle\sum_{\alpha}\frac{1}{2}\tilde{\rho}_{\alpha}\|\mathbf{w}_{\alpha}\|^{2},$
(2.35a) $\displaystyle\bar{\mathscr{K}}=$
$\displaystyle~{}\frac{1}{2}\rho\|\mathbf{v}\|^{2},$ (2.35b)

where $\bar{\mathscr{K}}$ represents the kinetic energy of the mixture
variables, and where the other term is the kinetic energy of the constituents
utilizing the peculiar velocity. The second terms may also be expressed in
terms of volume-averaged quantities:

$\displaystyle\displaystyle\sum_{\alpha}\frac{1}{2}\tilde{\rho}_{\alpha}\|\mathbf{w}_{\alpha}\|^{2}=\displaystyle\sum_{\alpha}\frac{1}{2}\tilde{\rho}_{\alpha}\|\boldsymbol{\omega}_{\alpha}-\rho^{-1}\sum_{\alpha}\mathbf{J}_{\alpha}^{u}\|^{2}.$
(2.36)

##  3 Constitutive modeling

This section details the development of constitutive models under the
constraints of an energy-dissipative postulate. First, Section 3.1 outlines
the fundamental assumptions and modeling choices. Next, Section 3.2
establishes the constitutive modeling restriction introduced in Section 3.1.
Finally, in Section 3.4, selects particular constitutive models that adhere to
these established restrictions.

###  3.1 Assumptions and modeling choices

Rather than using the complete set of balance laws as given in (2.18a),
(2.18b), and (2.18c), we limit our focus to the simplified subset:

$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+{\rm div}\mathbf{h}_{\alpha}$
$\displaystyle=~{}\rho_{\alpha}^{-1}\gamma_{\alpha},$ (3.1a)
$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}$ $\displaystyle=~{}0,$ (3.1b)
$\displaystyle\mathbf{T}-\mathbf{T}^{T}$ $\displaystyle=~{}0,$ (3.1c)

where (3.1a) holds for constituents $\alpha=1,...,N$. At this point, the
system is comprised of the unknown quantities: volume fractions
$\phi_{\alpha}$ ($\alpha=1,...,N$), where we recall the identity (2.4b), mass-
averaged mixture velocity $\mathbf{v}$, peculiar velocities
$\mathbf{J}_{\alpha}$ ($\alpha=1,...,N$), mass transfer terms
$\gamma_{\alpha}$ ($\alpha=1,...,N$), and mixture stress $\mathbf{T}$. In
order to close the system we seek for constitutive models for the peculiar
velocities $\mathbf{J}_{\alpha}$ ($\alpha=1,...,N$), the constituent mass
transfer terms $\gamma_{\alpha}$ ($\alpha=1,...,N$) and the mixture stress
$\mathbf{T}$. Seeking for constitutive models for the peculiar velocities
$\mathbf{J}_{\alpha}$ could be perceived as a simplification procedure.
Namely, substituting a constitutive model (in Section 3.4), in general,
violates the continuum mixture theory definitions (2.12). We discard these
definitions (2.12). Thus, instead of working with $N$ velocities quantities
$\mathbf{v}_{\alpha}$, the simplified system contains a single unknown
velocity quantity $\mathbf{v}$ and constitutive models for peculiar velocities
$\mathbf{J}_{\alpha}$. This is compatible with the structure of the system:
the full system is composed of $N$ linear momentum (mixture) balance law
whereas the simplified system contains a single linear momentum balance law.
Additionally, we enforce the balance condition for the peculiar velocities
(2.13a) and the mass transfer terms (2.27a) as follows:

$\displaystyle\sum_{\alpha}\mathbf{J}_{\alpha}=$ $\displaystyle~{}0,$ (3.2a)
$\displaystyle\sum_{\alpha}\gamma_{\alpha}=$ $\displaystyle~{}0.$ (3.2b)

The system (3.1) contains the unknown variables $\mathbf{v}$ and
$\phi_{\alpha}$ ($\alpha=1,...,N$). We emphasize that the set
$\left\\{\phi_{\alpha}\right\\}_{\alpha=1,...,N}$ is comprised of $N-1$
independent variables due to the summation condition (2.5b). As such, system
(3.1) has a degenerate nature; it contains $N+1$ equations for $N$ variables
(we preclude (2.25c) in this count). Section 3.2 restores the balance by means
of a Lagrange multiplier construction.

######  Remark 3.1 (Classification).

The above assumptions lead to a model that includes $N$ constituent mass
balance laws along with a single momentum balance law. According to the
classification by [41], this configuration aligns best with a class-I model.

We adopt the well-known Coleman-Noll procedure [47] as a guiding principle to
design constitutive models. For this purpose we postulate the energy-
dissipation law:

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=\mathscr{W}-\mathscr{D},$
(3.3)

satisfying $\mathscr{D}\geq 0$. The total energy is comprised of the Helmholtz
free energy, the kinetic energy and the gravitational energy:

$\displaystyle\mathscr{E}=\displaystyle\int_{\mathcal{R}(t)}(\Psi+\bar{\mathscr{K}}+\mathscr{G})~{}{\rm
d}v.$ (3.4)

In this context, $\mathcal{R}=\mathcal{R}(t)\subset\Omega$ refers to a time-
dependent control volume with volume element ${\rm d}v$ and a unit outward
normal $\boldsymbol{\nu}$ that is transported by the velocity field
$\mathbf{v}$. Additionally, $\mathscr{W}$ represents a work rate term on the
boundary $\partial\mathcal{R}(t)$ (with boundary element ${\rm d}a$), and
$\mathscr{D}$ denotes the dissipation within the interior of $\mathcal{R}(t)$.

######  Remark 3.2 (Energy-dissipation postulate).

As mentioned in [23, 45], the energy-dissipation statement (3.3) does not
match the second law of thermodynamics for mixtures; it could be perceived as
approximation of it.

We postulate that the free energy to pertain to the constitutive class:

$\displaystyle\Psi=$
$\displaystyle~{}\hat{\Psi}\left(\left\\{\phi_{\alpha}\right\\}_{\alpha=1,...,N},\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha=1,...,N}\right).$
(3.5)

The volume fractions $\left\\{\phi_{\alpha}\right\\}_{\alpha=1,...,N}$ and
their gradients $\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha=1,...,N}$ are
not independent quantities due to the saturation constraint (2.5b). Thus, we
consider (3.5) only subject to the constraint (2.5b). Additionally, we
introduce the chemical potential quantities ($\alpha=1,...,N$):

$\displaystyle\hat{\mu}_{\alpha}=$
$\displaystyle~{}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}-{\rm
div}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}.$ (3.6)

Observe that the chemical potentials are individually arbitrary; for example,
addition of the term $(1-\sum_{\alpha}\phi_{\alpha})$ to $\hat{\Psi}$ does not
alter it, but it modifies the chemical potentials $\mu_{\alpha}$. However,
once a free energy within the class (3.5) has been specified, the associated
chemical potential is well-defined.

######  Remark 3.3 (Reduced free energy class).

Instead of utilizing the class (3.5), one can also directly enforce the
summation constraint (2.5b) to arrive at a class with reduced dependency. In
general, this breaks the symmetry of the approach, and therefore we do not
adopt this alternative here. We discuss this option in Appendix A.

######  Remark 3.4 (Concentration-dependent free energy class).

One can also work with a constituent class that depends on concentration
quantities. We discuss this option in Section 3.3.

###  3.2 Modeling restriction

Moving forward, we study in detail the restriction (3.3). First, we analyze
the evolution of the energy (3.4). Through the application of the Reynolds
transport theorem to the free energy $\hat{\Psi}$, we have:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm
d}v=\displaystyle\int_{\mathcal{R}(t)}\partial_{t}\hat{\Psi}~{}{\rm
d}v+\displaystyle\int_{\partial\mathcal{R}(t)}\hat{\Psi}\mathbf{v}\cdot\boldsymbol{\nu}~{}{\rm
d}a.$ (3.7)

We notice that the summation constraint (2.5b) does not alter the derivative
of the free energy class (3.5).

######  Lemma 3.5 (Derivative of the free energy).

The derivative of the free energy class (3.5) is given by:

$\displaystyle{\rm
d}\hat{\Psi}=\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha}),$ (3.8)

where ${\rm d}$ is the derivative operator.

###### Proof.

See Lemma A.2. ∎

Invoking Lemma 3.5 and the divergence theorem yields:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm
d}v=\displaystyle\int_{\mathcal{R}(t)}$ $\displaystyle~{}\hat{\Psi}~{}{\rm
div}\mathbf{v}+\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}\dot{\phi}_{\alpha}+\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\cdot\left(\nabla\phi_{\alpha}\right)\dot{}~{}{\rm
d}v.$ (3.9)

Integrating by parts provides:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm
div}\mathbf{v}+\sum_{\alpha}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha}-\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}:\nabla\mathbf{v}~{}{\rm
d}v$
$\displaystyle~{}+\displaystyle\int_{\partial\mathcal{R}(t)}\sum_{\alpha}\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\cdot\boldsymbol{\nu}~{}{\rm
d}a,$ (3.10)

where we have substituted the identity

$\displaystyle(\nabla\uppsi)\dot{}=\nabla(\dot{\uppsi})-(\nabla\uppsi)^{T}\nabla\mathbf{v}$
(3.11)

for $\uppsi=\phi_{\alpha}$. Examination of the various free energy terms shows
that these terms are well-defined.

######  Lemma 3.6 (Well-defined free energy terms).

The following free energy terms in (3.2) are well-defined:

$\displaystyle\sum_{\alpha}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha};\quad\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}};\quad\sum_{\alpha}\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}};$
(3.12)

###### Proof.

See Lemma A.3. ∎

Substituting the constituent mass balance laws (3.1a) provides:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm
div}\mathbf{v}+\sum_{\alpha}\hat{\mu}_{\alpha}\left(-\phi_{\alpha}{\rm
div}\mathbf{v}-\rho_{\alpha}^{-1}{\rm
div}\mathbf{H}_{\alpha}+\rho_{\alpha}^{-1}\zeta_{\alpha}\right)$
$\displaystyle-\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}:\nabla\mathbf{v}~{}{\rm
d}v+\displaystyle\int_{\partial\mathcal{R}(t)}\sum_{\alpha}\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\cdot\boldsymbol{\nu}~{}{\rm
d}a.$ (3.13)

By again applying integration by parts one can infer that:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\sum_{\alpha}\hat{\Psi}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\hat{\Psi}~{}{\rm
div}\mathbf{v}-\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}{\rm
div}\mathbf{v}+\sum_{\alpha}\nabla(\rho_{\alpha}^{-1}\hat{\mu}_{\alpha})\cdot\mathbf{H}_{\alpha}$
$\displaystyle~{}~{}~{}-\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}:\nabla\mathbf{v}+\sum_{\alpha}\rho_{\alpha}^{-1}\hat{\mu}_{\alpha}\zeta_{\alpha}~{}{\rm
d}v$
$\displaystyle~{}+\displaystyle\int_{\partial\mathcal{R}(t)}\left(\sum_{\alpha}\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\rho_{\alpha}^{-1}\hat{\mu}_{\alpha}\mathbf{H}_{\alpha}\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a.$ (3.14)

At this point additional constraints are demanded in order for the chemical
potential terms to be well-defined. Next, the evolution of the kinetic and
gravitational energies take the form (see [23] for details):

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\mathscr{K}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}-\nabla\mathbf{v}:\mathbf{T}+\rho\mathbf{v}\cdot\mathbf{g}~{}{\rm
d}v+\displaystyle\int_{\partial\mathcal{R}(t)}\mathbf{v}\cdot\mathbf{T}\boldsymbol{\nu}~{}{\rm
d}a$ (3.15a) $\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\mathscr{G}~{}{\rm d}v=$
$\displaystyle~{}-\displaystyle\int_{\mathcal{R}(t)}\rho\mathbf{v}\cdot\mathbf{g}~{}{\rm
d}v.$ (3.15b)

The superposition of (3.2) and (3.15) provides the evolution of the total
energy:

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(\rho_{\alpha}^{-1}\hat{\mu}_{\alpha}\mathbf{H}_{\alpha}-\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a$
$\displaystyle~{}-\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla(\rho_{\alpha}^{-1}\hat{\mu}_{\alpha})\cdot\mathbf{H}_{\alpha}-\rho_{\alpha}^{-1}\hat{\mu}_{\alpha}\zeta_{\alpha}\right)~{}{\rm
d}v.$ (3.16)

As aforementioned in Section 3.1, the system of balance laws (3.1) subject the
balance conditions (3.2) is degenerate. Namely, the terms $\nabla\mathbf{v}$,
$\mathbf{H}_{\alpha}$ and $\zeta_{\alpha}$ are connected via (2.29b). This
manifests itself in the energy dissipation statement (3.2). The degeneracy
needs to be eliminated in order to exploit the energy-dissipation condition as
a guiding principle for constitutive modeling. To this purpose we enforce
(2.29b) with the Lagrange multiplier construction:

$\displaystyle 0=$ $\displaystyle~{}\lambda\left({\rm
div}\mathbf{v}+\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\nabla\cdot\mathbf{H}_{\alpha}-\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\zeta_{\alpha}\right),$
(3.17)

where $\lambda$ is the scalar Lagrange multiplier.

######  Remark 3.7 (Lagrange multiplier constraint).

Recalling Lemma 2.4, we observe that the Lagrange multiplier $\lambda$
enforces the constraint (2.29c). As such, in absence of mass transfer
($\gamma_{\alpha}=0$, $\alpha=1,...,N$), it constrains ${\rm
div}\mathbf{u}=0$.

Integrating (3.17) over $\mathcal{R}(t)$ and subtracting the result from (3.2)
provides:

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(g_{\alpha}\mathbf{H}_{\alpha}-\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a$
$\displaystyle~{}-\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\lambda\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla
g_{\alpha}\cdot\mathbf{H}_{\alpha}-g_{\alpha}\zeta_{\alpha}\right)~{}{\rm
d}v,$ (3.18)

where we have utilized Gauß divergence theorem, and where we have defined the
(generalized) chemical potential quantities:

$\displaystyle g_{\alpha}:=$
$\displaystyle~{}\rho_{\alpha}^{-1}\hat{\mu}_{\alpha,\lambda},$ (3.19a)
$\displaystyle\hat{\mu}_{\alpha,\lambda}:=$
$\displaystyle~{}\hat{\mu}_{\alpha}+\lambda.$ (3.19b)

We identify the rate of work and the dissipation respectively as:

$\displaystyle\mathscr{W}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(g_{\alpha}\mathbf{H}_{\alpha}-\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a,$ (3.20a) $\displaystyle\mathscr{D}=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\lambda\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla
g_{\alpha}\cdot\mathbf{H}_{\alpha}-g_{\alpha}\zeta_{\alpha}\right)~{}{\rm
d}v.$ (3.20b)

Given the arbitrary nature of the control volume $\mathcal{R}=\mathcal{R}(t)$,
the fulfillment of the energy dissipation law is contingent upon satisfying
the local inequality:

$\displaystyle\left(\mathbf{T}+\lambda\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle-\sum_{\alpha}\nabla
g_{\alpha}\cdot\mathbf{H}_{\alpha}-\sum_{\alpha}g_{\alpha}\zeta_{\alpha}$
$\displaystyle\geq 0.$ (3.21)

######  Remark 3.8 (Compatibility with continuum mixture theory).

This section has demonstrated that the energy dissipation postulate (3.3) is
fulfilled when the local inequality (3.2) is satisfied. As mentioned in Remark
3.2, the energy dissipation postulate is an approximation of the second law of
mixture theory. However, we emphasize that the presented derivations are fully
compatible with continuum mixture theory.

We finalize this section with a remark on the connection between the chemical
potentials and the Lagrange multiplier. Recalling the saturation constraint
(2.5b), we recognize:

$\displaystyle\lambda+\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}=\sum_{\alpha}\hat{\mu}_{\alpha,\lambda}\phi_{\alpha}.$
(3.22)

This observation reveals that chemical potentials in (3.2) solely occur in the
form $\hat{\mu}_{\alpha,\lambda}$. In other words, the chemical potentials
$\hat{\mu}_{\alpha}$ are tightly connected with the Lagrange multiplier
$\lambda$. This is consistent with the examination that the addition of
$\sum_{\alpha}\phi_{\alpha}-1$ does not alter the free energy. Indeed, we
have:

$\displaystyle\Psi=$
$\displaystyle~{}\hat{\Psi}\left(\left\\{\phi_{\alpha}\right\\}_{\alpha=1,...,N},\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha=1,...,N}\right)+\lambda\left(\sum_{\alpha}\phi_{\alpha}-1\right),$
(3.23)

and the associated chemical potential quantities ($\alpha=1,...,N$) naturally
include the Lagrange multiplier $\lambda$:

$\displaystyle\hat{\mu}_{\alpha,\lambda}=\dfrac{\partial\Psi}{\partial\phi_{\alpha}}-{\rm
div}\dfrac{\partial\Psi}{\partial\nabla\phi_{\alpha}},$ (3.24)

where we recall (3.19b).

###  3.3 Alternative free energy classes

As mentioned in Remark 3.4, as an alternative for (3.5), we explore the
approach of working with a class that depends on concentration. This
exploration is motivated by its occurrence in the literature on two-phase
models (e.g. [16]). We consider the following constitutive class:

$\displaystyle\Psi=$
$\displaystyle~{}\check{\Psi}\left(\left\\{c_{\alpha}\right\\},\left\\{\nabla
c_{\alpha}\right\\}\right),$ (3.25)

subject to the summation constraint (2.5a). Alongside with free energy class
(3.25) we introduce the chemical potential quantities ($\alpha=1,...,N$):

$\displaystyle\check{\mu}_{\alpha}=$
$\displaystyle~{}\dfrac{\partial\hat{\Psi}}{\partial c_{\alpha}}-{\rm
div}\dfrac{\partial\hat{\Psi}}{\partial\nabla c_{\alpha}}.$ (3.26)

In Appendix C we provide the derivation of the modeling restriction that
emerges from the constitutive class (3.25). The modeling restriction takes the
form:

$\displaystyle\left(\mathbf{T}+\check{\lambda}\mathbf{I}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle-\sum_{\alpha}\nabla\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\cdot\mathbf{H}_{\alpha}-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\zeta_{\alpha}$
$\displaystyle\geq 0,$ (3.27)

where $\check{\lambda}$ is the Lagrange multiplier associated with the
constraint (2.29b). Noting that the volume fractions and concentrations are
connected via (2.7):

$\displaystyle\phi_{\alpha}=$
$\displaystyle~{}\phi_{\alpha}\left(\left\\{c_{\beta}\right\\}\right),$
(3.28a) $\displaystyle c_{\beta}=$
$\displaystyle~{}c_{\beta}\left(\left\\{\phi_{\alpha}\right\\}\right),$
(3.28b)

the identification

$\displaystyle\hat{\Psi}\left(\left\\{\phi_{\alpha}\right\\}\right)=\check{\Psi}\left(\left\\{c_{\beta}\right\\}\right)$
(3.29)

reveals that the free energy classes coincide. Given that the initial modeling
restriction is the same for both classes we conclude that the resulting
modeling restrictions must coincide as well. In other words, the modeling
restriction is independent of the choice of order parameters.

######  Theorem 3.9 (Equivalence of modeling restrictions).

The modeling restrictions (3.2) and (3.3) are equivalent.

An alternative path to show equivalence of the modeling restrictions, one
could apply the variable transformation (3.28) defined in (2.7) to show that
(3.3) coincides with (3.2). The challenge of this approach lies in the
combination of working with $N$ order parameters that satisfy the saturation
constraints (2.5). We discuss this approach in Appendix B.

Guided by Theorem 3.9, we proceed with the formulation of the modeling
restriction presented in (3.2).

###  3.4 Selection of constitutive models

By means of the Colemann-Noll concept, we utilize (3.2) as a guiding principle
to design constitutive models. Inspired by the specific form of the constraint
(3.2), we restrict ourselves to mixture stress tensors $\mathbf{T}$,
constituent peculiar velocities $\mathbf{J}_{\alpha}$, and constituent mass
transfer terms $\gamma_{\alpha}$ that belong to the constitutive classes:

$\displaystyle\mathbf{T}=$
$\displaystyle~{}\hat{\mathbf{T}}\left(\left\\{\phi_{\alpha}\right\\},\left\\{\nabla\phi_{\alpha}\right\\},\left\\{\hat{\mu}_{\alpha,\lambda}\right\\}\right),$
(3.30a) $\displaystyle\mathbf{J}_{\alpha}=$
$\displaystyle~{}\hat{\mathbf{J}}_{\alpha}\left(\left\\{\phi_{\alpha}\right\\},\left\\{\nabla
g_{\alpha}\right\\}\right),$ (3.30b) $\displaystyle\mathbf{j}_{\alpha}=$
$\displaystyle~{}\hat{\mathbf{j}}_{\alpha}\left(\left\\{\phi_{\alpha}\right\\},\left\\{\nabla
g_{\alpha}\right\\}\right),$ (3.30c) $\displaystyle\zeta_{\alpha}=$
$\displaystyle~{}\hat{\zeta}_{\alpha}\left(\left\\{\phi_{\alpha}\right\\},\left\\{g_{\alpha}\right\\}\right),$
(3.30d)

and define
$\hat{\mathbf{H}}_{\alpha}=\hat{\mathbf{J}}_{\alpha}+\hat{\mathbf{j}}_{\alpha}$.
Generally speaking, the introduction of the class (3.30b) deviates from
continuum mixture theory. Arguably, this is not the case when
$\hat{\mathbf{J}}_{\alpha}=0$, which models the situation of matching
velocities $\mathbf{v}_{\alpha}=\mathbf{v}_{\beta}$. We return to this case in
Section 4.2.

We do not seek the most complete constitutive theory, rather our goal is to
find a set of practical constitutive models compatible with (3.2). To this
end, we aim to identify constitutive models (3.30) so that all three terms in
(3.2) are positive:

$\displaystyle\left(\mathbf{T}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha,\lambda}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle\geq 0,$ (3.31a) $\displaystyle-\sum_{\alpha}\nabla
g_{\alpha}\cdot\mathbf{J}_{\alpha}$ $\displaystyle\geq 0,$ (3.31b)
$\displaystyle-\sum_{\alpha}\nabla g_{\alpha}\cdot\mathbf{j}_{\alpha}$
$\displaystyle\geq 0,$ (3.31c)
$\displaystyle-\sum_{\alpha}g_{\alpha}\zeta_{\alpha}$ $\displaystyle\geq 0.$
(3.31d)

######  Remark 3.10 (Onsager reciprocal relations).

As mentioned above, our objective is to find a set of practical constitutive
models. As more complete theory follows from working with the original
constraint (3.2), and extending the dependency of the classes (3.30). In
particular, the classes may be interconnected. The well-known Onsager
reciprocal relations take a central place in this framework. We refer to [48,
49].

In the following we provide constitutive models for the mixture stress tensor,
constituent peculiar velocities, and constituent mass transfer, respectively.

Mixture stress tensor. We select the following constitutive model for the
stress tensor:

$\displaystyle\hat{\mathbf{T}}=-\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\left(\sum_{\alpha}\hat{\mu}_{\alpha,\lambda}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}+\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I}),$ (3.32)

subject to the symmetry condition:

$\displaystyle\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}=\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\otimes\nabla\phi_{\alpha},$
(3.33)

where the scalar field $\nu\geq 0$ is the mixture dynamic viscosity,
$\bar{\lambda}\geq-2/d$ is a scalar, and $d$ is the number of dimensions. The
condition (3.33) ensures compatibility with the angular momentum constraint
(2.25c).

######  Lemma 3.11 (Compatibility mixture stress tensor).

The mixture stress tensor (3.32) adheres to the constraint (3.31a).

###### Proof.

An elementary calculation gives:

$\displaystyle\left(\hat{\mathbf{T}}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\sum_{\alpha}\hat{\mu}_{\alpha,\lambda}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle=$ $\displaystyle
2\nu\left(\nabla^{s}\mathbf{v}-\frac{1}{d}({\rm
div}\mathbf{v})\mathbf{I}\right):\left(\nabla^{s}\mathbf{v}-\frac{1}{d}({\rm
div}\mathbf{v})\mathbf{I}\right)+\nu\left(\bar{\lambda}+\frac{2}{d}\right)\left({\rm
div}\mathbf{v}\right)^{2}$ $\displaystyle\geq 0.$ (3.34)

∎

Constituent peculiar velocities. We remark that the balance condition (2.13)
allows us to convert the constraint (3.31b) into:

$\displaystyle-\sum_{\alpha}\nabla
g_{\alpha}\cdot\mathbf{J}_{\alpha}=-\sum_{\alpha}\nabla(g_{\alpha}-g_{\beta})\cdot\mathbf{J}_{\alpha}$
$\displaystyle\geq 0,$ (3.35)

for some $\beta\in\left\\{1,...,N\right\\}$. Based on this formulation, we
choose the following constituent peculiar velocities:

$\displaystyle\hat{\mathbf{J}}_{\alpha}=-\sum_{\beta}\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right),$
(3.36)

for some constitutive quantity
$\mathbf{M}_{\alpha{\beta}}=\mathbf{M}_{{\beta}\alpha}$ that is a positive
semi-definite tensor for given $\alpha$ and ${\beta}$, with the same
dependencies as (3.31b). We refer to it as a mobility tensor.

######  Remark 3.12 (Lagrange multiplier in constituent peculiar velocities).

In most incompressible $N$-component models, the Lagrange multiplier $\lambda$
does not appear in the constituent peculiar velocities
$\hat{\mathbf{J}}_{\alpha}$. In contrast, the proposed framework indicates the
presence of the Lagrange multiplier (as part of $g_{\alpha}$). In the case of
matching constituent densities ($\rho_{\alpha}=\rho$ for $\alpha=1,...,N$) the
Lagrange multiplier vanishes:
$g_{\alpha}-g_{\beta}=\rho^{-1}(\hat{\mu}_{\alpha}-\hat{\mu}_{\beta}$).

######  Remark 3.13 (Related constitutive models).

The model (3.36) matches (for the isotropic case
$\mathbf{M}_{\alpha{\beta}}=M_{\alpha{\beta}}\mathbf{I}$) that of [45]. The
model (3.36) also closely resembles the form adopted in [26]. Both closure
models involve the Lagrange multiplier $\lambda$; a difference lies in fact
that the model proposed by [26] depends on the numbering of the constituents.

######  Lemma 3.14 (Compatibility constituent peculiar velocities).

The choice (3.36) aligns with both the balance (2.13), and the restriction
(3.31b).

###### Proof.

The symmetry property $\mathbf{M}_{\alpha{\beta}}=\mathbf{M}_{{\beta}\alpha}$
ensures conformity with (2.13), and adherence to (3.31b) follows from the
sequence of identities:

$\displaystyle-\sum_{\alpha}\hat{\mathbf{J}}_{\alpha}\cdot\nabla g_{\alpha}=$
$\displaystyle~{}\sum_{\alpha,{\beta}}\left(\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)\right)\cdot\nabla
g_{\alpha}$ $\displaystyle=$
$\displaystyle~{}\tfrac{1}{2}\sum_{\alpha,{\beta}}\left(\nabla\left(g_{\alpha}-g_{\beta}\right)\right)^{T}\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle\geq$ $\displaystyle~{}0.$ (3.37)

∎

Constituent diffusive flux. Analogously to the constituent peculiar
velocities, we select:

$\displaystyle\hat{\mathbf{j}}_{\alpha}=-\sum_{\beta}\mathbf{K}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right),$
(3.38)

for some constitutive quantity
$\mathbf{K}_{\alpha{\beta}}=\mathbf{K}_{{\beta}\alpha}$ that is a positive
semi-definite tensor for given $\alpha$ and ${\beta}$, with the same
dependencies as (3.31c).

######  Lemma 3.15 (Compatibility constituent diffusive fluxes).

The choice (3.38) aligns with both the balance (2.28b), and the restriction
(3.31c).

Constituent mass transfer. We select the constituent mass transfer terms
analogously to the constituent peculiar velocities:

$\displaystyle\hat{\zeta}_{\alpha}=-\sum_{\beta}m_{\alpha{\beta}}\left(g_{\alpha}-g_{\beta}\right),$
(3.39)

where the scalar mobility $m_{\alpha{\beta}}=m_{{\beta}\alpha}\geq 0$ vanishes
in the single fluid region $m_{\alpha{\beta}}(\phi_{\alpha}=0)=0$.

######  Lemma 3.16 (Compatibility mass transfer).

The choice (3.39) is compatible with the balance of mass supply (2.27b), and
the constraint (3.31d).

###### Proof.

Analogously to the proof of Lemma 3.14, compatibility with (2.27a) is
immediate from the symmetry property: $m_{\alpha{\beta}}=m_{{\beta}\alpha}$,
and consistency with (3.31d) follows from:

$\displaystyle-\sum_{\alpha}\hat{\zeta}_{\alpha}g_{\alpha}=\sum_{\alpha,{\beta}}m_{\alpha{\beta}}(g_{\alpha}-g_{\beta})g_{\alpha}=\tfrac{1}{2}\sum_{\alpha,{\beta}}m_{\alpha{\beta}}(g_{\alpha}-g_{\beta})^{2}\geq
0.$ (3.40)

∎

This finalizes the construction of constitutive models compatible with the
imposed energy-dissipative postulate. Substitution of the models (3.32),
(3.36), and (3.39) yields the class of incompressible $N$-component models:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda$
$\displaystyle+{\rm
div}\left(\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$
(3.41a) $\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{J}}_{\alpha}+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{j}}_{\alpha}-\hat{\zeta}_{\alpha}$ $\displaystyle=~{}0,$
(3.41b)
$\displaystyle\hat{\mu}_{\alpha}-\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}+{\rm
div}\left(\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)$
$\displaystyle=~{}0,$ (3.41c)
$\displaystyle\hat{\mathbf{J}}_{\alpha}+\sum_{\beta}\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle~{}=0,$ (3.41d)
$\displaystyle\hat{\mathbf{j}}_{\alpha}+\sum_{\beta}\mathbf{K}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle~{}=0,$ (3.41e)
$\displaystyle\hat{\zeta}_{\alpha}+\sum_{\beta}m_{\alpha{\beta}}\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle~{}=0.$ (3.41f)

Formulation (3.41) constitutes a class of models in the sense that particular
closure relations ($\mathbf{M}_{\alpha{\beta}}$, $\mathbf{K}_{\alpha{\beta}}$
and $m_{\alpha{\beta}}$) need to be specified. Given these relations, for each
specification of the free energy $\hat{\Psi}$, (3.41) is a well-defined closed
model; this model is invariant to the renumbering of the constituents,
invariant to the set of independent variables, reduction-consistent.
Additionally, it exhibits energy-dissipation, which we explicitly state in the
following theorem.

######  Theorem 3.17 (Compatibility energy dissipation).

The model (3.41) is compatible with the energy dissipation condition (3.3).

###### Proof.

This follows from Lemma 3.11, Lemma 3.14, and Lemma 3.16. In particular, the
dissipation takes the form:

$\displaystyle\mathscr{D}=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}2\nu\left(\nabla^{s}\mathbf{v}-\frac{1}{d}({\rm
div}\mathbf{v})\mathbf{I}\right):\left(\nabla^{s}\mathbf{v}-\frac{1}{d}({\rm
div}\mathbf{v})\mathbf{I}\right)+\nu\left(\lambda+\frac{2}{d}\right)\left({\rm
div}\mathbf{v}\right)^{2}$
$\displaystyle~{}\quad\quad+\frac{1}{2}\displaystyle\sum_{\alpha,{\beta}}\left(\nabla(g_{\alpha}-g_{\beta})\right)^{T}(\mathbf{M}_{\alpha{\beta}}+\mathbf{K}_{\alpha{\beta}})\nabla(g_{\alpha}-g_{\beta})$
$\displaystyle~{}\quad\quad+\frac{1}{2}\sum_{\alpha,{\beta}}m_{\alpha{\beta}}(g_{\alpha}-g_{\beta})^{2}~{}{\rm
d}v\geq 0.$ (3.42)

∎

##  4 Model characteristics

In this section, we explore the characteristics of the modeling framework
outlined in Section 3. To this end, we discuss alternative –equivalent –
formulations in Section 4.1. We present the case of matching velocities in
Section 4.2. Subsequently, Section 4.3 details the equilibrium
characteristics.

###  4.1 Alternative formulations

As discussed in Sections 2 and 3, the unified modeling framework outlined in
these sections is invariant to the choice of variables. However, it is
worthwhile to discuss some of the formulations that are associated with
particular variables.

First, we note that one can identify a pressure quantity in the model as:

$\displaystyle p:=\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}.$
(4.1)

With this choice, the model takes the more compact form:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla(\lambda+p)+{\rm
div}\left(\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$ (4.2a)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{H}}_{\alpha}-\hat{\zeta}_{\alpha}$ $\displaystyle=~{}0.$
(4.2b)

In accordance with the first metaphysical principle of continuum mixture
theory, the mixture free energy is comprised of constituent free energies:

$\displaystyle\hat{\Psi}=\sum_{\alpha}\hat{\Psi}_{\alpha},$ (4.3)

where $\hat{\Psi}_{\alpha}$ are the volume-measure constituent free energies.
Utilizing (4.3) we observe that the pressure satisfies Dalton’s law:

$\displaystyle p=$ $\displaystyle~{}\sum_{\alpha}p_{\alpha},$ (4.4a)
$\displaystyle p_{\alpha}=$
$\displaystyle~{}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}_{\alpha},$ (4.4b)

where $p_{\alpha}$ is the partial pressure of constituent $\alpha$. Thus, the
split (4.3) reveals that the system may be written as:

$\displaystyle\sum_{\beta}\left(\partial_{t}(\tilde{\rho}_{\beta}\mathbf{v})+{\rm
div}\left(\tilde{\rho}_{\beta}\mathbf{v}\otimes\mathbf{v}\right)+\phi_{\beta}\nabla\lambda+\nabla
p_{\beta}-\tilde{\rho}_{\beta}\mathbf{g}\right.$ $\displaystyle\left.+{\rm
div}\left(\nabla\phi_{\beta}\otimes\dfrac{\partial\hat{\Psi}_{\beta}}{\partial\nabla\phi_{\beta}}\right)\right)-{\rm
div}\left(\nu(2\nabla^{s}\mathbf{v}+\lambda({\rm
div}\mathbf{v})\mathbf{I})\right)$ $\displaystyle=~{}0,$ (4.5a)
$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{v})-\sum_{{\beta}=1}^{N}{\rm
div}\left(\mathbf{M}_{\alpha{\beta}}\nabla\left(g_{\alpha}-g_{\beta}\right)\right)+\sum_{{\beta}=1}^{N}m_{\alpha{\beta}}\left(g_{\alpha}-g_{\beta}\right)$
$\displaystyle=~{}0.$ (4.5b)

An alternative compact form is obtained with the aid of the following lemma.

######  Lemma 4.1 (Identity free energy).

The free energy contributions collapse into:

$\displaystyle{\rm
div}\left(\left(\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}-\hat{\Psi}\right)\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)=\sum_{\alpha}\phi_{\alpha}\nabla\hat{\mu}_{\alpha}.$
(4.6)

###### Proof.

See Lemma A.4. ∎

Invoking Lemma 4.1, model (3.41) takes a more compact form:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda+\sum_{\alpha}\phi_{\alpha}\nabla\hat{\mu}_{\alpha}$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$ (4.7a)
$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{H}}_{\alpha}-\hat{\zeta}_{\alpha}$ $\displaystyle=~{}0.$
(4.7b)

Considering the third and fourth term in the momentum equation in isolation,
this can be written as:

$\displaystyle\nabla\lambda+\sum_{\alpha}\phi_{\alpha}\nabla\hat{\mu}_{\alpha}=\sum_{\alpha}\phi_{\alpha}\nabla(\lambda+\hat{\mu}_{\alpha})=\sum_{\alpha}\tilde{\rho}_{\alpha}\nabla
g_{\alpha}.$ (4.8)

Similarly, in the mass balance (4.7b), we observe that the chemical potentials
$\hat{\mu}_{\alpha}$ and the Lagrange multiplier $\lambda$ appear solely as a
sum via $g_{\alpha}$.

Additionally, we note that the model can alternatively be written in a form
that more closely links to existing phase-field models:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda+\sum_{{\beta}}\phi_{\beta}\nabla\hat{\mu}_{\beta}$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$ (4.9a)
$\displaystyle\partial_{t}\rho+{\rm div}(\rho\mathbf{v})$
$\displaystyle=~{}0,$ (4.9b) $\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{H}}_{\alpha}-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0,$ (4.9c)

for $\alpha=1,...,N-1$. We note that the combination of the mixture mass
balance (4.9b) and the $N-1$ constituent balance laws (4.9c) are equivalent to
the $N$ balance laws (4.5b) or $N$ balance laws (4.9c).

While we refrain from discussing formulations that adopt concentration
variables, we discuss a formulation in terms of the volume-averaged velocity
$\mathbf{u}$. Inserting the constitutive model for the diffusive flux (3.36)
into the relation Lemma 2.4 we obtain:

$\displaystyle\mathbf{v}=\mathbf{u}+\rho^{-1}\displaystyle\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}.$
(4.10)

By substituting this identity, we express the model using the volume-averaged
velocity:

$\displaystyle\partial_{t}\left(\rho\mathbf{u}+\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\right)+{\rm
div}\left(\rho\mathbf{u}\otimes\mathbf{u}+\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\otimes\mathbf{u}\right.$
$\displaystyle+\left.\mathbf{u}\otimes\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}+\rho^{-1}\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\otimes\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\right)+\nabla\lambda+\sum_{{\beta}}\phi_{\beta}\nabla\hat{\mu}_{\beta}$
$\displaystyle-{\rm
div}\left(\nu\left(2\nabla^{s}\left(\mathbf{u}+\rho^{-1}\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\right)+\bar{\lambda}{\rm
div}\left(\mathbf{u}+\rho^{-1}\sum_{\beta}\hat{\mathbf{J}}_{\beta}^{u}\right)\mathbf{I}\right)\right)$
$\displaystyle=~{}0,$ (4.11a) $\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}\left(\phi_{\alpha}\mathbf{u}\right)+{\rm
div}\hat{\mathbf{h}}_{\alpha}^{u}+\rho_{\alpha}^{-1}{\rm
div}\hat{\mathbf{j}}_{\alpha}-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0.$ (4.11b)

Recalling (2.29c), we have:

$\displaystyle{\rm
div}\mathbf{u}=\sum_{\alpha}\rho_{\alpha}^{-1}(\hat{\zeta}_{\alpha}-{\rm
div}\hat{\mathbf{j}}_{\alpha}),$ (4.12)

where the right-hand side vanishes when densities match or mass transfer is
absent. Arguably, the formulation (4.11) is rather involved. We discuss a
simplification in the next subsection.

###  4.2 Matching velocities

We consider the case in which the peculiar velocities are zero; taking
$\mathbf{M}_{\alpha{\beta}}=0$, we find:

$\displaystyle\hat{\mathbf{J}}_{\alpha}=\hat{\mathbf{h}}_{\alpha}=\hat{\mathbf{J}}_{\alpha}^{u}=\hat{\mathbf{h}}_{\alpha}^{u}=0.$
(4.13)

As a consequence the mass-averaged and volume-averaged velocities are equal:

$\displaystyle\mathbf{v}=\mathbf{u}.$ (4.14)

This choice models the situation where the constituent velocities are
matching. We explicitly state the simplified formulations of the model:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda+\sum_{\beta}\phi_{\beta}\nabla\hat{\mu}_{\beta}$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$
(4.15a) $\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\mathbf{j}_{\alpha}-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0.$ (4.15b)

and obviously:

$\displaystyle\partial_{t}(\rho\mathbf{u})+{\rm
div}\left(\rho\mathbf{u}\otimes\mathbf{u}\right)+\nabla\lambda+\sum_{\beta}\phi_{\beta}\nabla\hat{\mu}_{\beta}$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{u}+\bar{\lambda}({\rm
div}\mathbf{u})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$
(4.16a) $\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{u})+\rho_{\alpha}^{-1}{\rm
div}\mathbf{j}_{\alpha}-\rho_{\alpha}^{-1}\hat{\zeta}_{\alpha}$
$\displaystyle=~{}0.$ (4.16b)

The formulation (4.16) demonstrates that a simplified – consistent – model in
terms of the volume-averaged velocity involves a straightforward momentum
equation. We emphasize that the volume-averaged velocity $\mathbf{u}$ is in
general not divergence-free (recall (2.29c)).

###  4.3 Equilibrium conditions

We utilize formulation (4.7) to study equilibrium properties. We characterize
the set equilibrium solutions
$\left\\{\mathbf{q}_{E}=(\mathbf{v}_{E},\phi_{\alpha,E},\lambda_{E},\mu_{\alpha,E})\right\\}$
of (4.7) as stationary solutions subject to boundary conditions for which the
dissipation vanishes: $\mathscr{D}(\mathbf{q}_{E})=0$. Invoking (3.4) yields
the conditions:

$\displaystyle\left(\nabla^{s}\mathbf{v}_{E}-\frac{1}{d}({\rm
div}\mathbf{v}_{E})\mathbf{I}\right):\left(\nabla^{s}\mathbf{v}_{E}-\frac{1}{d}({\rm
div}\mathbf{v}_{E})\mathbf{I}\right)=$ $\displaystyle~{}0,$ (4.17a)
$\displaystyle\left(\lambda_{E}+\frac{2}{d}\right)\left({\rm
div}\mathbf{v}_{E}\right)^{2}=$ $\displaystyle~{}0,$ (4.17b)
$\displaystyle\displaystyle\sum_{\alpha,{\beta}}\left(\nabla(g_{\alpha,E}-g_{{\beta},E})\right)^{T}(\mathbf{M}_{\alpha{\beta},E}+\mathbf{K}_{\alpha{\beta},E})\nabla(g_{\alpha,E}-g_{{\beta},E})=$
$\displaystyle~{}0,$ (4.17c)
$\displaystyle\sum_{\alpha,{\beta}}m_{\alpha{\beta},E}(g_{\alpha,E}-g_{{\beta},E})^{2}=$
$\displaystyle~{}0,$ (4.17d)

in $\Omega$, where the subscript $E$ denotes the equilibrium configuration of
the quantity. We deduce from (4.17a)-(4.17b) that $\mathbf{v}_{E}$ are rigid
body motions. Simplifying the analysis, we take $\mathbf{v}_{E}=0$, which
causes the inertia terms to vanish. Additionally, we assume the absence of
gravitational forces ($\mathbf{g}=0$). Substituting into the momentum equation
(4.7a) provides:

$\displaystyle\nabla\lambda_{E}+\sum_{\alpha}\phi_{\alpha,E}\nabla\hat{\mu}_{\alpha,E}=0.$
(4.18)

Recalling (4.8) we deduce:

$\displaystyle\sum_{\alpha}\tilde{\rho}_{\alpha,E}\nabla g_{\alpha,E}=0.$
(4.19)

Next, we distinguish the case that includes mass transfer, and the case that
does not ($m_{\alpha{\beta},E}=0$). Starting with the first case, we recall
that $\mathbf{M}_{\alpha{\beta},E}$ is degenerate. Consider the non-trivial
situation $\mathbf{M}_{\alpha{\beta},E}\neq\mathbf{0}$, (4.17c) implies that

$\displaystyle g_{\alpha,E}=g_{{\beta},E}+c_{\alpha,{\beta}},$ (4.20)

in the interfacial region ($-1<\phi_{\alpha,E}<1$). Here $c_{\alpha,{\beta}}$
are constants that satisfy
$\sum_{\alpha=1,...,N-1}c_{\alpha,\alpha+1}+c_{N,1}=0$. As a consequence,
$\hat{\mathbf{J}}_{\alpha,E}=0$ in equilibrium. The constraints (4.19) and
(4.20) form $N$ equilibrium conditions for $N$ variables. In case with mass
transfer we find $c_{\alpha,{\beta}}=0$, i.e.

$\displaystyle g_{\alpha,E}=g_{{\beta},E},$ (4.21)

and hence $\hat{\mathbf{J}}_{\alpha,E}=0$, and $\hat{\zeta}_{\alpha,E}=0$.

##  5 Connections to existing models

This section provides connections with existing models. First, we discuss the
binary-phase situation in Section 5.1. Next, in Section 5.3 discusses the link
to a model with $N$-momentum equations.

###  5.1 Binary-phase case

In this section we restrict to binary mixtures ($\alpha=1,2)$, and compare
with the framework presented in [23]. A formulation of this two-phase modeling
framework is:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda+\phi\nabla\breve{\mu}-{\rm
div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$ (5.1a)
$\displaystyle\partial_{t}\rho+{\rm div}(\rho\mathbf{v})$
$\displaystyle=~{}0.$ (5.1b) $\displaystyle\partial_{t}\phi+{\rm
div}(\phi\mathbf{v})-{\rm
div}\left(\breve{\mathbf{M}}\nabla\left(\breve{\mu}+\omega\lambda\right)\right)+{\rm
div}\left(\breve{m}\left(\breve{\mu}+\omega\breve{\lambda}\right)\right)$
$\displaystyle=~{}0.$ (5.1c)

Here $\phi$ is the phase-field quantity defined as the difference between
volume fractions:

$\displaystyle\phi(\phi_{1},\phi_{2})=\phi_{1}-\phi_{2},$ (5.2)

where we recall $\phi_{1}+\phi_{2}=1$. The chemical potential quantity is
defined as:

$\displaystyle\breve{\mu}=\dfrac{\partial\Psi}{\partial\phi}-{\rm
div}\dfrac{\partial\Psi}{\partial\nabla\phi}.$ (5.3)

Finally, the quantity $\omega$ is
$\omega=(\rho_{1}^{-1}-\rho_{2}^{-1})/(\rho_{1}^{-1}+\rho_{2}^{-1})$. The
model (4.7) takes for binary mixtures the following form:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\nabla\lambda+\phi_{1}\nabla\hat{\mu}_{1}+\phi_{2}\nabla\hat{\mu}_{2}$
$\displaystyle-{\rm div}\left(\nu(2\nabla^{s}\mathbf{v}+\bar{\lambda}({\rm
div}\mathbf{v})\mathbf{I})\right)-\rho\mathbf{g}$ $\displaystyle=~{}0,$ (5.4a)
$\displaystyle\partial_{t}\phi_{1}+{\rm
div}(\phi_{1}\mathbf{v})-\rho_{1}^{-1}{\rm
div}\left(\mathbf{M}\nabla\left(g_{1}-g_{2}\right)\right)+\rho_{1}^{-1}{\rm
div}\left(m\left(g_{1}-g_{2}\right)\right)$ $\displaystyle=~{}0,$ (5.4b)
$\displaystyle\partial_{t}\phi_{2}+{\rm
div}(\phi_{2}\mathbf{v})-\rho_{2}^{-1}{\rm
div}\left(\mathbf{M}\nabla\left(g_{2}-g_{1}\right)\right)+\rho_{2}^{-1}{\rm
div}\left(m\left(g_{2}-g_{1}\right)\right)$ $\displaystyle=~{}0,$ (5.4c)

where $\mathbf{M}=\mathbf{M}_{12}=\mathbf{M}_{21}$ and $m=m_{12}=m_{21}$. By
means of variable transformation, we aim to express the model in terms of the
quantities of model (5.1). The mass balance equations (5.4b) and (5.4c) can be
written as:

$\displaystyle\partial_{t}\phi+{\rm div}(\phi\mathbf{v})-{\rm
div}\left((\rho_{1}^{-1}+\rho_{2}^{-1})\mathbf{M}\nabla\left(g_{1}-g_{2}\right)\right)$
$\displaystyle+{\rm
div}\left((\rho_{1}^{-1}+\rho_{2}^{-1})m\left(g_{1}-g_{2}\right)\right)$
$\displaystyle=~{}0,$ (5.5a) $\displaystyle\partial_{t}\rho+{\rm
div}(\rho\mathbf{v})$ $\displaystyle=~{}0.$ (5.5b)

With the aim of identifying the two models, we select the relations
$\breve{\mathbf{M}}=(\rho_{1}^{-1}+\rho_{2}^{-1})^{2}\mathbf{M}$ and
$\breve{m}=(\rho_{1}^{-1}+\rho_{2}^{-1})^{2}m$. As a consequence, in case the
following identities hold:

$\displaystyle\phi_{1}\nabla\hat{\mu}_{1}+\phi_{2}\nabla\hat{\mu}_{2}=$
$\displaystyle~{}\phi\nabla\breve{\mu},$ (5.6a) $\displaystyle g_{1}-g_{2}=$
$\displaystyle~{}(\rho_{1}^{-1}+\rho_{2}^{-1})\left(\mu+\omega\lambda\right),$
(5.6b)

we find that (5.4) coincides with (5.1), i.e. the two models match. Recalling
that the chemical potentials $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$ depend on
precise formulation of the free energy, we impose the conditions:
$\mu_{1}+\mu_{2}=0$ and $\breve{\mu}=\mu_{1}=-\mu_{2}$. These conditions are
inspired by the chain rule for chemical potentials, where we note the
derivatives $\partial\phi/\partial\phi_{1}=1$,
$\partial\phi/\partial\phi_{2}=-1$ of (5.2). This shows that a particular
formulation of the binary constituent case coincides with the unified modeling
framework proposed in [23].

###  5.2 $N$-phase model Dong (2018)

The $N$-phase incompressible model proposed by Dong [28] is given by

$\displaystyle\rho\left(\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}\right)+\mathbf{J}^{\prime}\cdot\nabla\mathbf{u}+\nabla\lambda^{\prime}-{\rm
div}\left(\nu^{\prime}\nabla^{s}\mathbf{u}\right)$
$\displaystyle+\sum_{\beta}{\rm
div}\left(\nabla\phi_{\alpha}\otimes\frac{\partial\Psi^{\prime}}{\partial(\nabla\phi_{\beta})}\right)=$
$\displaystyle~{}0,$ (5.7a) $\displaystyle{\rm div}\mathbf{u}=$
$\displaystyle~{}0,$ (5.7b)
$\displaystyle\partial_{t}\phi_{\alpha}+\mathbf{u}\cdot\nabla\phi_{\alpha}-\sum_{\beta}{\rm
div}\left(m_{\alpha{\beta}}^{\prime}\nabla\left(\frac{\partial\Psi^{\prime}}{\partial\phi_{\beta}}-{\rm
div}\left(\frac{\partial\Psi^{\prime}}{\partial(\nabla\phi_{\beta})}\right)\right)\right)=$
$\displaystyle~{}0,$ (5.7c)

for $\alpha=1,...,N$, where $\lambda^{\prime}$ is the Lagrange multiplier
(pressure), $\nu^{\prime}$ is the dynamic viscosity, $\Psi^{\prime}$ is the
free energy, and $m_{\alpha{\beta}}^{\prime}$ is the mobility. For the purpose
of comparing the model (5.7) to the proposed framework, we define a chemical
potential quantity:

$\displaystyle\mu_{\beta}^{\prime}=\frac{\partial\Psi^{\prime}}{\partial\phi_{\beta}}-{\rm
div}\left(\frac{\partial\Psi^{\prime}}{\partial(\nabla\phi_{\beta})}\right).$
(5.8)

Invoking Lemma 4.1, we rewrite the model (5.7) as:

$\displaystyle\rho\left(\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}\right)+\mathbf{J}^{\prime}\cdot\nabla\mathbf{u}+\nabla\tilde{\lambda}^{\prime}+\sum_{\beta}\phi_{\alpha}\nabla\mu_{\alpha}^{\prime}-{\rm
div}\left(\nu^{\prime}\nabla^{s}\mathbf{u}\right)=$ $\displaystyle~{}0,$
(5.9a) $\displaystyle{\rm div}\mathbf{u}=$ $\displaystyle~{}0,$ (5.9b)
$\displaystyle\partial_{t}\phi_{\alpha}+\mathbf{u}\cdot\nabla\phi_{\alpha}-\sum_{\beta}{\rm
div}\left(m_{\alpha{\beta}}^{\prime}\nabla\mu_{\beta}^{\prime}\right)=$
$\displaystyle~{}0,$ (5.9c)

with

$\displaystyle\tilde{\lambda}^{\prime}=\lambda^{\prime}+\Psi^{\prime}-\sum_{\beta}\phi_{\alpha}\mu_{\alpha}^{\prime}.$
(5.10)

This model (5.9) is not compatible with the framework proposed in the current
paper. In particular, comparing (5.9) with (4.11), we observe that:

  * •

model (5.9) does not contain each of the peculiar velocity terms in the
momentum equation; this applies to both inertia and viscous terms;

  * •

model (5.9) does not include mass transfer terms;

  * •

the constitutive model for the diffusive flux in (5.9) is different; in
particular the terms do not scale with $g_{\alpha}-g_{\beta}$, and the
Lagrange multiplier is absent. As a consequence, the equilibrium conditions
are different.

###  5.3 Class-II mixture model

We compare the proposed unified modeling framework with an incompressible
mixture model presented in [45]:

$\displaystyle\partial_{t}\tilde{\rho}_{\alpha}+{\rm
div}(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha})+\sum_{\beta}\breve{m}_{\alpha{\beta}}(\breve{g}_{\alpha}-\breve{g}_{\beta})$
$\displaystyle=~{}0,$ (5.11a)
$\displaystyle\partial_{t}(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha})+{\rm
div}\left(\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha}\otimes\mathbf{v}_{\alpha}\right)+\phi_{\alpha}\nabla\left(\breve{\lambda}+\breve{\mu}_{\alpha}\right)$
$\displaystyle-{\rm
div}\left(\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{v}_{\alpha}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{v}_{\alpha}\right)\right)-\tilde{\rho}_{\alpha}\mathbf{b}$
$\displaystyle+\displaystyle\sum_{{\beta}}R_{\alpha{\beta}}(\mathbf{v}_{\alpha}-\mathbf{v}_{\beta})+\tfrac{1}{2}\sum_{\beta}\breve{m}_{\alpha{\beta}}(\breve{g}_{\alpha}-\breve{g}_{\beta})(\mathbf{v}_{\alpha}+\mathbf{v}_{\beta})$
$\displaystyle=~{}0,$ (5.11b)

for constituents $\alpha=1,...,N$. Here $\mathbf{v}_{\alpha}$ is the
constituent velocity, $\breve{\lambda}$ is a Lagrange
multiplier,$\tilde{\nu}_{\alpha}$ the constituent dynamical viscosity,
$\breve{\bar{\lambda}}_{\alpha}\geq 2/d$, $\nabla^{s}\mathbf{v}_{\alpha}$ the
constituent symmetric velocity gradient, and $\breve{m}_{\alpha{\beta}}$ and
$R_{\alpha{\beta}}$ are symmetric matrices (for the properties see [45]). This
model considers the free energy class:

$\displaystyle\Psi=$ $\displaystyle~{}\sum_{\alpha}\breve{\Psi}_{\alpha},$
(5.12a) $\displaystyle\breve{\Psi}_{\alpha}=$
$\displaystyle~{}\breve{\Psi}_{\alpha}\left(\phi_{\alpha},\nabla\phi_{\alpha}\right).$
(5.12b)

The associated constituent chemical potentials are defined as:

$\displaystyle\breve{\mu}_{\alpha}=$
$\displaystyle~{}\dfrac{\partial\breve{\Psi}_{\alpha}}{\partial\phi_{\alpha}}-{\rm
div}\dfrac{\partial\breve{\Psi}_{\alpha}}{\partial\nabla\phi_{\alpha}},$
(5.13)

and
$\breve{g}_{\alpha}=\rho_{\alpha}^{-1}(\breve{\mu}_{\alpha}+\breve{\lambda})$.

Inserting the class (5.12) into the proposed modeling framework, we find
$\breve{\mu}_{\alpha}=\hat{\mu}_{\alpha}$. Additionally, we identify
$\breve{\lambda}=\lambda$; consequently $\breve{g}_{\alpha}=g_{\alpha}$. In
contrast to the unified modeling framework presented in the current paper,
this model is comprised of $N$ mass balance equations, and $N$ momentum
balance equations. As such, we compare the $N$ mass balance laws, and the
single mixture momentum balance law of the models. Starting with the mass
balance laws, (5.11a) can be written as:

$\displaystyle\partial_{t}\phi_{\alpha}+{\rm
div}(\phi_{\alpha}\mathbf{v})+\rho_{\alpha}^{-1}{\rm
div}\mathbf{J}_{\alpha}+\rho_{\alpha}^{-1}\breve{\gamma}_{\alpha}$
$\displaystyle=~{}0,$ (5.14a)
$\displaystyle\breve{\gamma}_{\alpha}-\sum_{\beta}\breve{m}_{\alpha{\beta}}(\breve{g}_{\alpha}-\breve{g}_{\beta})$
$\displaystyle=~{}0.$ (5.14b)

This form is very similar to (4.7b); the key difference is that the peculiar
velocity $\mathbf{J}_{\alpha}$ governed by a constitutive model
$\mathbf{J}_{\alpha}=\hat{\mathbf{J}}_{\alpha}$ in the current paper, whereas
in (5.11) it follows from the constitutive velocities. With the identification
$\breve{m}_{\alpha{\beta}}=m_{\alpha{\beta}}$ the mass transfer terms match
(except for the difference $\hat{\mathbf{j}}_{\alpha}=0$ in (4.7b)). Focusing
on the momentum balance laws, addition of (5.11b) provides:

$\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)+\sum\phi_{\alpha}\nabla\left(\breve{\lambda}+\breve{\mu}_{\alpha}\right)$
$\displaystyle-{\rm
div}\left(\sum_{\alpha}\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{v}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{v}\right)\right)-\rho\mathbf{b}$ $\displaystyle-{\rm
div}\left(\sum_{\alpha}\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{w}_{\alpha}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{w}_{\alpha}\right)-\sum_{\alpha}\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha}\right)$
$\displaystyle~{}=0,$ (5.15)

where we have adopted the identities:

$\displaystyle\sum_{\alpha}\tilde{\rho}_{\alpha}\mathbf{v}_{\alpha}\otimes\mathbf{v}_{\alpha}=$
$\displaystyle~{}\rho\mathbf{v}\otimes\mathbf{v}+\sum_{\alpha}\tilde{\rho}_{\alpha}\mathbf{w}_{\alpha}\otimes\mathbf{w}_{\alpha},$
(5.16a)
$\displaystyle\sum_{\alpha}\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{v}_{\alpha}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{v}_{\alpha}\right)=$
$\displaystyle~{}\sum_{\alpha}\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{v}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{v}\right)$
$\displaystyle~{}+\sum_{\alpha}\breve{\nu}_{\alpha}\left(2\nabla^{s}\mathbf{w}_{\alpha}+\breve{\bar{\lambda}}_{\alpha}{\rm
div}\mathbf{w}_{\alpha}\right).$ (5.16b)

With the identifications $\nu=\sum_{\alpha}\breve{\nu}_{\alpha}$ and
$\bar{\lambda}=\breve{\bar{\lambda}}_{\alpha}$ the first two lines match the
momentum equation (4.7a). The last line in (5.3) consists of terms that we
absent in (4.7a). This is a direct consequence of energy-dissipation law (3.3)
and the introduction of the model
$\mathbf{J}_{\alpha}=\hat{\mathbf{J}}_{\alpha}$ in (3.30b). In the case of
matching constitutive velocities, as described in Section 4.2, these terms
vanish.

##  6 Conclusion and outlook

This paper presents a unified framework for $N$-component Navier-Stokes Cahn-
Hilliard Allen-Cahn mixture models with non-matching densities. The framework
finds its roots in continuum mixture theory, which serves as a fundamental
guiding principle for designing multi-physics models at large. The unified
framework proposes a (phase-field) system of $N$ mass balance laws, and $1$
momentum balance law, that is invariant to the set of fundamental variables,
has an energy-dissipative structure, is reduction-consistent, symmetric with
respect to the numbering of the phases, and provides well-defined equilibrium
solutions. More specifically, we draw the following conclusions:

  * •

The form of the balance laws is invariant to the set of fundamental variables;
at both the constituent and mixture levels (Sections 2.2 and 2.3).

  * •

The free energy class depends on all volume fractions (and their gradients)
(Sections 3.1 and 3.2); this provides symmetry with respect to the numbering
of the constituents.

  * •

Chemical potentials are tightly connected to the Lagrange multiplier that
enforces volume conservation; these quantities occur only as superposition
(Section 3.2).

  * •

The unified framework is invariant to the set of independent variables; both
before and after constitutive modeling (Sections 3.3 and 3.4).

  * •

Constitutive quantities are such that the resulting model exhibits energy-
dissipation (Section 3.4).

  * •

Consistency with the single-phase equations requires mobility quantities to be
degenerate (Section 3.4).

  * •

Equilibrium solutions are determined by the balance of (generalized) chemical
potentials (see Section 4.3).

  * •

The framework reduces to the binary case for a particular form of the free
energy (see Section 5.1), is closely connected to a class-II model (see
Section 5.3). The model of [28] does not fit into the framework (see Section
5.2).

While the proposed unified framework offers insight into the modeling of
$N$-phase flows, we do not claim that it is complete. Therefore, we delineate
potential future research directions. First, it is important to study the
implications of the particular form of the free energy model, such as
equilibrium characteristics. To this purpose, we acknowledge the existence of
numerous $N$-phase free energy closure models (see e.g. [37]. Second, it is
essential to investigate the sharp interface asymptotic behavior (e.g. jump
conditions at interfaces) of the $N$-phase model. The last point concerns the
design of (property-preserving) numerical schemes. While the proposed system
is fully symmetric with respect to the set of variables, in the numerical
solution there are at least two roads one can pursue; (i) work with $N-1$
volume fractions and compute the $N$-th volume fraction from the others, or
(ii) work with the full set of volume fractions and enforce the saturation
constraint.

##  Appendix A Reduced free energy class, and proofs

We briefly discuss the free energy class with reduced dependency:

$\displaystyle\Psi=$
$\displaystyle~{}\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}\left(\left\\{\phi_{\alpha}\right\\}_{\alpha\neq{\beta}},\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha\neq{\beta}}\right),$
(A.1)

where the constituent number $\beta\in\left\\{1,...,N\right\\}$ is fixed, and
where both $\left\\{\phi_{\alpha}\right\\}_{\alpha\neq{\beta}}$ and
$\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha\neq{\beta}}$ consist of
independent variables. The class (A.1) is connected to (3.5) via the
identification:

$\displaystyle\Psi=$
$\displaystyle~{}\hat{\Psi}\left(\left\\{\phi_{\alpha}\right\\}_{\alpha\neq{\beta}},1-\sum_{\alpha\neq{\beta}}\phi_{\alpha},\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha\neq{\beta}},-\sum_{\alpha\neq{\beta}}\nabla\phi_{\alpha}\right)$
$\displaystyle=$
$\displaystyle~{}\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}\left(\left\\{\phi_{\alpha}\right\\}_{\alpha\neq{\beta}},\left\\{\nabla\phi_{\alpha}\right\\}_{\alpha\neq{\beta}}\right).$
(A.2)

The associated chemical potentials take the form:

$\displaystyle\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{\alpha}^{(\beta)}=\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\phi_{\alpha}}-{\rm
div}\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\nabla\phi_{\alpha}}.$
(A.3)

######  Lemma A.1 (Chemical potentials reduced class).

The chemical potentials of the reduced class may be expressed as

$\displaystyle\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{\alpha}^{(\beta)}=\hat{\mu}_{\alpha}-\hat{\mu}_{\beta}.$
(A.4)

###### Proof.

Direct evaluation of the partial derivatives provides:

$\displaystyle\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\phi_{\alpha}}=$
$\displaystyle~{}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}+\dfrac{\partial\hat{\Psi}}{\partial\phi_{\beta}}\dfrac{\partial\phi_{\beta}}{\partial\phi_{\alpha}}=\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}-\dfrac{\partial\hat{\Psi}}{\partial\phi_{\beta}},$
(A.5a)
$\displaystyle\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\nabla\phi_{\alpha}}=$
$\displaystyle~{}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\beta}}\dfrac{\partial\phi_{\beta}}{\partial\phi_{\alpha}}=\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\beta}}.$
(A.5b) The linearity of the divergence operator concludes the proof.

∎

######  Lemma A.2 (Derivative of the free energy).

The derivative of the free energy class (3.5) is given by:

$\displaystyle{\rm
d}\hat{\Psi}=\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha}),$ (A.6)

where ${\rm d}$ is the derivative operator.

###### Proof.

Inserting (A.5), the derivative of $\Psi$ takes the form

$\displaystyle{\rm d}\hat{\Psi}={\rm
d}\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{(\beta)}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha})$ $\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\left(\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}-\dfrac{\partial\hat{\Psi}}{\partial\phi_{\beta}}\right){\rm
d}\phi_{\alpha}+\sum_{\alpha\neq{\beta}}\left(\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\beta}}\right){\rm
d}(\nabla\phi_{\alpha})$ $\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha})$
$\displaystyle~{}-\dfrac{\partial\hat{\Psi}}{\partial\phi_{\beta}}\sum_{\alpha\neq{\beta}}{\rm
d}\phi_{\alpha}-\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\beta}}\sum_{\alpha\neq{\beta}}{\rm
d}(\nabla\phi_{\alpha})$ $\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha\neq{\beta}}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha})+\dfrac{\partial\hat{\Psi}}{\partial\phi_{\beta}}{\rm
d}\phi_{\beta}+\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\beta}}{\rm
d}(\nabla\phi_{\beta})$ $\displaystyle=$
$\displaystyle~{}\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}{\rm
d}\phi_{\alpha}+\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}{\rm
d}(\nabla\phi_{\alpha}),$ (A.7)

where we have invoked $\sum_{\alpha}{\rm d}\phi_{\alpha}=0$ and
$\sum_{\alpha}{\rm d}(\nabla\phi_{\alpha})=0$. The latter expression matches
the unconstrained derivative.

∎

######  Lemma A.3 (Well-defined free energy terms).

The following free energy terms in (3.2) are well-defined:

$\displaystyle\sum_{\alpha}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha};\quad\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}};\quad\sum_{\alpha}\dot{\phi}_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}};$
(A.8)

###### Proof.

We show that the first term is well-defined; the others follow similarly.
Utilizing an argumentation analogously to that of the proof of Lemma A.2, we
have the sequence of identities:

$\displaystyle\sum_{\alpha}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha}=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha}+\hat{\mu}_{\beta}\dot{\phi}_{\beta}$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\hat{\mu}_{\alpha}\dot{\phi}_{\alpha}-\hat{\mu}_{\beta}\sum_{\alpha\neq{\beta}}\dot{\phi}_{\alpha}$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\left(\hat{\mu}_{\alpha}-\hat{\mu}_{\beta}\right)\dot{\phi}_{\alpha}$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha\neq{\beta}}\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{\alpha}^{(\beta)}\dot{\phi}_{\alpha},$
(A.9)

where we have utilized Lemma A.1 in the last identity, where we have invoked
$\sum_{\alpha}\phi_{\alpha}=1$. Since the latter expression is well-defined,
so is the initial one. ∎

######  Lemma A.4 (Free energy identity).

The following identity holds:

$\displaystyle\sum_{\alpha}\phi_{\alpha}\nabla\mu_{\alpha}=\nabla\left(\sum_{\alpha}\phi_{\alpha}\hat{\mu}_{\alpha}-\hat{\Psi}\right)+{\rm
div}\left(\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right).$
(A.10)

###### Proof.

Expanding the derivatives of the right-hand size term yields:

$\displaystyle\nabla\left(\sum_{\alpha}\phi_{\alpha}\mu_{\alpha}-\hat{\Psi}\right)+{\rm
div}\left(\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)=$
$\displaystyle\sum_{\alpha}\phi_{\alpha}\nabla\mu_{\alpha}+\sum_{\alpha}\nabla\phi_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}-\sum_{\alpha}\nabla\phi_{\alpha}{\rm
div}\left(\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)-\nabla\hat{\Psi}$
$\displaystyle+\sum_{\alpha}\nabla\phi_{\alpha}{\rm
div}\left(\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}\right)+\sum_{\alpha}\left(\mathbf{H}\phi_{\alpha}\right)\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}=$
$\displaystyle\sum_{\alpha}\phi_{\alpha}\nabla\mu_{\alpha}-\nabla\hat{\Psi}+\sum_{\alpha}\nabla\phi_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\phi_{\alpha}}+\sum_{\alpha}\left(\mathbf{H}\phi_{\alpha}\right)\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}},$
(A.11)

where $\mathbf{H}\phi_{\alpha}$ is the hessian of $\phi_{\alpha}$. Observing
that the latter three terms in the final expression in (A) vanish completes
the proof. ∎

##  Appendix B Equivalence of modeling restrictions

This section discusses the approach of demonstrating equivalence of the (3.2)
and (3.3) via variable transformation.

First, we recall the variable transformation (2.7)

$\displaystyle\phi_{\alpha}=$
$\displaystyle~{}\frac{c_{\alpha}}{\rho_{\alpha}}\left(\sum_{\beta}\frac{c_{\beta}}{\rho_{\beta}}\right)^{-1}=\frac{c_{\alpha}}{\rho_{\alpha}}\rho,$
(B.1a) $\displaystyle c_{\alpha}=$
$\displaystyle~{}\rho_{\alpha}\phi_{\alpha}\left(\sum_{\beta}\rho_{\beta}\phi_{\beta}\right)^{-1}=\rho_{\alpha}\phi_{\alpha}\rho^{-1},$
(B.1b)

where we note:

$\displaystyle\rho=\hat{\rho}\left(\left\\{c_{\beta}\right\\}\right)$
$\displaystyle~{}\left(\sum_{\beta}\frac{c_{\beta}}{\rho_{\beta}}\right)^{-1},$
(B.2a)
$\displaystyle\rho^{-1}=\check{\rho}^{-1}\left(\left\\{\phi_{\beta}\right\\}\right)$
$\displaystyle~{}=\left(\sum_{\beta}\rho_{\beta}\phi_{\beta}\right)^{-1},$
(B.2b)

for $\alpha=1,...,N$.

######  Lemma B.1 (Transformation maps are not invertible).

The maps (B.1) are not invertible.

###### Proof.

A straightforward evaluation provides the elements of the Jacobian mappings:

$\displaystyle\dfrac{\partial\phi_{\beta}}{\partial c_{\alpha}}=$
$\displaystyle~{}\rho\rho_{\beta}^{-1}\left(\delta_{\alpha{\beta}}-\rho
c_{\beta}\rho_{\alpha}^{-1}\right),$ (B.3a) $\displaystyle\dfrac{\partial
c_{\beta}}{\partial\phi_{\alpha}}=$
$\displaystyle~{}\rho^{-1}\rho_{\beta}\left(\delta_{\alpha{\beta}}-\rho^{-1}\phi_{\beta}\rho_{\alpha}\right),$
(B.3b)

where $\delta_{\alpha{\beta}}$ is the kroneckerdelta. Summation over
${\beta}=1,...,N$ yields:

$\displaystyle\sum_{\beta}\dfrac{\partial\phi_{\beta}}{\partial c_{\alpha}}=$
$\displaystyle~{}0,$ (B.4a) $\displaystyle\sum_{\beta}\dfrac{\partial
c_{\beta}}{\partial\phi_{\alpha}}=$ $\displaystyle~{}0.$ (B.4b)

Hence, each of the columns of the Jacobian sums to zero. Thus the columns are
linearly dependent, and consequently the determinants of the both mappings
vanish:

$\displaystyle{\rm det}\dfrac{\partial\phi_{\beta}}{\partial c_{\alpha}}=$
$\displaystyle~{}0,$ (B.5a) $\displaystyle{\rm det}\dfrac{\partial
c_{\beta}}{\partial\phi_{\alpha}}=$ $\displaystyle~{}0.$ (B.5b)

∎

Even though the transformation maps are not invertible, we can still apply the
variable transformation. We recall the chain rule for the chemical potential.

######  Lemma B.2 (Chain rule chemical potentials).

Provided that the free energies are specified (well-defined), we have the
chain rule for chemical potentials:

$\displaystyle\check{\mu}_{\alpha}=$
$\displaystyle~{}\sum_{\beta}\hat{\mu}_{\beta}\dfrac{\partial\phi_{\beta}}{\partial
c_{\alpha}},$ (B.6a) $\displaystyle\hat{\mu}_{\alpha}=$
$\displaystyle~{}\sum_{\beta}\check{\mu}_{\beta}\dfrac{\partial
c_{\beta}}{\partial\phi_{\alpha}}.$ (B.6b)

###### Proof.

We show (B.6a) and note that (B.6b) follows similarly. A direct computation
yields:

$\displaystyle\check{\mu}_{\alpha}=$
$\displaystyle~{}\dfrac{\partial\hat{\Psi}\left(\phi_{{\beta}}\left(\left\\{c_{\gamma}\right\\}\right),\sum_{\gamma}\frac{\partial\phi_{\beta}}{\partial
c_{\gamma}}\nabla c_{\gamma}\right)}{\partial c_{\alpha}}-{\rm
div}\dfrac{\partial\hat{\Psi}\left(\phi_{{\beta}}\left(\left\\{c_{\gamma}\right\\}\right),\sum_{\gamma}\frac{\partial\phi_{\beta}}{\partial
c_{\gamma}}\nabla c_{\gamma}\right)}{\partial\nabla c_{\alpha}}$
$\displaystyle=$
$\displaystyle~{}\sum_{\beta}\dfrac{\partial\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)}{\partial\phi_{\beta}}\dfrac{\partial\phi_{\beta}}{\partial
c_{\alpha}}$
$\displaystyle~{}+\sum_{\beta}\dfrac{\partial\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)}{\partial\nabla\phi_{\beta}}\cdot\left(\sum_{\gamma}\nabla
c_{\gamma}~{}\frac{\partial^{2}\phi_{\beta}}{\partial c_{\alpha}\partial
c_{\gamma}}\right)$ $\displaystyle~{}-\sum_{\beta}{\rm
div}\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)\frac{\partial\phi_{\beta}}{\partial
c_{\alpha}}$
$\displaystyle~{}-\sum_{\beta}\dfrac{\partial\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)}{\partial\nabla\phi_{\beta}}\cdot\nabla\left(\frac{\partial\phi_{\beta}}{\partial
c_{\alpha}}\right)$ $\displaystyle=$
$\displaystyle~{}\sum_{\beta}\left(\dfrac{\partial\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)}{\partial\phi_{\beta}}-{\rm
div}\left(\dfrac{\partial\hat{\Psi}\left(\left\\{\phi_{{\beta}}\right\\},\left\\{\nabla\phi_{\beta}\right\\}\right)}{\partial\nabla\phi_{\beta}}\right)\right)\dfrac{\partial\phi_{\beta}}{\partial
c_{\alpha}}$ $\displaystyle=$
$\displaystyle~{}\sum_{\beta}\hat{\mu}_{\beta}\dfrac{\partial\phi_{\beta}}{\partial
c_{\alpha}}.$ (B.7)

∎

######  Lemma B.3 (Relations between chemical quantities).

The chemical potential quantities are related via the following identities:

$\displaystyle\hat{\mu}_{\alpha}=$
$\displaystyle~{}\rho^{-1}\rho_{\alpha}\left(\check{\mu}_{\alpha}-\sum_{\beta}\check{\mu}_{\beta}c_{\beta}\right),$
(B.8a) $\displaystyle\check{\mu}_{\alpha}=$
$\displaystyle~{}\rho\rho_{\alpha}^{-1}\left(\hat{\mu}_{\alpha}-\sum_{\beta}\hat{\mu}_{\beta}\phi_{\beta}\right).$
(B.8b)

###### Proof.

This follows from substituting (B.3) into (B.2). ∎

######  Lemma B.4 (Matching Korteweg tensors).

Provided that the free energies are specified (well-defined), the Korteweg
stress tensors of the both modeling choices are identical:

$\displaystyle\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}=\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla c_{\alpha}}.$
(B.9)

###### Proof.

This follows from (B.3) and (2.5b):

$\displaystyle\sum_{\beta}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla c_{\alpha}}=$
$\displaystyle~{}\sum_{\alpha}\left(\sum_{\beta}\dfrac{\partial
c_{\alpha}}{\partial\phi_{\beta}}\nabla\phi_{\beta}\right)\otimes\left(\sum_{\gamma}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\gamma}}\dfrac{\partial\phi_{\gamma}}{\partial
c_{\alpha}}\right)$ $\displaystyle=$
$\displaystyle~{}\sum_{\alpha}\left(\sum_{\beta}\rho^{-1}\rho_{\alpha}\left(\delta_{\alpha{\beta}}-\rho^{-1}\rho_{\beta}\phi_{\alpha}\right)\nabla\phi_{\beta}\right)$
$\displaystyle\quad\otimes\left(\sum_{\gamma}\rho\rho_{\gamma}^{-1}\left(\delta_{\alpha{\gamma}}-\rho\rho_{\alpha}^{-1}c_{\gamma}\right)\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\gamma}}\right)$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha}\left(\rho^{-1}\rho_{\beta}\nabla\phi_{\alpha}-\rho^{-2}\rho_{\alpha}\phi_{\alpha}\nabla\rho\right)$
$\displaystyle\quad\otimes\left(\rho\rho_{\alpha}^{-1}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\rho\rho_{\alpha}^{-1}\sum_{\gamma}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\gamma}}\phi_{\gamma}\right)$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}-\rho^{-1}\nabla\rho\otimes\sum_{\alpha}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}$
$\displaystyle~{}-\sum_{\alpha}\nabla\phi_{\alpha}\otimes\sum_{\gamma}\phi_{\gamma}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\gamma}}+\rho^{-1}\nabla\rho\otimes\sum_{\gamma}\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\gamma}}$
$\displaystyle=$
$\displaystyle~{}\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}.$
(B.10)

∎

######  Theorem B.5 (Equivalence modeling restrictions).

Provided that the free energies are specified (well-defined), the modeling
restrictions (3.2) and (3.3) are equivalent.

###### Proof.

We select the following relations between the Lagrange multipliers of the two
modeling choices:

$\displaystyle\check{\lambda}=$
$\displaystyle~{}\hat{\lambda}+\sum_{\beta}\hat{\mu}_{\beta}\phi_{\beta}.$
(B.11)

Invoking Lemma B.4 and substituting the relation (B.11) provides:

$\displaystyle\check{\lambda}\mathbf{I}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}=\hat{\lambda}\mathbf{I}+\sum_{\alpha}\nabla\phi_{\alpha}\otimes\dfrac{\partial\hat{\Psi}}{\partial\nabla\phi_{\alpha}}+\left(\hat{\mu}_{\alpha}\phi_{\alpha}-\check{\Psi}\right)\mathbf{I}.$
(B.12)

In a similar fashion we find

$\displaystyle\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}=$
$\displaystyle~{}\rho_{\alpha}^{-1}\hat{\lambda}+\rho_{\alpha}^{-1}\left(\hat{\mu}_{\alpha}-\sum_{\beta}\hat{\mu}_{\beta}\phi_{\beta}\right)+\rho_{\alpha}^{-1}\sum_{\alpha}\hat{\mu}_{\beta}\phi_{\beta}=g_{\alpha},$
(B.13)

and conclude:

$\displaystyle-\sum_{\alpha}\nabla\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\cdot\mathbf{J}_{\alpha}=$
$\displaystyle~{}-\sum_{\alpha}\nabla g_{\alpha}\cdot\mathbf{J}_{\alpha}$
(B.14a)
$\displaystyle-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\gamma_{\alpha}=$
$\displaystyle~{}-\sum_{\alpha}g_{\alpha}\gamma_{\alpha}.$ (B.14b)

∎

Finally, we note that the non-invertibility of the transformation maps (B.1)
and the arbitrariness of the chemical potentials (Section 3) can lead to
cumbersome properties. For example, we can deduce

$\displaystyle\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}=$
$\displaystyle~{}0,$ (B.15a)
$\displaystyle\sum_{\alpha}\check{\mu}_{\alpha}c_{\alpha}=$
$\displaystyle~{}0,$ (B.15b)

from (B.8). In particular, it follows from

$\displaystyle\sum_{\alpha}\hat{\mu}_{\alpha}\phi_{\alpha}=$
$\displaystyle~{}\sum_{\alpha}\left(\sum_{\beta}\check{\mu}_{\beta}\frac{\partial
c_{\beta}}{\partial\phi_{\alpha}}\right)\phi_{\alpha},$ (B.16a)
$\displaystyle\sum_{\alpha}\check{\mu}_{\alpha}c_{\alpha}=$
$\displaystyle~{}\sum_{\alpha}\left(\sum_{\beta}\hat{\mu}_{\beta}\frac{\partial\phi_{\beta}}{\partial
c_{\alpha}}\right)c_{\alpha},$ (B.16b)

alongside with the identities:

$\displaystyle\sum_{\alpha}\dfrac{\partial
c_{\beta}}{\partial\phi_{\alpha}}\phi_{\alpha}=$ $\displaystyle~{}0,$ (B.17a)
$\displaystyle\sum_{\alpha}\dfrac{\partial\phi_{\beta}}{\partial
c_{\alpha}}c_{\alpha}=$ $\displaystyle~{}0.$ (B.17b)

In this situation we have $\hat{\lambda}=\check{\lambda}$. We note that the
identities (B.15) do not hold in general; however there exist free energies
for which it applies. We recall (3.22)-(3.6).

##  Appendix C Alternative constitutive modeling

In this section we provide some brief details on the constitutive modeling
based on concentration variables. Section C.1 outlines the modeling
assumptions, and Section C.2 derives the modeling restriction.

###  C.1 Assumptions and modeling choices

We use the balance laws (3.1), where the mass balance laws are now written in
terms of concentration variables:

$\displaystyle\rho(\partial_{t}c_{\alpha}+\mathbf{v}\cdot\nabla
c_{\alpha})+{\rm div}\mathbf{J}_{\alpha}$ $\displaystyle=~{}\gamma_{\alpha},$
(C.1a) $\displaystyle\partial_{t}(\rho\mathbf{v})+{\rm
div}\left(\rho\mathbf{v}\otimes\mathbf{v}\right)-{\rm
div}\mathbf{T}-\rho\mathbf{b}$ $\displaystyle=~{}0,$ (C.1b)
$\displaystyle\mathbf{T}-\mathbf{T}^{T}$ $\displaystyle=~{}0,$ (C.1c)

where (C.1a) holds for constituents $\alpha=1,...,N$. We use the energy-
dissipation law (3.3):

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=\mathscr{W}-\mathscr{D},$
(C.2)

with dissipation $\mathscr{D}\geq 0$, and recall (3.4). We postulate the free
energy to pertain to the constitutive class:

$\displaystyle\Psi=$
$\displaystyle~{}\check{\Psi}\left(\left\\{c_{\alpha}\right\\}_{\alpha=1,...,N},\left\\{\nabla
c_{\alpha}\right\\}_{\alpha=1,...,N}\right),$ (C.3)

subject to the summation constraint (2.5a), and introduce the chemical
potential quantities ($\alpha=1,...,N$):

$\displaystyle\check{\mu}_{\alpha}=$
$\displaystyle~{}\dfrac{\partial\check{\Psi}}{\partial c_{\alpha}}-{\rm
div}\dfrac{\partial\check{\Psi}}{\partial\nabla c_{\alpha}}.$ (C.4)

###  C.2 Modeling restriction

By applying Reynolds transport theorem, the divergence theorem and integration
by parts, and identity (3.11), the evolution of the free energy $\check{\Psi}$
takes the form:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\check{\Psi}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\check{\Psi}~{}{\rm
div}\mathbf{v}+\sum_{\alpha}\check{\mu}_{\alpha}\dot{c}_{\alpha}-\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}:\nabla\mathbf{v}~{}{\rm d}v$
$\displaystyle~{}+\displaystyle\int_{\partial\mathcal{R}(t)}\sum_{\alpha}\dot{c}_{\alpha}\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}\cdot\boldsymbol{\nu}~{}{\rm d}a.$ (C.5)

Substituting the constituent mass balance laws (3.1a), and again applying
integration by parts, yields:

$\displaystyle\dfrac{{\rm d}}{{\rm
d}t}\displaystyle\int_{\mathcal{R}(t)}\sum_{\alpha}\check{\Psi}~{}{\rm d}v=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\check{\Psi}~{}{\rm
div}\mathbf{v}+\sum_{\alpha}\nabla(\rho^{-1}\check{\mu}_{\alpha})\cdot\mathbf{J}_{\alpha}$
$\displaystyle~{}~{}~{}-\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}:\nabla\mathbf{v}+\sum_{\alpha}\rho^{-1}\check{\mu}_{\alpha}\gamma_{\alpha}~{}{\rm
d}v$
$\displaystyle~{}+\displaystyle\int_{\partial\mathcal{R}(t)}\left(\sum_{\alpha}\dot{c}_{\alpha}\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\rho^{-1}\check{\mu}_{\alpha}\mathbf{J}_{\alpha}\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a.$ (C.6)

Addition of (3.15) and (C.2) provides the evolution of the total energy:

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}\mathbf{J}_{\alpha}-\dot{c}_{\alpha}\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}\right)\right)\cdot\boldsymbol{\nu}~{}{\rm d}a$
$\displaystyle~{}-\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla(\rho^{-1}\check{\mu}_{\alpha})\cdot\mathbf{J}_{\alpha}-\rho^{-1}\check{\mu}_{\alpha}\gamma_{\alpha}\right)~{}{\rm
d}v.$ (C.7)

Analogously to Section 3.2, we restore the degenerate nature of (C.2) via a
Lagrange multiplier construction:

$\displaystyle 0=$ $\displaystyle~{}\check{\lambda}{\rm
div}\mathbf{v}+\nabla\left(\check{\lambda}\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}\right)-\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}\cdot\nabla\check{\lambda}-\check{\lambda}\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\gamma_{\alpha},$
(C.8)

where $\check{\lambda}$ is the scalar Lagrange multiplier. Integrating (C.8)
over $\mathcal{R}(t)$ and subtracting the result from (C.2) provides:

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\mathscr{E}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}\mathbf{J}_{\alpha}-\dot{c}_{\alpha}\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}\right)-\check{\lambda}\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a$
$\displaystyle~{}-\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\check{\lambda}\mathbf{I}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\cdot\mathbf{J}_{\alpha}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.-\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\gamma_{\alpha}\right)~{}{\rm
d}v.$ (C.9)

The rate of work and the dissipation take the forms:

$\displaystyle\mathscr{W}=$
$\displaystyle~{}\displaystyle\int_{\partial\mathcal{R}(t)}\left(\mathbf{v}^{T}\mathbf{T}-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}\mathbf{J}_{\alpha}-\dot{c}_{\alpha}\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}\right)-\check{\lambda}\displaystyle\sum_{\alpha}\rho_{\alpha}^{-1}\mathbf{J}_{\alpha}\right)\cdot\boldsymbol{\nu}~{}{\rm
d}a,$ (C.10a) $\displaystyle\mathscr{D}=$
$\displaystyle~{}\displaystyle\int_{\mathcal{R}(t)}\left(\mathbf{T}+\check{\lambda}\mathbf{I}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{\alpha}\left(-\nabla\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\cdot\mathbf{J}_{\alpha}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.-\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\gamma_{\alpha}\right)~{}{\rm
d}v.$ (C.10b)

Given that the control volume $\mathcal{R}=\mathcal{R}(t)$, can be chosen
arbitrarily, adhering to the energy dissipation law requires that the
following local inequality is satisfied:

$\displaystyle\left(\mathbf{T}+\check{\lambda}\mathbf{I}+\sum_{\alpha}\nabla
c_{\alpha}\otimes\dfrac{\partial\check{\Psi}}{\partial\nabla
c_{\alpha}}-\check{\Psi}\mathbf{I}\right):\nabla\mathbf{v}$
$\displaystyle-\sum_{\alpha}\nabla\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\cdot\mathbf{J}_{\alpha}-\sum_{\alpha}\left(\rho^{-1}\check{\mu}_{\alpha}+\rho_{\alpha}^{-1}\check{\lambda}\right)\gamma_{\alpha}$
$\displaystyle\geq 0.$ (C.11)

## References

  * Hirt and Nichols [1981] C.W. Hirt and B.D. Nichols.  Volume of fluid (VOF) Method for the Dynamics of Free Boundaries.  _Journal of Computational Physics_ , 39:201–225, 1981.
  * Sethian [2001] J.A. Sethian.  Evolution, implementation, and application of level set and fast marching methods for advancing fronts.  _Journal of computational physics_ , 169(2):503–555, 2001.
  * ten Eikelder and Akkerman [2021] M.F.P. ten Eikelder and I. Akkerman.  A novel diffuse-interface model and a fully-discrete maximum-principle-preserving energy-stable method for two-phase flow with surface tension and non-matching densities.  _Computer Methods in Applied Mechanics and Engineering_ , 379:113751, 2021.
  * Bothe [2022] D. Bothe.  Sharp-interface continuum thermodynamics of multicomponent fluid systems with interfacial mass.  _International Journal of Engineering Science_ , 179:103731, 2022.
  * Cahn and Hilliard [1958] J.W. Cahn and J.E. Hilliard.  Free energy of a nonuniform system. I. Interfacial free energy.  _The Journal of chemical physics_ , 28:258–267, 1958.
  * Cahn [1959] J.W. Cahn.  Free energy of a nonuniform system. II. Thermodynamic basis.  _The Journal of chemical physics_ , 30:1121–1124, 1959.
  * Anderson et al. [1998] D.M. Anderson, G.B. McFadden, and A.A. Wheeler.  Diffuse-interface methods in fluid mechanics.  _Annual review of fluid mechanics_ , 30:139–165, 1998.
  * Gomez and van der Zee [2018] H. Gomez and K.G. van der Zee.  Computational phase-field modeling.  _Encyclopedia of Computational Mechanics Second Edition_ , pages 1–35, 2018.
  * Steinbach [2009] I. Steinbach.  Phase-field models in materials science.  _Modelling and simulation in materials science and engineering_ , 17:073001, 2009.
  * Boettinger et al. [2002] W.J. Boettinger, J.A. Warren, C. Beckermann, and A. Karma.  Phase-field simulation of solidification.  _Annual review of materials research_ , 32:163–194, 2002.
  * Ambati et al. [2015] M. Ambati, T. Gerasimov, and L. De Lorenzis.  Phase-field modeling of ductile fracture.  _Computational Mechanics_ , 55:1017–1040, 2015.
  * Yue et al. [2004] P. Yue, J.J. Feng, C. Liu, and J. Shen.  A diffuse-interface method for simulating two-phase flows of complex fluids.  _Journal of Fluid Mechanics_ , 515:293–317, 2004.
  * Hohenberg and Halperin [1977] P.C. Hohenberg and B.I. Halperin.  Theory of dynamic critical phenomena.  _Reviews of Modern Physics_ , 49:435, 1977.
  * Gurtin [1996] M.E. Gurtin.  Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance.  _Physica D: Nonlinear Phenomena_ , 92:178–192, 1996.
  * Gurtin et al. [1996] M.E. Gurtin, D. Polignone, and Vinals J.  Two-phase binary fluids and immiscible fluids described by an order parameter.  _Mathematical Models and Methods in Applied Sciences_ , 6:815–831, 1996.
  * Lowengrub and Truskinovsky [1998] J. Lowengrub and L. Truskinovsky.  Quasi–incompressible Cahn–Hilliard fluids and topological transitions.  _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , 454:2617–2654, 1998.
  * Abels et al. [2012] H. Abels, H. Garcke, and G. Grün.  Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities.  _Mathematical Models and Methods in Applied Sciences_ , 22:1150013, 2012.
  * Aki et al. [2014] G.L. Aki, W. Dreyer, J. Giesselmann, and C. Kraus.  A quasi-incompressible diffuse interface model with phase transition.  _Mathematical Models and Methods in Applied Sciences_ , 24:827–861, 2014.
  * Kay and Welford [2007] D. Kay and R. Welford.  Efficient numerical solution of Cahn–Hilliard–Navier–Stokes fluids in 2D.  _SIAM Journal on Scientific Computing_ , 29(6):2241–2257, 2007.
  * Guo et al. [2014] Z. Guo, P. Lin, and J.S. Lowengrub.  A numerical method for the quasi-incompressible Cahn–Hilliard–Navier–Stokes equations for variable density flows with a discrete energy law.  _Journal of Computational Physics_ , 276:486–507, 2014.
  * Khanwale et al. [2022] M.A. Khanwale, K. Saurabh, M. Fernando, V.M. Calo, H. Sundar, J.A. Rossmanith, and B. Ganapathysubramanian.  A fully-coupled framework for solving Cahn-Hilliard Navier-stokes equations: Second-order, energy-stable numerical methods on adaptive octree based meshes.  _Computer Physics Communications_ , 280:108501, 2022.
  * ten Eikelder and Schillinger [2024] M.F.P. ten Eikelder and D. Schillinger.  The divergence-free velocity formulation of the consistent Navier-Stokes Cahn-hilliard model with non-matching densities, divergence-conforming discretization, and benchmarks.  _Journal of Computational Physics_ , 513:113148, 2024.
  * ten Eikelder et al. [2023] M.F.P. ten Eikelder, K.G. van der Zee, I. Akkerman, and D. Schillinger.  A unified framework for Navier-Stokes Cahn-Hilliard models with non-matching densities.  _Mathematical Models and Methods in Applied Sciences_ , 33:175–221, 2023.
  * Kim and Lowengrub [2005] J. Kim and J. Lowengrub.  Phase field modeling and simulation of three-phase flows.  _Interfaces and free boundaries_ , 7:435–466, 2005.
  * Heida et al. [2012] M. Heida, J. Málek, and K.R. Rajagopal.  On the development and generalizations of Cahn–Hilliard equations within a thermodynamic framework.  _Zeitschrift für angewandte Mathematik und Physik_ , 63:145–169, 2012.
  * Li and Wang [2014] J. Li and Q. Wang.  A class of conservative phase field models for multiphase fluid flows.  _Journal of Applied Mechanics_ , 81:021004, 2014.
  * Dong [2015] S. Dong.  Physical formulation and numerical algorithm for simulating N immiscible incompressible fluids involving general order parameters.  _Journal of Computational Physics_ , 283:98–128, 2015.
  * Dong [2018] S. Dong.  Multiphase flows of N immiscible incompressible fluids: a reduction-consistent and thermodynamically-consistent formulation and associated algorithm.  _Journal of Computational Physics_ , 361:1–49, 2018.
  * Huang et al. [2021] Z. Huang, G. Lin, and A.M. Ardekani.  A consistent and conservative model and its scheme for n-phase-m-component incompressible flows.  _Journal of Computational Physics_ , 434:110229, 2021.
  * Boyer and Lapuerta [2006] F. Boyer and C. Lapuerta.  Study of a three component Cahn-Hilliard flow model.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 40:653–687, 2006.
  * Boyer et al. [2010] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar, and M. Quintard.  Cahn–Hilliard/Navier–Stokes model for the simulation of three-phase flows.  _Transport in Porous Media_ , 82:463–483, 2010.
  * Tóth et al. [2016] G.I. Tóth, M. Zarifi, and B. Kvamme.  Phase-field theory of multicomponent incompressible Cahn-Hilliard liquids.  _Physical Review E_ , 93:013126, 2016.
  * Zhang and Wang [2016] Q. Zhang and X.-P. Wang.  Phase field modeling and simulation of three-phase flow on solid surfaces.  _Journal of Computational Physics_ , 319:79–107, 2016.
  * Baňas and Nürnberg [2017] L. Baňas and R. Nürnberg.  Numerical approximation of a non-smooth phase-field model for multicomponent incompressible flow.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 51:1089–1117, 2017.