2.1. Preliminaries

In the continuum theory of mixtures, the material body $\mathscr{B}$ comprises $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 $\boldsymbol{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:

(2.1) \begin{align} \boldsymbol{x} := \boldsymbol{\chi }_{\alpha }(\boldsymbol{X}_{\alpha },t), \end{align}

where $\boldsymbol{x} \in \Omega$ , with $\Omega \in \mathbb{R}^d$ the domain (dimension $d$ ). We refer for more details on continuum mixture theory to Truesdell & Toupin (Reference Truesdell and Toupin1960), and sketch the situation in figure 2.

Figure 2.

Situation sketch continuum mixture theory.

We introduce the constituent partial mass density $\tilde {\rho }_{\alpha }$ and specific mass density $\rho _{\alpha }\gt 0$ , respectively, as

(2.2a) \begin{align} \tilde {\rho }_{\alpha }(\boldsymbol{x},t) :=&\, \displaystyle \lim _{ \vert V \vert \rightarrow 0} \frac {M_{\alpha }(V)}{\vert V \vert }, \end{align}

(2.2b) \begin{align} \rho _{\alpha }(\boldsymbol{x},t) :=&\, \displaystyle \lim _{\vert V_{\alpha }\vert \rightarrow 0} \frac {M_{\alpha }(V)}{\vert V_{\alpha }\vert }, \end{align}

where $V \subset \Omega$ (measure $\vert V \vert$ ) is an arbitrary control volume around $\boldsymbol{x}$ , $V_{\alpha } \subset V$ (measure $\vert V_{\alpha }\vert$ ) is the volume of constituent $\alpha$ so that $V =\cup _{\alpha }V_{\alpha }$ . 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:

(2.3) \begin{align} \rho (\boldsymbol{x},t):=&\, \displaystyle \lim _{ \vert V \vert \rightarrow 0} \frac {M(V)}{\vert V \vert }=\displaystyle \sum _{\alpha }\tilde {\rho }_\alpha (\boldsymbol{x},t). \end{align}

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

(2.4a) \begin{align} c_\alpha (\boldsymbol{x},t):=&\, \displaystyle \lim _{ \vert V \vert \rightarrow 0} \frac {M_{\alpha }(V)}{M(V)}=\frac {\tilde {\rho }_\alpha }{\rho }, \end{align}

(2.4b) \begin{align} \phi _\alpha (\boldsymbol{x},t):=&\, \displaystyle \lim _{ \vert V \vert \rightarrow 0} \frac {|V_{\alpha }|}{|V|}=\frac {\tilde {\rho }_\alpha }{\rho _\alpha }, \end{align}

which sum up to one:

(2.5a) \begin{align} \displaystyle \sum _{\alpha } c_{\alpha }(\boldsymbol{x},t)=&\, 1, \end{align}

(2.5b) \begin{align} \displaystyle \sum _{\alpha } \phi _{\alpha }(\boldsymbol{x},t)=&\, 1. \end{align}

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

(2.6) \begin{align} \rho _\alpha (\boldsymbol{x},t) = \rho _\alpha . \end{align}

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

(2.7a) \begin{align} \phi _\alpha =&\, \frac {c_\alpha }{\rho _\alpha }\left (\sum _{\beta }\frac {c_{\beta }}{\rho _{\beta }}\right )^{-1}, \end{align}

(2.7b) \begin{align} c_\alpha =&\, \rho _\alpha \phi _\alpha \left (\sum _{\beta }\rho _{\beta }\phi _{\beta }\right )^{-1}, \end{align}

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

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

(2.8) \begin{align} \grave {\psi }_\alpha =\partial _t\psi _{\alpha }(\boldsymbol{X}_{\alpha },t) \vert _{\boldsymbol{X}_\alpha }. \end{align}

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

(2.9) \begin{align} \boldsymbol{v}_{\alpha }(\boldsymbol{x},t)=\partial _t\boldsymbol{\chi }_{\alpha }(\boldsymbol{X}_{\alpha },t) \vert _{\boldsymbol{X}_\alpha } = \grave {\boldsymbol{\chi }}_\alpha . \end{align}

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

(2.10a) \begin{align} \boldsymbol{v}(\boldsymbol{x},t) =&\, \displaystyle \sum _{\alpha } c_\alpha (\boldsymbol{x},t) \boldsymbol{v}_\alpha (\boldsymbol{x},t), \end{align}

(2.10b) \begin{align} \boldsymbol{u}(\boldsymbol{x},t) =&\, \displaystyle \sum _{\alpha } \phi _\alpha (\boldsymbol{x},t) \boldsymbol{v}_\alpha (\boldsymbol{x},t). \end{align}

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

(2.11a) \begin{align} \boldsymbol{w}_{\alpha }(\boldsymbol{x},t):=&\,\boldsymbol{v}_{\alpha }(\boldsymbol{x},t)-\boldsymbol{v}(\boldsymbol{x},t), \end{align}

(2.11b) \begin{align} \boldsymbol{\omega }_{\alpha }(\boldsymbol{x},t):=&\,\boldsymbol{v}_{\alpha }(\boldsymbol{x},t)-\boldsymbol{u}(\boldsymbol{x},t). \end{align}

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

(2.12a) \begin{align} \boldsymbol{J}_\alpha :=&\, \tilde {\rho }_\alpha \boldsymbol{w}_\alpha , \end{align}

(2.12b) \begin{align} \boldsymbol{h}_\alpha = &\, \phi _\alpha \boldsymbol{w}_\alpha , \end{align}

(2.12c) \begin{align} \boldsymbol{J}_\alpha ^u = &\, \tilde {\rho }_\alpha \boldsymbol{\omega }_\alpha, \end{align}

(2.12d) \begin{align} \boldsymbol{h}_\alpha ^u = &\, \phi _\alpha \boldsymbol{\omega }_\alpha . \end{align}

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

(2.13a) \begin{align} \displaystyle \sum _{\alpha } \boldsymbol{J}_\alpha =&\, 0, \end{align}

(2.13b) \begin{align} \displaystyle \sum _{\alpha } \boldsymbol{h}_\alpha ^u =&\, 0. \end{align}

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

(2.14a) \begin{align} \boldsymbol{u} =&\, \boldsymbol{v} + \displaystyle \sum _\alpha \rho _\alpha ^{-1} \boldsymbol{J}_\alpha = \boldsymbol{v} + \sum _\alpha \boldsymbol{h}_\alpha , \end{align}

(2.14b) \begin{align} \boldsymbol{v} =&\, \boldsymbol{u} + \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u. \end{align}

Proof. These relations result from the following sequences of identities:

(2.15a) \begin{align} \boldsymbol{u} &= \sum _\alpha \phi _\alpha \boldsymbol{v}_\alpha = \sum _\alpha \phi _\alpha \boldsymbol{w}_\alpha + \sum _\alpha \phi _\alpha \boldsymbol{v} = \sum _\alpha \rho _\alpha ^{-1} \boldsymbol{J}_\alpha + \boldsymbol{v}, \end{align}

(2.15b) \begin{align} 0 &= \sum _\alpha \boldsymbol{J}_\alpha = \sum _\alpha \tilde {\rho }_\alpha (\boldsymbol{v}_\alpha -\boldsymbol{v}) = \sum _\alpha \tilde {\rho }_\alpha (\boldsymbol{v}_\alpha -\boldsymbol{u} + \boldsymbol{u}-\boldsymbol{v}) \nonumber \\ & = \sum _\alpha \boldsymbol{J}_\alpha ^u + \rho (\boldsymbol{u}-\boldsymbol{v}).\\ \nonumber \end{align}

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

(2.16a) \begin{align} \boldsymbol{J}_\alpha =&\, \boldsymbol{J}_\alpha ^u - c_\alpha \sum _{\beta } \boldsymbol{J}_{\beta }^u, \end{align}

(2.16b) \begin{align} \boldsymbol{J}_\alpha ^u =&\, \boldsymbol{J}_\alpha - \tilde {\rho }_\alpha \sum _{\beta } \rho _{\beta }^{-1} \boldsymbol{J}_{\beta }, \end{align}

(2.16c) \begin{align} \boldsymbol{h}_\alpha =&\, \boldsymbol{h}_\alpha ^u - \phi _\alpha \rho ^{-1} \sum _{\beta } \rho _{\beta }\boldsymbol{h}_{\beta }^u, \end{align}

(2.16d) \begin{align} \boldsymbol{h}_\alpha ^u =&\, \boldsymbol{h}_\alpha - \phi _\alpha \sum _{\beta } \boldsymbol{h}_{\beta }. \end{align}

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

(2.17) \begin{align} \dot {\psi }(\boldsymbol{x},t) =&\, \partial _t \psi (\boldsymbol{x},t) + \boldsymbol{v}(\boldsymbol{x},t)\ {\boldsymbol\cdot}\ \nabla \psi (\boldsymbol{x},t). \end{align}

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 $\boldsymbol{x}\in \Omega$ and $t \gt 0$ :

(2.18a) \begin{align} \partial _t \tilde {\rho }_\alpha + \textrm{div}(\tilde {\rho }_\alpha \boldsymbol{v}_\alpha ) &=\, \gamma _\alpha , \end{align}

(2.18b) \begin{align} \partial _t (\tilde {\rho }_\alpha \boldsymbol{v}_\alpha ) + \textrm{div} \left ( \tilde {\rho }_\alpha \boldsymbol{v}_\alpha \otimes \boldsymbol{v}_\alpha \right ) - \textrm{div} \boldsymbol{T}_\alpha - \tilde {\rho }_\alpha \boldsymbol{b}_\alpha &=\, \boldsymbol{\pi }_\alpha , \end{align}

(2.18c) \begin{align} \boldsymbol{T}_\alpha -\boldsymbol{T}_\alpha ^{T} &=\,\boldsymbol{N}_\alpha . \end{align}

Equations (2.18a ) describe the local constituent mass balance laws, where the interaction terms $\gamma _{\alpha }$ denote the mass supply of constituent $\alpha$ due to chemical reactions with the other constituents. Then, (2.18b ) represent the local constituent linear momentum balance laws, where $\boldsymbol{T}_\alpha$ is the Cauchy stress tensor of constituent $\alpha$ , $\boldsymbol{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 ( $\boldsymbol{b}_\alpha = \boldsymbol{b}$ for $\alpha = 1, \dots , N$ ) throughout the article. Additionally, we restrict to gravitational body forces: $\boldsymbol{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 $\boldsymbol{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

(2.19) \begin{align} \gamma _\alpha = \zeta _\alpha - \textrm{div} \boldsymbol{j}_\alpha . \end{align}

The mass balance laws (3.1a ) take the form

(2.20) \begin{align} \partial _t \tilde {\rho }_\alpha + \textrm{div}(\tilde {\rho }_\alpha \boldsymbol{v}_\alpha ) + \textrm{div} \boldsymbol{j}_\alpha = \zeta _\alpha . \end{align}

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

(2.21a) \begin{align} \partial _t \tilde {\rho }_\alpha + \textrm{div}(\tilde {\rho }_\alpha \boldsymbol{v}) + \textrm{div} (\boldsymbol{J}_\alpha + \boldsymbol{j}_\alpha ) &=\, \zeta _\alpha , \end{align}

(2.21b) \begin{align} \partial _t \tilde {\rho }_\alpha + \textrm{div}(\tilde {\rho }_\alpha \boldsymbol{u}) + \textrm{div} (\boldsymbol{J}_\alpha ^u + \boldsymbol{j}_\alpha ) &=\, \zeta _\alpha , \end{align}

(2.21c) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \textrm{div} \boldsymbol{h}_\alpha + \rho _\alpha ^{-1}\textrm{div} \boldsymbol{j}_\alpha &=\, \rho _\alpha ^{-1}\zeta _\alpha , \end{align}

(2.21d) \begin{align} \partial _t \phi _\alpha + \textrm{div}\left (\phi _\alpha \boldsymbol{u}\right )+ \textrm{div}\boldsymbol{h}_\alpha ^u + \rho _\alpha ^{-1}\textrm{div} \boldsymbol{j}_\alpha &=\, \rho _\alpha ^{-1}\zeta _\alpha , \end{align}

(2.21e) \begin{align} \rho \partial _t c_\alpha + \rho \boldsymbol{v} {\boldsymbol\cdot} \nabla c_\alpha + \textrm{div} (\boldsymbol{J}_\alpha + \boldsymbol{j}_\alpha ) &=\, \zeta _\alpha , \end{align}

(2.21f) \begin{align} \rho \partial _t c_\alpha + \rho \boldsymbol{u} {\boldsymbol\cdot} \nabla c_\alpha + \textrm{div} (\boldsymbol{J}_\alpha ^u + \boldsymbol{j}_\alpha ) - c_\alpha \textrm{div}\!\left ( \sum _{\beta } \boldsymbol{J}_{\beta }^u \right ) &=\, \zeta _\alpha . \end{align}

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:

(2.22) \begin{align} \partial _t \left (\frac {c_\alpha }{\rho _\alpha }\left (\sum _{\beta }\frac {c_{\beta }}{\rho _{\beta }}\right )^{-1}\right ) + \textrm{ div}\left (\frac {c_\alpha }{\rho _\alpha }\left (\sum _{\beta }\frac {c_{\beta }}{\rho _{\beta }}\right )^{-1} \boldsymbol{v}\right ) &\nonumber \\ +\, \textrm{div} \boldsymbol{h}_\alpha + \rho _\alpha ^{-1}\textrm{div} \boldsymbol{j}_\alpha & = \rho _\alpha ^{-1}\zeta _\alpha . \end{align}

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

(2.23a) \begin{align} \partial _t (\tilde {\rho }_\alpha \boldsymbol{v} + \boldsymbol{J}_\alpha ) + \textrm{div} \left ( \tilde {\rho }_\alpha \boldsymbol{v}\otimes \boldsymbol{v} +\boldsymbol{J}_\alpha \otimes \boldsymbol{v} + \boldsymbol{v}\otimes \boldsymbol{J}_\alpha \right )&\nonumber \\ - \, \textrm{div} \left (\boldsymbol{T}_\alpha - \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \right ) - \tilde {\rho }_\alpha \boldsymbol{b}_\alpha &=\, \boldsymbol{\pi }_\alpha , \end{align}

(2.23b) \begin{align} \partial _t (\tilde {\rho }_\alpha \boldsymbol{u} + \boldsymbol{J}_\alpha ^u) + \textrm{div} \left ( \tilde {\rho }_\alpha \boldsymbol{u}\otimes \boldsymbol{u} +\boldsymbol{J}_\alpha ^u \otimes \boldsymbol{u} + \boldsymbol{u}\otimes \boldsymbol{J}_\alpha ^u \right )&\nonumber \\ + \, \textrm{div} \left (\tilde {\rho }_\alpha \boldsymbol{\omega }_\alpha \otimes \boldsymbol{\omega }_\alpha - \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \right )&\nonumber \\ -\, \textrm{div} \left (\boldsymbol{T}_\alpha - \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \right ) - \tilde {\rho }_\alpha \boldsymbol{b}_\alpha &=\, \boldsymbol{\pi }_\alpha . \end{align}

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

(2.24a) \begin{align} \mathscr{K}_\alpha =&\,\tilde {\rho }_\alpha \|\boldsymbol{v}_\alpha \|^2/2, \end{align}

(2.24b) \begin{align} \mathscr{G}_\alpha =&\,\tilde {\rho }_\alpha b y, \end{align}

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

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, utilise the formulations (2.21a ) and (2.23a ) and invoke the identity (2.13a ) to obtain

(2.25a) \begin{align} \partial _t \rho + \textrm{div}(\rho \boldsymbol{v}) &=\, 0, \end{align}

(2.25b) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v} \otimes \boldsymbol{v} \right ) - \textrm{div} \boldsymbol{T} - \rho \boldsymbol{b} &=\,0, \end{align}

(2.25c) \begin{align} \boldsymbol{T}-\boldsymbol{T}^{T} &=\,0, \end{align}

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

(2.26a) \begin{align} \boldsymbol{T} =&\, \sum _\alpha \boldsymbol{T}_\alpha -\tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha , \end{align}

(2.26b) \begin{align} \boldsymbol{b} =&\,\frac {1}{\rho }\sum _\alpha \tilde {\rho }_\alpha \boldsymbol{b}_\alpha , \end{align}

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

(2.27a) \begin{align} \displaystyle \sum _\alpha \gamma _\alpha =&\, 0, \end{align}

(2.27b) \begin{align} \displaystyle \sum _\alpha \boldsymbol{\pi }_\alpha =&\, 0, \end{align}

(2.27c) \begin{align} \displaystyle \sum _\alpha \boldsymbol{N}_\alpha =&\, 0, \end{align}

and where we invoke (2.27a ) via

(2.28a) \begin{align} \sum _\alpha \zeta _\alpha =&\, 0, \end{align}

(2.28b) \begin{align} \sum _\alpha \boldsymbol{j}_\alpha =&\, 0. \end{align}

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$ , $\boldsymbol{v}$ , $\boldsymbol{T}$ and $\boldsymbol{b}$ , as defined earlier. 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

(2.29a) \begin{align} \partial _t \rho + \textrm{div}\bigg(\rho \bigg(\boldsymbol{u} + \rho ^{-1}\sum _\alpha \boldsymbol{J}_\alpha ^u \bigg)\bigg) &=\, \sum _\alpha \gamma _\alpha = 0, \end{align}

(2.29b) \begin{align} \textrm{div}\bigg(\boldsymbol{v} + \sum _\alpha \boldsymbol{h}_\alpha \bigg) &=\, \sum _\alpha \rho _\alpha ^{-1}\gamma _\alpha , \end{align}

(2.29c) \begin{align} \textrm{div}\boldsymbol{u} &=\, \sum _\alpha \rho _\alpha ^{-1}\gamma _\alpha , \end{align}

(2.29d) \begin{align} 0 &=\, \sum _\alpha \gamma _\alpha , \end{align}

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 Lemma2.4. Next, note that (2.29b ) also follows from the summation over the constituents of (2.22). With the aid of Lemma2.4, one can infer that (2.29b ) and (2.29c ) are identical. Furthermore, $\boldsymbol{v} + \sum _\alpha \boldsymbol{h}_\alpha = \boldsymbol{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 }$ ). Next, (2.29d ) complies with the balance condition (2.27a ) and shows that no velocity divergence equation results from using concentration variables. Finally, the summation of (2.23) yields

(2.30) \begin{align} \partial _t \bigg(\rho \bigg(\boldsymbol{u} + \rho ^{-1}\sum _\alpha \boldsymbol{J}_\alpha ^u\bigg)\bigg) + \textrm{div} \bigg( \rho \bigg(\boldsymbol{u}\otimes \boldsymbol{u} + \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u \otimes \boldsymbol{u} + \rho ^{-1}\boldsymbol{u}\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u \bigg)\bigg)\nonumber \\ + \, \textrm{div} \bigg( \sum _\alpha \tilde {\rho }_\alpha \boldsymbol{\omega }_\alpha \otimes \boldsymbol{\omega }_\alpha - \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \bigg)- \textrm{div} \boldsymbol{T} - \rho \boldsymbol{b} =\,0. \end{align}

Invoking Lemma2.4 and Lemma2.5, this may be written as

(2.31) \begin{align}& \partial _t \bigg(\rho \bigg(\boldsymbol{u} + \rho ^{-1}\sum _\alpha \boldsymbol{J}_\alpha ^u\bigg)\bigg) + \textrm{div} \bigg( \rho \bigg(\boldsymbol{u}\otimes \boldsymbol{u} + \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u \otimes \boldsymbol{u} + \rho ^{-1}\boldsymbol{u}\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u \bigg)\bigg)\nonumber \\& \qquad + \, \textrm{div} \bigg( \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u\bigg) - \textrm{div} \boldsymbol{T} - \rho \boldsymbol{b} =\,0. \end{align}

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

(2.32a) \begin{align} \textrm{div} \bigg ( \rho \bigg (\boldsymbol{u}\otimes \boldsymbol{u} + \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u \otimes \boldsymbol{u} + \rho ^{-1}\boldsymbol{u}\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u \bigg )\bigg ) = \nonumber \\ \textrm{div} \bigg ( \rho \boldsymbol{v} \otimes \boldsymbol{v} -\rho (\boldsymbol{v}-\boldsymbol{u})\otimes (\boldsymbol{v}-\boldsymbol{u})\bigg ), \end{align}

(2.32b) \begin{align} \textrm{div} \bigg ( \sum _\alpha \tilde {\rho }_\alpha \boldsymbol{\omega }_\alpha \otimes \boldsymbol{\omega }_\alpha - \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \bigg ) = \textrm{div} \bigg (\rho (\boldsymbol{v}-\boldsymbol{u})\otimes (\boldsymbol{v}-\boldsymbol{u})\bigg ). \end{align}

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

(2.33a) \begin{align} \partial _t \rho + \textrm{div}\bigg (\rho \boldsymbol{u} + \sum _\alpha \boldsymbol{J}_\alpha ^u \bigg ) &=\, 0, \end{align}

(2.33b) \begin{align} \partial _t \bigg (\rho \boldsymbol{u} + \sum _\alpha \boldsymbol{J}_\alpha ^u\bigg ) + \textrm{div} \bigg ( \rho \boldsymbol{u}\otimes \boldsymbol{u} + \sum _\alpha \boldsymbol{J}_\alpha ^u \otimes \boldsymbol{u} \bigg .& \nonumber \\ + \bigg . \boldsymbol{u}\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u + \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u\otimes \sum _\alpha \boldsymbol{J}_\alpha ^u\bigg ) - \textrm{div} \boldsymbol{T} - \rho \boldsymbol{b} &=\,0, \end{align}

(2.33c) \begin{align} \boldsymbol{T}-\boldsymbol{T}^{T} &=\,0. \end{align}

The various forms presented in this section show that the set of balance laws, on both constituent level (§ 2.2) and mixture level (§ 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:

(2.34a) \begin{align} \mathscr{K} =&\, \displaystyle \sum _{\alpha } \mathscr{K}_\alpha , \end{align}

(2.34b) \begin{align} \mathscr{G} =&\, \displaystyle \sum _{\alpha } \mathscr{G}_\alpha . \end{align}

The kinetic energy of the mixture can be decomposed as

(2.35a) \begin{align} \mathscr{K} =&\, \bar {\mathscr{K}} + \displaystyle \sum _{\alpha } \tfrac {1}{2} \tilde {\rho }_\alpha \|\boldsymbol{w}_\alpha \|^2, \end{align}

(2.35b) \begin{align} \bar {\mathscr{K}} =&\, \tfrac {1}{2} \rho \|\boldsymbol{v}\|^2, \end{align}

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

(2.36) \begin{align} \displaystyle \sum _{\alpha } \frac {1}{2} \tilde {\rho }_\alpha \|\boldsymbol{w}_\alpha \|^2 = \displaystyle \sum _{\alpha } \frac {1}{2} \tilde {\rho }_\alpha \Bigg\|\boldsymbol{\omega }_\alpha - \rho ^{-1} \sum _\alpha \boldsymbol{J}_\alpha ^u \Bigg\|^2. \end{align}

3.1. Assumptions and modelling 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:

(3.1a) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div} \boldsymbol{H}_\alpha &=\, \rho _\alpha ^{-1}\zeta _\alpha , \end{align}

(3.1b) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) - \textrm{div} \boldsymbol{T} - \rho \boldsymbol{b} &=\,0, \end{align}

(3.1c) \begin{align} \boldsymbol{T}-\boldsymbol{T}^{T} &=\,0, \end{align}

with $\boldsymbol{H}_\alpha := \boldsymbol{J}_\alpha + \boldsymbol{j}_\alpha$ , where (3.1a ) holds for constituents $\alpha =1,\ldots ,N$ . At this point, the system comprises the unknown quantities: volume fractions $\phi _\alpha$ ( $\alpha = 1,\ldots ,N$ ), where we recall the identity (2.4b ), mass-averaged mixture velocity $\boldsymbol{v}$ , peculiar velocities $\boldsymbol{J}_\alpha$ ( $\alpha = 1,\ldots ,N$ ), mass transfer terms $\zeta _\alpha , \boldsymbol{j}_\alpha$ ( $\alpha = 1,\ldots ,N$ ) and mixture stress $\boldsymbol{T}$ . In order to close the system, we seek for constitutive models for $\boldsymbol{J}_\alpha$ , $\boldsymbol{j}_\alpha$ , $\zeta _\alpha$ and $\boldsymbol{T}$ . Seeking for constitutive models for the peculiar velocities $\boldsymbol{J}_\alpha$ could be perceived as a simplification procedure. Namely, substituting a constitutive model (in § 3.4), in general, violates the continuum mixture theory definitions (2.12). We discard these definitions (2.12) in the following, but design models compatible with Lemma2.4 and Lemma2.5 to ensure invariance with respect to the set of fundamental variables. Instead of working with $N$ velocities quantities $\boldsymbol{v}_\alpha$ , the simplified system contains a single unknown velocity quantity $\boldsymbol{v}$ and constitutive models for peculiar velocities $\boldsymbol{J}_\alpha$ . This is compatible with the structure of the system: the full system is composed of $N$ linear momentum (mixture) balance laws 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:

(3.2a) \begin{align} \sum _\alpha \boldsymbol{J}_\alpha =&\, 0, \end{align}

(3.2b) \begin{align} \sum _\alpha \boldsymbol{j}_\alpha =&\, 0, \end{align}

(3.2c) \begin{align} \sum _\alpha \zeta _\alpha =&\, 0, \end{align}

where we recall the decomposition (2.19). The system (3.1) contains the unknown variables $\boldsymbol{v}$ and $\phi _\alpha$ ( $\alpha = 1,\ldots ,N$ ). We emphasise that directly enforcing the summation condition (2.5b ) at this point would imply that the set $\left \{\phi _\alpha \right \}_{\alpha =1,\ldots ,N}$ comprises $N-1$ independent variables. As such, system (3.1) would have a degenerate nature; it contains $N+1$ equations for $N$ variables (we preclude (2.25c ) in this count). Instead, a natural approach to restore the balance is by enforcing the constraint with a Lagrange multiplier construction, see § 3.2.

We adopt the well-known Coleman–Noll procedure (Coleman & Noll Reference Coleman and Noll1974) as a guiding principle to design constitutive models. For this purpose, we postulate the energy-dissipation law:

(3.3) \begin{align} \frac {\textrm{d}}{\textrm{d}t} \mathscr{E} = \mathscr{W} - \mathscr{D}, \end{align}

satisfying $\mathscr{D}\geqslant 0$ . The total energy comprises the Helmholtz free energy, the kinetic energy and the gravitational energy:

(3.4) \begin{align} \mathscr{E} = \displaystyle \int _{\mathcal{R}(t)}(\Psi + \bar {\mathscr{K}} + \mathscr{G})\,\textrm{d}v. \end{align}

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

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

(3.5) \begin{align} \Psi =&\, \hat {\Psi }\left (\left \{\phi _\alpha \right \}_{\alpha =1,\ldots ,N},\left \{\nabla \phi _\alpha \right \}_{\alpha =1,\ldots ,N}\right ), \end{align}

and introduce the chemical potential quantities ( $\alpha = 1,\ldots ,N$ )

(3.6) \begin{align} \hat {\mu }_\alpha =&\, \frac { \partial \hat {\Psi }}{\partial \phi _\alpha } - \textrm{div}\frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }. \end{align}

At this point the volume fractions $\left \{\phi _\alpha \right \}_{\alpha =1,\ldots ,N}$ (and their gradients $\left \{\nabla \phi _\alpha \right \}_{\alpha =1,\ldots ,N}$ ) are independent quantities and (3.5) and (3.6) are obviously well-defined. However, this is no longer the case when the saturation constraint (2.5b ) would be directly enforced, which would make the chemical potentials individually arbitrary. Namely, addition of the term $(1-\sum _\alpha \phi _\alpha )$ to $\hat {\Psi }$ does not alter it, but it modifies the chemical potentials $\mu _\alpha$ . We return to this point at the end of this subsection.

3.2. Modelling restriction

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

(3.7) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\Psi } \,\textrm{d}v = \displaystyle \int _{\mathcal{R}(t)} \partial _t \hat {\Psi } \,\textrm{d}v + \displaystyle \int _{\partial \mathcal{R}(t)} \hat {\Psi } \boldsymbol{v} {\boldsymbol\cdot} \boldsymbol{\nu } \,\textrm{d}a. \end{align}

We notice that directly enforcing the summation constraint (2.5b ) would 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) , i.e.

(3.8) \begin{align} \textrm{d}\hat {\Psi } = \sum _{\alpha } \frac {\partial \hat {\Psi }}{\partial \phi _\alpha } \textrm{d}\phi _\alpha + \sum _{\alpha } \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } \textrm{d}(\nabla \phi _\alpha ), \end{align}

is not altered by enforcing the summation constraint (2.5b), where $\textrm{d}$ is the derivative operator.

Invoking Lemma3.5 and the divergence theorem yields

(3.9) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\int _{\mathcal{R}(t)} \hat {\Psi } \,\textrm{d}v = \int _{\mathcal{R}(t)} &\ \hat {\Psi } \,\textrm{div} \boldsymbol{v} + \sum _{\alpha }\frac {\partial \hat {\Psi }}{\partial \phi _\alpha } \dot {\phi }_\alpha +\sum _{\alpha } \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }{\boldsymbol{\cdot}} \left (\nabla \phi _\alpha \right )^{\boldsymbol{\cdot}}\,\textrm{d}v. \end{align}

Integrating by parts provides

(3.10) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\Psi } \,\textrm{d}v =&\, \displaystyle \int _{\mathcal{R}(t)} \hat {\Psi }\,\textrm{div} \boldsymbol{v} + \sum _{\alpha } \hat {\mu }_\alpha \dot {\phi }_\alpha - \sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }: \nabla \boldsymbol{v} \,\textrm{d}v \nonumber \\ &\,+ \displaystyle \int _{\partial \mathcal{R}(t)}\sum _{\alpha } \dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }{\boldsymbol\cdot} \boldsymbol{\nu } \,\textrm{d}a, \end{align}

where we have substituted the identity

(3.11) \begin{align} (\nabla \unicode{x03C8} )^{\boldsymbol{\cdot}} = \nabla (\dot {\psi }) - (\nabla \psi )^{T}\nabla \boldsymbol{v} \end{align}

for $\psi = \phi _\alpha$ . We note that the free energy terms are well-defined.

Lemma 3.6 (Well-defined free energy terms). The following free energy terms in(3.10)are well-defined when enforcing the summation constraint (2.5b):

(3.12) \begin{align} \sum _{\alpha } \hat {\mu }_\alpha \dot {\phi }_\alpha ; \quad \sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }; \quad \sum _{\alpha } \dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }. \end{align}

Substituting the constituent mass balance laws (3.1a ) provides

(3.13) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \hat {\Psi } \,\textrm{d}v =&\, \displaystyle \int _{\mathcal{R}(t)} \hat {\Psi }\,\textrm{div} \boldsymbol{v} + \sum _{\alpha } \hat {\mu }_\alpha \left(-\phi _\alpha \textrm{ div}\boldsymbol{v} - \rho _\alpha ^{-1}\textrm{div} \boldsymbol{H}_\alpha + \rho _\alpha ^{-1}\zeta _\alpha \right) \nonumber \\ &- \sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }: \nabla \boldsymbol{v} \,\textrm{d}v + \displaystyle \int _{\partial \mathcal{R}(t)}\sum _{\alpha } \dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }{\boldsymbol\cdot} \boldsymbol{\nu } \,\textrm{d}a, \end{align}

where we recall $\boldsymbol{H}_\alpha = \boldsymbol{J}_\alpha + \boldsymbol{j}_\alpha$ . By again applying integration by parts one can infer that

(3.14) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \sum _{\alpha } \hat {\Psi } \,\textrm{d}v = &\, \displaystyle \int _{\mathcal{R}(t)} \hat {\Psi }\,\textrm{div} \boldsymbol{v} -\sum _{\alpha } \hat {\mu }_\alpha \phi _\alpha \textrm{div}\boldsymbol{v} + \sum _{\alpha } \nabla (\rho _\alpha ^{-1}\hat {\mu }_\alpha )\ {\boldsymbol\cdot}\ \boldsymbol{H}_\alpha \nonumber \\ &\ \ \ - \sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }: \nabla \boldsymbol{v} + \sum _{\alpha } \rho _\alpha ^{-1}\hat {\mu }_\alpha \zeta _\alpha \,\textrm{d}v\nonumber \\ &\,+ \displaystyle \int _{\partial \mathcal{R}(t)}\sum _{\alpha } \left (\dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }-\rho _\alpha ^{-1}\hat {\mu }_\alpha \boldsymbol{H}_\alpha \right )\,{\boldsymbol\cdot}\, \boldsymbol{\nu } \,\textrm{d}a. \end{align}

Next, the evolution of the kinetic and gravitational energies take the form (see ten Eikelder et al. (Reference ten Eikelder, van der Zee, Akkerman and Schillinger2023) for details)

(3.15a) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \mathscr{K} \,\textrm{d}v =&\, \displaystyle \int _{\mathcal{R}(t)} - \nabla \boldsymbol{v} : \boldsymbol{T}+\rho \boldsymbol{v}\ {\boldsymbol\cdot}\ \boldsymbol{b} \,\textrm{d}v+ \displaystyle \int _{\partial \mathcal{R}(t)} \boldsymbol{v}\ {\boldsymbol\cdot}\ \boldsymbol{T} \boldsymbol{\nu } \,\textrm{d}a, \end{align}

(3.15b) \begin{align} \frac {\textrm{d}}{\textrm{d}t}\displaystyle \int _{\mathcal{R}(t)} \mathscr{G} \,\textrm{d}v =&\,- \displaystyle \int _{\mathcal{R}(t)} \rho \boldsymbol{v}\ {\boldsymbol\cdot}\ \boldsymbol{b}\,\textrm{d}v. \end{align}

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

(3.16) \begin{align} \frac {\textrm{d}}{\textrm{d}t} \mathscr{E} = &\, \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^{T}\boldsymbol{T}-\sum _{\alpha } \left (\rho _\alpha ^{-1}\hat {\mu }_\alpha \boldsymbol{H}_\alpha -\dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }\right )\right ){\boldsymbol\cdot}\, \boldsymbol{\nu } \,\textrm{d}a \nonumber \\ &\,- \displaystyle \int _{\mathcal{R}(t)} \left (\boldsymbol{T} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_\alpha \phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v}\nonumber \\ &\,\,\,\,\,\,\,\,\,\,\,+ \sum _{\alpha }\Big(-\nabla \big(\rho _\alpha ^{-1}\hat {\mu }_\alpha \big)\ {\boldsymbol\cdot}\ \boldsymbol{H}_\alpha - \rho _\alpha ^{-1}\hat {\mu }_\alpha \zeta _\alpha \Big)\,\textrm{d}v. \end{align}

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

(3.17) \begin{align} 0=&\, \lambda \left ( \textrm{div}\boldsymbol{v} + \displaystyle \sum _{\alpha } \rho _\alpha ^{-1} \textrm{div} \boldsymbol{H}_\alpha - \displaystyle \sum _{\alpha }\rho _\alpha ^{-1}\zeta _\alpha \right ), \end{align}

where $\lambda$ is the scalar Lagrange multiplier.

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

(3.18) \begin{align} \frac {\textrm{d}}{\textrm{d}t} \mathscr{E} = &\, \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^{T}\boldsymbol{T}-\sum _{\alpha } \left (g_\alpha \boldsymbol{H}_\alpha -\dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }\right )\right ){\boldsymbol\cdot}\, \boldsymbol{\nu } \,\textrm{d}a \nonumber \\ &\,- \displaystyle \int _{\mathcal{R}(t)} \left (\boldsymbol{T} + \lambda \boldsymbol{I} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_\alpha \phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v}\nonumber \\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \sum _{\alpha }\left (-\nabla g_\alpha\ {\boldsymbol\cdot}\ \boldsymbol{H}_\alpha - g_\alpha \zeta _\alpha \right )\,\textrm{d}v, \end{align}

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

(3.19a) \begin{align} g_\alpha :=&\, \rho _\alpha ^{-1}\hat {\mu }_{\alpha ,\lambda }, \end{align}

(3.19b) \begin{align} \hat {\mu }_{\alpha ,\lambda } :=&\, \hat {\mu }_\alpha +\lambda . \end{align}

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

(3.20a) \begin{align} \mathscr{W} =&\, \displaystyle \int _{\partial \mathcal{R}(t)}\left (\boldsymbol{v}^{T}\boldsymbol{T}-\sum _{\alpha } \left (g_\alpha \boldsymbol{H}_\alpha -\dot {\phi }_\alpha \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }\right )\right ){\boldsymbol\cdot}\, \boldsymbol{\nu } \,\textrm{ d}a, \end{align}

(3.20b) \begin{align} \mathscr{D} =&\,\displaystyle \int _{\mathcal{R}(t)} \left (\boldsymbol{T} + \lambda \boldsymbol{I} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_\alpha \phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v}\nonumber \\ &\qquad + \sum _{\alpha }\left (-\nabla g_\alpha\ {\boldsymbol\cdot}\ \boldsymbol{H}_\alpha -g_\alpha \zeta _\alpha \right )\,\textrm{d}v. \end{align}

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

(3.21) \begin{align} \left (\boldsymbol{T} + \lambda \boldsymbol{I} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_\alpha \phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v}&\nonumber \\ - \sum _{\alpha }\nabla g_\alpha {\boldsymbol\cdot}\ \boldsymbol{H}_\alpha -\sum _{\alpha } g_\alpha \zeta _\alpha &\geqslant 0. \end{align}

We finalise this section with a remark on the connection between the chemical potentials and the Lagrange multiplier. Directly enforcing the saturation constraint (2.5b ) provides

(3.22) \begin{align} \lambda + \sum _{\alpha }\hat {\mu }_\alpha \phi _\alpha = \sum _\alpha \hat {\mu }_{\alpha ,\lambda } \phi _\alpha . \end{align}

This observation reveals that chemical potentials in (3.21) occur solely 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$ should not alter the free energy. Indeed, we have

(3.23) \begin{align} \Psi =&\, \hat {\Psi }\left (\left \{\phi _\alpha \right \}_{\alpha =1,\ldots ,N},\left \{\nabla \phi _\alpha \right \}_{\alpha =1,\ldots ,N}\right ) + \lambda \left (\sum _\alpha \phi _\alpha -1 \right ), \end{align}

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

(3.24) \begin{align} \hat {\mu }_{\alpha ,\lambda } = \frac { \partial \Psi }{\partial \phi _\alpha } - \textrm{div}\frac {\partial \Psi }{\partial \nabla \phi _\alpha }, \end{align}

where we recall (3.19b ). We do not directly enforce (2.5b ) and continue with the left-hand side of (3.22).

3.3. Alternative free energy classes

As mentioned in Remark3.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. Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998). We consider the following constitutive class:

(3.25) \begin{align} \Psi =&\, \check {\Psi }\left (\left \{c_\alpha \right \},\left \{\nabla c_\alpha \right \}\right ), \end{align}

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

(3.26) \begin{align} \check {\mu }_\alpha =&\, \frac { \partial \hat {\Psi }}{\partial c_\alpha } - \textrm{div}\frac {\partial \hat {\Psi }}{\partial \nabla c_\alpha }. \end{align}

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

(3.27) \begin{align} \left (\boldsymbol{T} + \check {\lambda }\boldsymbol{I} +\sum _{\alpha } \nabla c_\alpha \otimes \frac {\partial \check {\Psi }}{\partial \nabla c_\alpha } - \check {\Psi }\boldsymbol{I}\right ):\nabla \boldsymbol{v}&\nonumber \\ - \sum _{\alpha }\nabla \left ( \rho ^{-1}\check {\mu }_\alpha + \rho _\alpha ^{-1}\check {\lambda } \right ){\boldsymbol\cdot}\ \boldsymbol{H}_\alpha -\sum _{\alpha } \left ( \rho ^{-1}\check {\mu }_\alpha + \rho _\alpha ^{-1}\check {\lambda } \right ) \zeta _\alpha &\geqslant 0, \end{align}

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):

(3.28a) \begin{align} \phi _\alpha =&\, \phi _\alpha ( \{c_{\beta } \}), \end{align}

(3.28b) \begin{align} c_{\beta } =&\, c_{\beta }(\{\phi _\alpha \}), \end{align}

the identification

(3.29) \begin{align} \hat {\Psi }(\{\phi _\alpha\}) =\check {\Psi }(\{c_{\beta }\}) \end{align}

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

For an alternative path to show equivalence of the modelling restrictions, one could apply the variable transformation (3.28) defined in (2.7) to show that (3.27) coincides with (3.21). We discuss this approach in Appendix B.

Guided by Theorem3.9, we proceed with the formulation of the modelling restriction presented in (3.21).

3.4. Selection of constitutive models

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

(3.30a) \begin{align} \boldsymbol{T} =&\, \hat {\boldsymbol{T}}\left (\nabla \boldsymbol{v}, \left \{\phi _\alpha \right \}, \left \{\nabla \phi _\alpha \right \}, \left \{g_{\alpha }\right \}, \left \{\nabla g_{\alpha }\right \}\right ), \end{align}

(3.30b) \begin{align} \boldsymbol{J}_\alpha =&\, \hat {\boldsymbol{J}}_\alpha \left (\left \{\phi _\alpha \right \}, \left \{\nabla g_\alpha \right \}\right ), \end{align}

(3.30c) \begin{align} \boldsymbol{j}_\alpha =&\, \hat {\boldsymbol{j}}_\alpha \left (\left \{\phi _\alpha \right \}, \left \{\nabla g_\alpha \right \}\right ), \end{align}

(3.30d) \begin{align} \zeta _\alpha =&\, \hat {\zeta }_\alpha \left (\left \{\phi _\alpha \right \}, \left \{g_\alpha \right \}\right ), \end{align}

and define $\hat {\boldsymbol{H}}_\alpha = \hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha$ . Generally speaking, the introduction of the class (3.30b ) deviates from continuum mixture theory. Arguably, a natural approximation is simply taking $\hat {\boldsymbol{J}}_\alpha = 0$ , which, for instance, models the situation of matching velocities $\boldsymbol{v}_\alpha =\boldsymbol{v}_{\beta }$ . We return to this case in § 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.21). To this end, we aim to identify constitutive models (3.30) so that all three terms in (3.21) are positive, which occurs when

(3.31a) \begin{align} \left (\boldsymbol{T} + \lambda \boldsymbol{I} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v}&\geqslant 0, \end{align}

(3.31b) \begin{align} -\sum _{\alpha }\nabla g_\alpha {\boldsymbol\cdot}\ \boldsymbol{J}_\alpha &\geqslant 0, \end{align}

(3.31c) \begin{align} -\sum _{\alpha }\nabla g_\alpha {\boldsymbol\cdot}\ \boldsymbol{j}_\alpha &\geqslant 0, \end{align}

(3.31d) \begin{align} -\sum _{\alpha } g_\alpha \zeta _\alpha &\geqslant 0. \end{align}

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:

(3.32) \begin{align} \hat {\boldsymbol{T}} = -\lambda \boldsymbol{I}-\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha }-\left (\sum _{\alpha } \hat {\mu }_{\alpha }\phi _\alpha - \hat {\Psi }\right )\boldsymbol{I} + \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}), \end{align}

subject to the symmetry condition

(3.33) \begin{align} \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } = \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } \otimes \nabla \phi _\alpha , \end{align}

where the scalar field $\nu \geqslant 0$ is the mixture dynamic viscosity, $\bar {\lambda } \geqslant -2/d$ is a scalar, and $d$ is the number of dimensions. Possible choices for the mixture viscosity include $\nu = \sum _\alpha \nu _\alpha \phi _\alpha$ and $\nu = \sum _\alpha \nu _\alpha c_\alpha$ , where $\nu _\alpha$ are constituent viscosities. The condition (3.33) ensures compatibility with the angular momentum constraint (2.25c ).

Proof. An elementary calculation gives

(3.34) \begin{align} \left (\hat {\boldsymbol{T}} +\lambda \boldsymbol{I} +\sum _{\alpha } \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } +\left (\sum _{\alpha } \hat {\mu }_{\alpha }\phi _\alpha - \hat {\Psi }\right )\boldsymbol{I}\right ):\nabla \boldsymbol{v} &= \nonumber \\ 2 \nu \left ( \nabla ^s \boldsymbol{v} - \frac {1}{d} (\textrm{div} \boldsymbol{v}) \boldsymbol{I}\right ):\left (\nabla ^s \boldsymbol{v} - \frac {1}{d} (\textrm{div} \boldsymbol{v}) \boldsymbol{I}\right ) + \nu \left (\bar {\lambda } + \frac {2}{d}\right )\left (\textrm{div} \boldsymbol{v}\right )^2 &\geqslant 0. \end{align}

Constituent peculiar velocities. We choose the peculiar velocities of the form

(3.35) \begin{align} \hat {\boldsymbol{J}}_\alpha = -\sum _{\beta } \boldsymbol{M}_{\alpha {\beta }} \nabla g_{\beta }, \end{align}

with mobility tensor $\boldsymbol{M}_{\alpha {\beta }}=\boldsymbol{M}_{{\beta }\alpha }$ . The mobility tensor is positive definite ( $\boldsymbol{y}_\alpha ^{T} \boldsymbol{M}_{\alpha {\beta }}\boldsymbol{y}_{\beta } \geqslant 0$ for all $\boldsymbol{y}_\alpha \in \mathbb{R}^d, \alpha = 1,\ldots ,N$ ), has the same dependencies as (3.31b ), is compatible with $\sum _\alpha \boldsymbol{M}_{\alpha {\beta }} = \sum _\alpha \boldsymbol{M}_{{\beta }\alpha } = 0$ for all ${\beta }=1,..,N$ , and vanishes in the single fluid region $\boldsymbol{M}_{\alpha {\beta }}|_{\phi _{\gamma }=1}=0, {\gamma } =1,\ldots ,N$ (thus is degenerate). We note that the symmetry requirement follows from the Onsager reciprocal relations, the positive definiteness from (3.31b ), and the zero sum of rows and columns from (3.2a ). A possible choice for the mobility tensor is $\boldsymbol{M}_{\alpha {\beta }}=-\boldsymbol{M}_0 \tilde {\rho }_\alpha \tilde {\rho }_{\beta }$ for $\alpha \neq {\beta }$ , and $\boldsymbol{M}_{\alpha \alpha } = \boldsymbol{M}_{0}\tilde {\rho }_\alpha \sum _{{\gamma }\neq \alpha } \tilde {\rho }_{\gamma }$ for some $\boldsymbol{M}_0$ that does not depend on the constituent number.

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

(3.36) \begin{align} \hat {\boldsymbol{j}}_\alpha = -\sum _{\beta } \boldsymbol{K}_{\alpha {\beta }} \nabla g_{\beta }, \end{align}

for some positive definite constitutive tensor $\boldsymbol{K}_{\alpha {\beta }}=\boldsymbol{K}_{{\beta }\alpha }$ compatible with $\sum _\alpha \boldsymbol{K}_{\alpha {\beta }}=\sum _\alpha \boldsymbol{K}_{{\beta }\alpha }=0$ , with the same dependencies as (3.30c ), and which vanishes in the single fluid region $\boldsymbol{K}_{\alpha {\beta }}|_{\phi _{\gamma }=1}=0, {\gamma } =1,\ldots ,N$ . Similarly as for the peculiar velocity, a possible choice for the mobility tensor is $\boldsymbol{K}_{\alpha {\beta }}=-\boldsymbol{K}_0 \tilde {\rho }_\alpha \tilde {\rho }_{\beta }$ for $\alpha \neq {\beta }$ , and $\boldsymbol{K}_{\alpha \alpha } = \boldsymbol{K}_{0}\tilde {\rho }_\alpha \sum _{{\gamma }\neq \alpha } \tilde {\rho }_{\gamma }$ for some $\boldsymbol{K}_0$ that is not dependent on the constituent number.

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

(3.37) \begin{align} \hat {\zeta }_\alpha = -\sum _{\beta } m_{\alpha {\beta }} g_{\beta }, \end{align}

where the positive definite scalar mobility $m_{\alpha {\beta }}=m_{{\beta }\alpha }$ has the same dependencies as (3.30d ), is compatible with $\sum _\alpha m_{\alpha {\beta }} = \sum _\alpha m_{{\beta }\alpha } = 0$ , and vanishes in the single fluid region $m_{\alpha {\beta }}|_{\phi _{\gamma }=1}=0, {\gamma } =1,\ldots ,N$ .

Remark 3.16 (Related constitutive models). In the case,

(3.38) \begin{align} \textbf{M}_{\alpha {\beta }} = \begin{cases} -\hat {\textbf{M}}_{\alpha {\beta }} & \text{ if }\alpha \neq {\beta },\\ \sum _{{\gamma }\neq \alpha } \hat {\textbf{M}}_{\alpha {\gamma }} & \text{ if }\alpha = {\beta }, \end{cases} \end{align}

for some symmetric $\hat {\textbf{M}}_{\alpha {\beta }}$ , we find

(3.39) \begin{align} \hat {\boldsymbol{J}}_\alpha =&\, - \sum _{{\beta }\neq \alpha } \textbf{M}_{\alpha {\beta }}\nabla \hat {\mu }_{\beta } - \textbf{M}_{\alpha \alpha }\nabla \hat {\mu }_\alpha \nonumber \\ =&\, \sum _{{\beta }\neq \alpha } \hat {\textbf{M}}_{\alpha {\beta }}\nabla \hat {\mu }_{\beta } - \sum _{{\gamma }\neq \alpha }\hat {\textbf{M}}_{\alpha {\gamma }}\nabla \hat {\mu }_\alpha \nonumber \\ =&\, -\sum _{{\beta }} \hat {\textbf{M}}_{\alpha {\beta }}\nabla (\hat {\mu }_\alpha -\hat {\mu }_{\beta }). \end{align}

This model matches (for the isotropic case $\textbf{M}_{\alpha {\beta }} = {\rm M}_{\alpha {\beta }}\textbf{I}$ ) that of ten Eikelder et al. (Reference ten Eikelder, van der Zee and Schillinger2024). It also closely resembles the form adopted in Li & Wang (Reference Li and Wang2014). Both closure models involve the Lagrange multiplier $\lambda$ ; a difference lies in the fact that the model proposed by Li & Wang (Reference Li and Wang2014) depends on the numbering of the constituents. Finally, we note that forms similar to (3.39) may be adopted for the diffusive fluxes and the mass transfer terms.

This finalises the construction of constitutive models compatible with the imposed energy-dissipative postulate. Substitution of the models (3.32), (3.35) and (3.37) yields the class of incompressible $N$ -phase models:

(3.40a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla \lambda - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )& \nonumber \\ + \textrm{div} \left (\left (\sum _\alpha \hat {\mu }_\alpha \phi _\alpha -\hat {\Psi }\right )\boldsymbol{I} + \sum _\alpha \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } \right )-\rho \boldsymbol{b}&=\, 0, \end{align}

(3.40b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\hat {\boldsymbol{J}}_\alpha + \rho _\alpha ^{-1}\textrm{ div}\hat {\boldsymbol{j}}_\alpha - \hat {\zeta }_\alpha &=\,0, \end{align}

(3.40c) \begin{align} \textrm{div} \boldsymbol{v} + \sum _\alpha \rho ^{-1}_\alpha \textrm{div}(\hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha ) - \rho ^{-1}_\alpha \hat {\zeta }_\alpha &=\,0, \end{align}

(3.40d) \begin{align} \hat {\mu }_\alpha - \frac {\partial \hat {\Psi }}{\partial \phi _\alpha } + \textrm{div} \left ( \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } \right )&=\,0, \end{align}

(3.40e) \begin{align} \hat {\boldsymbol{J}}_\alpha + \sum _{\beta } \boldsymbol{M}_{\alpha {\beta }} \nabla g_{\beta }&\,=0, \end{align}

(3.40f) \begin{align} \hat {\boldsymbol{j}}_\alpha + \sum _{\beta } \boldsymbol{K}_{\alpha {\beta }}\nabla g_{\beta } &\,=0, \end{align}

(3.40g) \begin{align} \hat {\zeta }_\alpha + \sum _{\beta } m_{\alpha {\beta }} g_{\beta } &\,=0, \end{align}

where (3.40b ) and (3.40c ) and the initial condition $\sum _{\beta } \phi _{{\beta }}(\boldsymbol{x},t=0) =1$ together ensure $\sum _{\beta } \phi _{{\beta }}(\boldsymbol{x},t) =1$ . Formulation (3.40) constitutes a class of models in the sense that particular closure relations ( $\boldsymbol{M}_{\alpha {\beta }}$ , $\boldsymbol{K}_{\alpha {\beta }}$ , $m_{\alpha {\beta }}$ and $\Psi$ ) need to be specified. A possible Cahn–Hilliard-type free energy is $\hat {\Psi }=\Psi _0(\left \{\phi _\alpha \right \})+{1}/{2}\sum _{\alpha {\beta }} \sigma _{\alpha {\beta }}\nabla \phi _\alpha {\boldsymbol\cdot} \nabla \phi _{\beta }$ , with $\sigma$ positive semidefinite. Given these relations, (3.40) is a well-defined closed model; this model is invariant to the renumbering of the constituents, invariant to the set of independent variables, and reduction-consistent. Additionally, it exhibits energy dissipation, which we state explicitly in the following theorem.

Proof. This follows from Lemma3.11, Lemma3.13 and Lemma3.15. In particular, the dissipation takes the form

(3.41) \begin{align} \mathscr{D} = &\, \displaystyle \int _{\mathcal{R}(t)} 2 \nu \left ( \nabla ^s \boldsymbol{v} - \frac {1}{d} (\textrm{div} \boldsymbol{v}) \boldsymbol{I}\right ):\left (\nabla ^s \boldsymbol{v} - \frac {1}{d} (\textrm{div} \boldsymbol{v}) \boldsymbol{I}\right )+ \nu \left (\bar\lambda + \frac {2}{d}\right )\left (\textrm{div} \boldsymbol{v}\right )^2 \nonumber \\ &\,\quad \quad + \displaystyle \sum _{\alpha ,{\beta }} \left (\nabla g_\alpha \right )^{T} \boldsymbol{B}_{\alpha {\beta }} \nabla g_{\beta }+ \sum _{\alpha ,{\beta }} m_{\alpha {\beta }}g_\alpha g_{\beta }\,\textrm{d}v \geqslant 0, \end{align}

with $\boldsymbol{B}_{\alpha {\beta }}=\boldsymbol{M}_{\alpha {\beta }}+\boldsymbol{K}_{\alpha {\beta }}$ .

4.1. Alternative formulations

As discussed in § 2 and § 3, the unified modelling 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

(4.1) \begin{align} p: = \sum _\alpha \hat {\mu }_\alpha \phi _\alpha -\hat {\Psi }. \end{align}

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

(4.2a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla (\lambda + p) + \textrm{div} \left ( \sum _{\beta } \nabla \phi _{\beta } \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _{\beta }} \right )& \nonumber \\ - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(4.2b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\hat {\boldsymbol{H}}_\alpha - \hat {\zeta }_\alpha &=\,0, \end{align}

(4.2c) \begin{align} \textrm{div} \boldsymbol{v} + \sum _\alpha \rho ^{-1}_\alpha \textrm{div}(\hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha ) - \rho ^{-1}_\alpha \hat {\zeta }_\alpha &=\,0. \end{align}

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

(4.3) \begin{align} \hat {\Psi } = \sum _\alpha \hat {\Psi }_\alpha , \end{align}

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

(4.4a) \begin{align} p = &\, \sum _\alpha p_\alpha , \end{align}

(4.4b) \begin{align} p_\alpha = &\, \sum _{{\beta }}\hat {\mu }_{\alpha {\beta }} \phi _{\beta } - \hat {\Psi }_\alpha , \end{align}

where $p_\alpha$ is the partial pressure of constituent $\alpha$ with $\hat {\mu }_{\alpha {\beta }}=\partial _{\phi _{\beta }}\hat {\Psi }_\alpha - \textrm{div}\partial _{\nabla \phi _{\beta }}\hat {\Psi }_\alpha$ . Thus, the split (4.3) reveals that the system may be written as

(4.5a) \begin{align} \sum _{\beta } \left ( \partial _t (\tilde {\rho }_{\beta } \boldsymbol{v}) + \textrm{div} \left ( \tilde {\rho }_{\beta } \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla p_{\beta } -\tilde {\rho }_{\beta } \boldsymbol{b} \right . &\nonumber \\ \left .+\, \textrm{div} \left ( \sum _{\gamma } \nabla \phi _{\gamma } \otimes \frac {\partial \hat {\Psi }_{\beta }}{\partial \nabla \phi _{\gamma }} \right )\right ) +\nabla \lambda - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar\lambda (\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right ) &=\, 0, \end{align}

(4.5b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\hat {\boldsymbol{H}}_\alpha - \hat {\zeta }_\alpha &=\,0, \end{align}

(4.5c) \begin{align} \textrm{div} \boldsymbol{v} + \sum _\alpha \rho ^{-1}_\alpha \textrm{div}(\hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha ) - \rho ^{-1}_\alpha \hat {\zeta }_\alpha &=\,0. \end{align}

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

(4.6) \begin{align} \textrm{div} \left (\left (\sum _\alpha \hat {\mu }_\alpha \phi _\alpha -\hat {\Psi }\right )\boldsymbol{I} + \sum _\alpha \nabla \phi _\alpha \otimes \frac {\partial \hat {\Psi }}{\partial \nabla \phi _\alpha } \right ) = \sum _\alpha \phi _\alpha \nabla \hat {\mu }_\alpha . \end{align}

Invoking Lemma4.1, model (3.40) takes a more compact form:

(4.7a) \begin{align} \partial _t (\rho \boldsymbol{v}) \!+\! \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) \!+\! \nabla \lambda \!+\! \sum _\alpha \phi _\alpha \nabla \hat {\mu }_\alpha - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}\!+\!\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &= 0, \end{align}

(4.7b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\hat {\boldsymbol{H}}_\alpha - \hat {\zeta }_\alpha &=\,0, \end{align}

(4.7c) \begin{align} \textrm{div} \boldsymbol{v} + \sum _\alpha \rho ^{-1}_\alpha \textrm{div}(\hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha ) - \rho ^{-1}_\alpha \hat {\zeta }_\alpha &=\,0. \end{align}

Considering the third and fourth terms in the momentum equation in isolation, when enforcing (2.5b ) these terms can be written as

(4.8) \begin{align} \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 . \end{align}

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:

(4.9a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla \lambda + \sum _{{\beta }} \phi _{\beta }\nabla \hat {\mu }_{\beta } &\nonumber \\ - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(4.9b) \begin{align} \partial _t \rho + \textrm{div}(\rho \boldsymbol{v}) &=\,0, \end{align}

(4.9c) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\hat {\boldsymbol{H}}_\alpha - \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

(4.9d) \begin{align} \textrm{div} \boldsymbol{v} + \sum _\alpha \rho ^{-1}_\alpha \textrm{div}(\hat {\boldsymbol{J}}_\alpha + \hat {\boldsymbol{j}}_\alpha ) - \rho ^{-1}_\alpha \hat {\zeta }_\alpha &=\,0. \end{align}

While we refrain from discussing formulations that adopt concentration variables, we discuss a formulation in terms of the volume-averaged velocity $\boldsymbol{u}$ . Inserting the constitutive model for the peculiar velocities (3.35) into Lemma2.4 we obtain

(4.10) \begin{align} \boldsymbol{v} = \boldsymbol{u} + \rho ^{-1}\displaystyle \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u. \end{align}

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

(4.11a) \begin{align} \partial _t \left (\rho \boldsymbol{u} + \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u\right ) + \textrm{div} \left ( \rho \boldsymbol{u}\otimes \boldsymbol{u} + \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u \otimes \boldsymbol{u} \right .& \nonumber \\+ \left . \boldsymbol{u}\otimes \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u + \rho ^{-1} \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u\otimes \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u\right ) + \nabla \lambda + \sum _{{\beta }} \phi _{\beta }\nabla \hat {\mu }_{\beta } & \nonumber \\- \textrm{div} \! \left ( \nu \! \left (2\nabla ^s \! \left (\boldsymbol{u} +\rho ^{-1} \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u \right )+\bar {\lambda }\textrm{ div} \! \left (\boldsymbol{u} +\rho ^{-1} \sum _{\beta } \hat {\boldsymbol{J}}_{\beta }^u \right ) \boldsymbol{I}\right ) \right )-\rho \boldsymbol{b}&=\, 0, \end{align}

(4.11b) \begin{align} \partial _t \phi _\alpha + \textrm{div}\left (\phi _\alpha \boldsymbol{u}\right )+ \textrm{div}\hat {\boldsymbol{h}}_\alpha ^u + \rho _\alpha ^{-1}\textrm{ div}\hat {\boldsymbol{j}}_\alpha -\rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

(4.11c) \begin{align} \textrm{div} \boldsymbol{u} +\sum _{\alpha } \rho _\alpha ^{-1} \hat {\boldsymbol{j}}_\alpha -\sum _{\alpha } \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

where the latter two terms in (4.11c ) vanish when densities match or mass transfer is absent (recall (2.29c )). 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 $\boldsymbol{M}_{\alpha {\beta }} = 0$ , we find

(4.12) \begin{align} \hat {\boldsymbol{J}}_\alpha =\hat {\boldsymbol{h}}_\alpha = \hat {\boldsymbol{J}}_\alpha ^u = \hat {\boldsymbol{h}}_\alpha ^u = 0. \end{align}

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

(4.13) \begin{align} \boldsymbol{v} = \boldsymbol{u}. \end{align}

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

(4.14a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla \lambda + \sum _{\beta }\phi _{\beta }\nabla \hat {\mu }_{\beta } &\nonumber \\ - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(4.14b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1} \textrm{div}\hat {\boldsymbol{j}}_\alpha - \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

(4.14c) \begin{align} \textrm{div} \boldsymbol{v} +\sum _{\alpha } \rho _\alpha ^{-1}\textrm{div} \hat {\boldsymbol{j}}_\alpha -\sum _{\alpha } \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

and, obviously,

(4.15a) \begin{align} \partial _t (\rho \boldsymbol{u}) + \textrm{div} \left ( \rho \boldsymbol{u}\otimes \boldsymbol{u} \right ) + \nabla \lambda + \sum _{\beta }\phi _{\beta }\nabla \hat {\mu }_{\beta } &\nonumber \\- \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{u}+\bar {\lambda }(\textrm{div}\boldsymbol{u}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(4.15b) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{u}) + \rho _\alpha ^{-1} \textrm{div}\hat {\boldsymbol{j}}_\alpha - \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0, \end{align}

(4.15c) \begin{align} \textrm{div} \boldsymbol{u} +\sum _{\alpha } \rho _\alpha ^{-1}\textrm{div} \hat {\boldsymbol{j}}_\alpha -\sum _{\alpha } \rho _\alpha ^{-1}\hat {\zeta }_\alpha &=\,0. \end{align}

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

4.3. Equilibrium conditions

We utilise formulation (4.7) to study equilibrium properties. We characterise the set equilibrium solutions $\left \{\textbf{q}_E=(\boldsymbol{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}(\textbf{q}_E) = 0$ . Invoking (3.41) yields the conditions

(4.16a) \begin{align} \left ( \nabla ^s \boldsymbol{v}_E - \frac {1}{d} (\textrm{div} \boldsymbol{v}_E) \boldsymbol{I}\right ):\left (\nabla ^s \boldsymbol{v}_E - \frac {1}{d} (\textrm{ div} \boldsymbol{v}_E) \boldsymbol{I}\right )=&\,0, \end{align}

(4.16b) \begin{align} \left (\bar\lambda _E + \frac {2}{d}\right )\left (\textrm{div} \boldsymbol{v}_E\right )^2 =&\,0, \end{align}

(4.16c) \begin{align} \displaystyle \sum _{\alpha ,{\beta }} \left (\nabla g_{\alpha ,E}\right )^T \boldsymbol{B}_{\alpha {\beta },E} \nabla g_{{\beta },E} =&\,0, \end{align}

(4.16d) \begin{align} \sum _{\alpha ,{\beta }} m_{\alpha {\beta },E}g_{\alpha ,E}g_{{\beta },E} =&\,0, \end{align}

in $\Omega$ , where $\boldsymbol{B}_{\alpha {\beta }}=\boldsymbol{M}_{\alpha {\beta }} + \boldsymbol{K}_{\alpha {\beta }}$ , and where the subscript $E$ denotes the equilibrium configuration of the quantity. We deduce from (4.16a )–(4.16b ) that $\boldsymbol{v}_E$ are rigid body motions. Simplifying the analysis, we take $\boldsymbol{v}_E = 0$ , which causes the inertia terms to vanish. Additionally, we assume the absence of gravitational forces ( $\boldsymbol{b} = 0$ ). Substituting into the momentum equation (4.7a ) provides

(4.17) \begin{align} \nabla \lambda _E + \sum _\alpha \phi _{\alpha ,E}\nabla \hat {\mu }_{\alpha ,E} = 0. \end{align}

Recalling (4.8), we deduce

(4.18) \begin{align} \sum _\alpha \tilde {\rho }_{\alpha ,E}\nabla g_{\alpha ,E} = 0. \end{align}

Next, from (4.16c ) we deduce that $\nabla g_{E}$ lies in the null space of $\boldsymbol{B}_{E}$ in the sense $\sum _{\beta } \boldsymbol{B}_{\alpha {\beta },E}\nabla g_{{\beta },E} = 0$ ( $\alpha = 1,\ldots ,N$ ), and hence $\hat {\boldsymbol{J}}_{\alpha ,E} + \hat {\boldsymbol{j}}_{\alpha ,E} = 0$ in equilibrium. In the special case $\boldsymbol{B}_{\alpha {\beta }}=-\boldsymbol{B}_0 \tilde {\rho }_\alpha \tilde {\rho }_{\beta }$ for $\alpha \neq {\beta }$ , and $\boldsymbol{B}_{\alpha \alpha } = \boldsymbol{B}_{0}\tilde {\rho }_\alpha \sum _{{\gamma }\neq \alpha } \tilde {\rho }_{\gamma }$ for some $\boldsymbol{B}_0$ that does not depend on $\tilde {\rho }_\alpha$ , $\alpha =1,\ldots ,N$ ; this coincides with (4.18). Similarly, (4.16d ) provides $\sum _{\beta } m_{\alpha {\beta },E} g_{{\beta },E} = 0$ , and hence $\hat {\zeta }_{\alpha ,E}=0$ ( $\alpha = 1,\ldots ,N$ ).

5.1. Binary-phase case

In this section, we restrict to binary mixtures ( $\alpha = 1,2)$ , and compare with the framework presented in ten Eikelder et al. (Reference ten Eikelder, van der Zee, Akkerman and Schillinger2023). A formulation of this two-phase modelling framework is

(5.1a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla \lambda + \phi \nabla \breve {\mu } - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(5.1b) \begin{align} \partial _t \rho + \textrm{div}(\rho \boldsymbol{v}) &=\,0. \end{align}

(5.1c) \begin{align} \partial _t \phi + \textrm{div}(\phi \boldsymbol{v}) - \textrm{div} \left ( \breve {\boldsymbol{M}} \nabla \left (\breve {\mu }+\omega \lambda \right )\right ) +\breve {m} \left (\breve {\mu }+\omega \breve {\lambda }\right ) &=\,0. \end{align}

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

(5.2) \begin{align} \phi = \phi _1 - \phi _2, \end{align}

where we recall $\phi _1+\phi _2=1$ . The chemical potential quantity is defined as

(5.3) \begin{align} \breve {\mu } = \frac {\partial \Psi }{\partial \phi } - \textrm{div}\frac {\partial \Psi }{\partial \nabla \phi }. \end{align}

Finally, the quantity $\omega$ is $\omega = (\rho _1^{-1}-\rho _2^{-1})/(\rho _1^{-1}+\rho _2^{-1})$ . On the other hand, the model (4.7) takes for binary mixtures the following form (where we directly enforce $\phi _1+\phi _2=1$ ):

(5.4a) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \nabla \lambda + \phi _1\nabla \hat {\mu }_1 + \phi _2\nabla \hat {\mu }_2 &\nonumber \\ - \textrm{div} \left ( \nu (2\nabla ^s \boldsymbol{v}+\bar {\lambda }(\textrm{div}\boldsymbol{v}) \boldsymbol{I}) \right )-\rho \boldsymbol{b} &=\, 0, \end{align}

(5.4b) \begin{align} \partial _t \phi _1 + \textrm{div}(\phi _1 \boldsymbol{v}) - \rho _1^{-1}\textrm{div} \left ( \boldsymbol{M} \nabla \left (g_1 -g_2\right )\right ) +\rho _1^{-1}m \left (g_1 -g_2\right ) &=\,0, \end{align}

(5.4c) \begin{align} \partial _t \phi _2 + \textrm{div}(\phi _2 \boldsymbol{v}) - \rho _2^{-1}\textrm{div} \left ( \boldsymbol{M} \nabla \left (g_2 -g_1\right )\right ) +\rho _2^{-1} m \left (g_2 -g_1\right ) &=\,0, \end{align}

where $\boldsymbol{M} = \boldsymbol{M}_{12}=\boldsymbol{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

(5.5a) \begin{align} \partial _t \phi \!+\! \textrm{div}(\phi \boldsymbol{v}) - \textrm{div} \left ( (\rho _1^{-1}\!+\!\rho _2^{-1})\boldsymbol{M} \nabla \left (g_1 -g_2\right )\right ) \!+\! (\rho _1^{-1}+\rho _2^{-1})m \left (g_1 -g_2\right ) &=0, \end{align}

(5.5b) \begin{align} \partial _t \rho \!+\! \textrm{div}(\rho \boldsymbol{v}) &=0. \end{align}

With the aim of comparing the two models, we select the relations $\breve {\boldsymbol{M}}=(\rho _1^{-1}+\rho _2^{-1})^2\boldsymbol{M}$ and $\breve {m}=(\rho _1^{-1}+\rho _2^{-1})^2m$ , which converts (5.5a ) into

(5.6) \begin{align} \partial _t \phi + \textrm{div}(\phi \boldsymbol{v}) - \textrm{div} \left ( \breve {\boldsymbol{M}} \nabla \left (\grave {\mu } +\omega \lambda \right )\right ) +\breve {m} \left (\grave {\mu }+\omega \lambda \right ) &=\,0, \end{align}

with $\grave {\mu } = (\rho _1^{-1}+\rho _2^{-1})^{-1}(\rho _1^{-1}\mu _1-\rho _2^{-1}\mu _2)$ . As a consequence, in case the identities

(5.7a) \begin{align} \phi _1 \nabla \hat {\mu }_1 + \phi _2 \nabla \hat {\mu }_2 =&\, \phi \nabla \breve {\mu }, \end{align}

(5.7b) \begin{align} \grave {\mu } =&\,\breve {\mu }, \end{align}

hold, we find that (5.4) coincides with (5.1). This is in general not the case, i.e. in general the two models do not match. (An $N$ -phase theory that reduces to existing two-phase models emerges when working with $N-1$ order parameters $\phi _\alpha$ , rather than the current case of $N$ order parameters $\phi _\alpha$ .) There are, however, specific situations in which the models coincide, for instance when $\mu _1 + \mu _2 = 0$ and $\breve {\mu } = \mu _1 = -\mu _2$ . These conditions are inspired by the chain rule for chemical potentials, where $\phi =\phi (\phi _1,\phi _2)=\phi _1-\phi _2$ so that $\partial \phi /\partial \phi _1 = 1$ , $\partial \phi /\partial \phi _2 = -1$ of (5.2).

5.2. $N$ -phase model Dong (Reference Dong2018)

The $N$ -phase incompressible model proposed by Dong (Reference Dong2018) is given by

(5.8a) \begin{align} \rho \left (\partial _t \boldsymbol{u} + \boldsymbol{u}\,{\boldsymbol\cdot}\, \nabla \boldsymbol{u} \right ) + \boldsymbol{J}' \,{\boldsymbol\cdot}\, \nabla \boldsymbol{u} + \nabla \lambda ' - \textrm{div}\left ( \nu ' \nabla ^s \boldsymbol{u} \right ) &\nonumber \\ + \sum _{\beta } \textrm{div}\left (\nabla \phi _\alpha \otimes \frac {\partial \Psi '}{\partial (\nabla \phi _{\beta })}\right ) =&\,0, \end{align}

(5.8b) \begin{align} \textrm{div}\boldsymbol{u} =&\,0, \end{align}

(5.8c) \begin{align} \partial _t \phi _\alpha + \boldsymbol{u}\,{\boldsymbol\cdot}\, \nabla \phi _\alpha - \sum _{\beta } \textrm{div}\left ( m_{\alpha {\beta }}' \nabla \left ( \frac {\partial \Psi '}{\partial \phi _{\beta }} - \textrm{div}\left (\frac {\partial \Psi '}{\partial (\nabla \phi _{\beta })}\right )\right ) \right ) =&\,0, \end{align}

for $\alpha = 1,\ldots , N$ , where $\lambda '$ is the Lagrange multiplier pressure, $\nu '$ is the dynamic viscosity, $\Psi '$ is the free energy, $\boldsymbol{J}'$ is the peculiar velocity, and $m_{\alpha {\beta }}'$ is the mobility. For the purpose of comparing the model (5.8) to the proposed framework, we define the chemical potential:

(5.9) \begin{align} \mu _{\beta }' = \frac {\partial \Psi '}{\partial \phi _{\beta }} - \textrm{div}\left (\frac {\partial \Psi '}{\partial (\nabla \phi _{\beta })}\right ). \end{align}

Invoking Lemma4.1, we rewrite the model (5.8) as

(5.10a) \begin{align} \rho \left (\partial _t \boldsymbol{u} + \boldsymbol{u}{\boldsymbol\cdot} \nabla \boldsymbol{u} \right ) + \boldsymbol{J}'{\boldsymbol\cdot} \nabla \boldsymbol{u} + \nabla \tilde {\lambda }' + \sum _{\beta } \phi _\alpha \nabla \mu _\alpha ' - \textrm{div}\left ( \nu ' \nabla ^s \boldsymbol{u} \right ) =&\,0, \end{align}

(5.10b) \begin{align} \textrm{div}\boldsymbol{u} =&\,0, \end{align}

(5.10c) \begin{align} \partial _t \phi _\alpha + \boldsymbol{u}{\boldsymbol\cdot} \nabla \phi _\alpha - \sum _{\beta } \textrm{div}\left ( m_{\alpha {\beta }}' \nabla \mu _{\beta }' \right ) =&\,0, \end{align}

with

(5.11) \begin{align} \tilde {\lambda }' = \lambda ' + \Psi ' - \sum _{\beta } \phi _\alpha \mu _\alpha '. \end{align}

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

1. (i) Model (5.10) does not contain each of the peculiar velocity terms in the momentum equation; this applies to both inertia and viscous terms.

2. (ii) Model (5.10) does not include mass transfer terms.

3. (iii) The constitutive model for the diffusive flux in (5.10) is different; in particular the Lagrange multiplier is absent. As a consequence, the equilibrium conditions are different.

5.3. Class-II mixture model

We compare the proposed unified modelling framework with an incompressible mixture model presented in ten Eikelder et al. (Reference ten Eikelder, van der Zee and Schillinger2024):

(5.12a) \begin{align} \partial _t \tilde {\rho }_\alpha + \textrm{div}(\tilde {\rho }_\alpha \boldsymbol{v}_\alpha ) + \sum _{\beta } \breve {m}_{\alpha {\beta }}(\breve {g}_\alpha -\breve {g}_{\beta }) &=\, 0, \end{align}

(5.12b) \begin{align} \partial _t (\tilde {\rho }_\alpha \boldsymbol{v}_\alpha ) + \textrm{div} \left ( \tilde {\rho }_\alpha \boldsymbol{v}_\alpha \otimes \boldsymbol{v}_\alpha \right ) + \phi _\alpha \nabla \left (\breve {\lambda } + \breve {\mu }_\alpha \right ) \nonumber \\ - \textrm{div}\left (\breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{v}_\alpha + \breve {\bar {\lambda }}_\alpha \textrm{ div}\boldsymbol{v}_\alpha \right )\right ) -\tilde {\rho }_\alpha \boldsymbol{b}&\nonumber \\ +\displaystyle \sum _{{\beta }} R_{\alpha {\beta }} (\boldsymbol{v}_\alpha -\boldsymbol{v}_{\beta }) + \tfrac {1}{2}\sum _{\beta } \breve {m}_{\alpha {\beta }}(\breve {g}_\alpha -\breve {g}_{\beta })(\boldsymbol{v}_\alpha +\boldsymbol{v}_{\beta }) &=\, 0, \end{align}

for constituents $\alpha = 1,\ldots ,N$ . Here $\boldsymbol{v}_\alpha$ is the constituent velocity, $\breve {\lambda }$ is a Lagrange multiplier, $\tilde {\nu }_\alpha$ the constituent dynamical viscosity, $\breve {\bar {\lambda }}_\alpha \geqslant 2/d$ , $\nabla ^s \boldsymbol{v}_\alpha$ the constituent symmetric velocity gradient, and $\breve {m}_{\alpha {\beta }}$ and $R_{\alpha {\beta }}$ are symmetric matrices (for the properties see ten Eikelder et al. Reference ten Eikelder, van der Zee and Schillinger2024). This model considers the free energy class:

(5.13a) \begin{align} \Psi =&\, \sum _\alpha \breve {\Psi }_\alpha , \end{align}

(5.13b) \begin{align} \breve {\Psi }_\alpha =&\, \breve {\Psi }_\alpha \left (\phi _\alpha , \nabla \phi _\alpha \right ). \end{align}

The associated constituent chemical potentials are defined as

(5.14) \begin{align} \breve {\mu }_\alpha =&\, \frac { \partial \breve {\Psi }_\alpha }{\partial \phi _\alpha } - \textrm{div}\frac {\partial \breve {\Psi }_\alpha }{\partial \nabla \phi _\alpha }, \end{align}

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

Inserting the class (5.13) into the proposed modelling 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 modelling framework presented in the current paper, this model comprises $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.12a ) can be written as

(5.15a) \begin{align} \partial _t \phi _\alpha + \textrm{div}(\phi _\alpha \boldsymbol{v}) + \rho _\alpha ^{-1}\textrm{div}\boldsymbol{J}_\alpha + \rho _\alpha ^{-1}\breve {\gamma }_\alpha &=\, 0, \end{align}

(5.15b) \begin{align} \breve {\gamma }_\alpha - \sum _{\beta } \breve {m}_{\alpha {\beta }}(\breve {g}_\alpha -\breve {g}_{\beta })&=\,0. \end{align}

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

(5.16) \begin{align} \partial _t (\rho \boldsymbol{v}) + \textrm{div} \left ( \rho \boldsymbol{v}\otimes \boldsymbol{v} \right ) + \sum \phi _\alpha \nabla \left (\breve {\lambda } + \breve {\mu }_\alpha \right )& \nonumber \\ - \textrm{div}\left (\sum _\alpha \breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{v} + \breve {\bar {\lambda }}_\alpha \textrm{div}\boldsymbol{v}\right )\right ) -\rho \boldsymbol{b}&\nonumber \\ - \textrm{div}\left (\sum _\alpha \breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{w}_\alpha + \breve {\bar {\lambda }}_\alpha \textrm{ div}\boldsymbol{w}_\alpha \right )- \sum _\alpha \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha \right )&\,=0, \end{align}

where we have adopted the identities

(5.17a) \begin{align} \sum _\alpha \tilde {\rho }_\alpha \boldsymbol{v}_\alpha \otimes \boldsymbol{v}_\alpha =&\, \rho \boldsymbol{v} \otimes \boldsymbol{v} + \sum _\alpha \tilde {\rho }_\alpha \boldsymbol{w}_\alpha \otimes \boldsymbol{w}_\alpha , \end{align}

(5.17b) \begin{align} \sum _\alpha \breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{v}_\alpha + \breve {\bar {\lambda }}_\alpha \textrm{div}\boldsymbol{v}_\alpha \right ) =&\, \sum _\alpha \breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{v} + \breve {\bar {\lambda }}_\alpha \textrm{div}\boldsymbol{v}\right ) \nonumber \\ &\,+ \sum _\alpha \breve {\nu }_\alpha \left (2\nabla ^s \boldsymbol{w}_\alpha + \breve {\bar {\lambda }}_\alpha \textrm{div}\boldsymbol{w}_\alpha \right ). \end{align}

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.16) consists of terms that are absent in (4.7a ). This is a direct consequence of energy-dissipation law (3.3) and the introduction of the model $\boldsymbol{J}_\alpha =\hat {\boldsymbol{J}}_\alpha$ in (3.30b ). In the case of matching constitutive velocities, as described in § 4.2, these terms vanish.