Finite element methods for multicomponent convection-diffusion

Abstract

We develop finite element methods for coupling the steady-state Onsager–Stefan–Maxwell (OSM) equations to compressible Stokes flow. These equations describe multicomponent flow at low Reynolds number, where a mixture of different chemical species within a common thermodynamic phase is transported by convection and molecular diffusion. Developing a variational formulation for discretizing these equations is challenging: the formulation must balance physical relevance of the variables and boundary data, regularity assumptions, tractability of the analysis, enforcement of thermodynamic constraints, ease of discretization and extensibility to the transient, anisothermal and nonideal settings. To resolve these competing goals, we employ two augmentations: the first enforces the definition of mass-average velocity in the OSM equations, while its dual modifies the Stokes momentum equation to enforce symmetry. Remarkably, with these augmentations we achieve a Picard linearization of symmetric saddle point type, despite the equations not possessing a Lagrangian structure. Exploiting structure mandated by linear irreversible thermodynamics, we prove the inf-sup condition for this linearization, and identify finite element function spaces that automatically inherit well-posedness. We verify our error estimates with a numerical example, and illustrate the application of the method to nonideal fluids with a simulation of the microfluidic mixing of hydrocarbons.

finite element methods, compressible Stokes equations, Stefan–Maxwell equations, multicomponent diffusion, linear irreversible thermodynamics, Stokes–Onsager–Stefan–Maxwell equations

1. Introduction

Many fluids consist of mixtures; for example, air is a mixture of nitrogen, oxygen, carbon dioxide and other species. In many situations, it is not necessary to resolve the motions of the individual species, such as when modelling the flow of air over an aircraft. However, in other contexts, detailed knowledge of the transport of individual species is required. Examples include biological applications, where one may be interested in the transport of oxygen and carbon dioxide in blood, in chemical engineering, where one may be interested in separating or combining the constituents of petroleum, or in electrochemistry, where the performance of a lithium-ion battery is often limited by the transport of lithium ions within an electrolyte. We describe this situation as a multicomponent flow, where a fluid is composed of |$ 2 \leq n \in \mathbb{N}_{+}$| distinct chemical species in a common thermodynamic phase. The primary contributions of our work are a novel variational formulation of a system of equations describing nonideal miscible isothermal multicomponent flow, a Picard linearization of these equations that possesses symmetric saddle point structure (despite the equations not arising from a Lagrangian) and the numerical analysis of a structure-preserving finite element discretization. In particular, we identify the structural relationships required of the finite element spaces for the different variables that allow for the continuous well-posedness to be inherited automatically in the discretization.

More specifically, this work considers multicomponent flow in the concentrated (i.e. general) solution regime, as opposed to the simpler and more commonly studied dilute solution regime. We now provide an overview of this distinction, for self-containment. The dilute approximation applies when a single species called the solvent (conventionally assigned index |$i = n$|⁠) is taken to have a concentration very far in excess of the remaining species (⁠|$i < n$|⁠), each of which is called a solute; for example, a classical problem in computational fluid dynamics concerns the tracking of tracers, present in small proportions in a solvent by which they are convected. Thus, if |$\varOmega \subset \mathbb{R}^{d}$| (⁠|$d\in \{2, 3\}$|⁠) is the domain containing the mixture, the fluid density |$\rho :\varOmega \to \mathbb{R}$| (⁠|$\textrm{kg}/\textrm{m}^{3}$|⁠) varies negligibly with solute content in a dilute solution and approximately coincides with the mass density of the pure solvent. This decouples the flow, the material’s bulk motion, from the mass transport, the motion of individual molecular constituents comprising the material. One can thus solve for the flow velocity, and then employ this velocity in a system of independent advection-diffusion equations for the mass transport of each species. A typical dilute solution problem—which, we emphasize, we do not consider in this paper—at low Mach number (where the density is assumed constant) is to solve:

$$ \begin{alignat}{3} \frac{\partial (\rho v)}{\partial t} -\operatorname{div}\left( 2 \eta \varepsilon(v) \right) + \operatorname{div}\left( \rho v \otimes v \right) + \nabla p &= \rho f, \end{alignat} $$

(1.1a)

$$ \begin{alignat}{3} \operatorname{div} v &= 0, \end{alignat} $$

(1.1b)

$$ \begin{alignat}{3} \frac{\partial c_{i}}{\partial t} + \operatorname{div}\left(c_{i} v\right) + \operatorname{div} J_{i} &= r_{i}, \quad && i = 1, \ldots, n - 1, \end{alignat} $$

(1.1c)

$$ \begin{alignat}{3} J_{i} &= -D_{i} \nabla c_{i}, \quad && i = 1, \ldots, n - 1, \end{alignat} $$

(1.1d)

where |$\varOmega $| is bounded Lipschitz, |$v: \varOmega \to \mathbb{R}^{d}$| is the flow velocity (⁠|$\textrm{m}/\textrm{s}$|⁠), |$\varepsilon $| the symmetric gradient operator |$\varepsilon (v) := \frac{1}{2} (\nabla v + (\nabla v)^{\top } )$||$(1/\textrm{s})$|⁠, |$\eta> 0$| the shear viscosity (⁠|$\textrm{Pa}\cdot \textrm{s}$|⁠), |$p: \varOmega \to \mathbb{R}$| the pressure (Pa), |$f: \varOmega \to \mathbb{R}^{d}$| the body acceleration (⁠|$\textrm{m}/\textrm{s}^{2}$|⁠) induced within |$\varOmega $| due to the action of external fields, |$c_{i}: \varOmega \to \mathbb{R}$| the molar concentration (⁠|$\textrm{mol}/\textrm{m}^{3}$|⁠) of solute |$i$| in the solution, |$J_{i}: \varOmega \to \mathbb{R}^{d}$| its diffusive flux (⁠|$\textrm{mol}/\textrm{m}^{2}\textrm{s}$|⁠), |$r_{i}: \varOmega \to \mathbb{R}$| its volumetric rate of generation or depletion (⁠|$\textrm{mol}/\textrm{m}^{3}\textrm{s})$|⁠, and |$D_{i}> 0$| its Fickian diffusion coefficient (⁠|$\textrm{m}^{2}/\textrm{s}$|⁠). The velocity |$v_{i}: \varOmega \to \mathbb{R}^{d}$| (⁠|$\textrm{m}/\textrm{s}$|⁠) of each individual species is given by

$$ \begin{align}& c_{i} v_{i} = c_{i} v + J_{i},\end{align} $$

(1.2)

decomposing the transport of each species into a convective and a diffusive contribution. The dilute solution regime is characterized by the approximation |$v \approx v_{n}$|⁠, that the bulk motion of the fluid coincides with the motion of the solvent.

While the dilute solution approximation has been applied to great effect (Levich, 1962; Bird et al., 2002; Cussler, 2009; Deen, 2016), it fails starkly when no particular species is present in great excess, the concentrated solution regime of interest in this work. Several problems arise when attempting to relax the dilute solution approximation and formulate and discretize models for concentrated solutions. In concentrated solutions the very notion of ‘flow velocity’ becomes ambiguous, because the bulk motion of the fluid need not coincide with any particular species velocity and these species velocities in general will be distinct. One can still identify a natural composition-dependent definition of |$v$| in the concentrated case, however. The density of the fluid is given by

$$ \begin{align}& \rho := \sum_{i=1}^{n} M_{i} c_{i},\end{align} $$

(1.3)

in which |$M_{i}> 0$| represents the molar mass (⁠|$\textrm{kg}/\textrm{mol}$|⁠) of species |$i$|⁠. Using (1.2), the continuity equations (1.1c) may be rephrased in terms of species velocities as

$$ \begin{align}& \frac{\partial c_{i}}{\partial t} + \operatorname{div}(c_{i} v_{i}) = r_{i}, \quad i = 1, \ldots, n.\end{align} $$

(1.4)

Time differentiation of (1.3), followed by elimination of the concentration derivatives with (1.4), yields

$$ \begin{align}& \frac{\partial \rho}{\partial t} = \sum_{i} M_{i} \left( r_{i} - \operatorname{div} \left( c_{i} v_{i} \right) \right) = -\operatorname{div} \left( \sum_{i} M_{i} c_{i} v_{i} \right).\end{align} $$

(1.5)

Here, the last equality incorporates the premise that homogeneous chemical reactions conserve atoms, which requires that |$\sum _{i} M_{i} r_{i} = 0$|⁠.¹ Equation (1.5) is consistent with the common understanding of mass continuity if the flow velocity |$v$| within a multicomponent fluid is identified as the so-called mass-average velocity, defined as (Hirschfelder et al., 1954, p. 454)

$$ \begin{align}& v := \sum_{i} \omega_{i} v_{i},\end{align} $$

(1.6)

a convex combination in which

$$ \begin{align}& \omega_{i} := \frac{M_{i} c_{i}}{\rho}\end{align} $$

(1.7)

defines the (dimensionless) mass fraction of species |$i$|⁠. Indeed, rewriting (1.5) in terms of the mass-average velocity yields

$$ \begin{align}& \frac{\partial \rho}{\partial t} = - \operatorname{div} (\rho v),\end{align} $$

(1.8)

thereby recovering the mass continuity equation familiar from fluid mechanics. In our formulation we will solve for both the mass-average velocity and the individual species velocities. The mass-average velocity is governed by a momentum balance, typically expressed in the form of the Cauchy equation

$$ \begin{align}& \frac{\partial (\rho v)}{\partial t} - \operatorname{div}\tau + \operatorname{div}\left(\rho v\otimes v\right) + \nabla p = \rho f,\end{align} $$

(1.9)

where |$\tau : \varOmega \to \mathbb{R}^{d \times d}_{\textrm{sym}}$| is the dissipative (viscous) internal stress tensor (Pa),² to be specified with a constitutive law. If Newton’s law of viscosity is used, then (1.9) reduces to the Navier–Stokes momentum equation (1.1a).

Another issue to address when moving to the concentrated solution regime is the choice of constitutive law for the diffusive fluxes. In dilute solutions each solute interacts at a molecular level almost solely with solvent molecules, and so the diffusive solute fluxes |$J_{i}$| can each be modelled by Fick’s law (Fick, 1855). In concentrated solutions the constitutive laws for mass transport become incomplete, because Fick’s law (1.1d) fails to take into account all possible species-species interactions. Even in the case of simple diffusion (where |$v = 0$| uniformly), the diffusive flux of a given species can generally be driven by a concentration gradient of any other species in the solution—a phenomenon known as cross-diffusion. The theory of linear irreversible thermodynamics, pioneered by Onsager (1931a,b, 1945), enables the thermodynamically consistent generalization of Fick’s law (1.1d) to the concentrated solution regime. This formalism is described in the next subsection.

1.1 Onsager–Stefan–Maxwell equations

Within a multi-species solution, the Onsager–Stefan–Maxwell (OSM) equations relate the diffusion driving forces|$d_{i}: \varOmega \to \mathbb{R}^{d}$| to the species velocities |$v_{i}$| via

$$ \begin{align}& d_{i} = \sum_{j} \textbf{M}_{ij} v_{j}, \quad i = 1,\ldots, n,\end{align} $$

(1.10)

where the diffusion driving forces |$d_{i}$| (⁠|$\textrm{Pa}/\textrm{m}$|⁠) incorporate the effects of various state variable gradients (Giovangigli, 1999, Eq. (2.5.4)), and where |$\textbf{M}$| is the Onsager transport coefficient matrix with entries

$$ \begin{align}& \textbf{M}_{ij} = \left\{ \renewcommand{\arraystretch}{1.2} \begin{array}{ l@{~} l} -\frac{RT c_{i} c_{j}}{\mathscr{D}_{ij}c_{\textrm{T}}} &\ \textrm{if}\ i \neq j, \\ \sum_{k=1, k \neq i}^{n}\frac{RTc_{i} c_{k}}{\mathscr{D}_{ik}c_{\textrm{T}}} &\ \textrm{if}\ i = j. \end{array} \right.\end{align} $$

(1.11)

Here |$R> 0$| is the ideal gas constant (⁠|$\textrm{J}/\textrm{mol}\cdot \textrm{K}$|⁠), |$T> 0$| the ambient temperature (K), |$c_{\textrm{T}}$| denotes the total concentration defined as

$$ \begin{align}& c_{\textrm{T}}: = \sum_{i} c_{i},\end{align} $$

(1.12)

and |$\mathscr{D}_{ij} \in \mathbb{R}$| represents the Stefan–Maxwell diffusivity (⁠|$\textrm{m}^{2}/\textrm{s}$|⁠) of species |$i$| through species |$j\neq i$|⁠. The Stefan–Maxwell diffusivities are symmetric in the species indices, |$\mathscr{D}_{ij} = \mathscr{D}_{ji}$|⁠, and |$\mathscr{D}_{ii}$| is not defined. In the present discussion we restrict attention to the case where every |$\mathscr{D}_{ij}$| is constant, which in turn requires each to be positive (Van-Brunt et al., 2021), but in general the Stefan–Maxwell diffusivities may depend on the species concentrations, temperature and pressure. (1.10) is often presented as

$$ \begin{align}& d_{i} = \sum_{\substack{j=1 \\ j\neq i}}^{n} \frac{RT c_{i} c_{j}}{ \mathscr{D}_{ij} c_{\textrm{T}}} ({v}_{i} - {v}_{j}),\end{align} $$

(1.13)

which follows from (1.10) and (1.11).

In an isothermal but nonisobaric fluid, the diffusion driving forces may be identified as (Bird et al., 2002, Eq. (24.1-8))

$$ \begin{align}& d_{i} := -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p, \quad i = 1, \ldots, n,\end{align} $$

(1.14)

in which |$\mu _{i}: \varOmega \to \mathbb{R}$| is the chemical potential (⁠|$\textrm{J}/\textrm{mol}$|⁠) of species |$i$|⁠. The chemical potential represents the partial derivative of the Gibbs free energy (a quantity describing the total amount of work a system can deliver to isothermal, isobaric surroundings) with respect to the number of moles of a given species |$i$| at constant temperature and pressure, and holding the molar contents of all other species fixed. Under isothermal conditions, chemical potentials are related to the concentrations and pressure via a constitutive law, discussed below in Section 1.4.

As a consequence of the statistical reciprocal relations developed by Onsager (1931a,b) and the second law of thermodynamics, the transport matrix |$\textbf{M}$| is everywhere symmetric positive semidefinite (Monroe et al., 2015). Thus, the Stefan–Maxwell diffusivities are symmetric in the species indices, |$\mathscr{D}_{ij} = \mathscr{D}_{ji}$|⁠. There is a single null eigenvalue, with eigenvector |$(1, \ldots , 1)^{\top }$|⁠:

$$ \begin{align}& \sum_{j} \textbf{M}_{ij} = 0.\end{align} $$

(1.15)

This nullspace is necessary to distinguish convection, which is nondissipative, from diffusion. A consequence of this nullspace is that one may shift each |$v_{i}$| in (1.10) by the mass-average velocity,

$$ \begin{align}& d_{i} = \sum_{j} \textbf{M}_{ij} \left(v_{j} - v\right), \quad i = 1, \ldots, n,\end{align} $$

(1.16)

so that the transport matrix acts on terms proportional to the diffusive flux |$J_{i}$|⁠. Thus (1.10) can be understood as an implicit constitutive relation for the diffusive fluxes (Bulíček et al., 2021).

The symmetry of |$\textbf{M}$| combined with (1.15) allows one to show that

$$ \begin{align}& \sum_{i} d_{i} = \sum_{i} \left(-c_{i} \nabla \mu_{i} + \omega_{i}\nabla p\right) = 0,\end{align} $$

(1.17)

an expression of the isothermal, nonisobaric Gibbs–Duhem equation from equilibrium thermodynamics. For further detail on the historical development and mathematical structure of the OSM equations, we refer the reader to Van-Brunt et al. (2021).

1.2 Stokes equation

A further consequence of the nullspace of |$\textbf{M}$| (1.15) is that the species velocities cannot be determined from the diffusion driving forces alone. They can be recovered, however, by incorporating the Cauchy momentum equation (1.9) to specify the mass-average velocity (1.6), which is required for a complete description of the overall transport problem; this in turn requires a constitutive law for the viscous stress. For isothermal Newtonian fluids, the viscous stress |$\tau $| relates to the linearized strain rate |$\varepsilon (v)$|⁠, the symmetric part of the velocity gradient, through

$$ \begin{align}& \tau = 2\eta\left(\varepsilon(v) - \frac{\operatorname{tr}\varepsilon(v)}{d}\mathbb{I}\right) + \zeta\operatorname{tr}(\varepsilon(v))\mathbb{I},\end{align} $$

(1.18)

where |$\zeta> 0$| is the bulk viscosity (⁠|$\textrm{Pa}\cdot \textrm{s}$|⁠), or equivalently

$$ \begin{align}& \varepsilon(v) = \frac{1}{2\eta}\tau + \left(\frac{1}{d^{2}\zeta} - \frac{1}{2\eta d}\right)(\operatorname{tr}\tau)\mathbb{I} =: \mathscr{A}\tau,\end{align} $$

(1.19)

where |$\mathscr{A}:\mathbb{R}^{d\times d}_{\textrm{sym}}\to \mathbb{R}^{d\times d}_{\textrm{sym}}$| denotes the compliance tensor. The full Cauchy stress |$\sigma : \varOmega \to \mathbb{R}^{d \times d}_{\textrm{sym}}$| (Pa) may then be decomposed as

$$ \begin{align}& \sigma = \tau - p\mathbb{I}.\end{align} $$

(1.20)

We further consider steady-state creeping flow, under which assumptions Stokes’ equation,

$$ \begin{align}& \operatorname{div}\tau - \nabla p = -\rho f,\end{align} $$

(1.21)

follows from the momentum balance (1.9).

The OSM equations (1.10) are written in force-explicit form: the equations express the species velocities (fluxes) as implicit variables. We choose also to write the Newtonian constitutive equation (1.19) in this manner, expressing the thermodynamic force (the linearized strain rate) in terms of the corresponding flux (the viscous stress) (Hirschfelder et al., 1954). Typically in computational fluid dynamics, a flux-explicit formulation is obtained by using an explicit constitutive relation such as (1.18) to eliminate the Cauchy and viscous stresses in the first instance. For our overall coupled system (stated later in (1.29)), we do not eliminate the viscous stress, but include it as an implicit variable to be solved for. While this choice increases the computational cost, it has substantial benefits; the viscous stress plays a fundamental role in the calculation of local entropy production, but more significantly, we show in Section 2 that the resulting system of equations can be cast as a symmetric perturbed saddle point-like system, which is conducive to both theoretical analysis and (we anticipate) efficient linear solvers.

1.3 Augmentation of the diffusion transport matrix and the Stokes momentum balance

The variational formulation of our equations must enforce the relation (1.6), between the bulk (mass-average) velocity and the species velocities. We employ the augmentation approach introduced by Helfand (1960) and later used by Giovangigli (1990); Ern & Giovangigli (1994); Giovangigli (1999). We augment each OSM equation (1.13) by adding the mass-average velocity constraint (1.6) to both sides with a prefactor |$\gamma> 0$|⁠:

$$ \begin{align}& d_{i} + \gamma\omega_{i} v = \sum_{j\neq i} \frac{RT c_{i} c_{j}}{ \mathscr{D}_{ij} c_{\textrm{T}}} ({v}_{i} - {v}_{j}) + \gamma\omega_{i} \sum_{j} \omega_{j} {v}_{j}.\end{align} $$

(1.22)

The prefactor |$\gamma $| formally has units of Pa|$\cdot $|s|$/\textrm{m}^{2}$|⁠. For convenience, we define the augmented transport matrix

$$ \begin{align}& \textbf{M}^{\gamma}_{ij} := \textbf{M}_{ij} + \gamma\omega_{i} \omega_{j},\end{align} $$

(1.23)

so that (1.22) can be stated as |$d_{i} + \gamma \omega _{i} v = \sum _{j} \textbf{M}_{ij}^{\gamma } v_{j}$|⁠. We may then compute directly

$$ \begin{align}& \sum_{i,j} {v}_{i} \cdot \textbf{M}^{\gamma}_{ij} {v}_{j} = \frac{1}{2} \sum_{i,j \atop j \neq i} \frac{ RT c_{i}c_{j}}{ \mathscr{D}_{ij}c_{\textrm{T}}} \vert{v}_{i} - {v}_{j}\vert^{2} + \gamma\left\vert \sum_{j} \omega_{j} {v}_{j} \right\vert^{2}\end{align} $$

(1.24)

to show that the augmented transport matrix is symmetric positive definite (Van-Brunt et al., 2021). This was used in Van-Brunt et al. (2021) to construct coercive bilinear forms for the pure Stefan–Maxwell diffusion problem. With this augmentation, although the transport matrix is a priori singular, one can nevertheless recover the species velocities from the driving forces by coupling with the mass-average velocity constraint (1.6).

The augmentation (1.22) modifies a constitutive law of the system, which will induce coercivity of a certain bilinear form below. However, this comes at the cost of symmetry. To recover symmetry, we add a ‘dual’ augmentation to the Stokes equation (1.21)

$$ \begin{align}& \operatorname{div}\tau - \nabla p = -\rho f + \gamma \sum_{j} \omega_{j}(v - v_{j}).\end{align} $$

(1.25)

With these two augmentations (1.22) and (1.25), an important bilinear form defined later in (3.3a) will be both symmetric and coercive on an appropriate kernel. This greatly aids the proofs of well-posedness for the continuous and discrete problems, as we will demonstrate in Section 3.2 and Section 4.1.

To close the equations, we must relate the chemical potentials to the concentrations via a thermodynamic equation of state, discussed next.

1.4 The chemical potential and equation of state

Our variational formulation will solve for the chemical potential |$\mu _{i}$| of each species |$i$|⁠. This has several advantages. First, this allows for a general formula for the diffusion driving forces (1.14), independent of the materials considered. If we were to make the (perhaps more obvious) choice of solving for concentrations |$c_{i}$| as the primary variables instead, the form of the diffusion driving forces would change in a material-dependent manner. Second, our choice allows for a decoupling in the linearization we employ: the primary mixed system to solve only depends on the material via the diffusion coefficients and viscosities, with any nonideality confined to the computation of concentrations and density postprocessed at every iteration using material-dependent thermodynamic constitutive relations discussed below. Third, together with the choice to solve for the viscous stress as described in Section 1.2, this decoupling endows the equations to solve with a symmetric perturbed saddle point-like structure.³

Generally each species concentration |$c_{i}$| can be inferred from |$\{\mu _{i}\}_{i=1}^{n}$| and |$p$|⁠, given thermodynamically consistent constitutive laws for the chemical potential, and an equation of state which relates |$c_{\textrm{T}}$| to pressure and composition. Within an isothermal ideal gas, the former relation is simply

$$ \begin{align}& c_{i} = \frac{p^{\ominus}}{RT}\exp\left(\frac{\mu_{i} - \mu^{\ominus}_{i}}{RT}\right),\end{align} $$

(1.26)

for some known reference pressure |$p^{\ominus }$| and a set of reference chemical potentials |$\{\mu ^{\ominus }_{i}\}_{i}$|⁠. A general relation for nonideal systems is

$$ \begin{align}& \mu_{i} = \mu_{i}^{\ominus} + RT \ln (\gamma_{i} x_{i}),\end{align} $$

(1.27)

where |$x_{i} := c_{i}/c_{\textrm{T}}$| is the mole fraction, and |$\gamma _{i}$| the (dimensionless) activity coefficient, of species |$i$|⁠. (Within a system made up of |$n$| species, specifying |$n - 1$| mole fractions determines the composition referred to earlier.) In nonideal solutions, the reference potentials |$\mu ^{\ominus }_{i}$| generally depend on the temperature and pressure (Guggenheim, 1985; Atkins & de Paula, 2010); they determine the value of the molar Gibbs free energy of pure species |$i$| at the |$T$| and |$p$| values of interest. Activity coefficients generally depend on temperature, pressure and composition; the definition of the reference state further requires that they approach unity at infinite dilution, i.e. |$\lim _{x_{i} \rightarrow 0} \gamma _{i} = 1$|⁠. Constitutive laws (1.27) suffice to determine the mole fractions within nonideal solutions. To obtain the concentrations, elementary thermodynamic principles can be used to derive a state equation for volume, which may be expressed as

$$ \begin{align}& c_{\textrm{T}} = \frac{1}{\sum_{i} V_{i} x_{i}},\end{align} $$

(1.28)

in which |$V_{i}> 0$| is the partial molar volume (⁠|$\textrm{m}^{3}/\textrm{mol}$|⁠) of species |$i$|⁠. Formally, the partial molar volume is a material property that quantifies the change in total fluid volume with respect to the number of moles of a species |$i$| at constant temperature and pressure, holding all other species contents fixed. Maxwell relations derived from the Gibbs free energy also require that |$V_{i}$| quantifies the partial derivative of |$\mu _{i}$| with respect to |$p$| (Goyal & Monroe, 2017, Eq. (28)). This dependence, part of which is embedded in the pressure dependence of |$\mu _{i}^{\ominus }$|⁠, may be regarded as given data that is experimentally measurable (Doyle & Newman, 1997). Note that for a definition of the Gibbs free energy in terms of chemical potentials to be thermodynamically consistent, it must imply the state equation (1.28).

Our linearization below is designed so that the concentrations are calculated from the chemical potentials and pressure. This trivially guarantees positivity of the concentrations, but more significantly, the model is able to incorporate nonideality by employing chemical potential constitutive laws more general than (1.26), such as (1.27). We intend for this choice to facilitate future research into convection-diffusion problems with alternative equations of state, for which |$c_{i}$| may (for example) depend on all|$\mu _{j}$|⁠, and on temperature.

1.5 Coupled problem statement

Our goal is to find and analyze a variational formulation and structure-preserving finite element discretization of the following problem: given data |$f$| and |$\{r_{i}\}_{i=1}^{n}$|⁠, find chemical potentials |$\{\mu _{i}\}_{i=1}^{n}$|⁠, viscous stress |$\tau $|⁠, pressure |$p$|⁠, species velocities |$\{v_{i}\}_{i=1}^{n}$| and convective velocity |$v$|⁠, satisfying

$$ \begin{align} -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p + \gamma\omega_{i} v &= \sum_{j} \textbf{M}^{\gamma}_{ij} v_{j}~\forall i, && \text{(augmented OSM equations)} \end{align} $$

(1.29a)

$$ \begin{align} \varepsilon(v) &= \mathscr{A} \tau, && \text{(viscous stress constitutive law)} \end{align} $$

(1.29b)

$$ \begin{align} \operatorname{div} \tau - \nabla p -\gamma \sum_{j} \omega_{j}(v - v_{j}) &= -\rho f, && \text{(augmented Stokes equation)} \end{align} $$

(1.29c)

$$ \begin{align} \operatorname{div} (c_{i} v_{i}) &= r_{i}~\forall i, && \text{(species continuity equation)} \end{align} $$

(1.29d)

$$ \begin{align} \operatorname{div}(v) &= \operatorname{div} \left( \sum_{j} \omega_{j} v_{j} \right)\!, && \text{(mass-average velocity constraint)} \end{align} $$

(1.29e)

for an augmentation parameter |$\gamma \geq 0$|⁠, where |$\{c_{i}, \omega _{i}\}_{i=1}^{n}, \rho $| are functions of the unknowns via thermodynamic constitutive laws such as (1.27) and (1.28), and algebraic relations (1.3), (1.7). We shall introduce appropriate boundary conditions in Section 2.2. We call the system (1.29) the (augmented) Stokes–Onsager–Stefan–Maxwell (SOSM) system.⁴ When the convection term |$\operatorname{div} \left ( \rho v \otimes v \right )$| is incorporated into (1.29c), we call this the Navier–Stokes–Onsager–Stefan–Maxwell (NSOSM) system.

Note that in the system we only directly enforce the divergence of the mass-average velocity constraint (1.29e), which may be interpreted as the compressible generalization of the standard divergence constraint (1.1b); this choice gives rise to a saddle point-like structure, as we show in the next section. Nevertheless, the full constraint (1.6) is incorporated via the augmentations (1.22) and (1.25), as discussed further in Remark 2.2. This constraint reduction, combined with the augmentations, may be regarded as a principal novelty of the system (1.29).

1.6 Relation to existing literature and outline

For dilute solutions with constant solvent concentration (and no volumetric generation or depletion of the solvent, |$r_{n} = 0$|⁠), the (N)SOSM equations reduce to the incompressible (Navier–)Stokes equations, as well as convection-diffusion equations constituted by Fick’s law for each solute. These equations have been studied for many decades, and effective numerical techniques are available. We do not attempt a systematic review here, but mention Thomée (2006); Hundsdorfer & Verwer (2013); Elman et al. (2014); Stynes & Stynes (2018) as gateways to this literature. In this regime, the momentum solve and the equation for the transport of concentration are decoupled using incompressibility.

Our formulation (1.29) solves for the viscous stress as an unknown variable. Of most relevance to this aspect of our approach is the work of Carstensen et al. (2012), who discretized the stress in an incompressible stress-velocity Stokes system using the same stress elements of Arnold & Winther (2002) that we shall employ.

Related systems of equations have been formulated and analyzed, including the coupling of the Stefan–Maxwell equations with the incompressible Navier–Stokes equations by Chen & Jüngel (2015), the compressible Navier–Stokes equations by Bothe & Druet (2021), the Darcy momentum equation by Ostrowski & Rohde (2020) and the Cahn–Hilliard equations by Huo et al. (2022). The coupling of an anisothermal NSOSM system to surface phenomena by sorption was formulated by Bothe & Dreyer (2015); Souček et al. (2019) proposed an extension to the Stefan–Maxwell equations in the presence of chemically reacting constituents.

Numerical methods for solving the NSOSM equations have received much less attention. The only works of which the authors are aware are those of Ern, Giovangigli and coauthors, including a monograph (Ern & Giovangigli, 1994) and a series of other works (Ern & Giovangigli, 1998, 1999; Giovangigli, 1999; Burman et al., 2004) which apply multicomponent transport to combustion modelling for ideal gas mixtures. These schemes use sophisticated finite difference methods, with the important exception of Burman et al. (2004), which uses a finite element method with additional least-squares terms to stabilize the formulation. The authors are unaware of any literature that addresses numerical methods for SOSM or NSOSM systems in the nonideal case.

For OSM models of isobaric ideal gases under simple diffusion, several recent papers have addressed numerical approaches, including a finite element method proposed by McLeod & Bourgault (2014), a finite volume method by Cancés et al. (2020) and a finite difference scheme by Bondesan et al. (2019). Such works typically employ a reference velocity as prescribed data.

Recently a finite element scheme for simple isobaric OSM diffusion in ideal gases was proposed by a subset of the current authors (Van-Brunt et al., 2021). They solved the augmented OSM equations (1.22) combined with the species continuity equations (1.29d). The present paper builds on the foundation established in Van-Brunt et al. (2021), but now fully incorporates momentum, nonideality and pressure-driven diffusion. In contrast to this prior work, we are able here to avoid a generalized saddle point formulation, and in Section 3 will derive a symmetric perturbed saddle point system to be solved at each nonlinear iteration—a more classical linear algebraic structure for which many solvers have been developed (Benzi et al., 2005). However, due to the more complex form (1.14) of the driving force, and since we solve for the chemical potentials to allow for nonideal fluids, we are not able to enforce the Gibbs–Duhem relation (1.17) to machine precision, as achieved in Van-Brunt et al. (2021).

 

Remark 1.1

Many cross-diffusion systems, such as those describing multiagent systems in mathematical biology (Carrillo et al., 2018), arise from a gradient flow of an associated entropy functional. Unfortunately, although the OSM system admits an associated thermodynamic energy—the Gibbs free energy—we are not able to show equivalence of the (S)OSM system to the Euler–Lagrange stationarity condition of any energy or Lagrangian functional, and hence cannot exploit any gradient flow structure. Instead, our mathematical line of attack will be to exploit the positive definiteness of the augmented transport matrix |$\textbf{M}^{\gamma }$| from (1.23). With our augmentations of the equations, the Picard scheme we propose below in Section 3 nevertheless gives rise to symmetric linearized problems to solve at each nonlinear iteration.

The remainder of this work is organized as follows. In Section 2, we derive a novel variational formulation of the fully coupled nonlinear SOSM problem, incorporating boundary conditions and augmentation terms, as a nonlinear perturbed saddle point-like system, using a novel solution-dependent test space relating to the thermodynamic driving force; our principal discovery is the duality between the diffusion driving forces, and the combination of species continuity equations with the divergence of the mass-average velocity constraint. Section 3 proposes a Picard-like linearization, which is proven to be well-posed under physically reasonable assumptions. Section 4 identifies appropriate finite element spaces, and the structural relations which should hold between them, for a well-posed and convergent discretization of this linearization. We then validate our convergence results numerically. Finally, we illustrate our method by simulating the steady mixing of liquid benzene and cyclohexane in a two-dimensional microfluidic laminar-flow device.

2. Variational formulation

We employ standard notation for the Sobolev space |$H^{k}(\varOmega ; \mathbb{X})$| (or |$L^{2}(\varOmega ; \mathbb{X})$| when |$k = 0$|⁠) with domain |$\varOmega \subset \mathbb{R}^{d}$| and codomain |$\mathbb{X}$|⁠, and associated norm |$\|\cdot \|_{k}$| and seminorm |$|\cdot |_{k}$|⁠. We denote by |$\mathbb{S} = \mathbb{R}^{d\times d}_{\textrm{sym}}$| the space of |$d\times d$| symmetric tensors. The symbol |$\lesssim $| denotes inequality up to a constant that may depend on mesh regularity, but not mesh spacing |$h$|⁠. Let |$L^{2}_{0}(\varOmega ) := \{z\in L^{2}(\varOmega ) \int{\!\!\!\!\!-}_{\varOmega } z~\textrm{d}x = 0\}$|⁠. Let |$\varGamma = \partial \varOmega $| and let |$\langle \cdot ,\cdot \rangle _{\varGamma }$| denote the |$(H^{-1/2}\times H^{1/2})(\varGamma ;\mathbb{R}\ \textrm{or}\ \mathbb{R}^{d})$| dual pairing.

2.1 Integrability of pressure gradients

For isothermal, isobaric multicomponent diffusion in an ideal gas mixture as originally considered by Maxwell (1867) and Stefan (1871), it suffices to work with driving forces of the form

$$ \begin{align}& d_{i} = -RT \nabla c_{i}.\end{align} $$

(2.1)

In a variational formulation of the nonisobaric case, one would like to integrate the pressure gradient term in our diffusion driving forces (1.14) by parts, to reduce the regularity requirement on |$p$|⁠. However, it is not obvious how to do so, since the mass fractions |$\omega _{i}$| are spatially varying. Formally, for a species velocity test function |$u$| (omitting surface terms for simplicity),

$$ \begin{align}& \int_{\varOmega} u\cdot(\omega_{i}\nabla p)~\textrm{d}x = -\int_{\varOmega}\operatorname{div}(u\omega_{i})p~\textrm{d}x = -\int_{\varOmega}(\omega_{i}\operatorname{div} u + u\cdot\nabla\omega_{i})p~\textrm{d}x,\end{align} $$

(2.2)

suggesting |$u$| be taken from a space requiring at least |$\operatorname{div} u\in L^{2}(\varOmega )$|⁠; however, we wish to avoid the introduction of the rightmost term into the bilinear forms and instead exploit the coercivity on |$L^{2}(\varOmega ;\mathbb{R}^{d})$| induced by the augmented Onsager transport matrix (1.24) (as in Van-Brunt et al. (2021)).

In order to rigorously incorporate the effect of pressure-driven diffusion, we are therefore led to consider the somewhat unorthodox possibility of formulating the Stokes subproblem with pressure lying in |$H^{1}(\varOmega )$|⁠. Typically, the condition that |$p \in H^{1}(\varOmega )$| may be provided by elliptic regularity results for the pressure field, but to the authors’ knowledge, the a priori square-integrability of pressure gradients (i.e. for which, we emphasize, pressure is defined to lie in |$H^{1}(\varOmega )$|⁠) has not been considered for the Stokes system, except in the incompressible case (e.g. Beirão et al. (2022), and in the incompressible case at the discrete level in Stenberg (1989)). This condition is also suggested by the case of pure Stefan–Maxwell diffusion for nonisobaric ideal gases. Here the driving forces are as follows:

$$ \begin{align}& d_{i} = -RT \nabla c_{i} + \omega_{i} \nabla p,\end{align} $$

(2.3)

which suggests considering each |$c_{i}$| (and hence |$c_{\textrm{T}}$|⁠) to lie in |$H^{1}(\varOmega )$|⁠, which forces the pressure to lie in the same space due to the ideal equation of state |$p = c_{\textrm{T}}RT$|⁠.

In general, one must distinguish between the thermodynamic pressure |$p$|⁠, which we use throughout this paper, and the mechanical pressure |$p_{\textrm{m}} := -{\operatorname{tr}\sigma }/{d}$|⁠. The mechanical pressure is related to the spherical Cauchy stress by |$\operatorname{sph}\sigma := \frac{\operatorname{tr}\sigma }{d}\mathbb{I} = -p_{\textrm{m}}\mathbb{I}$|⁠, and to |$p$| by

$$ \begin{align}& p = p_{\textrm{m}} + \zeta\operatorname{div} v.\end{align} $$

(2.4)

In the context of multicomponent flow, this decomposition is discussed in further detail by Bothe & Dreyer (2015). Even in the simpler incompressible limit where |$p = p_{\textrm{m}}$|⁠, we cannot expect extra regularity of |$\nabla p = -\operatorname{div}(\operatorname{sph}\sigma )$| because |$H(\textrm{div};{\mathbb{S}})$|⁠, the natural space for |$\sigma $|⁠, is not closed under taking spherical parts.⁵ Consequently, we do not take |$p\in H^{1}(\varOmega )$|⁠, but as a compromise consider a weaker condition defined by the combined viscous stress-pressure space

$$ \begin{align}& \{(\tau, p)\in L^{2}(\varOmega;\mathbb{S})\times L^{2}(\varOmega)~\vert~\operatorname{div} \tau - \nabla p\in L^{2}(\varOmega;\mathbb{R}^{d})\} \left(\supsetneq H(\textrm{div};{\mathbb{S}})\times H^{1}(\varOmega)\right),\end{align} $$

(2.5)

and assign to it the weaker norm |$\|\tau \|_{0}^{2} + \|p\|^{2}_{0} + \|\!\operatorname{div}\tau - \nabla p\|_{0}^{2}$|⁠. This space and norm were previously employed by Manouzi & Farhloul (2001) in an analysis of a non-Newtonian incompressible Stokes flow where |$\tau $| was taken to be the deviatoric shear stress. Membership of the space (2.5) is naturally interpretable as the square-integrability of the divergence of the full Cauchy stress, i.e. that |$\sigma = \tau - p\mathbb{I}\in H(\textrm{div};{\mathbb{S}})$|⁠. Together with an analogous condition for the chemical potential gradient to be detailed next, this weaker condition will account for the pressure gradient in the driving forces.⁶

2.2 Fully coupled variational formulation

In this subsection, we derive a variational formulation for the stationary problem as a nonlinear perturbed saddle point-like system. We have found the following statement of the problem to be a feasible tradeoff between the (competing) goals of: physical relevance of variables and boundary data, regularity assumptions, numerical implementability and effectiveness, analytic tractability of continuous and discrete well-posedness, enforcement of fundamental thermodynamic relations and extensibility to the anisothermal and nonideal settings.

For boundary data, we prescribe mass flux and molar fluxes:

$$ \begin{alignat}{2} \rho v &= g_{\textrm{v}} \in H^{1/2}(\varGamma;\mathbb{R}^{d}) &&\ \textrm{on}\ \varGamma, \end{alignat} $$

(2.6a)

$$ \begin{alignat}{2} c_{i} v_{i} \cdot \textbf{n} &= g_{i} \in H^{-1/2}(\varGamma) &&\ \textrm{on}\ \varGamma, \quad i = 1, \ldots, n. \end{alignat} $$

(2.6b)

For consistency with the mass-average velocity constraint (1.6), we require

$$ \begin{align}& \sum_{i} M_{i} g_{i} = g_{\textrm{v}} \cdot \textbf{n},\end{align} $$

(2.7)

with equality in |$H^{-1/2}(\varGamma )$|⁠. We further impose conditions

$$ \begin{align} \int_{\varOmega} p~\textrm{d}x = \int_{\varOmega} \mu_{i}~\textrm{d}x = 0, \quad i = 1, \ldots, n,\end{align} $$

(2.8)

on the pressure and chemical potentials. Typically, the equation of state will require or imply strict positivity of |$p$| everywhere, in which case this condition should be understood as |$\int{\!\!\!\!\!-} _{\varOmega } p~\textrm{d}x = p^{\ominus }> 0$| and that |$p$| be shifted by the known value |$p^{\ominus }$| as a postprocessing step.

Denote an |$n$|-tuple of functions by the notation |$\tilde{x} := \{x_{i} \}_{i=1}^{n}$|⁠. Let |$Q = L^{2}(\varOmega ;\mathbb{R}^{d})^{n}\times L^{2}(\varOmega ;\mathbb{R}^{d})$| with norm |$\|(\tilde{v}, v)\|_{Q}^{2} := \|\tilde{v}\|_{0}^{2} + \|v\|_{0}^{2}$|⁠. For formal derivation of the weak form, we assume the solution tuple |$(\tilde{\mu }, \tau , p, \tilde{v}, v)$| to be smooth on |$\overline{\varOmega }$|⁠, and consider choosing |$(\tilde{w}, s, q)$| from the solution-dependent potential-stress-pressure test space

$$ \begin{align}& \varTheta := \left\{ \left(\tilde{w}, s, q \right) \in L^{2}_{0}(\varOmega)^{n} \times L^{2}(\varOmega;\mathbb{S})\times L^{2}_{0}(\varOmega)~\middle|\begin{array}{rl}\operatorname{div} s - \nabla q &\in L^{2}(\varOmega;\mathbb{R}^{d}),\\ - c_{i} \nabla w_{i} + \omega_{i}\nabla q &\in L^{2}(\varOmega;\mathbb{R}^{d})~\forall i\\ \end{array} \right\}.\end{align} $$

(2.9)

Here it is understood that the |$\{c_{i}, \omega _{i}\}_{i}$| are computed from the solution tuple. Multiplying the |$i$|th continuity equation (1.29d) by |$w_{i}$|⁠, the divergence of the mass-average velocity constraint (1.29e) by |$q$|⁠, and contracting the stress constitutive law (1.29b) with |$s$|⁠, we obtain

$$ \begin{align}& \sum_{i}\left(\operatorname{div}(c_{i}v_{i}) - r_{i}\right)w_{i} + \operatorname{div}\left(v - \sum_{i}\omega_{i}v_{i}\right)q + (\mathscr{A}\tau - \varepsilon(v)):s = 0,\end{align} $$

(2.10)

and hence

$$ \begin{align}&\int_{\Omega} \sum_{i} (\operatorname{div}(c_{i}v_{i})w_{i} - \operatorname{div}(\omega_{i}v_{i})q) + \mathscr{A}\tau:s - (s - q\mathbb{I}):\varepsilon(v)~\textrm{d}x = \int_{\Omega}\sum_{i} r_{i} w_{i}~\textrm{d}x.\end{align} $$

(2.11)

Integrating by parts yields

$$ \begin{align}\int_{\Omega} &\mathscr{A}\tau:s + \sum_{i}(-c_{i} \nabla w_{i} + \omega_{i} \nabla q)\cdot v_{i} + (\operatorname{div} s - \nabla q)\cdot v~\textrm{d}x\nonumber \\ &= \left\langle (s - q\mathbb{I})\textbf{n}, v\right\rangle_{\Gamma} + \sum_{i}\left\langle c_{i}v_{i}\cdot\textbf{n}, -w_{i} + \frac{\omega_{i}}{c_{i}}q\right\rangle_{\Gamma} + \int_{\Omega} \tilde{r}\cdot\tilde{w}~\textrm{d}x\nonumber\\ &= \left\langle (s - q\mathbb{I})\textbf{n}, \frac{g_{\textrm{v}}}{\rho}\right\rangle_{\Gamma} + \sum_{i}\left\langle g_{i}, -w_{i} + \frac{M_{i}}{\rho}q\right\rangle_{\Gamma} + \int_{\Omega} \tilde{r}\cdot\tilde{w}~\textrm{d}x.\end{align} $$

(2.12)

Now for each |$i = 1, \ldots , n$|⁠, we take the scalar product of |${u}_{i} \in L^{2}(\varOmega ;\mathbb{R}^{d})$| with the augmented OSM equation (1.29a) and integrate over |$\varOmega $| to obtain

$$ \begin{align}& \int_{\varOmega}\left(-c_{i} \nabla \mu_{i} + \omega_{i} \nabla p\right) \cdot{u}_{i} - {u}_{i} \cdot \sum_{j} \textbf{M}_{ij} {v}_{j} - \gamma\omega_{i} \sum_{j} \omega_{j} (v_{j} - v) \cdot u_{i}~\textrm{d}x = 0.\end{align} $$

(2.13)

Taking the inner product of the augmented Cauchy momentum balance (1.29c) with |$u\in L^{2}(\varOmega ;\mathbb{R}^{d})$| yields

$$ \begin{align}& \int_{\varOmega} (\operatorname{div} \tau - \nabla p) \cdot u - \gamma \left(\sum_{j} \omega_{j} (v - v_{j}) \right)\cdot u~\textrm{d}x = -\int_{\varOmega} \rho f \cdot u~\textrm{d}x.\end{align} $$

(2.14)

We sum (2.13) over |$i$| and add (2.14) to derive (using the symmetry of |$\textbf{M}_{ij}$|⁠)

$$ \begin{align}\begin{aligned} & \int_{\varOmega} \sum_{i} (- c_{i}\nabla \mu_{i}+\omega_{i}\nabla p) \cdot u_{i}+(\operatorname{div}\tau - \nabla p)\cdot u \\[-6pt] & - \sum_{i,j} {v}_{i}\cdot\textbf{M}_{ij} {u}_{j} - \gamma \left(\sum_{j} \omega_{j} (v_{j} - v) \right) \cdot \left( \sum_{i}\omega_{i} (u_{i} - u) \right)~\textrm{d}x = \int_{\Omega} -\rho f \cdot u~\textrm{d}x.\end{aligned}\end{align} $$

(2.15)

Note that both augmentations (1.22) and (1.25) were involved in deriving this expression.

Finally, we observe that by definition, we have |$\omega _{i}\in L^{\infty }(\varOmega )$| with |$\|\omega _{i}\|_{L^{\infty }(\varOmega )}\leq 1$|⁠. Moreover, we make the physically reasonable assumptions that the concentrations associated with the solution are uniformly bounded, |$c_{i} \in L^{\infty }(\varOmega )$|⁠, with |$c_{i} \geq \kappa> 0$| almost everywhere (a.e.), as in Van-Brunt et al. (2021) (which in turn implies |$\textbf{M}^{\gamma }_{ij}, \rho \in L^{\infty }(\varOmega )$|⁠, and |$\rho \geq \kappa \sum _{i} M_{i}> 0 $| a.e.), and that the density gradient is uniformly bounded, |$\nabla \rho \in L^{\infty }(\varOmega ;\mathbb{R}^{d})$|⁠.⁷

 

Definition 2.1

(Weak solution to the SOSM equations.) We define a weak solution to the augmented Stokes–Onsager–Stefan–Maxwell system to be a |$(2n+3)$|-tuple

$$ \begin{align}& \left(\{\mu_{i}\}_{i=1}^{n}, \tau, p, \{v_{i}\}_{i=1}^{n}, v\right) \in L^{2}_{0}(\varOmega)^{n} \times L^{2}(\varOmega;\mathbb{S}) \times L^{2}_{0}(\varOmega)\times \underbrace{L^{2}(\varOmega;\mathbb{R}^{d})^{n}\times L^{2}(\varOmega;\mathbb{R}^{d})}_{Q}\end{align} $$

(2.16)

inducing concentrations |$\{c_{i}\}_{i=1}^{n}$| through a constitutive law (such as (1.26)) implicitly defining |$c_{i} = c_{i}(\{\mu _{i}\}_{i=1}^{n}, p) \geq \kappa> 0$| a.e. for |$i = 1, \ldots , n$|⁠, such that

$$ \begin{alignat}{2} \|c_{i}\|_{L^{\infty}(\varOmega)} &<\infty,~i = 1, \ldots, n, \end{alignat} $$

(2.17a)

$$ \begin{alignat}{2} \|\nabla \rho\|_{L^{\infty}(\varOmega;\mathbb{R}^{d})} &<\infty, \end{alignat} $$

(2.17b)

$$ \begin{alignat}{2} \|\! \operatorname{div} \tau - \nabla p \|^{2}_{0} &<\infty, \end{alignat} $$

(2.17c)

$$ \begin{alignat}{2} \| - c_{i} \nabla \mu_{i} + \omega_{i} \nabla p \|^{2}_{0} &< \infty,~i = 1, \ldots, n, \end{alignat} $$

(2.17d)

and satisfying (2.12), (2.15) for all test tuples |$(\{w_{i}\}_{i=1}^{n}, s, q, \{u_{i}\}_{i=1}^{n}, u)\in \varTheta \times Q$|⁠, where |$\varTheta $| is defined in (2.9).

Observe that the solution tuple does not reside in any standard Sobolev space, but that the regularity assumptions placed on the solution tuple and test spaces ensure that the surface terms in (2.12) are well-defined. Recall that condition (2.17c) is the square-integrability of the Cauchy stress |$\sigma $| (as in Manouzi & Farhloul, 2001). The nonlinear integrability condition (2.17d) is to our knowledge a novel requirement, but also has a natural interpretation, namely the square-integrability of the diffusion driving forces:⁸

$$ \begin{align}& d_{i} \in L^{2}(\varOmega;\mathbb{R}^{d}).\end{align} $$

(2.18)

Moreover, we emphasize that this unorthodox formulation allows the rigorous incorporation of pressure diffusion via the pressure gradient on the left side of (1.29a), despite the fact that the pressure field is not a priori|$H^{1}$|-regular in the Stokes subsystem. Later in Section 4.4 we observe convergence of the diffusion driving forces in |$L^{2}$| and of the pressure in |$H^{1}$|⁠, but otherwise leave this consideration, and further investigation into the optimal nonlinear formulation of the SOSM system, as intriguing open questions.

 

Remark 2.2

In the derivation of (2.10), we used the divergence of the mass-average velocity constraint (1.29e), which ignores the |$\operatorname{curl}$| component in the Helmholtz decomposition of the mass-average velocity relationship (1.6). This choice ensures that the number of equations matches the number of unknown variables. The full constraint is weakly incorporated, however, via the augmentations (1.22) and (1.25). In the case of an ideal isobaric isothermal gas, a proof that the augmentations will enforce the full constraint (1.6) was given in Van-Brunt et al. (2021). By means of the first law of thermodynamics and the extensivity of the Gibbs free energy, we now show this is formally true in general, provided the constitutive equations relating concentrations and chemical potentials are thermodynamically rigorous in the sense that they arise from a Gibbs free energy functional.

Specifically, the first law and definition of the chemical potential imply that (see, for example, Moran et al., 2018, eq. (11.112))

$$ \begin{align}& \nabla p + \sum_{i} \mu_{i} \nabla c_{i} = \nabla \tilde{G}\end{align} $$

(2.19)

for the volumetric Gibbs free energy functional |$\tilde{G}$| (⁠|$\textrm{J}/\textrm{m}^{3}$|⁠), assumed to lie in |$H^{1}(\varOmega )$|⁠.

Choosing |$u = 0$| and |$u_{i} = \hat{u}\in L^{2}(\varOmega ;\mathbb{R}^{d})$| for each |$i$| in (2.15), we derive, by the nullspace (1.15),

$$ \begin{align}& \int_{\varOmega} \sum_{i} \left( -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p \right) \cdot \hat{u} + \gamma \left(v - \sum_{j} \omega_{j} v_{j}\right) \cdot \hat{u}~\textrm{d}x = 0.\end{align} $$

(2.20)

Choosing |$\hat{u} = \nabla q$| for |$q \in H^{1}(\varOmega )$|⁠, and |$w_{i} = 0$| for each |$i$| and |$s = 0$| in (2.12), and noting that |$\left \langle g_{i}, \frac{M_{i}}{\rho }q\right \rangle _{\varGamma } = \left \langle M_{i} g_{i} q, \frac{1}{\rho }\right \rangle _{\varGamma }$| for each |$i$|⁠, the final term may be written as

$$ \begin{align}& \int_{\varOmega}\left(v - \sum_{j}\omega_{j} v_{j}\right)\cdot\nabla q~\textrm{d}x = \left\langle qg_{\textrm{v}}\cdot\textbf{n}, \frac{1}{\rho}\right\rangle_{\varGamma} - \left\langle\sum_{i} M_{i} g_{i} q, \frac{1}{\rho}\right\rangle_{\varGamma} = 0\end{align} $$

(2.21)

by the compatibility condition (2.7). Thus,

$$ \begin{align}& \int_{\varOmega} \sum_{i} \left( -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p \right) \cdot \nabla q~\textrm{d}x = 0\end{align} $$

(2.22)

for each |$q \in H^{1}(\varOmega )$|⁠. Further rearrangement of the integrand yields

$$ \begin{align}& \sum_{i} -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p = \nabla p + \sum_{i} -\nabla \left( c_{i} \mu_{i} \right ) + \sum_{i} \mu_{i} \nabla c_{i} = -\sum_{i} \nabla \left ( c_{i} \mu_{i} \right ) + \nabla \tilde{G}.\end{align} $$

(2.23)

Then choosing |$\nabla q = -\sum _{i}\nabla (c_{i} \mu _{i} ) + \nabla \tilde{G}$|⁠, we derive the Gibbs–Duhem equation

$$ \begin{align}& \int_{\varOmega} \left\vert\sum_{i} -c_{i} \nabla \mu_{i} + \omega_{i} \nabla p \right\vert^{2}\,\textrm{d}x = 0.\end{align} $$

(2.24)

Then (2.20) reduces to

$$ \begin{align}& v = \sum_{i} \omega_{i} v_{i},\end{align} $$

(2.25)

showing the enforcement of the mass-average velocity constraint. Although this shows the preservation of the Gibbs–Duhem and mass-average velocity constraints in the nonlinear infinite dimensional case, it is significantly harder in the discretized nonideal system to enforce these exactly. This is a price we pay for solving for the chemical potential.

3. Linearization and well-posedness

3.1 Variational formulation of a generalized Picard scheme

In this section we derive a variational formulation of a generalized Picard linearization. Given a previous estimate for the potentials |$\tilde{\mu }^{k}$| and pressure |$p^{k}$| for |$k\geq 0$|⁠, we regard these as fixed quantities which determine the concentrations |$\tilde{c}^{k}$| via chemical potential constitutive laws and an appropriate equation of state such as (1.26). This in turn determines the density |$\rho ^{k}$|⁠, mass fractions |$\tilde{\omega }^{k}$|⁠, total concentration |$c^{k}_{\textrm{T}}$| and transport matrix |$\textbf{M}^{k}$| defined via (1.3), (1.7), (1.12) and (1.11), respectively. We then construct a linear system to solve for the next iterate |$((\tilde{\mu }^{k+1}, \tau ^{k+1}, p^{k+1}), (\tilde{v}^{k+1}, v^{k+1}))$|⁠. This update strategy is expected to be effective because the gradients of chemical potential, pressure and mass-average velocity primarily drive the dynamics of multicomponent flow; the role of the species concentrations is mostly confined to the effect of altering the drag coefficients in the transport matrix. We make the following physically reasonable assumptions about each iterate, in analogy to Definition 2.1.

 

Assumption 3.1

(Uniform positivity of concentrations.) For each |$k\geq 0$|⁠, we assume |$c^{k}_{i} \in L^{\infty }(\varOmega ), \rho ^{k} \in W^{1, \infty }(\varOmega )$|⁠, and that |$c^{k}_{i}\geq \kappa> 0$| a.e. for each |$i$|⁠.

This again implies |$\rho ^{k} \geq \kappa \sum _{i}M_{i}> 0$| a.e. We also assume henceforth that |$\gamma> 0$|⁠.

Given |$\tilde{c}^{k}$| and the corresponding |$\tilde{\omega }^{k}$|⁠, we define the iteration-dependent weighted function space

$$ \begin{align}& \varTheta^{k} := \left\{ (\tilde{w}, s, q) \in L^{2}_{0}(\varOmega)^{n} \times L^{2}(\varOmega;\mathbb{S})\times L^{2}_{0}(\varOmega)\middle|\begin{array}{rl}\operatorname{div} s - \nabla q&\in L^{2}(\varOmega;\mathbb{R}^{d}),\\ -c_{i}^{k} \nabla w_{i} + \omega_{i}^{k} \nabla q&\in L^{2}(\varOmega;\mathbb{R}^{d})~\forall i \end{array} \right\},\end{align} $$

(3.1)

whose defining conditions linearize those in (2.9). This mixed space is Hilbert with graph norm

$$ \begin{align}&\|(\tilde{w}, s,q)\|_{\Theta^{k}}^{2} := \sum_{i}\|w_{i}\|_{0}^{2} + \|s\|_{0}^{2} + \|q\|_{0}^{2} + \|\!\operatorname{div} s - \nabla q\|_{0}^{2} + \sum_{i} \left\|-c_{i}^{k} \nabla w_{i} + \omega_{i}^{k} \nabla q \right\|_{0}^{2}.\end{align} $$

(3.2)

We now formulate our linearized problem as a symmetric perturbed saddle point problem. Define |$A_{k}:Q\to Q^{*}, \varLambda :\varTheta ^{k}\to (\varTheta ^{k})^{*}, B_{k}:\varTheta ^{k}\to Q^{*}$| by

$$ \begin{alignat}{2} A_{k}(\tilde{v}, v; \tilde{u}, u) &:= \int_{\varOmega} \sum_{i,j} v_{i}\cdot\textbf{M}^{k}_{ij}u_{j}~\textrm{d}x + \gamma \int_{\varOmega} \left(\sum_{j}\omega^{k}_{j}(v_{j} - v)\right)\cdot\left(\sum_{i}\omega^{k}_{i}(u_{i} - u)\right)\textrm{d}x, \end{alignat} $$

(3.3a)

$$ \begin{alignat}{2} \varLambda(\tilde{\mu}, \tau, p; \tilde{w}, s, q) &:= \int_{\varOmega}\mathscr{A}\tau:s~\textrm{d}x, \end{alignat} $$

(3.3b)

$$ \begin{alignat}{2} B_{k}(\tilde{\mu}, \tau, p; \tilde{u}, u) &:= \int_{\varOmega}\sum_{i}\left(-c^{k}_{i}\nabla\mu_{i} + \omega^{k}_{i}\nabla p \right)\cdot u_{i} + (\operatorname{div}\tau - \nabla p)\cdot u~\textrm{d}x, \end{alignat} $$

(3.3c)

and the functionals

$$ \begin{align}\ell^{1}_{k}(\tilde{w}, s, q) &:= \left\langle(s - q\mathbb{I})\textbf{n}, \frac{g_{\textrm{v}}}{\rho^{k}}\right\rangle_{\Gamma} + \sum_{i}\left\langle g_{i}, - w_{i} + \frac{M_{i}}{\rho^{k}}q\right\rangle_{\Gamma} + \int_{\Omega}\tilde{r}\cdot\tilde{w}~\textrm{d}x,\nonumber\\ \ell^{2}_{k}(\tilde{u}, u) &:= -\int_{\Omega}\rho^{k} f\cdot u~\textrm{d}x.\end{align} $$

(3.4)

Note that under Assumption 3.1, each of the bilinear functionals is continuous; we will denote their norms as |$\|A_{k}\|$|⁠, |$\| \varLambda \|$| and |$\|B_{k}\|$|⁠, respectively. Our linearized problem is posed as follows: find |$((\tilde{\mu }^{k+1}, \tau ^{k+1}, p^{k+1}), (\tilde{v}^{k+1},v^{k+1})) \in \varTheta ^{k} \times Q$| such that

$$ \begin{align}\begin{aligned}\Lambda(\tilde{\mu}^{k+1}, \tau^{k+1}, p^{k+1}; \tilde{w}, s, q) + B_{k}(\tilde{w}, s, q; \tilde{v}^{k+1}, v^{k+1}) &= \ell^{1}_{k}(\tilde{w}, s, q)&&~\forall~(\tilde{w}, s, q)\in \Theta^{k}, \\ B_{k}(\tilde{\mu}^{k+1}, \tau^{k+1}, p^{k+1}; \tilde{u}, u) - A_{k}(\tilde{v}^{k+1}, v^{k+1}; \tilde{u}, u) &= \ell^{2}_{k}(\tilde{u}, u)&&~\forall~(\tilde{u}, u)\in Q,\end{aligned}\end{align} $$

(3.5)

i.e., defining the transpose |$B_{k}^{\top }:Q\to (\varTheta ^{k})^{*}$| in the canonical way,

$$\ $$

(3.6)

We note that the variational terms involving chemical potential and pressure gradients are precisely of the same variational form as the species continuity equations and the divergence of the mass-average velocity constraint, which can be seen by inspecting (2.12) and (2.15). This key insight is what leads to a symmetric system.

Our nonlinear iteration scheme is as follows: for an initial estimate of the concentrations |$\tilde{c}^{0}$|⁠, we solve the system (3.5) for the updated variables |$((\tilde{\mu }^{k+1}, \tau ^{k+1}, p^{k+1}), (\tilde{v}^{k+1},v^{k+1})) \in \varTheta ^{k} \times Q$|⁠, for |$k = 0, 1, 2, \ldots $|⁠. By the relations detailed in Section 1.4, these variables are used to calculate the updated concentrations |$\tilde{c}^{k+1}$|⁠. This is iterated until for some set tolerance |$\varepsilon> 0$|⁠,

$$ \begin{align}& \left( \|(\tilde{\mu}^{k+1}, \tau^{k+1}, p^{k+1} ) - (\tilde{\mu}^{k}, \tau^{k}, p^{k} ) \|^{2}_{\varTheta^{k}} + \| (\tilde{v}^{k+1}, v^{k+1}) - (\tilde{v}^{k}, v^{k}) \|^{2}_{Q} \right)^{1/2} \leq \varepsilon.\end{align} $$

(3.7)

3.2 Well-posedness of the linearized system

We will now prove that the saddle point system (3.5) is well-posed under Assumption 3.1. This will require a Poincaré-type inequality for the following seminorm on |$\varTheta ^{k}$|⁠:

$$ \begin{align}& |(\tilde{w}, s, q) |_{\varTheta^{k}}^{2}: =\|s\|_{0}^{2} +\|\! \operatorname{div} s - \nabla q\|_{0}^{2} + \sum_{i} \left\|-c_{i}^{k} \nabla w_{i} + \omega_{i}^{k} \nabla q \right\|_{0}^{2}.\end{align} $$

(3.8)

 

Lemma 3.2

(Poincaré inequality for the driving force.) Let |$\varOmega $| be a Lipschitz domain. Under Assumption 3.1, there exists |$K> 0$| such that for each |$(\tilde{\mu }, \tau , p) \in \varTheta ^{k}$|⁠,

$$ \begin{align}& \|(\tilde{\mu}, \tau, p) \|_{\varTheta^{k}} \leq K|(\tilde{\mu}, \tau, p) |_{\varTheta^{k}}.\end{align} $$

(3.9)

 

Proof of Lemma 3.2.

The first stage of the proof is to show that

$$ \begin{align}& \|p \|_{0} \lesssim \|\tau\|_{0} + \|\! \operatorname{div} \tau - \nabla p \|_{0},\end{align} $$

(3.10)

following and mildly generalizing Manouzi & Farhloul (2001, Lemma 4). Set |$\theta = \tau - p\mathbb{I} - r\mathbb{I}$| where |$r = \frac{1}{d | \varOmega |} \int _{\varOmega } \operatorname{tr}\tau ~\textrm{d}x$|⁠. Then

$$ \begin{align}& \| \tau - p\mathbb{I} \|_{0} \leq \| \theta \|_{0} + \|r \mathbb{I}\|_{0}.\end{align} $$

(3.11)

As |$\int _{\varOmega } \operatorname{tr}\theta ~\textrm{d}x = 0$|⁠, we can use Boffi et al. (2013, Proposition 9.1.1) to derive

$$ \begin{align}\|\tau - p\mathbb{I} \|_{0} &\lesssim \|\!\operatorname{dev}\theta \|_{0} + \|\! \operatorname{div} \theta\|_{0} + \| r\mathbb{I} \|_{0}\nonumber \\& \lesssim \|\tau \|_{0} + \|\! \operatorname{div} \tau - \nabla p \|_{0},\end{align} $$

(3.12)

where the deviator is defined as |$\operatorname{dev}\theta := \theta - \frac{\operatorname{tr}\theta }{d}\mathbb{I} = \operatorname{dev}\tau $|⁠. Now using

$$ \begin{align}& \sqrt{d}\|p\|_{0} \leq \|\tau - p \mathbb{I}\|_{0} + \|\tau \|_{0},\end{align} $$

(3.13)

the result (3.10) follows. For the second stage of the proof, we will show that

$$ \begin{align}& \|\mu_{i} \|_{0} \lesssim \|p\|_{0} + \left\|-c^{k}_{i}\nabla\mu_{i} + \omega^{k}_{i}\nabla p \right\|_{0}.\end{align} $$

(3.14)

This combined with (3.10) gives (3.9). To prove this second inequality, for each |$i$| we take the unique |$z_{i} \in H^{1}_{0}(\varOmega ;\mathbb{R}^{d})/\ker (\operatorname{div})$| such that |$\operatorname{div} z_{i} = \mu _{i}$|⁠. Then |$u_{i} := z_{i}/c^{k}_{i} \in L^{2}(\varOmega ;\mathbb{R}^{d})$| with |$\operatorname{div}(c^{k}_{i} u_{i}) = \mu _{i}$|⁠. With integration by parts we deduce

$$ \begin{align}& \int_{\varOmega} \left(-c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k} \nabla p\right) \cdot u_{i}~\textrm{d}x = \int_{\varOmega} |\mu_{i}|^{2} - M_{i}p\left( \frac{\mu_{i}}{\rho^{k}} - \frac{\nabla\rho^{k}}{(\rho^{k})^{2}}\cdot c^{k}_{i} u_{i}\right)\textrm{d}x.\end{align} $$

(3.15)

Upon rearrangement, we can derive the inequality

$$ \begin{align}\|\mu_{i}\|^{2}_{0} &\leq M_{i}\|p\|_{0}\left(\frac{\|\mu_{i}\|_{0}}{\kappa\sum_{j} M_{j}} + \|u_{i}\|_{0}\|c^{k}_{i}\|_{0}\left\|\frac{\nabla\rho^{k}}{(\rho^{k})^{2}}\right\|_{L^{\infty}(\Omega;\mathbb{R}^{d})}\right) + \left\|-c_{i}^{k} \nabla \mu_{i} + \omega_{i} \nabla p \right\|_{0} \| u_{i} \|_{0}\nonumber\\ &\leq \kappa^{-1} \|p\|_{0}\left(\|\mu_{i}\|_{0} + \|u_{i} \|_{0} \left\|c^{k}_{i}\right\|_{L^{\infty}(\Omega)}\|\nabla \ln \rho^{k} \|_{L^{\infty}(\Omega;\mathbb{R}^{d})}\right) + \left\|-c_{i}^{k} \nabla \mu_{i} + \omega_{i} \nabla p \right\|_{0} \| u_{i} \|_{0}.\end{align} $$

(3.16)

By the bounded inverse theorem, |$\operatorname{div}$| admits a bounded left inverse, so |$\|z_{i}\|_{1}\lesssim \|\mu _{i}\|_{0}$| and thus

$$ \begin{align}& \| u_{i} \|_{0} \leq \kappa^{-1}\|z_{i}\|_{0} \lesssim \kappa^{-1}\|\mu_{i} \|_{0}.\end{align} $$

(3.17)

Combining this with (3.16), we can divide through by |$\| \mu _{i} \|_{0}$| to derive

$$ \begin{align}& \| \mu_{i} \|_{0} \lesssim \kappa^{-1}\|p\|_{0}\left(1 + \kappa^{-1}\left\|c^{k}_{i}\right\|_{L^{\infty}(\varOmega)}\|\nabla \ln \rho^{k} \|_{L^{\infty}(\varOmega)}\right) + \kappa^{-1} \left\|-c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k}\nabla p \right\|_{0}.\end{align} $$

(3.18)

Note that in particular, the two steps of this proof imply that

$$ \begin{align}& \begin{array}{rl} \operatorname{div}\tau - \nabla p &= 0\\ -c^{k}_{i}\nabla\mu_{i} + \omega_{i}^{k}\nabla p &= 0 \end{array} \bigg\} \Longrightarrow \bigg\{ \begin{array}{rl} \|p\|_{0}&\lesssim\|\tau\|_{0},\\ \|\mu_{i}\|_{0}&\lesssim\|p\|_{0}. \end{array}\end{align} $$

(3.19)

Physically, this implies that in the absence of (at least this linearization of) the driving force, one can recover the chemical potentials from the pressure. In this sense it is a generalization to the OSM framework of (Manouzi & Farhloul, 2001, Lemma 4), which is exactly the first line of (3.19): that in the absence of external forces, one can recover the pressure from the viscous stress. Since also the constant in (3.18) depends unfavourably on |$\kappa $| (the uniform lower bound on concentrations), we see also that such ‘recovery’ of the potentials becomes more unstable near the singular regime in which concentrations approach zero. Provided |$\kappa $| and the relative variation of the density are well-behaved across iterations, so too will be the resulting constant.

A further intermediate lemma we need to prove well-posedness is the following.

 

Lemma 3.3

(Coercivity of a perturbation to |$A_{k}$|⁠). Under Assumption 3.1 there exists |$\lambda _{\kappa }^{\gamma }> 0$| depending on |$\kappa $| and |$\gamma $| such that for all |$(\tilde{v},v) \in Q$|⁠,

$$ \begin{align}& \left(\frac{n}{2} + \frac{1}{4} \right) \lambda_{\kappa}^{\gamma} \| v \|_{0}^{2} + A_{k}( \tilde{v}, v; \tilde{v}, v ) \geq \frac{ \lambda_{\kappa}^{\gamma}}{4} \|(\tilde{v}, v)\|_{Q}^{2}.\end{align} $$

(3.20)

 

Proof.

Let us define |$\delta _{i} = v_{i} - v$| for each |$i$|⁠. Then we can explicitly compute

$$ \begin{align}& A_{k}(\tilde{v},v; \tilde{v},v) = \int_{\varOmega} \sum_{i,j} \delta_{i} \cdot \textbf{M}_{ij}^{k,\gamma} \delta_{j}~\textrm{d}x,\end{align} $$

(3.21)

where |$\textbf{M}^{k,\gamma }$| is defined using |$\tilde{c}^{k}$| via (1.23). It follows from Van-Brunt et al. (2021, Lemma 4.1) that this is a coercive bilinear form in |$\tilde{\delta }$|⁠: for some |$\lambda _{\kappa }^{\gamma }$|⁠,

$$ \begin{align}& A_{k}(\tilde{v}, v; \tilde{v}, v) \geq \frac{\lambda_{\kappa}^{\gamma}}{2} \sum_{i} \| \delta_{i} \|^{2}_{0} = \frac{\lambda_{\kappa}^{\gamma}}{2} \sum_{i} \| v_{i} -v\|^{2}_{0};\end{align} $$

(3.22)

hence,

$$ \begin{align}\left(\frac{n}{2} + \frac{1}{4} \right) \lambda_{\kappa}^{\gamma} \| v \|_{0}^{2} + A_{k}(\tilde{v}, v; \tilde{v}, v) &\geq \frac{\lambda_{\kappa}^{\gamma}}{2} \sum_{i} \left(\| v_{i} -v\|^{2}_{0} + \| v\|^{2}_{0} \right) + \frac{\lambda_{\kappa}^{\gamma}}{4} \|v\|^{2}_{0} \nonumber \\ & \geq \frac{\lambda_{\kappa}^{\gamma}}{4} \sum_{i} \| v_{i} \|_{0}^{2} + \frac{\lambda_{\kappa}^{\gamma}}{4} \|v\|^{2}_{0} = \frac{ \lambda_{\kappa}^{\gamma}}{4} \|(\tilde{v}, v)\|_{Q}^{2}.\end{align} $$

(3.23)

Despite the complexity of the original fully coupled physics problem, our constructed formulation allows us to invoke standard theory for well-posedness of Babuška (1971).

 

Theorem 3.4

(Well-posedness of the Picard linearization). Under Assumption 3.1, there exists a unique solution to the perturbed saddle point system (3.5).

 

Proof.

For a given |$(\mathfrak{p}, \mathfrak{q}) := \left ( (\tilde{\mu }, \tau , p), (\tilde{v}, v) \right ) \in \varTheta ^{k} \times Q$| and |$(\mathfrak{s}, \mathfrak{u}) := \left ( (\tilde{w}, s, q), (\tilde{u},u) \right ) \in \varTheta ^{k} \times Q$| we will define the bounded bilinear form |$\mathscr{G}: (\varTheta ^{k} \times Q) \times (\varTheta ^{k} \times Q) \rightarrow \mathbb{R}$| as

$$ \begin{align}& \mathscr{G}\left(\mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u} \right) := \varLambda ( \mathfrak{p}; \mathfrak{s} ) + B_{k} (\mathfrak{s}; \mathfrak{q}) + B_{k}(\mathfrak{p}; \mathfrak{u}) - A_{k}(\mathfrak{q}; \mathfrak{u}).\end{align} $$

(3.24)

We will prove Babuška’s inf-sup condition, namely that there exists a constant |$c> 0$| such that for each |$(\mathfrak{p}, \mathfrak{q}) \in \varTheta ^{k}\times Q$| there is |$(\mathfrak{s},\mathfrak{u}) \in \varTheta ^{k}\times Q$| such that

$$ \begin{align}& \frac{\mathscr{G}\left(\mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u}\right)}{\|(\mathfrak{s}, \mathfrak{u}) \|_{\varTheta^{k} \times Q}} \geq c \|(\mathfrak{p}, \mathfrak{q} )\|_{\varTheta^{k} \times Q},\end{align} $$

(3.25)

with product norm |$\| (\mathfrak{p}, \mathfrak{q}) \|_{\varTheta ^{k} \times Q}^{2} := \| \mathfrak{p}\|_{\varTheta ^{k}}^{2} + \| \mathfrak{q}\|_{Q}^{2}$|⁠. Note that |$\mathscr{G}$| is defined on the product of a space with itself and is symmetric, and so only the one inf-sup condition (3.25) needs to be verified. Proving (3.25) will be accomplished by showing that for a constant |$c> 0$|⁠, for each |$(\mathfrak{p}, \mathfrak{q}) \in \varTheta ^{k} \times Q$| there is |$(\mathfrak{s}, \mathfrak{u}) \in \varTheta ^{k} \times Q$| such that |$ \mathscr{G}\left ( \mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u} \right ) \geq c \| (\mathfrak{p}, \mathfrak{q}) \|_{\varTheta ^{k} \times Q}^{2}, $| and for a |$C> 0$| independent of |$(\mathfrak{p},\mathfrak{q})$|⁠,

$$ \begin{align}& C \| (\mathfrak{p}, \mathfrak{q} ) \|_{\varTheta^{k} \times Q} \geq \| (\mathfrak{s}, \mathfrak{u} ) \|_{\varTheta^{k} \times Q}.\end{align} $$

(3.26)

This, combined with our Poincaré-type inequality, will imply well-posedness. We do this by fixing |$(\mathfrak{s}, \mathfrak{u})$| as the ansatz

$$ \begin{align} &w_{i} = C_{1} \mu_{i}, \qquad s = C_{1} \tau + C_{2} s_{v}, \qquad &&q = C_{1}p, \nonumber\\ &u_{i} = C_{3} \left(-c^{k}_{i} \nabla \mu_{i} + \omega^{k}_{i} \nabla p\right) - C_{1} v_{i}, \qquad &&u = -C_{1} v + C_{4} (\operatorname{div} \tau - \nabla p).\end{align} $$

(3.27)

Here |$C_{1}, \ldots> 0$| are constants to be set later and |$s_{v} \in H(\textrm{div};{\mathbb{S}})$| is chosen to satisfy |$\operatorname{div} s_{v} = v$| and |$\|s_{v}\|_{H(\textrm{div};{\mathbb{S}})} \leq C_{\varSigma } \|v\|_{0}$| for a constant |$C_{\varSigma }$| independent of |$v$| (its existence is guaranteed by Arnold (2018, Theorem 8.1)). It is clear that (3.26) holds. With these choices of test functions we may compute

$$ \begin{align}\mathscr{G}\left( \mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u} \right) = \int_{\Omega}& \mathscr{A} \tau : (C_{1} \tau + C_{2} s_{v})~\textrm{d}x + C_{3} \sum_{i} \left\| - c^{k}_{i} \nabla \mu_{i} + \omega^{k}_{i} \nabla p \right\|_{0}^{2}\nonumber \\ &+ C_{4} \|\! \operatorname{div} \tau - \nabla p \|^{2}_{0} + C_{2} \|v\|^{2}_{0} - A_{k} \left( \tilde{v}, v; \tilde{u}, u \right)\!,\end{align} $$

(3.28)

and observe that the final term on the right may be written equivalently as

$$ \begin{align}& C_{1} A_{k} \left( \tilde{v}, v; \tilde{v}, v \right) - A_{k} \left( \tilde{v}, v; \left\{C_{3}\left(-c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k} \nabla p\right)\right\}_{i}, C_{4}(\operatorname{div} \tau - \nabla p) \right).\end{align} $$

(3.29)

With |$\lambda _{\kappa }^{\gamma }$| given by Lemma 3.3 we now choose |$C_{2} = \|A_{k}\|^{2}(2n + 1) \lambda _{\kappa }^{\gamma }$| and assume that |$C_{1} \geq 4 \|A_{k} \|^{2}$|⁠. Then

$$ \begin{align}& C_{2}\|v\|^{2}_{0} + C_{1} A_{k} \left( \tilde{v},v; \tilde{v}, v \right) \geq \|A_{k}\|^{2} \lambda_{\kappa}^{\gamma} \| (\tilde{v}, v)\|_{Q}^{2}.\end{align} $$

(3.30)

Using uniform positive definiteness of the compliance tensor, the bound on |$s_{v}$| and boundedness of the operators |$\varLambda , A_{k}$|⁠, we can proceed to bound (3.28) from below as

$$ \begin{align}\mathscr{G}\left( \mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u} \right) &\geq \alpha C_{1} \| \tau \|^{2}_{0} + \|A_{k}\|^{2} \lambda_{\kappa}^{\gamma} \| (\tilde{v}, v)\|_{Q}^{2} + C_{3} \sum_{i} \left\| -c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k} \nabla p \right\|_{0}^{2} \nonumber\\ & + C_{4} \|\! \operatorname{div} \tau - \nabla p \|_{0}^{2} - \lambda_{\kappa}^{\gamma} (n+1)C_{\Sigma} \|A_{k} \|^{2} \| \Lambda \| \|v \|_{0} \| \tau \|_{0} \nonumber\\ & - \| A_{k} \| \| (\tilde{v}, v)\|_{Q} \left( C^{2}_{3}\sum_{i} \left\| -c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k} \nabla p \right\|^{2}_{0} + C^{2}_{4}\|\! \operatorname{div} \tau - \nabla p \|^{2}_{0} \right)^{1/2}.\end{align} $$

(3.31)

Here |$\alpha> 0$| is such that |$\int _{\varOmega } \mathscr{A}\tau : \tau ~\textrm{d}x \geq \alpha \| \tau \|_{0}^{2}$| for all |$\tau \in L^{2}(\varOmega ;\mathbb{S})$|⁠. The desired inequality may now be derived by judiciously selecting the constants |$C_{1}, C_{3}, C_{4}$| (typically by choosing |$C_{1} \gg C_{2} \gg \textrm{max} (C_{3},C_{4})$|⁠) and using the weighted Young inequality. For concreteness, constants we might pick are

$$ \begin{align} C_{1} = \left(\frac{\lambda_{\kappa}^{\gamma} C^{2}_{\Sigma} \| \Lambda\|^{2} (n+1)^{2} }{\alpha} + 2\right)\| A_{k}\|^{2}, \qquad C_{3} = C_{4} = \lambda_{\kappa}^{\gamma}.\end{align} $$

(3.32)

With this choice and our Poincaré-type inequality from Lemma 3.2, combined with the inequality (3.26) we may derive

$$ \begin{align}\mathscr{G}\left( \mathfrak{p}, \mathfrak{q}; \mathfrak{s}, \mathfrak{u} \right) &\geq 2\alpha\|A_{k}\|^{2}\|\tau\|_{0}^{2} + \frac{\lambda_{\kappa}^{\gamma}}{6}\|A_{k}\|^{2} \| (\tilde{v}, v)\|_{Q}^{2} \nonumber\\&\quad+ \frac{ \lambda_{\kappa}^{\gamma}}{4} \left( \sum_{i} \left\| -c_{i}^{k} \nabla \mu_{i} + \omega_{i}^{k} \nabla p \right\|_{0}^{2} + \|\! \operatorname{div} \tau - \nabla p \|^{2}_{0} \right) \nonumber\\ & \gtrsim | \mathfrak{p} |_{\Theta^{k}}^{2} + \|\mathfrak{q}\|_{Q}^{2} \gtrsim \| (\mathfrak{p}, \mathfrak{q}) \|_{\Theta^{k} \times Q}^{2} \gtrsim \| (\mathfrak{p}, \mathfrak{q}) \|_{\Theta^{k} \times Q} \| (\mathfrak{s}, \mathfrak{u}) \|_{\Theta^{k} \times Q},\end{align} $$

(3.33)

which is the statement of the Babuška condition (3.25).

 

Remark 3.5

The analysis presented here relies heavily on Assumption 3.1, that |$c_{i} \geq \kappa \geq 0$|⁠. This assumption may be brought to the edge of validity (or |$\kappa $| may become very small) in many situations, such as near the intakes of two mixing substances as is explored in the second numerical example, Section 4.5. In the case of ideal gas diffusion decoupled from convection, the influence of |$\kappa $| was explored in Van-Brunt et al. (2021) where it was found that |$\lambda _{\kappa } = \textrm{O}(\kappa ^{-1})$| and that numerical robustness was maintained for very small values of |$\kappa $|⁠.

4. Discretization and numerical experiments

We now assume that |$\varOmega $| is polytopal, and admits a quasi-uniform triangulation |$\mathscr{T}_{h}$| with simplicial elements of maximal diameter |$h$|⁠. Denote conforming finite element spaces for the discrete solution tuple by

$$ \begin{align}\underbrace{\left(X_{h}^{n}\times \Sigma_{h}\times P_{h}\right)}_{=:~\Theta^{k}_{h}}\times\underbrace{\left(W_{h}^{n}\times V_{h}\right)}_{=:~Q_{h}} \subset \underbrace{\left(L^{2}_{0}(\Omega)^{n} \times L^{2}(\Omega;\mathbb{S})\times L^{2}_{0}(\Omega)\right)}_{\supset~\Theta^{k}}\times \underbrace{(L^{2}(\Omega;\mathbb{R}^{d})^{n}\times L^{2}(\Omega; \mathbb{R}^{d}))}_{=~Q}.\end{align} $$

(4.1)

Here |$\varTheta ^{k}_{h}$| is independent of |$k$| as a set, but inherits an iteration-dependent norm described below; |$Q_{h}$| inherits the norm of |$Q$|⁠. Our discretized linear problem after |$k\geq 0$| nonlinear iterations therefore reads: seek |$((\tilde{\mu }_{h}, \tau _{h}, p_{h}), (\tilde{v}_{h}, v_{h})) \in \varTheta ^{k}_{h}\times Q_{h}$| such that

$$ \begin{align}\begin{aligned}\Lambda(\tilde{\mu}_{h}, \tau_{h}, p_{h}; \tilde{w}_{h}, s_{h}, q_{h}) + B_{k,h}(\tilde{w}_{h}, s_{h}, q_{h}; \tilde{v}_{h}, v_{h}) &= \ell^{1}_{k,h}(\tilde{w}_{h}, s_{h}, q_{h})&&\forall(\tilde{w}_{h}, s_{h}, q_{h})\in \Theta^{k}_{h}, \\ B_{k,h}(\tilde{\mu}_{h}, \tau_{h}, p_{h}; \tilde{u}_{h}, u_{h}) - A_{k,h}(\tilde{v}_{h}, v_{h}; \tilde{u}_{h}, u_{h}) &= \ell^{2}_{k,h}(\tilde{u}_{h}, u_{h})&&\forall(\tilde{u}_{h}, u_{h})\in Q_{h}.\end{aligned}\end{align} $$

(4.2)

The terms where |$A_{k,h}, B_{k,h}$| are obtained from |$A_{k}, B_{k}$| and |$\ell ^{1}_{k,h}, \ell ^{2}_{k,h}$| from |$\ell ^{1}_{k}, \ell ^{2}_{k}$|⁠, respectively, by replacing the discretely computed concentrations |$c^{k}_{i}$| and inverse density |$(\rho ^{k})^{-1}$| with discrete approximations; we use these to define a norm |$\|\cdot \|_{\varTheta ^{k}_{h}}$| for |$\varTheta ^{k}_{h}$| in analogy to (3.2).⁹ In block form, the linearized discrete problem reads

$$\ $$

(4.3)

4.1 Structure-preservation and well-posedness

We have already argued the need for pressure regularity greater than |$L^{2}$|⁠. We therefore employ the continuous Lagrange element of degree |$t\geq 1$|⁠, |$P_{h} = \textrm{CG}_{t}(\mathscr{T}_{h})$|⁠, for the pressure. It appears natural to take the chemical potential space |$X_{h}$| to be CG elements of at least the same degree |$r \geq t$|⁠, |$X_{h} = \textrm{CG}_{r}(\mathscr{T}_{h})$|⁠, from the diffusion driving forces (1.14), and since the chemical potentials are used to calculate the concentrations, which one would like to approximate to high order due to their physical importance.

The mass-average velocity constraint (1.6) suggests that the species velocity space be contained in the space used for convective velocity, |$W_{h}\subset V_{h}$|⁠. For the Stokes subsystem, it is desirable that the pair |$(\varSigma _{h}\times P_{h}, V_{h})$| be inf-sup compatible, for which it is sufficient to have that the full ‘divergence of the Cauchy stress’, |$\sigma = (\tau , p)\mapsto \operatorname{div} \tau - \nabla p$|⁠, is surjective onto |$V_{h}$|⁠. By the regularity choice (2.5) for the pressure, it is thus natural to apply |$\operatorname{div}$|-conforming tensor elements to discretize the viscous stress. The conservation of angular momentum is equivalent to the symmetry of the Cauchy stress, which by the decomposition (1.20) is equivalent to that of the viscous stress; for now, we consider exactly symmetric elements for the viscous stress (such as those proposed in Arnold & Winther, 2002; and Arnold et al., 2008) since this obviates the need for a symmetry-enforcing Lagrange multiplier which would add a further field to our |$(2n + 3)$|-field formulation.

However, if one would like to repeat at the discrete level the proof of Theorem 3.4, it is necessary for |$\operatorname{div}:\varSigma _{h}\to V_{h}$| to be surjective, allowing us to construct the discrete analogue of the tensor field |$s_{v}$| in the ansatz (3.27). This is stronger than surjectivity of |$(\tau , p)\mapsto \operatorname{div}\tau - \nabla p, \varSigma _{h}\times P_{h}\to V_{h}$|⁠, but in practice is equivalent because many appropriate choices of |$\varSigma _{h}$| are designed to discretize the full Cauchy stress. Furthermore, the discrete analogue of the constant |$C_{\varSigma }$| (and hence the resulting inf-sup constant) will a priori depend on |$h$|⁠; it is therefore necessary to assume that such |$s_{v}$| can be constructed stably.

 

Assumption 4.1

(Stable right inverse for the divergence.) There exists |$C^{\varSigma }$| independent of |$h$| such that for each |$u_{h}\in V_{h}$| and for the unique |$s_{h}\in \varSigma _{h}/\ker (\operatorname{div})$| with |$\operatorname{div} s_{h} = u_{h}$|⁠, there holds |$\|s_{h}\|_{H(\operatorname{div};\mathbb{S})}\leq C^{\varSigma }\|u_{h}\|_{0}$|⁠.

This is true for (for example) stress elements discretizing an elasticity complex which admits bounded commuting projections to the subcomplex, as is for example the case for the Arnold–Winther and Arnold–Awanou–Winther elements (products of the finite element exterior calculus (Arnold, 2018, Sec. 8.8)), as proved in Arnold & Winther (2002) and Arnold et al. (2008) by explicit construction of such projections. The other relations are summarized below as,¹⁰

$$\ $$

(4.4)

where |$\twoheadrightarrow $| indicates surjectivity. The bottom row corresponds to the final segment of a discrete stress elasticity complex, with stress space refined for Stokes flow due to the decomposition (1.20). We will need the conditions of Lemma 3.2 to hold exactly in the discretization. This will in general require that we approximate the concentrations, |$c^{k}_{i}$|⁠, and density reciprocal, |$(\rho ^{k})^{-1}$|⁠, in specific discrete function spaces. The interpolation of these terms will be denoted by |$c^{k}_{i,h}$| and |$\rho ^{k,-1}_{h}$|⁠, respectively.

Finally, to show well-posedness of the discrete problem, we require an additional condition which does not fit neatly onto (4.4).

 

Assumption 4.2

(Discrete driving force.) The operator given by

$$ \begin{align}& \textrm{d}_{h}^{i,k}(w_{h}, q_{h}) := -c_{i,h}^{k}\nabla w_{h} + \omega^{k}_{i,h}\nabla q_{h},\end{align} $$

(4.5)

where |$\omega ^{k}_{i,h} := M_{i} c^{k}_{i,h} \rho ^{k,-1}_{h}$|⁠, is well-defined as a map |$\textrm{d}_{h}^{i,k}:X_{h}\times P_{h}\to W_{h}$|⁠, i.e., it takes values in |$W_{h}$|⁠.

 

Remark 4.3

Note that Lemma 3.2 required the gradient of |$\rho ^{-1}$|⁠, and so |$\rho ^{k,-1}_{h}$| should at least be a continuous piecewise linear function. In order to minimize the polynomial degree for |$W_{h}$| arising from Assumption 4.2, it is advantageous to interpolate |$c^{k}_{i,h}$| onto the space |$\textrm{DG}_{0}$|⁠. These coefficients do not affect the accuracy of the discretization, only the quality of the linearization, and nonlinear convergence appears robust regardless of this approximation.

 

Theorem 4.4

(Discrete inheritance of well-posedness.) Under Assumptions 4.1 and 4.2 and the relations specified in (4.4), the system (4.2) is well-posed, uniformly in |$h$|⁠.

 

Proof.

Due to the structural conditions demanded in the Assumptions, by inspection the choices of test functions (3.27) are valid. As a consequence we may completely replicate the argument presented in the infinite-dimensional case, and derive the analogue of condition (3.25) with constant independent of |$h$|⁠.

4.2 Error estimates

Following Xu & Zikatanov (2003, Theorem 2), for fixed |$k$| we infer the abstract error estimate

$$ \begin{align}& \left\|(\underline{\tilde{\mu}}^{k+1}, \underline{\tau}^{k+1}, \underline{p}^{k+1}) - (\tilde{\mu}_{h}, \tau_{h}, p_{h})\right\|_{\varTheta^{k}_{h}} + \left\|(\underline{\tilde{v}}^{k+1}, \underline{v}^{k+1}) - (\tilde{v}_{h}, v_{h})\right\|_{Q}\lesssim \mathbb{E}_{\varTheta^{k}_{h}} + \mathbb{E}_{Q_{h}},\end{align} $$

(4.6)

where

$$ \begin{align}\begin{aligned}\mathbb{E}_{\Theta^{k}_{h}} &:= \inf_{(\tilde{w}_{h},s_{h},q_{h})\in\Theta^{k}_{h}}\left\|(\underline{\tilde{\mu}}^{k+1}, \underline{\tau}^{k+1}, \underline{p}^{k+1}) - (\tilde{w}_{h},s_{h},q_{h})\right\|_{\Theta^{k}_{h}},\\ \mathbb{E}_{Q_{h}} &:= \inf_{(\underline{\tilde{u}}_{h}, u_{h})\in Q_{h}}\left\|(\underline{\tilde{v}}^{k+1}, \underline{v}^{k+1}) - (\tilde{u}_{h}, u_{h})\right\|_{Q}.\end{aligned}\end{align} $$

(4.7)

Here the tuple |$((\underline{\tilde{\mu }}^{k+1}, \underline{\tau }^{k+1}, \underline{p}^{k+1}), (\underline{\tilde{v}}^{k+1}, \underline{v}^{k+1}))$| is defined as the exact solution to (3.5), but with |$B_{k}, A_{k}, \ell ^{1}_{k}, \ell ^{2}_{k}$| replaced with |$B_{k, h}, A_{k, h}, \ell ^{1}_{k, h}, \ell ^{2}_{k, h}$|⁠, respectively—that is, the solution of the system (3.5) with the estimates of the concentrations and inverse density replaced by |$c_{i,h}^{k}$| and |$\rho ^{k,-1}_{h}$|⁠, respectively.

To derive a practical error estimate, we will need to bound the quantities |$\mathbb{E}_{\varTheta ^{k}_{h}}$| and |$\mathbb{E}_{Q_{h}}$| by interpolation estimates specific to the choice of finite element spaces, by estimating |$\|\cdot \|_{\varTheta ^{k}_{h}}, \|\cdot \|_{Q}$| in terms of standard Sobolev norms. To this end we can readily check that

$$ \begin{align}\begin{aligned}\mathbb{E}_{\Theta^{k}_{h}}&\lesssim \max\left(1, \sum_{i} \left\|c^{k}_{i,h} \right\|_{L^{\infty}(\Omega)}\right)\inf_{\tilde{w}_{h} \in X_{h}^{n}} \left\| \underline{\tilde{\mu}}^{k+1} - \tilde{w}_{h} \right\|_{1} \\ &+ \max\left(1, \sum_{i} \left\|\omega^{k}_{i,h}\right\|_{L^{\infty}(\Omega)}\right) \inf_{q_{h} \in P_{h}} \left\| \underline{p}^{k+1} - q_{h} \right\|_{1} + \inf_{s_{h} \in \Sigma_{h}} \left\|\underline{\tau}^{k+1} - s_{h} \right\|_{H(\operatorname{div};\mathbb{S})}\!, \\ \mathbb{E}_{Q_{h}} &\lesssim \inf_{\tilde{u}_{h} \in W_{h}^{n}} \left\| \underline{\tilde{v}}^{k+1} - \tilde{u}_{h} \right\|_{0} + \inf_{u_{h} \in V_{h} } \left\| \underline{v}^{k+1} - u_{h} \right\|_{0}.\end{aligned}\end{align} $$

(4.8)

4.3 Examples of suitable finite elements

Having derived abstract error estimates, we now consider two simple combinations of finite elements satisfying the structural conditions (4.4) and Assumptions 4.1 and 4.2.

The design and implementation of stress elements which exactly enforce symmetry and |$\operatorname{div}$|-conformity is notoriously difficult; in 2D, one choice of such elements is the conforming Arnold–Winther element (Arnold & Winther, 2002), recently incorporated into the Firedrake finite element library (Rathgeber et al., 2016; Aznaran et al., 2022a). In the lowest-order case we denote this element by |$\textrm{AW}^{c}_{3}$|⁠. Specifying

$$ \begin{alignat}{2} X_{h} = P_{h} &= \textrm{CG}_{1}(\mathscr{T}_{h})\cap L^{2}_{0}(\varOmega), \end{alignat} $$

(4.9a)

$$ \begin{alignat}{2} \varSigma_{h} &= \textrm{AW}^{c}_{3}(\mathscr{T}_{h};\mathbb{S}), \end{alignat} $$

(4.9b)

$$ \begin{alignat}{2} W_{h} = V_{h} &= \textrm{DG}_{1}(\mathscr{T}_{h};\mathbb{R}^{d}), \end{alignat} $$

(4.9c)

and further assuming that the discretely computed |$c^{k}_{i}$| and |$(\rho ^{k})^{-1}$| have been interpolated into |$\textrm{DG}_{0}$| and |$\textrm{CG}_{1}$|⁠, respectively, then this discretization satisfies the structural properties (4.4) and Assumptions 4.1 and 4.2, hence we deduce the error estimate (4.6).

Let |$\varPi _{h}^{\textrm{CG}_{1}}:H^{2}(\varOmega )\to \textrm{CG}^{1}(\mathscr{T}_{h}), \varPi ^{\textrm{AW}^{c}_{3}}_{h}:H^{1}(\varOmega ;\mathbb{S})\to \textrm{AW}^{c}_{3}(\mathscr{T}_{h};\mathbb{S})$|⁠, and |$\varPi _{h}^{\textrm{DG}_{1}^{d}}:H^{1}(\varOmega ;\mathbb{R}^{d})\to \textrm{DG}_{1}(\mathscr{T}_{h};\mathbb{R}^{d})$| be the associated interpolation operators. We then have the following estimates under sufficient regularity assumptions (for details we refer to Arnold & Winther, 2002; Logg et al., 2012, Ch. 3; and Boffi et al., 2013, p. 72):

$$ \begin{alignat}{2} &\left\|\tilde{\mu} - \varPi^{\textrm{CG}_{1}}_{h} \tilde{\mu}\right\|_{1} \lesssim h |\tilde{\mu} |_{2}, \end{alignat} $$

(4.10a)

$$ \begin{alignat}{2} &\left\| p - \varPi_{h}^{\textrm{CG}_{1}}p \right\|_{1} \lesssim h |p|_{2}, \end{alignat} $$

(4.10b)

$$ \begin{alignat}{2} &\left\| \tau - \varPi^{\textrm{AW}^{c}_{3}}_{h} \tau \right\|_{0} + h \left\| \operatorname{div} \left(\tau - \varPi_{h}^{\textrm{AW}^{c}_{3}} \tau\right) \right\|_{0} \lesssim h^{2} |\tau|_{2}, \end{alignat} $$

(4.10c)

$$ \begin{alignat}{2} &\left\| (\tilde{v}, v) - \varPi^{\textrm{DG}_{1}^{d}}_{h}(\tilde{v}, v) \right\|_{Q} \lesssim h^{2} | (\tilde{v}, v) |_{1}, \end{alignat} $$

(4.10d)

where |$\varPi ^{\textrm{CG}_{1}}_{h}, \varPi ^{\textrm{DG}_{1}^{d}}_{h}$| have been applied component-wise. Using these estimates for the interpolation operators and the error estimate (4.6), we can derive the error bound

$$ \begin{align}& \left\|(\underline{\tilde{\mu}}^{k+1}, \underline{\tau}^{k+1}, \underline{p}^{k+1}) - (\tilde{\mu}_{h}, \tau_{h}, p_{h})\right \|_{\varTheta^{k}_{h}} + \left\| (\underline{\tilde{v}}^{k+1}, \underline{v}^{k+1}) - (\tilde{v}_{h}, v_{h}) \right\|_{Q} \lesssim h.\end{align} $$

(4.11)

We will see in practice that the order of convergence for several fields is actually higher, but the error of the species velocities and the driving forces remains |$\textrm{O} (h)$|⁠.

A second class of finite elements may be found by replacing (4.9a) with

$$ \begin{align}& X_{h} = \textrm{CG}_{2}(\mathscr{T}_{h})\cap L^{2}_{0}(\varOmega), P_{h} = \textrm{CG}_{1}(\mathscr{T}_{h})\cap L^{2}_{0}(\varOmega),\end{align} $$

(4.12)

and leaving the others unchanged. Again, the structural conditions are satisfied if |$c^{k}_{i}$| and |$(\rho ^{k})^{-1}$| are interpolated into |$\textrm{DG}_{0}$| and |$\textrm{CG}_{1}$|⁠, respectively. A similar error analysis again confers an error bound of |$\textrm{O} (h)$|⁠, though shortly we will see that this order is actually higher in practice.

 

Remark 4.5

These estimates bound the error between the discrete solutions at iteration |$k + 1, ((\tilde{\mu }_{h}, \tau _{h}, p_{h}), (\tilde{v}_{h}, v_{h}))$| and the continuous solution of the nonlinear scheme |$((\underline{\tilde{\mu }}^{k+1}, \underline{\tau }^{k+1}, \underline{p}^{k+1}), (\underline{\tilde{v}}^{k+1}, \underline{v}^{k+1}))$| with the same (discrete) coefficients. In principle this is incomplete, as ideally we would derive error estimates against the continuous solution |$((\tilde{\mu }^{k+1}, \tau ^{k+1}, p^{k+1}), (\tilde{v}^{k+1}, v^{k+1}))$| at iteration |$k + 1$| with the exact (continuous) coefficients. Estimates on such consistency errors were analysed for a simpler system in Van-Brunt et al. (2021) and some rationale was provided as to why in practice this is not an issue, based on the formal order of the consistency error being strictly less than the discretization error. We expect a similar (if laborious) analysis could be performed, following the strategy in Van-Brunt et al. (2021). In general, the consistency errors are expected to be |$\textrm{O}(h^{2})$|⁠, which will be borne out in the numerical examples.

 

Remark 4.6

We emphasize that we have aimed to identify appropriate structural conditions between finite element spaces in order to preserve desirable properties of the SOSM system—in particular, conditions which allow mimicry of well-posedness proofs from the infinite-dimensional problem—rather than to advocate specifically for the elements used here. We expect it is possible to use Lagrange multipliers to weakly enforce the symmetry of the viscous stress, which would allow for the choice of higher polynomial degrees.

The system matrix of our discrete linearized system (4.3) has symmetric perturbed saddle point structure, and although indefinite, is such that both the blocks |$\varLambda , A_{k,h}$| are symmetric positive semidefinite. These matrix properties hold independently of the particular material considered. We expect that this structure should be conducive to the development of fast preconditioners.

4.4 Validation with manufactured solutions

We now verify our scheme, implemented in Firedrake (Rathgeber et al., 2016). Firedrake currently only supports symmetry-enforcing stress elements in 2D, and we thus restrict ourselves to the case |$d=2$|⁠. Throughout these experiments we chose the penalty parameter |$\gamma = 0.1$| and constant functions for the initial guesses, and the linear systems were solved using the sparse LU factorization of MUMPS (Amestoy et al., 2001) via PETSc (Balay et al., 2019).

To validate our error estimates, we construct a manufactured solution for an ideal gas on the unit square |$\varOmega = (0,1)^{2}$|⁠. If |$RT = 1$|⁠, the diffusion coefficients are given by |$\mathscr{D}_{ij} = D_{i}D_{j}$| for |$D_{j}> 0$|⁠, and |$g:\mathbb{R}^{2}\to \mathbb{R}$| is smooth, then one can check that an analytical solution to the OSM subsystem (1.10) is given by

$$ \begin{align}& c_{i} = \exp\left(\frac{g}{D_{i}}\right), \quad v_{i} = D_{i}\nabla g,\end{align} $$

(4.13)

which together implicitly define every other quantity (for given shear and bulk viscosities) apart from the chemical potentials. We compute the latter by inverting the ideal gas relation (1.26) with |$p^{\ominus } = \int{\!\!\!\!\!-} _{\varOmega } p~\textrm{d}x, \mu _{i}^{\ominus } = \int{\!\!\!\!\!-} _{\varOmega } c_{i}~\textrm{d}x~\forall i$|⁠.¹¹ The molar mass of each species was set to |$1$|⁠, and |$\zeta , \eta $| to |$0.1$|⁠. The initial guesses for the concentrations |$\tilde{c}^{0}$| were set as |$c^{0}_{i} = \int{\!\!\!\!\!-} _{\varOmega } c_{i}~\textrm{d}x$|⁠, i.e. as the average of the exact concentrations.

We used |$D_{i} = \frac{1}{2} + \frac{i}{4}, i = 1, 2, 3$|⁠, and |$g(x,y) = \frac{xy}{5}(1 - x)(1 - y)$| to generate Figure 1, the log-log error plot for the overall algorithm, which demonstrates the negligible effect of the consistency error induced by the discrete interpolations |$c^{k}_{i,h}, \rho ^{k,-1}_{h}$|⁠, and verifies the error estimate (4.11) (for the 0^th-order terms).

Error plots for two finite element families: (4.9) (left) and (4.12) (right).

Fig. 1.

The tolerance in the outer solver was |$10^{-7}$| in the |$\ell ^{2}$| norm, and took |$6$| nonlinear iterations on the coarsest mesh of |$4 \times 4$|⁠, to |$7$| iterations on finest mesh of |$32 \times 32$|⁠. We have denoted in Figure 1|$d_{i,h}$| as the discrete driving force defined by (4.5) at the final iteration, and |$\sigma _{h} := \tau _{h} - p_{h} \mathbb{I}$|⁠. Note that we observe |$\textrm{O}(h^{2})$| convergence in the |$L^{2}$| norms of the chemical potential, stress and pressure. As in Van-Brunt et al. (2021), this suggests that the error estimates could be improved, for example by duality arguments.

We now turn to the higher-order terms in (4.11). Due to our construction of the function spaces (3.1) for the linearized problem, it is the norm |$\|\cdot \|_{\varTheta ^{k}_{h}}$| with respect to which we have proved convergence of the solution tuple. It is natural to wonder whether this is an artefact of our constructed function spaces. To answer this, we measure convergence of the chemical potential gradients |$\nabla \mu _{i}$|⁠, pressure gradient |$\nabla p$| and divergence |$\operatorname{div}\tau $| of the nonequilibrium stress to their true values, compared to the convergence of the nonlinear diffusion driving forces and the divergence of the full Cauchy stress. For the former quantities, there is a priori no reason to expect any convergence.

Remarkably, we observe in Figure 2 that (at least for the higher-order family) not only do the components |$\nabla \mu _{i}, \nabla p, \operatorname{div}\tau $| converge, but in fact there is convergence of the nonlinear diffusion driving forces |$d_{i}$| and of the divergence of the full Cauchy stress |$\operatorname{div}\sigma $|⁠, and at a rate one order higher than these individual components |$\nabla \mu _{i}, \nabla p, \operatorname{div}\tau $|⁠; this suggests that, rather than being a mathematical artefact of our formulation, the conditions defining the |$\varTheta ^{k}$| space capture the underlying thermodynamic quantities of interest. This also provides intriguing (if circumstantial) evidence towards the physical relevance of our nonlinear formulation in Definition 2.1.

Higher-order convergence in |$L^{2}$| of the divergence of the full Cauchy stress, and driving forces, for the finite element family (4.12).

Fig. 2.

4.5 Benzene-cyclohexane mixture

Cyclohexane (CH) is important in the petrochemical industry as it is used to synthesize a variety of products, such as nylon. It is primarily produced through the hydrogenation of benzene (CH), resulting in a benzene-cyclohexane mixture. Separation of cyclohexane from this mixture is difficult due to their similar vaporization temperatures (Villaluenga & Tabe-Mohammadi, 2000). Since liquid mixtures of these components are important in the chemical industry, most of the required thermodynamic and dynamical property data are readily available in the literature. Because it provides a tractable nonideal example for which a complete set of material properties is known, we consider here a microfluidic chamber in which Stokes flows of liquid benzene and cyclohexane mix.

The required transport parameters (measured at |$298.15 \: \textrm{K}$|⁠) may be found in Guevara-Carrion et al. (2016). We observe from these data that the Stefan–Maxwell diffusivity and the shear viscosity are both approximately constant with respect to composition and pressure, and will be approximated as |$\mathscr{D}_{12} = 2.1 \times 10^{-9}~\textrm{m}^{2}/\textrm{s}$| and |$6 \times 10^{-4}~\textrm{Pa} \cdot \textrm{s}$|⁠, respectively. Lacking accurate data for the bulk viscosities of either benzene or cyclohexane, we set them to be effectively zero, choosing |$\zeta = 10^{-7}~\textrm{Pa} \cdot \textrm{s}$|⁠. (Numerical experiments confirmed that a value of this order has no discernible impact on the output of the simulation.) The molar masses used in the simulation are |$0.078~\textrm{kg}/\textrm{mol}$| for benzene and |$0.084~{\textrm{kg}/\textrm{mol}}$| for cyclohexane. The ambient pressure was taken as |$p^{\ominus } = 10^{5}~\textrm{Pa}$|⁠.

Although benzene and cyclohexane are fully miscible liquids, they form a nonideal solution. Information relating the chemical potentials to the mole fractions is therefore required. This is accomplished using a Margules model (Green & Perry, 2007) for activity coefficients, the parameters of which were reported by Tasić et al. (1978) as |$A_{{C_{6}H_{6}}, {C_{6}H_{12}}} = 0.4498$| and |$A_{{C_{6}H_{12}},{C_{6}H_{6}}} = 0.4952$|⁠. This well-known model parameterizes activity coefficients in terms of a minimal set of functions which maintain thermodynamic rigour.

To aid convergence, we use an under-relaxation scheme with respect to the concentrations in our nonlinear solver, with a relaxation parameter of |$0.1$|⁠. That is, we update the concentration as |$c^{*,k+1}_{i}$| where

$$ \begin{align}& c_{i}^{*,k+1} := 0.9c_{i}^{k} + 0.1c_{i}^{k+1}.\end{align} $$

(4.14)

This is necessary due to stiffness of the problem, which owes to the fact that the mixtures are essentially fully separated at the inlets to the apparatus.

To calculate the total concentration of the mixture we use

$$ \begin{align}& c_{\textrm{T}} = \frac{c_{{C_{6}H_{6}}}^{\textrm{ref}} c_{{C_{6}H_{12}}}^{\textrm{ref}}}{ x_{{C_{6}H_{6}}}c_{{C_{6}H_{12}}}^{\textrm{ref}} +x_{{C_{6}H_{12}}} c_{{C_{6}H_{6}}}^{\textrm{ref}}},\end{align} $$

(4.15)

where |$c_{-}^{\textrm{ref}}$| denotes the concentration (inverse molar volume) of the pure species at |$10^{5} \: \textrm{Pa}$| and |$298.15 \: \textrm{K}$|⁠, approximately |$9.20 \:\: \textrm{mol} \:\textrm{L}^{-1}$| and |$11.23 \:\: \textrm{mol} \: \textrm{L}^{-1}$| for cyclohexane and benzene, respectively. (4.15) is derived from (1.28) under the assumption that the partial molar volumes of the two components are independent of the solution’s composition.

We consider a two-dimensional channel configuration where two inlets converge into a single outlet in a truncated ‘Y’ configuration. A diagram of the mesh is shown in Figure 3.

Mesh used in the benzene-cyclohexane mixing simulation. Each unit along the axes corresponds to a physical distance of |$2$|mm.

Fig. 3.

At the top inlet, pure benzene enters and at the bottom pure cyclohexane, at speeds |$v_{{C_{6}H_{6}}}^{\textrm{ref}}$| and |$v_{{C_{6}H_{12}}}^{\textrm{ref}}$|⁠, respectively. Rather than symmetrize these speeds, superior mixing results are obtained by symmetrizing the molar fluxes at the inlets. In other words, we impose the condition

$$ \begin{align}& c_{{C_{6}H_{6}}}^{\textrm{ref}} v_{{C_{6}H_{6}}}^{\textrm{ref}} = c_{{C_{6}H_{12}}}^{\textrm{ref}} v_{{C_{6}H_{12}}}^{\textrm{ref}}.\end{align} $$

(4.16)

We will specify that |$v_{{C_{6}H_{6}}}$| enters at a speed of |$4 \mu \textrm{m} \: \textrm{s}^{-1}$|⁠. This prescribes an inlet speed for cyclohexane of approximately |$3.28 \mu \textrm{m} \: \textrm{s}^{-1}$|⁠. A parabolic profile across each inlet and outlet is imposed, as this is consistent with the no-slip condition and the characteristics of a plane Poiseuille flow. Concretely, this means the boundary conditions at the inlets are the normal component of

$$ \begin{align}& c_{i} v_{i} = - 2x(x - 1) c_{i}^{\textrm{ref}} \begin{bmatrix} 2 \\ \pm 1 \end{bmatrix}\!,\end{align} $$

(4.17)

where |$i = {C_{6}H_{6}}$| or |${C_{6}H_{12}}$| and the sign in the 2nd entry is negative for benzene and positive for cyclohexane. For both species, the normal conditions at the outlet are

$$ \begin{align}& c_{i} v_{i} = - 2(y-0.5)(y + 0.5)c_{i}^{\textrm{ref}} \begin{bmatrix} 1 \\ 0 \end{bmatrix}.\end{align} $$

(4.18)

Results for the fields computed by the simulation are visualized in Figs 4 and 5. We may observe that the pressure profile is smooth. We also note that, although the mass-average velocity exhibits rather simple, predictable behaviour, the flow fields of the individual species are significantly more complex and that these three flow profiles are cleanly distinguished. We see that both species develop convective rolls—behaviour markedly different from the bulk velocity.

Plot of change in pressure in the mixing chamber, with streamlines computed from the mass-average velocity.

Fig. 4.

Concentrations of benzene (left) and cyclohexane (right), with streamlines computed from their velocities.

Fig. 5.

 

Remark 4.7

The concentrations of each substance become very small near the inlet of the other, with benzene reaching a concentration of less than |$10$| mmol/L. Although convergence was achieved despite the value of |$\kappa $| being consequently very small, further investigation is needed to ensure robust convergence for vanishing |$\kappa $| at large Péclet numbers. This stands in contrast to the pure diffusion problem considered in Van-Brunt et al. (2021).

4.6 Code availability

For reproducibility, the exact software versions used to generate the numerical results in this paper are archived on Zenodo (Aznaran et al., 2022b); our implementation employs a nondimensionalization of the SOSM system. All code containing further implementation details of the microfluidic example as well as the mesh used and scripts for the associated plots, are available at https://bitbucket.org/FAznaran/sosm-numerics/.

5. Conclusions and outlook

We have proposed a formulation and numerical method for the Stokes–Onsager–Stefan–Maxwell equations of multicomponent flow, proving continuous and discrete inf-sup conditions for a linearization of the system with saddle point structure. This structure arises from the duality between the diffusion driving forces, and the combination of species continuity equations with the divergence of the mass-average velocity constraint. We obtained error estimates in a norm corresponding to a space requiring square-integrable diffusion driving forces and total stress divergence, and verified these with numerical experiments.

This work represents a first step towards the simulation of nonideal mixtures; further physical extensions will be required for realistic applications in chemical engineering. Of particular interest are the analysis of the transient problem, the incorporation of thermal effects based on the framework proposed in Van-Brunt et al. (2022), the weak enforcement of symmetry for the viscous stress tensor to ease the extension of the method to three dimensions (Boffi et al., 2013, Section 9.2), the consideration of the case of vanishing bulk viscosity as encountered in dilute monatomic gases (Hirschfelder et al., 1954), and the extension to NSOSM transport (i.e. non-negligible Reynolds number).

Rigorous investigation into a notion of weak solution more refined than Definition 2.1 incorporating (for example) integrability of thermodynamic pressure gradients, and appropriate Picard or Newton linearization, would also be of significant interest. We also remark that a proof of convergence of the Picard iteration could be used to prove the existence of a solution tuple for the infinite-dimensional nonlinear SOSM system.

Funding

This work was supported by the Engineering and Physical Sciences Research Council Centre for Doctoral Training in Partial Differential Equations: Analysis and Applications (EP/L015811/1); the Engineering and Physical Sciences Research Council (grants EP/R029423/1 and EP/W026163/1); the MathWorks, Inc., the Clarendon fund scholarship and the Faraday institution SOLBAT project; and Multiscale Modelling projects (subawards FIRG007 and FIRG003 under grant EP/P003532/1).

Acknowledgements

The authors are also grateful to J. Málek and M. P. Juniper for useful comments, and to the two anonymous referees for their valuable suggestions.

Footnotes

1

In multiphase flows, heterogeneous material exchange may occur, leaving a nonzero generation term in the mass continuity equation for a given phase. We limit the present discussion to single-phase flows.

2

In general, an internal stress |$\tau $| is characterized by |$\int _{\partial M} \tau ~\textrm{d}s$| expressing the net force exerted on the surroundings by a volume |$M \subset \varOmega $| on the closed surface |$\partial M$| that bounds it.

3

Extensive (but unreported) investigations into alternative tuples of fields for which to solve, and ways to weakly formulate the resulting fully coupled system, gave rise to ill-posed or analytically intractable Picard linearizations when the concentrations |$c_{i}$|⁠, or their normalizations the mole fractions |$x_{i} = c_{i}/c_{\textrm{T}}$|⁠, were solved for as primary unknowns.

4

The possibility of |$\gamma = 0$| is allowed as a means to recover the original equations, but our analysis will rely on taking |$\gamma> 0$|⁠.

5

In any case, the incompressible regime for which |$\rho $| is constant is physically irrelevant to the OSM framework for mass diffusion, which exhibits spatial heterogeneity of the density. We also remark that, viewing the pressure as a component of the full Cauchy stress, appealing to the Hodge decomposition of the stress space |$H(\operatorname{div};\mathbb{S})$| (Arnold, 2018, Theorem 4.5) does not endow that component with any extra regularity.

6

One alternative approach is provided by attempting to construct a smoother analogue of the stress elasticity complex associated with the Cauchy stress space (2.5), for which (2.5) is replaced by some superspace of |$H(\textrm{div};{\mathbb{S}})\times H^{1}(\varOmega )$|⁠, just as the Stokes complex is precisely a smoothing of the de Rham complex. We do not pursue this.

7

A stronger condition, that |$\rho \in W^{1,\infty }(\varOmega )$| is bounded below with |$\frac{\nabla \rho }{\rho }\in L^{\infty }(\varOmega ;\mathbb{R}^{d})$|⁠, was used to analyze a compressible Stokes flow in Caucao et al. (2016).

8

We conjecture that one could alternatively derive a formulation of the SOSM system dual to ours, which takes |$d_{i}\in L^{2}(\varOmega ;\mathbb{R}^{d})$| as a primary unknown. We also conjecture that the integrability assumptions in Definition 2.1 could potentially be relaxed, for example via Sobolev embeddings.

9

Note that this norm is mesh-dependent only in the sense of depending on the previous discrete solution.

10

Here |$\pi _{i}$| denotes projection onto the |$i$|th component, i.e. we require |$P_{h}\subset X_{h}$|⁠.

11

Physically, a normalization factor should be placed in front of the integral to ensure consistency of units. However, since this example is purely to verify the proposed scheme, units are not emphasized in this section.