Thermodynamically consistent phase-field modelling of activated solute transport in binary solvent fluids

Abstract

Active diffusion of substances in binary immiscible and incompressible fluids with different densities occurs universally in nature and industry, but relevant mathematical models and numerical simulation have been studied scarcely so far. In this paper, a thermodynamically consistent phase-field model is established to describe the activated solute transport in binary fluids with different densities. A mixed free-energy function for multiple solutes is proposed, which leads to different solute chemical potentials in binary solvent fluids, thus it has the ability to characterize the solubility difference of the solutes in two solvents. The two-phase flow is governed by a general hydrodynamic phase-field model that can account for general average velocity and different densities. The proposed model is derived rigorously using the second law of thermodynamics. Moreover, a general multi-component solute diffusion model is established using the Maxwell–Stefan approach, which involves the crossing influences between different solutes. To solve the model effectively, an efficient, linearized and decoupled numerical method is proposed for the model as well. The proposed numerical method can preserve the thermodynamical consistency, i.e. obeying an energy dissipation law at the discrete level, as well as guarantee the mass conservation law for the solutes and solvent fluids. Numerical experiments are carried out to show that the proposed model and numerical method can simulate various processes of the solute active diffusion in two-phase solvent fluids.

1. Introduction

In the normal diffusion process, solutes move from a region of higher concentration to a region of lower concentration if the solution is dilute (Wesselingh & Krishna 2000). The normal diffusion obeys Fick's law of diffusion, so it is also called Fickian diffusion. Fick's law is usually applicable for the solute diffusion in the single-phase fluid that may consist of multiple solvents. When the solvent consists of two-phase immiscible fluids, the solute movement from one fluid with lower solute concentration to the other with higher solute concentration may take place due to different solubility characteristics of the solute in two solvent fluids. This is called active diffusion, which is a universal phenomenon in nature and industry (Hass & Fine 2010; Drechsler et al. 2017; von den Driesch et al. 2020; Guzmán-Lastra, Löwen & Mathijssen 2021; Watson et al. 2020). For example, as an important greenhouse gas, CO$_2$ released through human activities has been raising the global temperature, and in turn the temperature increase makes more CO$_2$ originally stored in oceans transfer into the atmosphere. This thermally activated diffusion process of CO$_2$ influences drastically the global climate change (Watson et al. 2020). In biology, mass transfer between cells is often accomplished via a delicate process of active diffusion and advection (Drechsler et al. 2017). As a common phenomenon occurring in biology and material science, mass transfer through semi-permeable or conducting interfaces (Johansson et al. 2005; Gong, Gong & Huang 2014) can be attributed to specific active diffusion; an example is the regulation of intracellular ion concentrations in living cells despite the changes of concentrations in extracellular spaces. In industry, the removal and extraction of one material can be achieved by the use of its preference for specific solvents (Hass & Fine 2010).

Due to its critical importance, mathematical models and numerical simulation of solute diffusion in multiple solvent fluids have been an attractive topic. In the works of Gong et al. (2014) and Wang et al. (2020a), the immersed boundary methods have been developed to model and simulate mass transfer through permeable interfaces interacting with the surrounding fluids under large moving deformations, and in particular, an additional equation has been introduced to describe the temporal evolution of advective and diffusive fluxes across the interfaces. A thermodynamically consistent phase-field model (Qin et al. 2022) has been proposed to describe mass transfer across permeable moving interfaces; it considers the restricted diffusion problems where the diffusive flux across the interfaces depends on the conductance and concentration difference on each side. The above works follow with interest the active diffusion of mass transfer through porous membranes, in which the membranous interfaces have a critical influence on forming the contrast of concentrations of the substances in internal and external parts of the interfaces. There have been only a few works that use the phase-field models to describe universal active diffusion processes of solutes in two immiscible and incompressible fluids with different densities. This work concerns modelling and numerical simulation of such active diffusion problems using a thermodynamically consistent phase-field model.

The phase-field approach (Anderson, McFadden & Wheeler 1998; Kim 2012) has been employed extensively to model a variety of multiphase fluid and fluid–structure problems through the use of diffuse interfaces (Cahn & Hilliard 1958; Cahn & Allen 1977; Boyer 2002; Shen & Yang 2010; Gomez & Hughes 2011; Shen, Yang & Wang 2013; Guo et al. 2017; Yu & Yang 2017; Zhu et al. 2019). A key ingredient of phase-field models is to introduce proper free-energy functionals to characterize phase properties and behaviours. This allows a phase-field model to satisfy a specific energy dissipation law that is a direct consequence of the second law of thermodynamics for the isothermal system (Kou & Sun 2018a,b). The phase-field approach has been applied successfully to model the solute transport in two-phase flow. One main focus of attention is to model surfactants soluble in two fluids (e.g. oil and water) (van der Sman & van der Graaf 2006; Engblom et al. 2013; Garcke, Lam & Stinner 2014; Zhu et al. 2019). Surfactants can reduce the interfacial tension, thereby greatly affecting the dynamics of binary fluids; for instance, the use of surfactants can enhance the oil recovery rate significantly. van der Sman & van der Graaf (2006) proposed a diffuse interface model with surfactant adsorption based on a free-energy functional. Garcke, Lam & Stinner (2014) derived the comprehensive phase-field models of two-phase flow incorporating a soluble surfactant in a thermodynamically consistent way, which can provide deep understanding of impacts of surfactants on coalescence and segregation of droplets. A thermodynamically consistent model of two-phase flows with soluble surfactants and moving contact line has been derived further by Zhu et al. (2019). Although the solute transport is usually described by the convection–diffusion equations, modelling of the solute active diffusion demands specific free-energy functions. In order to characterize the active diffusion of solutes, we propose a mixed free-energy function for the solutes, which takes different energy parameters in different solvent fluids. This leads to a contrast between the solute chemical potentials in binary fluids despite the same solute concentrations. In other words, at the equilibrium state, the solute chemical potentials are equal in binary solvent fluids, whereas the solute concentration in one solvent fluid may be higher than that in the other fluid. Thus the proposed mixed free-energy has the ability to reflect the solubility difference of the solutes in two fluids. Moreover, we consider the solute composed of multiple substances and establish a general multi-component solute diffusion model using the Maxwell–Stefan approach (Wesselingh & Krishna 2000), which takes into consideration the crossing influences between different solutes.

For real two-phase flow, there is generally a density contrast between two different fluids, so general phase-field models should take this property into consideration. Several phase-field models with different densities have been developed in the literature (Anderson et al. 1998; Lowengrub & Truskinovsky 1998; Boyer 2002; Ding, Spelt & Shu 2007; Shen & Yang 2010; Shen et al. 2013; Liu, Shen & Yang 2015; Roudbari et al. 2018; Haddada & Tierra 2022). There are two kinds of models used popularly in the literature. The first is the model (Abels, Garcke & Grün 2012) that uses the volume-averaged velocity and involves a mass diffusion flux in the momentum balance equation to guarantee the thermodynamical consistency. The others are the models based on the mass-averaged velocity (Shen et al. 2013; Aki et al. 2014; Li & Wang 2014), which are the models used most commonly due to the feature that the standard mass conservation equation can be derived naturally using the mass-averaged velocity. The above models have different formulations due to different choices of averaged velocities. Recently, a general phase-field model has been proposed by Kou et al. (2020b), which features a unified formulation for different averaged velocities, but the solute transport has never been considered. In this paper, the proposed model is established combining the solute transport equations with the general phase-field model (Kou et al. 2020b) to account for two fluids with different densities. Moreover, the proposed model is proved rigorously to follow an energy dissipation law, thus it is thermodynamically consistent.

It is quite challenging to develop efficient energy-stable numerical methods for the phase-field models with different densities. There have been great efforts made for this topic in recent years (Shen & Yang 2010; Grün & Klingbeil 2014; Chen & Shen 2016; Guo et al. 2017; Yu & Yang 2017; Gong et al. 2018; Roudbari et al. 2018; Yang & Zhao 2019). It is recognized universally that numerical methods need to preserve the energy dissipation law at the discrete level so as to remove any restriction on the time step size as well as reduce spurious numerical solutions in practical simulations. In addition, mass conservation is required essentially for numerical methods as well since the large density contrast between binary fluids may cause serious mass loss (Yue, Zhou & Feng 2007; Aland & Voigt 2012; Ding & Yuan 2014; Wang et al. 2015; Guo et al. 2017). The energy-stable and mass-conservative numerical methods (Guo et al. 2017) have been developed for the quasi-incompressible model (Lowengrub & Truskinovsky 1998). An energy-stable finite element numerical scheme (Grün & Klingbeil 2014) has been designed for the model proposed by Abels et al. (2012). Several energy-stable numerical approximations of two classes of hydrodynamic phase-field models (Shen & Yang 2010; Abels et al. 2012) have been proposed by Gong et al. (2018). In the proposed model, the solute concentrations are fully coupled with the phase variable and velocity. Moreover, the model consists of a highly nonlinear system of partial differential equations. Hence it is a challenge to design efficient numerical methods for the proposed model. In this paper, we propose an energy-stable and mass-conservative numerical method for the proposed model through decoupling the tight relationships between the solute concentrations, the phase variable and velocity. Moreover, the proposed scheme is linearized and decoupled, thus it is easy to implement.

We now highlight the main contributions of this paper as follows.

1. (i) A thermodynamically consistent phase-field model for the activated solute transport in binary fluids is derived using the second law of thermodynamics. The solute transport equations are coupled with the general model of binary immiscible solvent fluids with different densities.

2. (ii) A mixed free-energy function for multiple solutes is proposed to characterize the solubility difference of the solutes in two solvents. Different from the solute free energies (van der Sman & van der Graaf 2006; Engblom et al. 2013; Garcke, Lam & Stinner 2014; Zhu et al. 2019), the proposed free energy yields different solute chemical potentials in binary solvent fluids, thus it has the ability to describe the active diffusion processes.

3. (iii) The modified Maxwell–Stefan equations are presented to take into consideration the crossing influences between different solutes, and the resulting diffusion coefficient matrix is symmetric positive definite and consistent with Onsager's reciprocal principle and the second law of thermodynamics.

4. (iv) An efficient, linearized, decoupled and energy-stable numerical method is developed to solve the proposed model. The method is proved to conserve the masses of both solutes and solvents as well as preserve the energy dissipation law at the semi-discrete and fully discrete levels.

The rest of this paper is organized as follows. In § 2, we introduce the free-energy functionals and derive the phase-field model with active diffusion using the second law of thermodynamics. In § 3, we design a linearized, decoupled, mass-conservative and energy-stable numerical method. A series of numerical experiments is carried out in § 4 to verify and validate the proposed model and numerical method. Finally, some concluding remarks are provided in § 5.

2. Mathematical model

In this section, we introduce the free energies and then derive a thermodynamically consistent model of active diffusion in two-phase immiscible and incompressible fluids. The immiscible and incompressible two-phase fluids are the solvents, and the solute is a multi-component mixture consisting of $N$ chemical substances. We use $\phi$ to denote the phase variable, which stands for the volume fraction of one phase, thereby taking values between 0 and 1. Let $c_l$ be the concentration of solute component $l$, where $1\leq l\leq N$. We also denote $\boldsymbol {c}=[c_1,\ldots,c_N]^{\rm T}$ and $c=\sum _{l=1}^Nc_l$.

2.1. Free energies

In order to describe the active diffusion of the solute mixture in binary fluids, we introduce the free-energy function

(2.1)\begin{equation} A(\boldsymbol{c},\phi)= \sum_{l=1}^N \phi\alpha_l c_l(\ln(c_l)-1-\gamma_l) +\sum_{l=1}^N(1-\phi)\beta_l c_l(\ln(c_l)-1-\delta_l) , \end{equation}

where $\alpha _l$, $\beta _l$, $\gamma _l$ and $\delta _l$ are the energy parameters. The solute chemical potentials are defined as the partial derivatives of the free-energy function:

(2.2)\begin{gather} \mu_{c,l}(\boldsymbol{c},\phi)=\frac{\partial A(\boldsymbol{c},\phi)}{\partial c_l} = \phi\alpha_l (\ln(c_l)-\gamma_l) +(1-\phi)\beta_l (\ln(c_l)-\delta_l) , \end{gather}

(2.3)\begin{gather}\mu_{c\phi}(\boldsymbol{c})=\frac{\partial A(\boldsymbol{c},\phi)}{\partial \phi} = \sum_{l=1}^N\alpha_l c_l(\ln(c_l)-1-\gamma_l) -\sum_{l=1}^N\beta_l c_l(\ln(c_l)-1-\delta_l) . \end{gather}

Physically, $\alpha _l c_l(\ln (c_l)-1-\gamma _l)$ and $\beta _l c_l(\ln (c_l)-1-\delta _l)$ represent the solute component free-energy functions in two different solvent fluids, each of which results from the ideal gas equation of state and is usually used to describe the conventional diffusion process in the respective bulk fluid. The parameters rely generally on the solubility of the component in two fluids under specific physical and chemical conditions, and they usually take different values for different solvent fluids. When the solute concentrations are homogeneous in binary fluids, there may still be contrasts of solute chemical potentials in two fluids, which make the solutes transfer from one solvent fluid to the other. This leads to the so-called active diffusion. At the equilibrium state, the chemical potentials of each solute component in binary fluids are the same, but the solute concentrations in two fluids may be different. This is the underlying mechanics that the mixed free energy is able to characterize the active diffusion of each solute in binary fluids.

Patankar (2016) established the equations for the solute equilibrium and the phase equilibrium. The logarithmic chemical potentials were introduced by Patankar (2016) to formulate the equilibrium equations, and the chemical potentials (2.2) can be viewed as their extensions with considerations of binary solvents and multiple solutes. Typically, the chemical potentials of a substance in two solvents must be equal at equilibrium (Patankar 2016). As a result, the equilibrium condition can be expressed as

(2.4)\begin{equation} \alpha_l (\ln(c_l^1)-\gamma_l) =\beta_l (\ln(c_l^2)-\delta_l),\quad l=1,\ldots,N, \end{equation}

where $c_l^1$ and $c_l^2$ are the equilibrium concentrations of solute $l$ in two bulk solvents. If $\alpha _l\neq \beta _l$ or $\gamma _l\neq \delta _l$, then at the equilibrium state, there will be a solute concentration contrast in two bulk solvents.

We now discuss several special cases of a single-component solute and also omit the subscripts in (2.1)–(2.3). Since the parameters $\alpha$ and $\beta$ usually rely on the temperature, the thermally activated diffusion can be modelled taking different values for $\alpha$ and $\beta$. The parameters $\gamma$ and $\delta$ reflect the equilibrium-state solute concentrations in two solvent fluids, so they can be adjusted to characterize the contrast of ion concentrations in intracellular and extracellular regions. In particular, if the solute has the same diffusion properties in two solvent fluids, then the proposed active diffusion can be reduced to Fick's diffusion. As a matter of fact, taking $\alpha =\beta =1$ and $\gamma =\delta =0$, the chemical potential (2.2) becomes $\mu _c= \ln (c)$, and we deduce from (2.32), which will be derived in § 2.2, that the diffusion flux $\boldsymbol {J}_c$ is proportional to the concentration gradient as

(2.5)\begin{equation} \boldsymbol{J}_c={-}{\mathcal{D}} c\,\boldsymbol{\nabla}\mu_c={-}{\mathcal{D}} c\,\boldsymbol{\nabla}\ln(c) ={-}{\mathcal{D}}\,\boldsymbol{\nabla} c, \end{equation}

which is Fick's diffusion law, and ${\mathcal {D}}$ becomes the Fickian diffusion coefficient. Thus the proposed active diffusion model is an extension of Fick's diffusion in the two immiscible fluid solvents.

The free-energy functional $F(\phi,\boldsymbol {\nabla } \phi )$ for the phase variable is composed of two terms:

(2.6)\begin{equation} F(\phi,\boldsymbol{\nabla} \phi)=\frac{\sigma}{\epsilon}\,f(\phi)+\frac{1}{2}\,\epsilon \sigma\,|\boldsymbol{\nabla} \phi|^2, \end{equation}

where $\sigma$ represents the interfacial tension between binary fluids, and $\epsilon$ is the interface thickness. The bulk free energy for the phase variable can be expressed commonly by the double-well potential

(2.7)\begin{equation} f(\phi)= \phi^2(1- \phi)^2, \end{equation}

or the logarithmic Flory–Huggins energy potential

(2.8)\begin{equation} f(\phi)= \phi\ln(\phi)+(1-\phi)\ln(1-\phi)+\theta(\phi- \phi^2), \end{equation}

where $\theta >2$ is the energy parameter (Wells, Kuhl & Garikipati 2006; Wang, Kou & Cai 2020b). The chemical potential of the phase variable is defined as the variational derivative of the free-energy functional $F$:

(2.9)\begin{equation} \mu_\phi=\frac{\sigma}{\epsilon}\,f'(\phi)-\epsilon\sigma\,\Delta \phi. \end{equation}

2.2. Model derivation

In what follows, we will derive the complete phase-field model using the second law of thermodynamics. Here, for a regular spatial domain $\varOmega$, we use the traditional notations to denote by $(\cdot,\cdot )$ and $\|\cdot \|$ the inner product and the norm of $L^2(\varOmega )$, $(L^2(\varOmega ))^d$ and $(L^2(\varOmega ))^{d\times d}$, where $d$ is the spatial dimension.

We start with the model derivation on the basis of the following mass conservation equations and incompressibility of fluids:

(2.10)\begin{gather} \frac{\partial c_l }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} c_l ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_{c,l}=0,\quad 1\leq l\leq N, \end{gather}

(2.11)\begin{gather}\frac{\partial \phi }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi=0, \end{gather}

(2.12)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} +\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi =0, \end{gather}

where $\boldsymbol {u}$ is the mixture-averaged velocity, and $\boldsymbol {J}_{c,l}$ and $\boldsymbol {J}_\phi$ are the diffusion fluxes that will be determined later. In (2.12), $\lambda =1-\varsigma$, where $\varsigma$ is a dimensionless parameter.

In practice, there are several measures to quantify the flow rate of multiphase multi-component flow. The measured velocity of multiphase multi-component flow is usually some kind of averaged velocity – for instance, mass-averaged velocity, volume-averaged velocity and molar-averaged velocity, depending upon the experimental methodology favoured in different scientific and industrial fields. We note that (2.12) can account for different averaged velocities. Indeed, $\phi$ stands for the volume fraction of one phase, and let $\phi _2$ denote the volume fraction of the other phase; then we have $\phi +\phi _2=1$ due to the immiscibility. Similar to (2.11), we have the mass balance equation

(2.13)\begin{equation} \frac{\partial \phi_2 }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi_2 ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_{\phi_2}=0. \end{equation}

If $\boldsymbol {u}$ is the volume-averaged velocity, then adding (2.11) and (2.13) together should yield

(2.14)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} =0, \end{equation}

which implies that $\boldsymbol {J}_\phi +\boldsymbol {J}_{\phi _2}=0$, namely, $\boldsymbol {J}_{\phi _2}=-\boldsymbol {J}_\phi$; in this case, $\lambda =0$ and $\varsigma =1$. Let $\varrho _1$ and $\varrho _2$ stand for the respective densities of two fluids, and $\rho =\phi \varrho _1+(1-\phi )\varrho _2$. If $\boldsymbol {u}$ is the mass-averaged velocity, then we have

(2.15)\begin{equation} \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u} ) =0, \end{equation}

which implies that $\varrho _1\boldsymbol {J}_\phi +\varrho _2\boldsymbol {J}_{\phi _2}=0$ or $\boldsymbol {J}_{\phi _2}=-({\varrho _1}/{\varrho _2})\boldsymbol {J}_\phi$; in this case, $\lambda =1-{\varrho _1}/{\varrho _2}$, $\varsigma ={\varrho _1}/{\varrho _2}$, and adding (2.11) and (2.13) together yields

(2.16)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} +\left(1-\frac{\varrho_1}{\varrho_2}\right)\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi =0. \end{equation}

From the above analysis, we can assume the relation $\boldsymbol {J}_{\phi _2}=-\varsigma \boldsymbol {J}_\phi$ in general for different averaged velocities. For instance, $\varsigma =1$ yields the volume-averaged velocity, while $\varsigma ={\varrho _1}/{\varrho _2}$ implies the mass-averaged velocity. Moreover, $\varsigma ={n_1}/{n_2}$ implies the mole-averaged velocity, where $n_1$ and $n_2$ stand for mole densities of two fluids. As a result, different choices of $\varsigma$ can represent different types of averaged velocities, and consequently, (2.12) is quite general.

We assume that the solution is dilute and the dissolved substances in binary fluids have no perceptible effect on the volume and density of the mixture. The mixture density is expressed as $\rho =\phi \varrho _1+(1-\phi )\varrho _2$, where $\varrho _1$ and $\varrho _2$ stand for the respective densities of the two fluids. The mass conservation equation of the mixture can be derived from (2.11) and (2.12) as

(2.17)\begin{equation} \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u} + \chi \boldsymbol{J}_\phi)=0, \end{equation}

where $\chi =\varrho _1-\varsigma \varrho _2$.

As pointed out by Brenner (2009a,b, 2013), the general physical principles will not be violated if the mass velocity $\boldsymbol {v}_m$ appearing in the continuity equation differs from the momentum velocity $\boldsymbol {v}_a$ appearing in the momentum balance equation, and the kinetic energy velocity $\boldsymbol {v}_k$ appearing in the energy equation. This finding yields the bi-velocity hydrodynamical models (Brenner 2009a,b, 2013). In this work, we take $\boldsymbol {v}_a=\boldsymbol {v}_k=\boldsymbol {u}$ and define a general mass velocity $\boldsymbol {v}_m=\boldsymbol {u} + \chi \rho ^{-1} \boldsymbol {J}_\phi$. Then (2.17) is rewritten as

(2.18)\begin{equation} \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{v}_m )=0, \end{equation}

which is the general continuity equation. We will show that the proposed model with the above choices of different velocities obeys the second law of thermodynamics.

We define the total solute free energy ($\mathcal {A}$), total phase free energy ($\mathcal {F}$), total kinetic energy ($\mathcal {U}$) and total gravitational potential energy ($\mathcal {G}$) in the domain $\varOmega$ as

(2.19)\begin{align} \mathcal{A} =\int_\varOmega A(\boldsymbol{c},\phi) \,{\rm d}\kern0.06em \boldsymbol{x},\quad \mathcal{F} =\int_\varOmega F(\phi,\boldsymbol{\nabla} \phi) \,{\rm d}\kern0.06em \boldsymbol{x},\quad \mathcal{U} =\frac{1}{2}\int_\varOmega \rho\,|\boldsymbol{u}|^2 \,{\rm d}\kern0.06em \boldsymbol{x},\quad \mathcal{G} =\int_\varOmega \rho gz \,{\rm d}\kern0.06em \boldsymbol{x} , \end{align}

where $g$ is the gravitational acceleration and $z$ is the height above a given reference plane.

The entropy equation for an isothermal fluid system (Kou & Sun 2018a,b) can be expressed as

(2.20)\begin{equation} T\,\frac{\partial S}{\partial t}={-}\frac{\partial(\mathcal{A}+ \mathcal{F}+ \mathcal{U}+ \mathcal{G})}{\partial t}-\int_{\varOmega} \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol\sigma)\,{\rm d}\kern0.06em \boldsymbol{x}, \end{equation}

where $S$ is the total entropy, $T$ is the temperature, and the last term represents the work done by the stress tensor ($\boldsymbol \sigma$).

For the boundary conditions, we assume that all boundary terms will vanish when integrating by parts is performed – for example, the periodic boundary conditions and homogeneous Neumann boundary conditions.

Taking into account (2.10) and (2.11), we derive the variance of total solute free energy with time as

(2.21)\begin{align} \frac{\partial\mathcal{A}}{\partial t} &= \sum_{l=1}^N\left(\mu_{c,l},\frac{\partial c_l }{\partial t}\right)+\left(\mu_{c\phi},\frac{\partial \phi }{\partial t}\right) \nonumber\\ &={-} \sum_{l=1}^N(\mu_{c,l},\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} c_l) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_{c,l})- (\mu_{c\phi},\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi)\nonumber\\ &= \sum_{l=1}^N(\boldsymbol{u}, c_l\,\boldsymbol{\nabla} \mu_{c,l}) +(\boldsymbol{u}, \phi\,\boldsymbol{\nabla} \mu_{c\phi}) + \sum_{l=1}^N (\boldsymbol{J}_{c,l},\boldsymbol{\nabla}\mu_{c,l})+ (\boldsymbol{J}_\phi,\boldsymbol{\nabla}\mu_{c\phi}) . \end{align}

The variance of total phase free energy with time is

(2.22)\begin{align} \frac{\partial\mathcal{F}}{\partial t} &= \left(\frac{\sigma}{\epsilon}\,f'(\phi)-\epsilon\sigma\,\Delta \phi,\frac{\partial \phi }{\partial t}\right) \nonumber\\ &= \left(\mu_\phi,\frac{\partial \phi }{\partial t}\right) \nonumber\\ &={-} (\mu_\phi,\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi)\nonumber\\ &= (\boldsymbol{u} , \phi\,\boldsymbol{\nabla}\mu_\phi) + (\boldsymbol{J}_\phi,\boldsymbol{\nabla}\mu_\phi). \end{align}

The variance equations of kinetic energy and gravitational potential energy with time are derived using (2.17) as

(2.23)\begin{align} \frac{\partial \mathcal{U}}{\partial t}&=\left(\rho\,\frac{\partial \boldsymbol{u}}{\partial t},\boldsymbol{u}\right)+\frac{1}{2}\left(\frac{\partial \rho}{\partial t},|\boldsymbol{u}|^2\right)\nonumber\\ &=\left(\rho\,\frac{\partial \boldsymbol{u}}{\partial t},\boldsymbol{u}\right)-\frac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u})+ \chi\,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_\phi,|\boldsymbol{u}|^2\right)\nonumber\\ &=\left(\rho\,\frac{\partial \boldsymbol{u}}{\partial t}+\rho\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}+ \chi\,\boldsymbol{J}_\phi \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u},\boldsymbol{u}\right) , \end{align}

(2.24)\begin{align} \frac{\partial \mathcal{G}}{\partial t} &=\left(\frac{\partial\rho}{\partial t},gz\right) ={-}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u} )+ \chi\,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_\phi,gz\right)\nonumber\\ &=(\rho\boldsymbol{u}+ \chi\,\boldsymbol{J}_\phi, g\,\boldsymbol{\nabla} z ) \nonumber\\ &={-}(\boldsymbol{u}, \rho\boldsymbol{\mathsf{g}})-( \boldsymbol{J}_\phi, \chi\boldsymbol{\mathsf{g}}) , \end{align}

where $\boldsymbol{\mathsf{g}} =-\boldsymbol {\nabla } z$. We split the total stress $\boldsymbol \sigma$ into two parts:

(2.25)\begin{equation} \boldsymbol\sigma= p\boldsymbol{I}+\boldsymbol\sigma_{irrev}, \end{equation}

where $p$ is the reversible part, i.e. the pressure, $\boldsymbol {I}$ is the identity tensor, and $\boldsymbol \sigma _{irrev}$ is the irreversible part. Thanks to (2.12), we derive

(2.26)\begin{align} (\boldsymbol\sigma,\boldsymbol{\nabla}\boldsymbol{u}) &=(p,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})+({\boldsymbol\sigma_{irrev}},\boldsymbol{\nabla}\boldsymbol{u})\nonumber\\ &={-}(p,\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi)+({\boldsymbol\sigma_{irrev}},\boldsymbol{\nabla}\boldsymbol{u}) \nonumber\\ &=(\lambda\,\boldsymbol{\nabla} p,\boldsymbol{J}_\phi)+({\boldsymbol\sigma_{irrev}},\boldsymbol{\nabla}\boldsymbol{u}) . \end{align}

Substituting (2.21)–(2.24) and (2.26) into (2.20), we can deduce that

(2.27)\begin{align} T\,\frac{\partial S}{\partial t}&={-}(\boldsymbol{J}_\phi,\boldsymbol{\nabla}\mu_\phi+\boldsymbol{\nabla}\mu_{c\phi}+\lambda\,\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}} ) - \sum_{l=1}^N(\boldsymbol{J}_{c,l},\boldsymbol{\nabla}\mu_{c,l})- ( \boldsymbol\sigma_{irrev},\boldsymbol{\nabla}\boldsymbol{u})\nonumber\\ &\quad -\left(\rho\,\frac{\partial \boldsymbol{u}}{\partial t}+(\rho\boldsymbol{u}+ \chi\boldsymbol{J}_\phi) \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol\sigma+\phi\,\boldsymbol{\nabla}(\mu_\phi+\mu_{c\phi})\right.\nonumber\\ &\quad \left.+\sum_{l=1}^N c_l\,\boldsymbol{\nabla} \mu_{c,l}-\rho\boldsymbol{\mathsf{g}},\boldsymbol{u}\right) . \end{align}

Galilean invariance yields from (2.27) that

(2.28)\begin{equation} \rho\,\frac{\partial \boldsymbol{u}}{\partial t}+(\rho\boldsymbol{u}+ \chi\boldsymbol{J}_\phi) \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol\sigma+\phi\, \boldsymbol{\nabla}(\mu_\phi+\mu_{c\phi})+ \sum_{l=1}^N c_l\,\boldsymbol{\nabla} \mu_{c,l}-\rho\boldsymbol{\mathsf{g}} =0. \end{equation}

For the irreversible system, in terms of Newtonian fluid theory, we have

(2.29)\begin{equation} \boldsymbol\sigma_{irrev}={-} \eta(\phi)\,(\boldsymbol{\nabla}\boldsymbol{u}+ \boldsymbol{\nabla}\boldsymbol{u}^{\rm T}) , \end{equation}

where $\eta$ is the fluid viscosity depending on $\phi$. The diffusion flux $\boldsymbol {J}_\phi$ should take the form

(2.30)\begin{equation} \boldsymbol{J}_\phi={-}M_\phi(\phi)\left(\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})+\lambda\,\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}}\right), \end{equation}

where $M_\phi (\phi )$ is the coefficient dependent on $\phi$. For the solute mixture, there exists an $N\times N$ phenomenological coefficient matrix

(2.31)\begin{equation} \boldsymbol{\mathcal{D}}=({\mathcal{D}}_{l,m})_{l,m=1}^N, \end{equation}

such that

(2.32)\begin{equation} \boldsymbol{J}_{c,l}={-}\sum_{m=1}^N{\mathcal{D}}_{l,m}\,\boldsymbol{\nabla}\mu_{c,m},\quad l=1,\ldots, N . \end{equation}

According to the famous Onsager's reciprocal principle, $\boldsymbol {\mathcal {D}}$ should be symmetric. Substituting (2.28), (2.29), (2.30) and (2.32) into (2.27), we deduce that

(2.33)\begin{align} T\,\frac{\partial S}{\partial t}&= \left\|M_\phi(\phi)^{1/2}\left(\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})+\lambda\,\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}}\right)\right\|^2 \nonumber\\ &\quad +\sum_{l=1}^N\sum_{m=1}^N({\mathcal{D}}_{l,m}\,\boldsymbol{\nabla}\mu_{c,m},\boldsymbol{\nabla}\mu_{c,l}) +\frac{1}{2}\left\|\eta(\phi)^{1/2}\,(\boldsymbol{\nabla}\boldsymbol{u}+ \boldsymbol{\nabla}\boldsymbol{u}^{\rm T})\right\|^2 . \end{align}

To obey the second law of thermodynamics, i.e. $T({\partial S}/{\partial t})\geq 0$, the matrix $\boldsymbol {\mathcal {D}}$ should be symmetric positive semi-definite.

2.3. Model equations

The proposed model can be summarized as follows:

(2.34a)\begin{gather} \frac{\partial c_l }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} c_l) +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_{c,l}=0,\quad l=1,\ldots, N, \end{gather}

(2.34b)\begin{gather}\boldsymbol{J}_{c,l}={-}\sum_{m=1}^N{\mathcal{D}}_{l,m}\,\boldsymbol{\nabla}\mu_{c,m},\quad l=1,\ldots, N, \end{gather}

(2.34c)\begin{gather}\frac{\partial \phi }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi ) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi=0, \end{gather}

(2.34d)\begin{gather}\boldsymbol{J}_\phi={-}M_\phi(\phi)\left(\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})+\lambda\,\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}}\right), \end{gather}

(2.34e)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} +\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi =0, \end{gather}

(2.34f)\begin{gather} \frac{\partial (\rho\boldsymbol{u})}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left((\rho\boldsymbol{u}+ \chi\boldsymbol{J}_\phi\right)\otimes \boldsymbol{u}) ={-}\boldsymbol{\nabla} p-\phi\,\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})\nonumber\\ {}- \sum_{l=1}^N c_l\,\boldsymbol{\nabla} \mu_{c,l}+\boldsymbol{\nabla}\boldsymbol{\cdot}\eta(\phi)\,(\boldsymbol{\nabla}\boldsymbol{u}+ \boldsymbol{\nabla}\boldsymbol{u}^{\rm T})+\rho\boldsymbol{\mathsf{g}} . \end{gather}

In (2.34), $\lambda =1-\varsigma$, $\chi =\varrho _1-\varsigma \varrho _2$, and $\varsigma$ is a given dimensionless parameter that represents the type of averaged velocity. The chemical potentials $\mu _c$, $\mu _{c\phi }$ and $\mu _\phi$ are given in (2.2), (2.3) and (2.9), respectively. Here, we take $\eta (\phi )=\phi \eta _1+(1-\phi )\eta _2$, where $\eta _1$ and $\eta _2$ are the viscosities of the two solvent fluids. As (2.34) shows, if $c_l=0$ in both phases, then the system (2.34) reduces to the regular phase-field model for which ample verification exercises can be found in our previous work (Kou et al. 2020b). Equation (2.34f) is the conservative form of the momentum balance equation, which is obtained by combining (2.28) and (2.17).

A simple form of $\boldsymbol {\mathcal {D}}$ is taken as a diagonal matrix with the principal diagonal elements ${\mathcal {D}}_{l,l}=D_{l}c_l$, where $D_{l}>0$. A general choice of $\boldsymbol {\mathcal {D}}$ can be obtained using the Maxwell–Stefan approach, which takes into consideration the crossing influences between different solutes.

The solute transport equations are generally expressed as

(2.35)\begin{equation} \frac{\partial c_l }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}_l c_l) =0,\quad l=1,\ldots, N, \end{equation}

where $\boldsymbol {u}_l$ is the velocity of solute $l$. Let us define $\boldsymbol {J}_{c,l}=c_l(\boldsymbol {u}_l-\boldsymbol {u})$; then we get (2.34a). The modified Maxwell–Stefan equations can be formulated as

(2.36)\begin{equation} \frac{c_l(\boldsymbol{u}_l-\boldsymbol{u})}{c D_{l,f}}+\sum_{m=1, m\neq l}^N\frac{c_lc_m(\boldsymbol{u}_l-\boldsymbol{u}_m)}{c^2 D_{l,m}}={-}c_l\,\boldsymbol{\nabla} \mu_{c,l},\quad l=1,\ldots, N, \end{equation}

where $c=\sum _{l=1}^Nc_l$, $D_{l,f}>0$ and $D_{l,m}=D_{m,l}>0$. In (2.36), the first term represents the friction between solute $l$ and the solvent fluids, which is usually ignored in the classical Maxwell–Stefan model (Wesselingh & Krishna 2000), the second term stands for the friction between solute $l$ and the other solutes, and the right-hand side is the driving force of solute $l$. Equation (2.36) describes the balance between the drag and driving forces for each solute.

We assume that $c_l>0$, $1\leq l\leq N$, and define a matrix $\boldsymbol {\mathcal {L}}=(\mathcal {L}_{l,m})_{l,m=1}^N$ as

(2.37)\begin{equation} \mathcal{L}_{l,l}= \frac{c_l}{c D_{l,f}} + \sum_{m=1, m\neq l}^N\frac{c_lc_m}{c^2 D_{l,m}},\quad \mathcal{L}_{l,m}={-}\frac{c_lc_m}{c^2 D_{l,m}},\quad m\neq l. \end{equation}

Due to $D_{l,m}=D_{m,l}$ and $c_l>0$, $\boldsymbol {\mathcal {L}}$ is symmetric positive definite, and so is its inverse matrix $\boldsymbol {\mathcal {L}}^{-1}$. Consequently, $\boldsymbol {u}_l-\boldsymbol {u}$, $l=1,\ldots, N$, can be determined uniquely by (2.36). Taking into account $\boldsymbol {J}_{c,l}=c_l(\boldsymbol {u}_l-\boldsymbol {u})$, we derive a full diffusion coefficient matrix from (2.36) as

(2.38)\begin{equation} \boldsymbol{\mathcal{D}}= \textrm{diag}(\boldsymbol{c})\,\boldsymbol{\mathcal{L}}^{{-}1}\,\textrm{diag}(\boldsymbol{c}), \end{equation}

which is symmetric positive definite and consistent with Onsager's reciprocal principle and the second law of thermodynamics.

If the interfacial effects on solutes are considered, then the diffusion coefficients $D_{l,f}$ and $D_{l,m}$ in (2.36) should be expressed as phase-dependent functions; for instance, $D_{l,f}$ (Qin et al. 2022) can be defined generally as

(2.39)\begin{equation} \frac{1}{D_{l,f}(\phi)}=\frac{\phi^2(1-\phi)^2}{\epsilon K} +\frac{1-\phi}{Q_{l,0}}+\frac{\phi}{Q_{l,1}}, \end{equation}

where $K$ is the interface permeability, and $Q_{l,0}$ and $Q_{l,1}$ are the diffusion coefficients of solute $l$ in two bulk solvents.

The derivations of the model imply the following energy dissipation law:

(2.40)\begin{align} \frac{\partial(\mathcal{A}+ \mathcal{F}+ \mathcal{U}+ \mathcal{G})}{\partial t}&={-}\sum_{l=1}^N\sum_{m=1}^N({\mathcal{D}}_{l,m}\,\boldsymbol{\nabla}\mu_{c,m},\boldsymbol{\nabla}\mu_{c,l})\nonumber\\ &\quad -\left\|M_\phi(\phi)^{1/2}\left(\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})+\lambda\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}}\right)\right\|^2\nonumber\\ &\quad -\frac{1}{2}\left\|\eta(\phi)^{1/2}\,(\boldsymbol{\nabla}\boldsymbol{u}+ \boldsymbol{\nabla}\boldsymbol{u}^{\rm T})\right\|^2 \nonumber\\ &\leq0, \end{align}

which shows that total free energy of a closed isothermal system will be decreasing with time.

We denote by $L_*$, $U_*$ and $M_*$ the characteristic length, the characteristic velocity and the characteristic diffusion constant, and denote by $c_*$ the reference concentration. Furthermore, we introduce the following dimensionless quantities:

(2.41a–g)\begin{align} \hat{\boldsymbol{x}}=\frac{\boldsymbol{x}}{L_*},\quad \hat{t} =\frac{tU_*}{L_*},\quad \hat{\boldsymbol{u}} =\frac{\boldsymbol{u}}{U_*},\quad \hat{\rho}=\frac{\rho}{\varrho_1},\quad \hat{\eta}=\frac{\eta}{\eta_1}, \quad \hat{p}=\frac{p}{\varrho_1U_*^2},\quad \hat{c}_l=\frac{c_l}{c_{*}}. \end{align}

Omitting the hats, we can write the dimensionless system of the model (2.34) as

(2.42a)\begin{gather} \frac{\partial c_l }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} c_l) +\frac{1}{Pe_c}\,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_{c,l}=0,\quad l=1,\ldots, N, \end{gather}

(2.42b)\begin{gather}\boldsymbol{J}_{c,l}={-}\sum_{m=1}^N{\mathcal{D}}_{l,m}\,\boldsymbol{\nabla}\mu_{c,m}, \quad l=1,\ldots, N, \end{gather}

(2.42c)\begin{gather}\frac{\partial \phi }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \phi) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi=0, \end{gather}

(2.42d)\begin{gather}\boldsymbol{J}_\phi={-}\frac{M_\phi(\phi)}{Pe_\phi}\left(\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})+\lambda\,\boldsymbol{\nabla} p- \chi \boldsymbol{\mathsf{g}}\right), \end{gather}

(2.42e)\begin{gather}\mu_\phi=\frac{1}{We}\left(\frac{1}{Cn}\,f(\phi)-Cn\,\Delta \phi\right), \end{gather}

(2.42f)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} +\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi =0, \end{gather}

(2.42g)\begin{gather}\frac{\partial (\rho\boldsymbol{u})}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left((\rho\boldsymbol{u}+ \chi\boldsymbol{J}_\phi)\otimes \boldsymbol{u}\right) =\frac{1}{Re}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\eta(\phi)\,(\boldsymbol{\nabla}\boldsymbol{u}+ \boldsymbol{\nabla}\boldsymbol{u}^{\rm T})\nonumber\\{}-\boldsymbol{\nabla} p-\phi\,\boldsymbol{\nabla}(\mu_\phi +\mu_{c\phi})-\sum_{l=1}^Nc_l\,\boldsymbol{\nabla} \mu_{c,l}+\frac{1}{Fr^2}\,\rho \boldsymbol{\mathsf{g}}, \end{gather}

where $\lambda =1-\varsigma$, $\chi =1-\varsigma ({\varrho _2}/{\varrho _1})$, $\boldsymbol{\mathsf{g}}$ is the unit vector with the gravitational direction, $\rho =\phi +(1-\phi )({\varrho _2}/{\varrho _1})$ and $\eta (\phi )=\phi +(1-\phi )({\eta _2}/{\eta _1})$. Here, $\varsigma$ remains a dimensionless parameter relating to $\boldsymbol {u}$, $c_l$ is the dimensionless concentration with the range $(0,1)$, $Pe_c=Pe_\phi = {L_*U_*}/{M_*}$ is the diffusional Péclet number, $We = {\varrho _1 U_*^2L_*}/{\sigma }$ is the Weber number, $Cn = {\epsilon }/{L_*}$ is the Cahn number as a measure of the interfacial thickness, $Re = {\varrho _1 U_*L_*}/{\eta _1}$ is the Reynolds number, and $Fr = {U_*}/{\sqrt {gL_*}}$ is the Froude number that describes the relative strength between the gravitational and inertial forces.

3. Numerical method

In this section, we propose an efficient, linear, decoupled, energy-stable and mass-conservative numerical method for the proposed model.

The total simulation time $[0, t_{f}]$ is divided as $0=t_0< t_1< t_2<\cdots < t_f$, and the time step size is denoted by $\tau _k=t_{k+1}-t_k$. The proposed scheme allows us to use non-uniform and adaptive time step sizes, and as a result, $\tau _k$ is not constant in general. The value of a function $v$ at time $t_k$ is denoted by $v^k$.

3.1. Discrete chemical potentials

An important ingredient of designing energy-stable numerical methods is to derive the discrete chemical potentials that ensure certain energy dissipation inequalities. The energy factorization approach (Kou, Sun & Wang 2020a; Wang et al. 2020b; Wang, Kou & Gao 2021) can produce efficient, linear and energy-stable numerical schemes.

We now derive the discrete chemical potentials for the active solute diffusion. The concavity of $\ln (c_l)$ implies that

(3.1)\begin{equation} \ln(c_l^{k+1})-\ln(c_l^k) \leq \frac{1}{c_l^k}\,(c_l^{k+1}- c_l^k). \end{equation}

Using (3.1), we deduce that

(3.2)\begin{align} c_l^{k+1}\ln(c_l^{k+1})- c_l^{k}\ln(c_l^{k})&= \ln(c_l^{k})\,(c_l^{k+1}- c_l^{k} )+c_l^{k+1}(\ln(c_l^{k+1})-\ln(c_l^k))\nonumber\\ &\leq\ln(c_l^{k})\,(c_l^{k+1}- c_l^{k})+\frac{c_l^{k+1}}{c_l^k}\,(c_l^{k+1}- c_l^k). \end{align}

Using (3.2), we derive the energy difference as

(3.3) \begin{align} A(\boldsymbol{c}^{k+1},\phi^{k+1})- A(\boldsymbol{c}^k,\phi^k) &=\sum_{l=1}^N \alpha_l \phi^k (c_l^{k+1}\ln(c_l^{k+1})- c_l^{k}\ln(c_l^{k}) )\nonumber\\ &\quad -\sum_{l=1}^N \alpha_l \phi^k(1+\gamma_l) ( c_l^{k+1}- c_l^k)\nonumber\\ &\quad +\sum_{l=1}^N\alpha_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\gamma_l\right) (\phi^{k+1}- \phi^k )\nonumber\\ &\quad +\sum_{l=1}^N\beta_l(1-\phi^k) \left(c_l^{k+1}\ln(c_l^{k+1})- c_l^{k}\ln(c_l^{k}) \right) \nonumber\\ &\quad -\sum_{l=1}^N\beta_l(1-\phi^k) (1+\delta_l)(c_l^{k+1}- c_l^{k})\nonumber\\ &\quad -\sum_{l=1}^N\beta_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\delta_l) (\phi^{k+1}- \phi^k \right)\nonumber\\ &\leq\sum_{l=1}^N\alpha_l \phi^k \left(\ln(c_l^{k})+\frac{c_l^{k+1}}{c_l^k}-1-\gamma_l\right)( c_l^{k+1}- c_l^k) \nonumber\\ &\quad +\sum_{l=1}^N\beta_l(1-\phi^k)\left(\ln(c_l^{k})+\frac{c_l^{k+1}}{c_l^k}-1-\delta_l\right)( c_l^{k+1}- c_l^k)\nonumber\\ &\quad +\sum_{l=1}^N\alpha_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\gamma_l\right) (\phi^{k+1}- \phi^k )\nonumber\\ &\quad -\sum_{l=1}^N\beta_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\delta_l\right) (\phi^{k+1}- \phi^k ). \end{align}

Based on (3.3), we define the discrete chemical potentials as

(3.4)\begin{align} \mu_{c,l}^{k+1}&= \alpha_l \phi^k \left(\ln(c_l^{k})+\frac{c_l^{k+1}}{c_l^k}-1-\gamma_l\right)\nonumber\\ &\quad +\beta_l(1-\phi^k)\left(\ln(c_l^{k})+\frac{c_l^{k+1}}{c_l^k}-1-\delta_l\right),\quad 1\leq l\leq N, \end{align}

(3.5)\begin{align} \mu_{c\phi}^{k+1}&= \sum_{l=1}^N\alpha_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\gamma_l\right)-\sum_{l=1}^N\beta_l c_l^{k+1}\left(\ln(c_l^{k+1})-1-\delta_l\right), \end{align}

which satisfy the inequality

(3.6)\begin{equation} A(\boldsymbol{c}^{k+1},\phi^{k+1})- A(\boldsymbol{c}^k,\phi^k) \leq \sum_{l=1}^N\mu_{c,l}^{k+1} (c_l^{k+1}-c_l^{k})+\mu_{c\phi}^{k+1} (\phi^{k+1}-\phi^{k}). \end{equation}

For the bulk phase free-energy function, we assume that the time discrete scheme is linear and satisfies the energy inequality

(3.7)\begin{equation} f(\phi^{k+1})-f(\phi^{k})\leq \ell(\phi^{k},\phi^{k+1})\,(\phi^{k+1}-\phi^{k}), \end{equation}

where $\ell (\phi ^{k},\phi ^{k+1})$ represents the discrete chemical potential function that is linear with respect to $\phi ^{k+1}$. For instance, the logarithmic Flory–Huggins energy potential has such a discrete chemical potential (Wang et al. 2020b) as

(3.8)\begin{equation} \ell(\phi^k,\phi^{k+1}) =\ln(\phi^{k})-\ln(1-\phi^{k})+ \frac{\phi^{k+1}}{\phi^{k}}-\frac{1-\phi^{k+1}}{1-\phi^{k}} +\theta(1-\phi^{k+1}-\phi^{k}). \end{equation}

A fine discrete chemical potential for the double-well potential can be found in Wang et al. (2021). For the gradient free energy, we have

(3.9)\begin{align} \tfrac{1}{2}\left(\|\boldsymbol{\nabla} \phi^{k+1}\|^2-\|\boldsymbol{\nabla} \phi^{k}\|^2\right) &=\left( \boldsymbol{\nabla} \phi^{k+1},\boldsymbol{\nabla}(\phi^{k+1}-\phi^{k})\right) -\tfrac{1}{2}\|\boldsymbol{\nabla}(\phi^{k+1}-\phi^k)\|^2\nonumber\\ &\leq\left( \boldsymbol{\nabla} \phi^{k+1},\boldsymbol{\nabla}(\phi^{k+1}-\phi^{k})\right)\nonumber\\ &={-} (\Delta \phi^{k+1},\phi^{k+1}-\phi^{k} ), \end{align}

which suggests the implicit discrete chemical potential $- \Delta \phi ^{k+1}$. Therefore, we define the discrete form of chemical potential $\mu _\phi$ as

(3.10)\begin{equation} \mu_\phi^{k+1}=\frac{\sigma}{\epsilon}\left(\ell(\phi^{k},\phi^{k+1})-\epsilon^2\,\Delta \phi^{k+1}\right). \end{equation}

3.2. Semi-implicit time marching scheme

On the basis of the chemical potentials given in (3.4), (3.5) and (3.10), we propose a semi-implicit scheme for the model (2.34) as

(3.11a)\begin{gather} \frac{ c_l^{k+1}- c_l^{k} }{\tau_k} + \boldsymbol{\nabla}\boldsymbol{\cdot}(c_l^{k}\boldsymbol{u}_{{\dagger}}^k) +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_{c,l}^{k+1}=0,\quad l=1,\ldots, N, \end{gather}

(3.11b)\begin{gather}\boldsymbol{J}_{c,l}^{k+1}={-}\sum_{m=1}^N{\mathcal{D}}_{l,m}^{k}\,\boldsymbol{\nabla}\mu_{c,m}^{k+1},\quad l=1,\ldots, N, \end{gather}

(3.11c)\begin{gather}\boldsymbol{u}_{{\dagger}}^k=\boldsymbol{u}^k-\sum_{m=1}^N\frac{\tau_k }{\rho^{k}}\,c_m^k\,\boldsymbol{\nabla} \mu_{c,m}^{k+1} , \end{gather}

(3.11d)\begin{gather}\frac{ \phi^{k+1}- \phi^{k} }{\tau_k} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\phi^{k}\boldsymbol{u}_\star^k) +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}=0, \end{gather}

(3.11e)\begin{gather}\boldsymbol{u}_\star^k=\boldsymbol{u}_{\dagger}^k-\frac{\tau_k}{\rho^k}\left( \phi^{k}\,\boldsymbol{\nabla}(\mu_\phi^{k+1} + \mu_{c\phi}^{k+1})+\boldsymbol{\nabla} p^{k+1}-\rho^k \boldsymbol{\mathsf{g}}\right), \end{gather}

(3.11f)\begin{gather}\boldsymbol{J}_\phi^{k+1}={-}M_\phi(\phi^{k})\left(\boldsymbol{\nabla}(\mu_\phi^{k+1} + \mu_{c\phi}^{k+1})+\lambda\,\boldsymbol{\nabla} p^{k+1}- \chi \boldsymbol{\mathsf{g}}\right), \end{gather}

(3.11g)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_\star^k +\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1} =0 , \end{gather}

(3.11h)\begin{gather} \frac{\rho^{k+1}\boldsymbol{u}^{k+1}-\rho^k\boldsymbol{u}_\star^k}{\tau_k}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left((\rho^k\boldsymbol{u}_\star^k+\chi\boldsymbol{J}_\phi^{k+1})\otimes \boldsymbol{u}^{k+1}\right) \nonumber\\{} -\boldsymbol{\nabla}\boldsymbol{\cdot}\eta(\phi^{k})\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ {\boldsymbol{\nabla}\boldsymbol{u}^{k+1}}^{\rm T})=0, \end{gather}

where $\rho ^{k}=\phi ^{k}\varrho _1+(1-\phi ^{k})\varrho _2$.

The resulting system can be solved sequentially using the decoupling relationships between unknowns. First, since the concentrations $\boldsymbol {c}^{k+1}$ are decoupled from $\phi ^{k+1}$, $p^{k+1}$ and $\boldsymbol {u}^{k+1}$, we solve all solute concentrations $\boldsymbol {c}^{k+1}$ together by substituting (3.11b) and (3.11c) into (3.11a). Second, taking into account the property that $\phi ^{k+1}$ and $p^{k+1}$ are decoupled from $\boldsymbol {u}^{k+1}$, we can substitute (3.11e) and (3.11f) into (3.11d) and (3.11g) to compute $\phi ^{k+1}$ and $p^{k+1}$ simultaneously. Finally, $\boldsymbol {u}^{k+1}$ is solved by (3.11h). The purpose of coupling the phase variable and pressure is to ensure the mass conservation of the solvent fluids. In addition to the decoupling relationships between different unknowns, the scheme (3.11) is totally linearized, so it is easy to implement in practice.

We consider the following homogeneous boundary conditions on the boundary $\partial \varOmega$:

(3.12)\begin{equation} \left.\begin{gathered} \boldsymbol{u} =0,\quad \boldsymbol{\nabla}(\mu_\phi + \mu_{c\phi} +\chi gh)\boldsymbol{\cdot}\boldsymbol{n} =0,\\ \boldsymbol{\nabla}(p+\varrho_2 gh)\boldsymbol{\cdot}\boldsymbol{n} =0,\quad \boldsymbol{\nabla}\mu_{c,l}\boldsymbol{\cdot}\boldsymbol{n} =0,\quad \boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n} =0, \end{gathered}\right\} \end{equation}

where $\boldsymbol {n}$ denotes an outward unit vector normal to $\partial \varOmega$.

The proposed scheme has the ability to guarantee the mass conservation law for both solutes and solvent fluids without any additional treatments.

Theorem 3.1 The scheme (3.11) conserves the total masses of both solutes and solvent fluids as

(3.13)\begin{gather} \int_\varOmega c_l^{k+1}\,{\rm d}\kern0.06em \boldsymbol{x}=\int_\varOmega c_l^{k}\,{\rm d}\kern0.06em \boldsymbol{x}=\cdots= \int_\varOmega c_l^{0}\,{\rm d}\kern0.06em \boldsymbol{x}, \quad k\geq0,\ l=1,\ldots, N, \end{gather}

(3.14)\begin{gather}\int_\varOmega \rho^{k+1}\,{\rm d}\kern0.06em \boldsymbol{x}=\int_\varOmega \rho^{k}\,{\rm d}\kern0.06em \boldsymbol{x}=\cdots= \int_\varOmega \rho^{0}\,{\rm d}\kern0.06em \boldsymbol{x},\quad k\geq0. \end{gather}

Proof. Integrating (3.11a) over $\varOmega$ and taking into account the boundary conditions, we obtain

(3.15)\begin{equation} \int_\varOmega\frac{ c_l^{k+1}- c_l^{k} }{\tau_k} \,{\rm d}\kern0.06em \boldsymbol{x} =0,\quad k\geq0,\ l=1,\ldots, N, \end{equation}

which yields (3.13) by induction. Let $\psi ^k=1-\phi ^{k}$. Combining (3.11d) and (3.11g) yields

(3.16)\begin{equation} \frac{ \psi^{k+1}-\psi^k }{\tau_k} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\psi^k\boldsymbol{u}_\star^k) - \varsigma\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}=0. \end{equation}

Multiplying (3.11d) by $\varrho _1$, and (3.16) by $\varrho _2$, and then summing them, we obtain the total mass balance equation

(3.17)\begin{equation} \frac{\rho^{k+1}-\rho^k }{\tau_k} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho^k\boldsymbol{u}_\star^k) +\chi\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}=0, \end{equation}

where $\chi =\varrho _1-\varsigma \varrho _2$. Integrating (3.17) over $\varOmega$ and taking into account the boundary conditions, we get (3.14) by induction.

We now prove that the proposed scheme follows a discrete energy dissipation law. The total discrete energy is defined as

(3.18)\begin{equation} \mathcal{E}^k = \mathcal{A}^k+\mathcal{F}^k+\mathcal{U}^k+\mathcal{G}^k, \end{equation}

where

(3.19)\begin{equation} \left.\begin{gathered} \mathcal{A}^k =\int_\varOmega A(\boldsymbol{c}^k,\phi^k) \,{\rm d}\kern0.06em \boldsymbol{x},\quad \mathcal{F}^k =\int_\varOmega F(\phi^{k},\boldsymbol{\nabla} \phi^{k}) \, {\rm d}\kern0.06em \boldsymbol{x},\\ \mathcal{U}^k =\frac{1}{2}\int_\varOmega \rho^k\,|\boldsymbol{u}^k|^2\, {\rm d}\kern0.06em \boldsymbol{x},\quad \mathcal{G}^k =\int_\varOmega \rho^k gz \,{\rm d}\kern0.06em \boldsymbol{x}. \end{gathered}\right\} \end{equation}

Theorem 3.2 Given that (3.7) holds, the scheme (3.11) satisfies the following discrete energy dissipation law:

(3.20)\begin{align} \frac{\mathcal{E}^{k+1}-\mathcal{E}^{k} }{\tau_k} &\leq{-} \sum_{l=1}^N\sum_{m=1}^N({\mathcal{D}}_{l,m}^k\,\boldsymbol{\nabla}\mu_{c,m}^{k+1},\boldsymbol{\nabla}\mu_{c,l}^{k+1})\nonumber\\ &\quad -\left\|M_\phi(\phi^{k})^{1/2}\left(\boldsymbol{\nabla}(\mu_\phi^{k+1} + \mu_{c\phi}^{k+1})+\lambda\,\boldsymbol{\nabla} p^{k+1}- \chi \boldsymbol{\mathsf{g}}\right)\right\|^2\nonumber\\ &\quad -\frac{1}{2}\left\|\eta(\phi^{k})^{1/2}\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ \boldsymbol{\nabla}{\boldsymbol{u}^{k+1}}^{\rm T})\right\|^2. \end{align}

Proof. Using (3.6), (3.11a) and (3.11d), we deduce that

(3.21)\begin{align} \frac{\mathcal{A}^{k+1}-\mathcal{A}^{k}}{\tau_k} &\leq\sum_{l=1}^N\left(\mu_{c,l}^{k+1},\frac{ c_l^{k+1}- c_l^{k} }{\tau_k}\right)+\left(\mu_{c\phi}^{k+1},\frac{ \phi^{k+1}- \phi^{k} }{\tau_k}\right)\nonumber\\ &\leq{-} \sum_{l=1}^N\left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}_{\dagger}^{k} c_l^{k} )+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_{c,l},\mu_{c,l}^{k+1}\right) \nonumber\\ &\quad - \left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}_*^{k} \phi^{k}),\mu_{c\phi}^{k+1}\right)-(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_\phi^{k+1}, \mu_{c\phi}^{k+1}) \nonumber\\ &\leq \sum_{l=1}^N ( \boldsymbol{u}_{\dagger}^{k} , c_l^{k}\,\boldsymbol{\nabla}\mu_{c,l}^{k+1} )+\sum_{l=1}^N(\boldsymbol{J}_{c,l}, \boldsymbol{\nabla}\mu_{c,l}^{k+1})\nonumber\\ &\quad +(\boldsymbol{u}_*^{k} , \phi^{k}\,\boldsymbol{\nabla}\mu_{c\phi}^{k+1}) +( \boldsymbol{J}_\phi^{k+1}, \boldsymbol{\nabla}\mu_{c\phi}^{k+1}) . \end{align}

It is deduced from (3.7) and (3.9) that

(3.22)\begin{equation} \mathcal{F}^{k+1}-\mathcal{F}^{k} \leq (\mu_\phi^{k+1}, \phi^{k+1}- \phi^{k}) . \end{equation}

Substituting (3.11d) into (3.22), we obtain

(3.23)\begin{align} \frac{\mathcal{F}^{k+1}-\mathcal{F}^{k}}{\tau_k} &\leq{-} \left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}_\star^k \phi^{k})+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}, \mu_\phi^{k+1}\right) \nonumber\\ &\leq(\boldsymbol{u}_\star^k , \phi^{k}\,\boldsymbol{\nabla}\mu_\phi^{k+1}) +(\boldsymbol{J}_\phi^{k+1}, \boldsymbol{\nabla}\mu_\phi^{k+1}) . \end{align}

Using (3.17), we derive the variation of gravitational potential energy as

(3.24)\begin{align} \frac{\mathcal{G}^{k+1}-\mathcal{G}^{k}}{\tau_k} &=\frac{1}{\tau_k}\,( \rho^{k+1}-\rho^k,gz)\nonumber\\ &={-}\left( \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho^k\boldsymbol{u}_\star^k ) +\chi\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}, gz\right) \nonumber\\ &=(\boldsymbol{u}_\star^k \rho^k, g\,\boldsymbol{\nabla} z)+(\boldsymbol{J}_\phi^{k+1}, \chi g\,\boldsymbol{\nabla} z)\nonumber\\ &={-}(\boldsymbol{u}_\star^k , \rho^k\boldsymbol{\mathsf{g}})-(\boldsymbol{J}_\phi^{k+1}, \chi \boldsymbol{\mathsf{g}}). \end{align}

Summing (3.21), (3.23) and (3.24) together and taking into account (3.11f), we obtain

(3.25)\begin{align} &\frac{\mathcal{A}^{k+1}+\mathcal{F}^{k+1}+\mathcal{G}^{k+1}-\mathcal{A}^{k}-\mathcal{F}^{k}-\mathcal{G}^{k}}{\tau_k}\nonumber\\ &\quad \leq \sum_{l=1}^N( \boldsymbol{u}_{\dagger}^{k} , c_l^{k}\,\boldsymbol{\nabla}\mu_{c,l}^{k+1} ) + \left(\boldsymbol{u}_\star^k , \phi^{k}\,\boldsymbol{\nabla}(\mu_\phi^{k+1}+\mu_{c\phi}^{k+1})-\rho^k \boldsymbol{\mathsf{g}}\right)\nonumber\\ &\qquad +\left(\boldsymbol{J}_\phi^{k+1}, \boldsymbol{\nabla}(\mu_\phi^{k+1} + \mu_{c\phi}^{k+1})- \chi \boldsymbol{\mathsf{g}}\right) \nonumber\\ &\qquad - \sum_{l=1}^N\sum_{m=1}^N({\mathcal{D}}_{l,m}^k\,\boldsymbol{\nabla}\mu_{c,m}^{k+1},\boldsymbol{\nabla}\mu_{c,l}^{k+1}) . \end{align}

In order to estimate the kinetic energy, we introduce the following intermediate kinetic energies:

(3.26)\begin{equation} \mathcal{U}_{\dagger}^{k}=\tfrac{1}{2}(\rho^k,|\boldsymbol{u}_{\dagger}^k|^2),\quad \mathcal{U}_\star^{k}=\tfrac{1}{2}(\rho^k,|\boldsymbol{u}_\star^k|^2). \end{equation}

The difference between $\mathcal {U}_{\dagger} ^{k}$ and $\mathcal {U}^{k}$ is estimated as

(3.27)\begin{align} \frac{\mathcal{U}_{\dagger}^{k}-\mathcal{U}^{k}}{\tau_k}&=\frac{1}{\tau_k}\left(\rho^{k}(\boldsymbol{u}_{\dagger}^{k} - \boldsymbol{u}^{k}),\boldsymbol{u}_{\dagger}^{k}\right)-\frac{1}{2\tau_k}\left(\rho^{k},|\boldsymbol{u}_{\dagger}^{k} - \boldsymbol{u}^{k}|^2\right)\nonumber\\ &\leq\frac{1}{\tau_k}\left(\rho^{k}(\boldsymbol{u}_{\dagger}^{k} - \boldsymbol{u}^{k}),\boldsymbol{u}_{\dagger}^{k}\right)\nonumber\\ &={-} \sum_{l=1}^N(c_l^{k}\,\boldsymbol{\nabla}\mu_{c,l}^{k+1},\boldsymbol{u}_{\dagger}^{k} ), \end{align}

where we have used (3.11c). Using (3.11e), we derive the difference between $\mathcal {U}_\star ^{k}$ and $\mathcal {U}_{\dagger} ^{k}$ as

(3.28)\begin{align} \frac{\mathcal{U}_\star^{k}-\mathcal{U}_{\dagger}^{k}}{\tau_k}&=\frac{1}{\tau_k}\left(\rho^{k}(\boldsymbol{u}_\star^{k} - \boldsymbol{u}_{\dagger}^{k}),\boldsymbol{u}_\star^{k}\right)-\frac{1}{2\tau_k}\left(\rho^{k},|\boldsymbol{u}_\star^{k} - \boldsymbol{u}_{\dagger}^{k}|^2\right)\nonumber\\ &\leq\frac{1}{\tau_k}\left(\rho^{k}(\boldsymbol{u}_\star^{k} - \boldsymbol{u}_{\dagger}^{k}),\boldsymbol{u}_\star^{k}\right)\nonumber\\ &={-} \left(\phi^{k}\,\boldsymbol{\nabla}(\mu_\phi^{k+1}+\mu_{c\phi}^{k+1})+\boldsymbol{\nabla} p^{k+1}-\rho^k \boldsymbol{\mathsf{g}},\boldsymbol{u}_\star^{k} \right). \end{align}

Using (3.11g), we deduce that

(3.29)\begin{align} (\boldsymbol{u}_\star^{k}, \boldsymbol{\nabla} p^{k+1})&={-}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_\star^{k}, p^{k+1})\nonumber\\ &=(\lambda\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1}, p^{k+1})\nonumber\\ &={-}(\boldsymbol{J}_\phi^{k+1}, \lambda\,\boldsymbol{\nabla} p^{k+1}). \end{align}

It is deduced from (3.28) and (3.29) that

(3.30)\begin{equation} \frac{\mathcal{U}_\star^{k}-\mathcal{U}_{\dagger}^{k}}{\tau_k}={-} \left(\phi^{k}\boldsymbol{\nabla}(\mu_\phi^{k+1}+\mu_{c\phi}^{k+1})-\rho^k \boldsymbol{\mathsf{g}},\boldsymbol{u}_\star^{k} \right) +(\boldsymbol{J}_\phi^{k+1}, \lambda\,\boldsymbol{\nabla} p^{k+1}) . \end{equation}

For $a,b,c,d\in \mathbb {R}$ and $c\geq 0$, we have

(3.31)\begin{equation} ab^2-cd^2=2b(ab-cd)-c(b-d)^2-b^2(a-c)\leq2b(ab-cd)-b^2(a-c). \end{equation}

The difference between $\mathcal {U}^{k+1}$ and $\mathcal {U}_\star ^k$ is estimated using (3.31) as

(3.32)\begin{align} \mathcal{U}^{k+1}-\mathcal{U}_\star^{k} &=\tfrac{1}{2}\left(\rho^{k+1},|\boldsymbol{u}^{k+1}|^2\right)-\tfrac{1}{2}\left(\rho^k,|\boldsymbol{u}_\star^{k}|^2\right)\nonumber\\ &\leq \left(\rho^{k+1}\boldsymbol{u}^{k+1}-\rho^k\boldsymbol{u}_\star^{k},\boldsymbol{u}^{k+1}\right) -\tfrac{1}{2}\left(\rho^{k+1}-\rho^k,|\boldsymbol{u}^{k+1}|^2\right). \end{align}

Substituting (3.11h) and (3.17) into (3.32) yields

(3.33)\begin{align} \frac{\mathcal{U}^{k+1}-\mathcal{U}_\star^{k}}{\tau_k} &\leq{-}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\left((\rho^k\boldsymbol{u}_\star^k+\chi\boldsymbol{J}_\phi^{k+1})\otimes \boldsymbol{u}^{k+1}\right),\boldsymbol{u}^{k+1}\right)\nonumber\\ &\quad +\left( \boldsymbol{\nabla}\boldsymbol{\cdot} \eta(\phi^{k})\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ {\boldsymbol{\nabla}\boldsymbol{u}^{k+1}}^{\rm T}),\boldsymbol{u}^{k+1}\right)\nonumber\\ &\quad +\tfrac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho^k\boldsymbol{u}_\star^k ) +\chi\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1},|\boldsymbol{u}^{k+1}|^2\right)\nonumber\\ &={-}\left((\rho^k\boldsymbol{u}_\star^k+\chi\boldsymbol{J}_\phi^{k+1})\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^{k+1} ,\boldsymbol{u}^{k+1}\right)\nonumber\\ &\quad +\left( \boldsymbol{\nabla}\boldsymbol{\cdot} \eta(\phi^{k})\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ {\boldsymbol{\nabla}\boldsymbol{u}^{k+1}}^{\rm T}),\boldsymbol{u}^{k+1}\right)\nonumber\\ &\quad -\tfrac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho^k\boldsymbol{u}_\star^k ) +\chi\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1},|\boldsymbol{u}^{k+1}|^2\right)\nonumber\\ &={-}\tfrac{1}{2}\left\|\eta(\phi^{k})^{1/2}\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ {\boldsymbol{\nabla}\boldsymbol{u}^{k+1}}^{\rm T})\right\|^2, \end{align}

where we have used the identities

(3.34)\begin{gather} (\rho^{k}\boldsymbol{u}_\star^{k}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^{k+1},\boldsymbol{u}^{k+1}) +\tfrac{1}{2}(\boldsymbol{\nabla}\boldsymbol{\cdot}\rho^{k}\boldsymbol{u}_\star^{k},|\boldsymbol{u}^{k+1}|^2)=0, \end{gather}

(3.35)\begin{gather}(\boldsymbol{J}_\phi^{k+1}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^{k+1},\boldsymbol{u}^{k+1})+\tfrac{1}{2}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{J}_\phi^{k+1},|\boldsymbol{u}^{k+1}|^2)=0. \end{gather}

Combining (3.27), (3.30) and (3.33) yields

(3.36)\begin{align} \frac{\mathcal{U}^{k+1}-\mathcal{U}^{k}}{\tau_k} &=\frac{\mathcal{U}^{k+1}-\mathcal{U}_\star^{k}+\mathcal{U}_\star^{k}-\mathcal{U}^{k}_{\dagger} +\mathcal{U}_{\dagger}^{k}-\mathcal{U}^{k}}{\tau_k}\nonumber\\ &\leq{-} \sum_{l=1}^N(c_l^{k}\,\boldsymbol{\nabla}\mu_{c,l}^{k+1},\boldsymbol{u}_{\dagger}^{k} )- \left(\phi^{k}\,\boldsymbol{\nabla}(\mu_\phi^{k+1}+\mu_{c\phi}^{k+1})-\rho^k \boldsymbol{\mathsf{g}},\boldsymbol{u}_\star^{k} \right) \nonumber\\ &\quad +(\boldsymbol{J}_\phi^{k+1}, \lambda\,\boldsymbol{\nabla} p^{k+1}) -\tfrac{1}{2}\left\|\eta(\phi^{k})^{1/2}\,(\boldsymbol{\nabla}\boldsymbol{u}^{k+1}+ \boldsymbol{\nabla}{\boldsymbol{u}^{k+1}}^{\rm T})\right\|^2. \end{align}

Finally, combining (3.25) and (3.36) yields (3.20).

3.3. Fully discrete scheme

For the spatial discretization, we employ the finite difference/volume methods on staggered grids (Tryggvason, Scardovelli & Zaleski 2011; Kou, Sun & Wang 2018) and treat the convection terms by the upwind strategy. We consider a rectangular domain $\varOmega =[0,L_x]\times [0,L_y]$, where $L_x>0$ and $L_y>0$. The velocity vector is expressed as $\boldsymbol {u}=[u,v]^\textrm {T}$. For simplification in presentation, we use a uniform division of $\varOmega$ as $0=x_0< x_1<\cdots < x_{n_x}=L_x$ and $0=y_0< y_1<\cdots < y_{n_y}=L_y,$ where $n_x$ and $n_y$ are positive integers. The mesh size is denoted by $h=x_{i+1}-x_i=y_{j+1}-y_j$, where the subscripts $0\leq i< n_x$ and $0\leq j< n_y$ are integers. The midpoints are defined as $x_{i+{1}/{2}} =x_i+\frac {1}{2}h$ and $y_{j+{1}/{2}} =y_j+\frac {1}{2}h$. The basic idea of the staggered grids is that different variables are discretized on different grids that are shifted by half a grid point. In order to represent different discrete variables, we define the following discrete function spaces:

(3.37a)\begin{gather} \mathcal{S}_c=\{c: (x_{i+{1}/{2}},y_{j+{1}/{2}})\mapsto\mathbb{R},\ 0\leq i\leq n_x-1,\ 0\leq j\leq n_y-1\}, \end{gather}

(3.37b)\begin{gather}\mathcal{S}_u=\{u: (x_{i},y_{j+{1}/{2}})\mapsto\mathbb{R},\ 0\leq i\leq n_x,\ 0\leq j\leq n_y-1\}, \end{gather}

(3.37c)\begin{gather}\mathcal{S}_v=\{v: (x_{i+{1}/{2}},y_{j})\mapsto\mathbb{R},\ 0\leq i\leq n_x-1,\ 0\leq j\leq n_y\}, \end{gather}

(3.37d)\begin{gather}\mathcal{S}_{pu}=\{\psi: (x_{i},y_{j})\mapsto\mathbb{R},\ 1\leq i\leq n_x-1,\ 0\leq j\leq n_y\}, \end{gather}

(3.37e)\begin{gather}\mathcal{S}_{pv}=\{\varphi: (x_{i},y_{j})\mapsto\mathbb{R},\ 0\leq i\leq n_x,\ 1\leq j\leq n_y-1\}. \end{gather}

For a discrete function $c\in \mathcal {S}_c$, we write $c_{i+{1}/{2},j+{1}/{2}}=c(x_{i+{1}/{2}},y_{j+{1}/{2}})$ in its component form. The other functions in the above function spaces are denoted likewise. Taking the homogeneous Neumann boundary conditions into consideration, we further define subsets of $\mathcal {S}_u$ and $\mathcal {S}_v$ as

(3.38a)\begin{gather} \mathcal{S}^0_u=\{u\in\mathcal{S}_u\mid u_{0,j+{1}/{2}}=u_{n_x,j+{1}/{2}}=0,\ 0\leq j\leq n_y-1\}, \end{gather}

(3.38b)\begin{gather}\mathcal{S}^0_v=\{v\in\mathcal{S}_v\mid v_{i+{1}/{2},0}=v_{i+{1}/{2},n_y}=0,\ 0\leq i\leq n_x-1\}. \end{gather}

For the purpose of treating conveniently the boundary conditions for the momentum balance equation, we introduce two extended discrete function spaces:

(3.39a)\begin{gather} \mathcal{S}_u^e=\left\{u: (x_{i},y_{j+{1}/{2}})\mapsto\begin{cases} \mathbb{R}, & 1\leq i\leq n_x-1,\ 0\leq j\leq n_y-1\\ 0, & i\in\{0,n_x\},\ 0\leq j\leq n_y-1\\ 0, & 1\leq i\leq n_x-1,j\in\{{-}1,n_y\}\\ \end{cases}\right\}, \end{gather}

(3.39b)\begin{gather}\mathcal{S}_v^e=\left\{v: (x_{i+{1}/{2}},y_{j})\mapsto\begin{cases} \mathbb{R}, & 0\leq i\leq n_x-1,\ 1\leq j\leq n_y-1\\ 0, & 0\leq i\leq n_x-1,\ j\in\{0,n_y\}\\ 0, & i\in\{{-}1,n_x\},\ 1\leq j\leq n_y-1 \end{cases}\right\}. \end{gather}

We now introduce the discrete operators, which will be used to formulate the fully discrete equations. For $c\in \mathcal {S}_c$, we define the difference operators

(3.40a)\begin{gather} d_x^cc_{i,j+{1}/{2}}= \frac{1}{h}\,(c_{i+{1}/{2},j+{1}/{2}}-c_{i-{1}/{2},j+{1}/{2}}),\quad 1\leq i\leq n_x-1,\ 0\leq j\leq n_y-1, \end{gather}

(3.40b)\begin{gather}d_y^cc_{i+{1}/{2},j}= \frac{1}{h}\,(c_{i+{1}/{2},j+{1}/{2}}-c_{i+{1}/{2},j-{1}/{2}}),\quad 0\leq i\leq n_x-1,\ 1\leq j\leq n_y-1, \end{gather}

where $d_x^cc\in \mathcal {S}_u^0$ and $d_y^cc\in \mathcal {S}_v^0$. For $c\in \mathcal {S}_c$, we define the average values of $c$ as

(3.41a)\begin{gather} \bar{c}_{i,j+{1}/{2}}=\begin{cases} \frac{1}{2}(c_{i+{1}/{2},j+{1}/{2}}+c_{i-{1}/{2},j+{1}/{2}}), & 1\leq i\leq n_x-1,\ 0\leq j\leq n_y-1,\\ c_{i+{1}/{2},j+{1}/{2}}, & i=0,\ 0\leq j\leq n_y-1,\\ c_{i-{1}/{2},j+{1}/{2}}, & i=n_x,\ 0\leq j\leq n_y-1, \end{cases} \end{gather}

(3.41b)\begin{gather}\bar{c}_{i+\frac{1}{2},j}=\begin{cases} \frac{1}{2}(c_{i+{1}/{2},j+{1}/{2}}+c_{i+{1}/{2},j-{1}/{2}}), & 0\leq i\leq n_x-1, \ 1\leq j\leq n_y-1,\\ c_{i+{1}/{2},j+{1}/{2}}, & j=0,\ 0\leq i\leq n_x-1.\\ c_{i+{1}/{2},j-{1}/{2}}, & j=n_y,\ 0\leq i\leq n_x-1. \end{cases} \end{gather}

The upwind strategy is applied to treat the convection terms, and for this purpose, we introduce the upwind values of $c\in \mathcal {S}_c$ in the presence of the velocity components $u\in \mathcal {S}_u^0$ and $v\in \mathcal {S}_v^0$ as

(3.42a)\begin{gather} \hat{c}_{i,j+{1}/{2}}=\begin{cases} c_{i-{1}/{2},j+{1}/{2}}, & u_{i,j+{1}/{2}}\geq 0,\\ c_{i+{1}/{2},j+{1}/{2}}, & u_{i,j+{1}/{2}}< 0, \end{cases} \end{gather}

(3.42b)\begin{gather}\hat{c}_{i+\frac{1}{2},j}=\begin{cases} c_{i+{1}/{2},j-{1}/{2}}, & v_{i+{1}/{2},j}\geq 0,\\ c_{i+{1}/{2},j+{1}/{2}}, & v_{i+{1}/{2},j}< 0. \end{cases} \end{gather}

For $u\in \mathcal {S}_u$ and $v\in \mathcal {S}_v$, we define the difference operators

(3.43a)\begin{gather} d_x^uu_{i+{1}/{2},j+{1}/{2}}=\frac{1}{h}\,(u_{i+1,j+{1}/{2}}-u_{i,j+{1}/{2}}), \end{gather}

(3.43b)\begin{gather}d_y^vv_{i+{1}/{2},j+{1}/{2}}=\frac{1}{h}\,(v_{i+{1}/{2},j+1}-v_{i+{1}/{2},j}). \end{gather}

For $u\in \mathcal {S}^e_u$ and $v\in \mathcal {S}^e_v$, we define the difference and average operators

(3.44a,b)\begin{gather} d_y^uu_{i,j}=\frac{1}{h}\,(u_{i,j+{1}/{2}}-u_{i,j-{1}/{2}}),\quad d_x^vv_{i,j}= \frac{1}{h}\,(v_{i+{1}/{2},j}-v_{i-{1}/{2},j}), \end{gather}

(3.45a,b)\begin{gather}a_x^u u_{i+{1}/{2},j+{1}/{2}}=\frac{1}{2}\,(u_{i,j+{1}/{2}}+u_{i+1,j+{1}/{2}}),\quad a_y^u u_{i,j}=\frac{1}{2}\,(u_{i,j-{1}/{2}}+u_{i,j+{1}/{2}}), \end{gather}

(3.46a,b)\begin{gather}a_y^v v_{i+{1}/{2},j+{1}/{2}}=\frac{1}{2}\,(v_{i+{1}/{2},j}+v_{i+{1}/{2},j+1}),\quad a_x^v v_{i,j}=\frac{1}{2}\,(v_{i-{1}/{2},j}+v_{i+{1}/{2},j}). \end{gather}

We can see that $d_x^uu\in \mathcal {S}_c$, $d_y^vv\in \mathcal {S}_c$, $d_y^uu\in \mathcal {S}_{pu}$ and $d_x^vv\in \mathcal {S}_{pv}$.

We are ready to present the fully discrete scheme. At each time step, we seek $c_l^{k+1}\in \mathcal {S}_c$, $\phi ^{k+1}\in \mathcal {S}_c$, $p^{k+1}\in \mathcal {S}_c$, $u^{k+1}\in \mathcal {S}^e_u$ and $v^{k+1}\in \mathcal {S}^e_v$ for given $c_l^{k}\in \mathcal {S}_c$, $\phi ^{k}\in \mathcal {S}_c$, $u^{k}\in \mathcal {S}^e_u$ and $v^{k}\in \mathcal {S}^e_v$. The fully discrete chemical potentials $\mu _{c,l}^{k+1},\mu _{c\phi }^{k+1}\in \mathcal {S}_c$ still take the pointwise forms given (3.4) and (3.5). Fully discrete intermediate velocity components $u^k_{{\dagger} }\in \mathcal {S}_u^0$ and $v^k_{{\dagger} }\in \mathcal {S}_v^0$ are expressed as

(3.47a,b)\begin{equation} u^k_{{\dagger}} = u^k-\sum_{l=1}^N\frac{\tau_k }{\bar{\rho}^k}\,\hat{c}_l^k d_x^c\mu_{c,l}^{k+1},\quad v^k_{{\dagger}} = v^k-\sum_{l=1}^N\frac{\tau_k }{\bar{\rho}^k}\,\hat{c}_l^k d_y^c\mu_{c,l}^{k+1}. \end{equation}

Let $\bar {{\mathcal {D}}}_{l,m}^k={\mathcal {D}}_{l,m}(\bar {\phi }^{k},\bar {\boldsymbol {c}}^{k})$. The solute concentration equations are fully discretized as

(3.48)\begin{equation} \frac{ c_l^{k+1}-c_l^{k}}{\tau_k } + d_x^u(u_{{\dagger}}^k\hat{c}_l^k+J_{c,l,x}^{k+1})+d_y^v(v_{{\dagger}}^k\hat{c}_l^k+J_{c,l,y}^{k+1}) =0,\quad l=1,\ldots,N, \end{equation}

where $J_{c,l,x}^{k+1}\in \mathcal {S}_u$ and $J_{c,l,y}^{k+1}\in \mathcal {S}_v$ are defined as

(3.49a,b)\begin{equation} J_{c,l,x}^{k+1}={-}\sum_{m=1}^N\bar{{\mathcal{D}}}_{l,m}^k d_x^c \mu_{c,m}^{k+1},\quad J_{c,l,y}^{k+1}={-}\sum_{m=1}^N\bar{{\mathcal{D}}}_{l,m}^k d_y^c \mu_{c,m}^{k+1}. \end{equation}

Let $\boldsymbol{\mathsf{g}}=[g_x,g_y]$. The fully discrete solvent transport equation is expressed as

(3.50)\begin{equation} \frac{\phi^{k+1}-\phi^{k}}{\tau_k } + d_x^u(u_{{\star}}^k\hat{\phi}^k+J_{\phi,x}^{k+1})+d_y^v(v_{{\star}}^k\hat{\phi}^k+J_{\phi,y}^{k+1}) =0, \end{equation}

where

(3.51)\begin{gather} u^k_{{\star}} = u_{\dagger}^k-\frac{\tau_k }{\bar{\rho}^k}\left(\hat{\phi}^k d_x^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})+d_x^cp^{k+1}-\hat{\rho}^k g_x\right), \end{gather}

(3.52)\begin{gather}v^k_{{\star}} = v_{\dagger}^k-\frac{\tau_k }{\bar{\rho}^k} \left(\hat{\phi}^k d_y^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})+d_y^cp^{k+1}-\hat{\rho}^k g_y\right), \end{gather}

(3.53)\begin{gather}J_{\phi,x}^{k+1}={-}M_\phi(\bar{\phi}^{k}) \left( d_x^c (\mu_\phi^{k+1}+\mu_{c\phi}^{k+1}+\lambda p^{k+1})- \chi g_x\right), \end{gather}

(3.54)\begin{gather}J_{\phi,y}^{k+1}={-}M_\phi(\bar{\phi}^{k})\left( d_y^c (\mu_\phi^{k+1}+\mu_{c\phi}^{k+1}+\lambda p^{k+1})- \chi g_y\right), \end{gather}

(3.55)\begin{gather}\mu_\phi^{k+1}=\frac{\sigma}{\epsilon}\left(\ell(\phi^{k},\phi^{k+1})-\epsilon^2\left(d_x^u(d_x^c \phi^{k+1})+d_y^v(d_y^c \phi^{k+1})\right)\right). \end{gather}

The fully discrete form of (3.11g) is expressed as

(3.56)\begin{equation} d_x^u (u_\star^k +\lambda J_{\phi,x}^{k+1}) +d_y^v (v_\star^k +\lambda J_{\phi,y}^{k+1})=0 . \end{equation}

In order to discretize the momentum balance equation, we define $U_u^{k}\in \mathcal {S}_c$, $U_v^{k}\in \mathcal {S}_{pv}$, $V_u^{k}\in \mathcal {S}_{pu}$ and $V_v^{k}\in \mathcal {S}_c$ as

(3.57a)\begin{gather} U_u^k= a_x^u(\hat{\rho}^k u_*^k +\chi J_{\phi,x}^{k+1}),\quad U^k_{v} = a_y^u(\hat{\rho}^k u_*^k +\chi J_{\phi,x}^{k+1}), \end{gather}

(3.57b)\begin{gather}V^k_{u}= a_x^v(\hat{\rho}^k v_*^k +\chi J_{\phi,y}^{k+1}),\quad V^k_{v} = a_y^v(\hat{\rho}^k v_*^k +\chi J_{\phi,y}^{k+1}). \end{gather}

For $u\in {\mathcal {S}}^e_u$ and $v\in {\mathcal {S}}^e_v$, we also define their upwind values as

(3.58a)\begin{gather} \hat{u}_{i+{1}/{2},j+{1}/{2}}=\left\{\begin{array}{@{}ll@{}} u_{i,j+{1}/{2}}, & U_{u,i+{1}/{2},j+{1}/{2}}\geq 0,\\ u_{i+1,j+{1}/{2}}, & U_{u,i+{1}/{2},j+{1}/{2}}< 0, \end{array}\right. \end{gather}

(3.58b)\begin{gather}\hat{u}_{i,j}=\left\{\begin{array}{@{}ll@{}} u_{i,j-{1}/{2}}, & V_{u,i,j}\geq 0,\\ u_{i,j+{1}/{2}}, & V_{u,i,j}< 0, \end{array}\right. \end{gather}

(3.58c)\begin{gather}\hat{v}_{i+{1}/{2},j+{1}/{2}}=\left\{\begin{array}{@{}ll@{}} v_{i+{1}/{2},j}, & V_{v,i+{1}/{2},j+{1}/{2}}\geq 0,\\ v_{i+{1}/{2},j+1}, & V_{v,i+{1}/{2},j+{1}/{2}}< 0, \end{array}\right. \end{gather}

(3.58d)\begin{gather}\hat{v}_{i,j}=\left\{\begin{array}{@{}ll@{}} v_{i-{1}/{2},j}, & U_{v,i,j}\geq 0,\\ v_{i+{1}/{2},j}, & U_{v,i,j}< 0. \end{array}\right. \end{gather}

For $\psi \in \mathcal {S}_{pu}\cup \mathcal {S}_{pv}$, we define the difference operators

(3.59a,b)\begin{equation} d_x^p \psi_{i+\frac{1}{2},j}=\frac{1}{h}\,(\psi_{i+1,j}-\psi_{i,j}), \quad d_y^p \psi_{i,j+{1}/{2}}=\frac{1}{h}\,(\psi_{i,j+1}-\psi_{i,j}). \end{equation}

Let $\eta ^k=\eta (\phi ^k)\in \mathcal {S}_c$, and let $\bar {\eta }_{i,j}^k$ be the average values of $\eta ^k$ at the vertices $(x_i,y_j)$. We denote

(3.60a–c)\begin{align} \varPsi^{k+1}=2\eta^kd_x^uu^{k+1}, \quad \varUpsilon^{k+1}= 2\eta^kd_y^vv^{k+1}, \quad \varTheta^{k+1}= \bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), \end{align}

where $\varPsi ^{k+1}\in \mathcal {S}_c$, $\varUpsilon ^{k+1}\in \mathcal {S}_c$ and $\varTheta ^{k+1}\in \mathcal {S}_{pu}\cup \mathcal {S}_{pv}$. The fully discrete scheme for the momentum balance equation (3.11h) is expressed as

(3.61)\begin{gather} \frac{\bar{\rho}^{k+1}u^{k+1}-\bar{\rho}^{k}u_{{\star}}^k}{\tau_k } + d_x^c(U_u^k\hat{u}^{k+1} )+d_y^p(V_u^k\hat{u}^{k+1}) = d_x^c\varPsi^{k+1}+d_y^p\varTheta^{k+1}, \end{gather}

(3.62)\begin{gather}\frac{\bar{\rho}^{k+1}v^{k+1}-\bar{\rho}^{k}v_{{\star}}^k}{\tau_k } +d_x^p(U_v^k\hat{v}^{k+1})+d_y^c(V_v^k \hat{v}^{k+1}) =d_x^p\varTheta^{k+1}+d_y^c\varUpsilon^{k+1}. \end{gather}

In what follows, we prove that the fully discrete scheme conserves the masses of both solutes and solvent fluids as well as following a discrete energy dissipation law. To this end, we define the following discrete inner products:

(3.63a)\begin{gather} (c,c')_c=h^2\sum_{i=0}^{n_x-1}\sum_{j=0}^{n_y-1}c_{i+{1}/{2},j+{1}/{2}}c'_{i+{1}/{2},j+{1}/{2}},\quad c,c' \in\mathcal{S}_c, \end{gather}

(3.63b)\begin{gather}( \psi,\psi')_{pu}=h^2\sum_{i=1}^{n_x-1}\sum_{j=0}^{n_y}\psi_{i,j}\psi'_{i,j},\quad \psi,\psi'\in\mathcal{S}_{pu}, \end{gather}

(3.63c)\begin{gather}(\varphi,\varphi')_{pv}=h^2\sum_{i=0}^{n_x}\sum_{j=1}^{n_y-1}\varphi_{i,j}\varphi'_{i,j},\quad \varphi,\varphi'\in\mathcal{S}_{pv}, \end{gather}

(3.63d)\begin{gather}( u,u')_u=h^2\sum_{i=1}^{n_x-1}\sum_{j=0}^{n_y-1}u_{i,j+{1}/{2}}u'_{i,j+{1}/{2}},\quad u,u'\in\mathcal{S}_u\cup\mathcal{S}_u^0\cup {\mathcal{S}}^e_u, \end{gather}

(3.63e)\begin{gather}(v,v')_v=h^2\sum_{i=0}^{n_x-1}\sum_{j=1}^{n_y-1}v_{i+{1}/{2},j}v'_{i+{1}/{2},j},\quad v,v'\in\mathcal{S}_v\cup\mathcal{S}_v^0\cup \mathcal{S}^e_v. \end{gather}

In addition, for $\psi,\psi '\in \mathcal {S}_{pu}\cup \mathcal {S}_{pv}$, we define

(3.64)\begin{equation} ( \psi,\psi')_{puv}=h^2\sum_{i=1}^{n_x-1}\sum_{j=1}^{n_y-1}\psi_{i,j}\psi'_{i,j}. \end{equation}

The following summation-by-parts formulas are obtained by direct calculations (Kou et al. 2018):

(3.65a)\begin{gather} (u,d_x^c c)_u={-}( d_x^u u,c)_c,\quad u\in \mathcal{S}_u^0\cup {\mathcal{S}}^e_u,\ c\in \mathcal{S}_c, \end{gather}

(3.65b)\begin{gather}( v,d_y^c c)_v={-}( d_y^v v,c)_c,\quad v\in \mathcal{S}_v^0\cup {\mathcal{S}}^e_v,\ c\in \mathcal{S}_c, \end{gather}

(3.65c)\begin{gather}( d_y^p \psi,u)_u={-}( \psi,d_y^u u)_{pu},\quad u\in {\mathcal{S}}^e_u,\ \psi\in \mathcal{S}_{pu}, \end{gather}

(3.65d)\begin{gather}(d_x^p \varphi,v)_v={-}( \varphi,d_x^v v)_{pv},\quad v\in {\mathcal{S}}^e_v,\ \varphi\in \mathcal{S}_{pv}. \end{gather}

Theorem 3.3 The fully discrete scheme conserves total masses of solutes and solvent fluids as

(3.66)\begin{gather} (c_l^{k+1},1)_c=(c_l^{k},1)_c=\cdots= (c_l^{0},1)_c,\quad k\geq0,\ l=1,\ldots, N, \end{gather}

(3.67)\begin{gather}(\rho^{k+1},1)_c=(\rho^{k},1)_c=\cdots= (\rho^{0},1)_c,\quad k\geq0. \end{gather}

Proof. Taking the discrete inner product of (3.48), we obtain

(3.68)\begin{equation} \frac{ 1}{\tau_k }\,(c_l^{k+1}-c_l^{k},1)_c + \left(d_x^u(u_{{\dagger}}^k\hat{c}_l^k+J_{c,l,x}^{k+1}),1\right)_c +\left(d_y^v(v_{{\dagger}}^k\hat{c}_l^k+J_{c,l,y}^{k+1}),1\right)_c=0. \end{equation}

Applying (3.65a) and (3.65b) to (3.68) yields

(3.69)\begin{equation} \frac{ 1}{\tau_k }\,(c_l^{k+1}-c_l^{k},1)_c =0, \end{equation}

which leads to (3.66) by induction. From (3.50), we can derive the mass balance equation of solvents as

(3.70)\begin{equation} \frac{ \rho^{k+1}-\rho^k }{\tau_k} + d_x^u(u_{{\star}}^k\hat{\rho}^k+\chi J_{\phi,x}^{k+1})+d_y^v(v_{{\star}}^k\hat{\rho}^k+\chi J_{\phi,y}^{k+1}) =0, \end{equation}

where $\chi =\varrho _1-\varsigma \varrho _2$. Taking the discrete inner product of (3.70) and then using (3.65a) and (3.65b), we can obtain (3.67) by induction.

We now prove that the fully discrete scheme follows a discrete energy dissipation law. The total discrete energy is defined as

(3.71)\begin{equation} \mathcal{E}_h^k = \mathcal{A}_h^k+\mathcal{F}_h^k+\mathcal{U}_h^k+\mathcal{G}_h^k, \end{equation}

where

(3.72)\begin{align} \left.\begin{gathered} \mathcal{A}_h^k =\left(\! A(\boldsymbol{c}^k,\phi^k) ,1\!\right)_c,\quad \mathcal{F}_h^k = \frac{\sigma}{\epsilon}\left(\! f(\phi^{k}),1\!\right)_c+\frac{\epsilon \sigma}{2}\left(\! (d_x^c\phi^{k},d_x^c\phi^{k})_u+ (d_y^c\phi^{k},d_y^c\phi^{k})_v\!\right),\\ \mathcal{U}_h^k =\frac{1}{2}\,( \bar{\rho}^ku^k,u^k)_u+\frac{1}{2}\,( \bar{\rho}^kv^k,v^k)_v,\quad \mathcal{G}_h^k =( \rho^k, gz)_c. \end{gathered}\right\} \end{align}

Theorem 3.4 Given that (3.7) holds, the fully discrete scheme satisfies the following discrete energy dissipation law:

(3.73)\begin{align} \frac{\mathcal{E}_h^{k+1}-\mathcal{E}_h^{k} }{\tau_k} &\leq{-} \sum_{l=1}^N\sum_{m=1}^N ( \bar{{\mathcal{D}}}_{l,m}^kd_x^c\mu_{c,m}^{k+1}, d_x^c\mu_{c,l}^{k+1})_u -\sum_{l=1}^N\sum_{m=1}^N(\bar{{\mathcal{D}}}_{l,m}^kd_y^c\mu_{c,m}^{k+1}, d_y^c\mu_{c,l}^{k+1})_v\nonumber\\ &\quad -\left(M_\phi(\bar{\phi}^{k})^{{-}1}J_{\phi,x}^{k+1},J_{\phi,x}^{k+1}\right)_u -\left(M_\phi(\bar{\phi}^{k})^{{-}1}J_{\phi,y}^{k+1}, J_{\phi,y}^{k+1}\right)_v \nonumber\\ &\quad - 2(\eta^kd_x^uu^{k+1}, d_x^u u^{k+1})_c- 2(\eta^kd_y^vv^{k+1}, d_y^v v^{k+1})_c\nonumber\\ &\quad - \left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_y^uu^{k+1}+d_x^vv^{k+1}\right)_{puv} . \end{align}

Proof. Taking into account (3.6), (3.48) and (3.50), and applying (3.65a) and (3.65b), we deduce that

(3.74)\begin{align} \frac{\mathcal{A}_h^{k+1}-\mathcal{A}_h^{k}}{\tau_k} &\leq\sum_{l=1}^N\left(\mu_{c,l}^{k+1},\frac{ c_l^{k+1}- c_l^{k} }{\tau_k}\right)_c+\left(\mu_{c\phi}^{k+1},\frac{ \phi^{k+1}- \phi^{k} }{\tau_k}\right)_c\nonumber\\ &={-} \sum_{l=1}^N(d_x^u(u_{{\dagger}}^k\hat{c}_l^k+J_{c,l,x}^{k+1})+d_y^v(v_{{\dagger}}^k\hat{c}_l^k+J_{c,l,y}^{k+1}) ,\mu_{c,l}^{k+1})_c \nonumber\\ &\quad - \left(d_x^u(u_{{\star}}^k\hat{\phi}^k+J_{\phi,x}^{k+1})+d_y^v(v_{{\star}}^k\hat{\phi}^k+J_{\phi,y}^{k+1}), \mu_{c\phi}^{k+1}\right)_c \nonumber\\ &= \sum_{l=1}^N ( u_{{\dagger}}^k\hat{c}_l^k+J_{c,l,x}^{k+1}, d_x^c\mu_{c,l}^{k+1})_u +\sum_{l=1}^N (v_{{\dagger}}^k\hat{c}_l^k+J_{c,l,y}^{k+1}, d_y^c\mu_{c,l}^{k+1})_v\nonumber\\ &\quad +( u_{{\star}}^k\hat{\phi}^k+J_{\phi,x}^{k+1}, d_x^c\mu_{c\phi}^{k+1})_u +( v_{{\star}}^k\hat{\phi}^k+J_{\phi,y}^{k+1}, d_y^c\mu_{c\phi}^{k+1})_v . \end{align}

Noticing that

(3.75)\begin{align} \tfrac{1}{2} (d_x^c\phi^{k+1},d_x^c\phi^{k+1})_u-\tfrac{1}{2} (d_x^c\phi^{k},d_x^c\phi^{k})_u &\leq\left(d_x^c\phi^{k+1},d_x^c(\phi^{k+1}-\phi^{k})\right)_u \nonumber\\ &={-}\left(d_x^u(d_x^c \phi^{k+1}),\phi^{k+1}-\phi^{k}\right)_c , \end{align}

(3.76)\begin{align} \frac{1}{2} (d_y^c\phi^{k+1},d_y^c\phi^{k+1})_v-\tfrac{1}{2} (d_y^c\phi^{k},d_y^c\phi^{k})_v &\leq\left(d_y^c\phi^{k+1},d_y^c(\phi^{k+1}-\phi^{k})\right)_v \nonumber\\&={-}\left(d_y^v(d_y^c \phi^{k+1}),\phi^{k+1}-\phi^{k}\right)_c, \end{align}

we deduce that

(3.77)\begin{align} \frac{\mathcal{F}_h^{k+1}-\mathcal{F}_h^{k}}{\tau_k} &\leq \left(\mu_\phi^{k+1}, \frac{\phi^{k+1}- \phi^{k}}{\tau_k}\right)_c \nonumber\\ &={-} \left(d_x^u(u_{{\star}}^k\hat{\phi}^k+J_{\phi,x}^{k+1})+d_y^v(v_{{\star}}^k\hat{\phi}^k+J_{\phi,y}^{k+1}) , \mu_\phi^{k+1}\right)_c \nonumber\\ &=(u_{{\star}}^k,\hat{\phi}^kd_x^c\mu_\phi^{k+1})_u+(J_{\phi,x}^{k+1},d_x^c\mu_\phi^{k+1})_u\nonumber\\ &\quad +(v_{{\star}}^k,\hat{\phi}^kd_y^c\mu_\phi^{k+1})_v+(J_{\phi,y}^{k+1}, d_y^c\mu_\phi^{k+1})_v . \end{align}

The variation of gravitational potential energy with time steps is estimated using (3.70) as

(3.78)\begin{align} \frac{\mathcal{G}_h^{k+1}-\mathcal{G}_h^{k}}{\tau_k} &=\frac{1}{\tau_k}\,( \rho^{k+1}-\rho^k,gz)_c\nonumber\\ &={-}\left( d_x^u(\hat{\rho}^ku_\star^k +\chi J_{\phi,x}^{k+1}), gz\right)_c-\left( d_y^v(\hat{\rho}^kv_\star^k +\chi J_{\phi,y}^{k+1}), gz\right)_c \nonumber\\ &={-}( u_\star^k, \hat{\rho}^k g_x)_u-( J_{\phi,x}^{k+1}, \chi g_x)_u-( v_\star^k, \hat{\rho}^k g_y)_v-( J_{\phi,y}^{k+1}, \chi g_y)_v. \end{align}

To estimate the variation of the kinetic energy with time steps, we introduce the following intermediate kinetic energies:

(3.79a)\begin{gather} \mathcal{U}_{h{\dagger}}^{k}=\tfrac{1}{2}( \bar{\rho}^ku_{\dagger}^k,u_{\dagger}^k)_u+\tfrac{1}{2}( \bar{\rho}^kv_{\dagger}^k,v_{\dagger}^k)_v, \end{gather}

(3.79b)\begin{gather}\mathcal{U}_{h\star}^{k}=\tfrac{1}{2}( \bar{\rho}^ku_\star^k,u_\star^k)_u+\tfrac{1}{2}( \bar{\rho}^kv_\star^k,v_\star^k)_v. \end{gather}

The difference between $\mathcal {U}_{h{\dagger} }^{k}$ and $\mathcal {U}^{k}$ is estimated using (3.47a,b) as

(3.80)\begin{align} \frac{\mathcal{U}_{h{\dagger}}^{k}-\mathcal{U}^{k}}{\tau_k} &\leq \frac{1}{\tau_k}\left(\bar{\rho}^{k}(u_{\dagger}^{k} - u^{k}),u_{\dagger}^{k}\right)_u +\frac{1}{\tau_k}\left(\bar{\rho}^{k}(v_{\dagger}^{k} - v^{k}),v_{\dagger}^{k}\right)_v\nonumber\\ &={-} \sum_{l=1}^N(\hat{c}_l^k d_x^c\mu_{c,l}^{k+1},u_{\dagger}^{k})_u -\sum_{l=1}^N(\hat{c}_l^k d_y^c\mu_{c,l}^{k+1},v_{\dagger}^{k})_v. \end{align}

Using (3.51) and (3.52), we derive the difference between $\mathcal {U}_{h\star }^{k}$ and $\mathcal {U}_{\dagger} ^{k}$ as

(3.81)\begin{align} \frac{\mathcal{U}_{h\star}^{k}-\mathcal{U}_{h{\dagger}}^{k}}{\tau_k}&\leq \frac{1}{\tau_k}\left(\bar{\rho}^{k}(u_\star^{k} - u_{\dagger}^{k}),u_\star^{k}\right)_u +\frac{1}{\tau_k}\left(\bar{\rho}^{k}(v_\star^{k} - v_{\dagger}^{k},v_\star^{k}\right)_v\nonumber\\ &={-}\left(\hat{\phi}^k d_x^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})+d_x^cp^{k+1}-\hat{\rho}^k g_x,u_\star^{k}\right)_u\nonumber\\ &\quad -\left(\hat{\phi}^k d_y^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})+d_y^cp^{k+1}-\hat{\rho}^k g_y,v_\star^{k}\right)_v. \end{align}

Using (3.56), we deduce that

(3.82)\begin{align} (d_x^cp^{k+1},u_\star^{k})_u+(d_y^cp^{k+1},v_\star^{k})_v&={-}(p^{k+1},d_x^u u_\star^{k}+d_y^v v_\star^k)_c\nonumber\\ &=(\lambda p^{k+1}, d_x^u J_{\phi,x}^{k+1}+d_y^v J_{\phi,y}^{k+1})_c\nonumber\\ &={-}(\lambda d_x^cp^{k+1},J_{\phi,x}^{k+1})_u-(\lambda d_y^cp^{k+1},J_{\phi,y}^{k+1})_v . \end{align}

It follows from (3.81) and (3.82) that

(3.83)\begin{align} \frac{\mathcal{U}_{h\star}^{k}-\mathcal{U}_{h{\dagger}}^{k}}{\tau_k}&={-} \left(\hat{\phi}^k d_x^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1}) -\hat{\rho}^k g_x,u_\star^{k}\right)_u+(\lambda d_x^cp^{k+1},J_{\phi,x}^{k+1})_u\nonumber\\ &\quad - \left(\hat{\phi}^k d_y^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})-\hat{\rho}^k g_y,v_\star^{k}\right)_v+(\lambda d_y^cp^{k+1},J_{\phi,y}^{k+1})_v. \end{align}

The difference between $\mathcal {U}_h^{k+1}$ and $\mathcal {U}_{h\star }^k$ is estimated using (3.31) as

(3.84)\begin{align} \mathcal{U}_h^{k+1}-\mathcal{U}_{h\star}^{k}&=\tfrac{1}{2}( \bar{\rho}^{k+1}u^{k+1},u^{k+1})_u-\tfrac{1}{2}( \bar{\rho}^{k}u_\star^{k},u_\star^{k})_u\nonumber\\ &\quad +\tfrac{1}{2}( \bar{\rho}^{k+1}v^{k+1},v^{k+1})_v-\tfrac{1}{2}( \bar{\rho}^{k}v_\star^{k},v_\star^{k})_v\nonumber\\ &\leq ( \bar{\rho}^{k+1}u^{k+1}-\bar{\rho}^{k}u_\star^{k},u^{k+1})_u-\tfrac{1}{2}( \bar{\rho}^{k+1}- \bar{\rho}^{k},|u^{k+1}|^2)_u\nonumber\\ &\quad + ( \bar{\rho}^{k+1}v^{k+1}-\bar{\rho}^{k}v_\star^{k},v^{k+1})_v-\tfrac{1}{2}( \bar{\rho}^{k+1}- \bar{\rho}^{k},|v^{k+1}|^2)_v. \end{align}

Taking into account (3.61) and (3.62), we deduce that

(3.85)\begin{align} \frac{1}{\tau_k}\,( \bar{\rho}^{k+1}u^{k+1}-\bar{\rho}^{k}u_\star^{k},u^{k+1})_u&={-}\left( d_x^c(U_u^k\hat{u}^{k+1} )+d_y^p(V_u^k\hat{u}^{k+1}) , u^{k+1}\right)_u\nonumber\\ &\quad + ( d_x^c\varPsi^{k+1}+d_y^p\varTheta^{k+1},u^{k+1})_u\nonumber\\ &=( U_u^k\hat{u}^{k+1} , d_x^u u^{k+1})_c+( V_u^k\hat{u}^{k+1}, d_y^u u^{k+1})_{pu}\nonumber\\ &\quad - (\varPsi^{k+1}, d_x^u u^{k+1})_c-( \varTheta^{k+1}, d_y^u u^{k+1})_{pu}, \end{align}

(3.86)\begin{align} \frac{1}{\tau_k}\,( \bar{\rho}^{k+1}v^{k+1}-\bar{\rho}^{k}v_\star^{k},v^{k+1})_v&={-}\left(d_x^p(U_v^k\hat{v}^{k+1})+d_y^c(V_v^k \hat{v}^{k+1}) , v^{k+1}\right)_v\nonumber\\ &\quad + ( d_x^p\varTheta^{k+1}+d_y^c\varUpsilon^{k+1},v^{k+1})_v\nonumber\\ &=( U_v^k\hat{v}^{k+1}, d_x^v v^{k+1})_{pv}+(V_v^k \hat{v}^{k+1}, d_y^v v^{k+1})_c\nonumber\\ &\quad - (\varUpsilon^{k+1}, d_y^v v^{k+1})_c- (\varTheta^{k+1}, d_x^v v^{k+1})_{pv}. \end{align}

Taking into account two averaged forms of the discrete mass balance equations

(3.87a)\begin{gather} \frac{ \bar{\rho}^{k+1} -\bar{\rho}^{k} }{\tau_k }+ d_x^c U_u^k + d_y^p V_u^k =0, \end{gather}

(3.87b)\begin{gather}\frac{ \bar{\rho}^{k+1} -\bar{\rho}^{k} }{\tau_k }+d_x^p U_v^k + d_y^c V_v^k =0, \end{gather}

we deduce that

(3.88)\begin{align} \frac{1}{2\tau_k}\left( \bar{\rho}^{k+1}- \bar{\rho}^{k},|u^{k+1}|^2\right)_u&={-}\frac{1}{2}\left( d_x^c U_u^k + d_y^p V_u^k,|u^{k+1}|^2\right)_u\nonumber\\ &=( U_u^k\bar{u}^{k+1}, d_x^uu^{k+1})_c+( V_u^k\bar{u}^{k+1}, d_y^u u^{k+1})_{pu}, \end{align}

(3.89)\begin{align} \frac{1}{2\tau_k}\left( \bar{\rho}^{k+1}- \bar{\rho}^{k},|v^{k+1}|^2\right)_v&={-}\frac{1}{2}\left( d_x^p U_v^k + d_y^c V_v^k, |v^{k+1}|^2\right)_v\nonumber\\ &=( U_v^k\bar{v}^{k+1}, d_x^v v^{k+1})_{pv}+( V_v^k\bar{v}^{k+1}, d_y^v v^{k+1})_{c}, \end{align}

where $\bar {u}$ and $\bar {v}$ are the average values of $u$ and $v$ on the boundaries of their respective grid cells. Substituting (3.85)–(3.89) into (3.84) yields

(3.90)\begin{align} \frac{\mathcal{U}_h^{k+1}-\mathcal{U}_{h\star}^{k}}{\tau_k}&\leq \left( U_u^k(\hat{u}^{k+1}-\bar{u}^{k+1}) , d_x^u u^{k+1}\right)_c+\left( V_u^k(\hat{u}^{k+1}-\bar{u}^{k+1}), d_y^u u^{k+1}\right)_{pu}\nonumber\\ &\quad +\left( U_v^k(\hat{v}^{k+1}-\bar{v}^{k+1}), d_x^v v^{k+1}\right)_{pv}+\left(V_v^k (\hat{v}^{k+1}-\bar{v}^{k+1}), d_y^v v^{k+1}\right)_c\nonumber\\ &\quad - (\varPsi^{k+1}, d_x^u u^{k+1})_c- ( \varTheta^{k+1}, d_y^u u^{k+1})_{pu}\nonumber\\ &\quad - (\varUpsilon^{k+1}, d_y^v v^{k+1})_c- (\varTheta^{k+1}, d_x^v v^{k+1})_{pv}. \end{align}

In terms of the upwind rules, we can estimate the first four terms of (3.90) as

(3.91)\begin{gather} \left( U_u^k(\hat{u}^{k+1}-\bar{u}^{k+1}) , d_x^u u^{k+1}\right)_c={-} (h\,|U_u^k|\,d_x^uu^{k+1}, d_x^u u^{k+1})_c\leq0, \end{gather}

(3.92)\begin{gather}\left( V_u^k(\hat{u}^{k+1}-\bar{u}^{k+1}), d_y^u u^{k+1}\right)_{pu}={-}(h\,|V_u^k|\,d_y^uu^{k+1}, d_y^u u^{k+1})_{pu}\leq0, \end{gather}

(3.93)\begin{gather}\left( U_v^k(\hat{v}^{k+1}-\bar{v}^{k+1}), d_x^v v^{k+1}\right)_{pv}={-}(h\,|U_v^k|\,d_x^v v^{k+1}, d_x^v v^{k+1})_{pv}\leq0, \end{gather}

(3.94)\begin{gather}\left(V_v^k (\hat{v}^{k+1}-\bar{v}^{k+1}), d_y^v v^{k+1}\right)_c={-}(h\,|V_v^k|\,d_y^v v^{k+1}, d_y^v v^{k+1})_c\leq0. \end{gather}

As a consequence of (3.91)–(3.94), the inequality (3.90) can be simplified as

(3.95)\begin{align} \frac{\mathcal{U}_h^{k+1}-\mathcal{U}_{h\star}^{k}}{\tau_k}&\leq{-} (\varPsi^{k+1}, d_x^u u^{k+1})_c- ( \varTheta^{k+1}, d_y^u u^{k+1})_{pu}\nonumber\\ &\quad - (\varUpsilon^{k+1}, d_y^v v^{k+1})_c- (\varTheta^{k+1}, d_x^v v^{k+1})_{pv}\nonumber\\ &={-} 2(\eta^kd_x^uu^{k+1}, d_x^u u^{k+1})_c- 2(\eta^kd_y^vv^{k+1}, d_y^v v^{k+1})_c\nonumber\\ &\quad - \left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_y^u u^{k+1}\right)_{pu}\nonumber\\ &\quad - \left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_x^v v^{k+1}\right)_{pv}, \end{align}

where the last equality is obtained using (3.60a–c). Taking into account

(3.96)\begin{align} &\left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_y^u u^{k+1}\right)_{pu}+\left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_x^v v^{k+1}\right)_{pv}\nonumber\\ &\quad \geq \left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_y^uu^{k+1}+d_x^vv^{k+1}\right)_{puv}, \end{align}

we combine (3.80), (3.83) and (3.95) to obtain

(3.97)\begin{align} \frac{\mathcal{U}_h^{k+1}-\mathcal{U}_h^{k}}{\tau_k} &=\frac{\mathcal{U}_h^{k+1}-\mathcal{U}_{h\star}^{k}}{\tau_k} +\frac{\mathcal{U}_{h\star}^{k}-\mathcal{U}^{k}_{h{\dagger}}}{\tau_k} +\frac{\mathcal{U}_{h{\dagger}}^{k}-\mathcal{U}_h^{k}}{\tau_k}\nonumber\\ &\leq{-} \sum_{l=1}^N(\hat{c}_l^k d_x^c\mu_{c,l}^{k+1},u_{\dagger}^{k})_u - \sum_{l=1}^N(\hat{c}_l^k d_y^c\mu_{c,l}^{k+1},v_{\dagger}^{k})_v\nonumber\\ &\quad - \left(\hat{\phi}^k d_x^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1}) -\hat{\rho}^k g_x,u_\star^{k}\right)_u+(\lambda d_x^cp^{k+1},J_{\phi,x}^{k+1})_u\nonumber\\ &\quad - \left(\hat{\phi}^k d_y^c(\mu_{\phi}^{k+1}+\mu_{c\phi}^{k+1})-\hat{\rho}^k g_y,v_\star^{k}\right)_v+(\lambda d_y^cp^{k+1},J_{\phi,y}^{k+1})_v\nonumber\\ &\quad - 2(\eta^kd_x^uu^{k+1}, d_x^u u^{k+1})_c- 2(\eta^kd_y^vv^{k+1}, d_y^v v^{k+1})_c\nonumber\\ &\quad - \left(\bar{\eta}^k (d_y^uu^{k+1}+d_x^vv^{k+1}), d_y^uu^{k+1}+d_x^vv^{k+1}\right)_{puv}. \end{align}

Finally, we deduce (3.73) by combining (3.74), (3.77), (3.78) and (3.97).

4. Numerical results

In this section, numerical experiments are performed to verify and validate the proposed model and numerical scheme.

4.1. Example 1

In this example, we consider the single-solute problems with the Dirichlet boundary conditions in the one-dimensional domain $\varOmega =[0,1]$. The velocity is set to be zero, and the solute transport equation is simplified as

(4.1)\begin{equation} \frac{c^{k+1}- c^{k} }{ \tau_k} -\boldsymbol{\nabla}\boldsymbol{\cdot} c^{k}\,\boldsymbol{\nabla} \mu_c^{k+1}=0. \end{equation}

We take $\tau _k=0.01$. Let $\epsilon$ be the interface width. The phase variable is time-independent and given by the formula

(4.2)\begin{equation} \phi(x) = \frac{1}{2}\left(\tanh\left( \frac{x-0.5}{\sqrt{2}\epsilon}\right)+1\right),\quad x\in[0,1], \end{equation}

which is also illustrated in figure 1.

Figure 1. Phase variable in Example 1.

In the first case, we take $\alpha =\beta =1$, $\gamma =\ln (0.1)$ and $\delta =\ln (0.5)$. The boundary condition $\mu _c=0$ is imposed on both endpoints, while the homogeneous initial concentration $c^0=0.3$ is taken in the entire domain. In terms of the boundary condition, the solute concentrations in the bulk solvent fluids are calculated as

(4.3)\begin{gather} c_0=\exp(\delta)=0.5,\quad \phi=0, \end{gather}

(4.4)\begin{gather}c_1 =\exp(\gamma )=0.1,\quad \phi=1. \end{gather}

For $\epsilon =0.01$, figure 2 illustrates dynamics of the solute concentration with time, showing that the system converges rapidly to a steady state. Figure 3 gives a comparison of the numerical solutions at the equilibrium state with the exact solutions, as well as the numerical errors that are the absolute values of the difference between numerical and exact solutions. The results demonstrate that the numerical solutions are in perfect agreement with the exact solutions. The sharp-interface limit of the proposed model can be derived via formal asymptotic analysis as in the literature (e.g. Garcke, Lam & Stinner 2014; Qin et al. 2022), and it will be our ongoing work. Here, we show the sharp-interface limit by numerical results in figure 4, which clearly indicates that the concentrations will approach the sharp-interface limit as the interface width $\epsilon$ approaches zero.

Figure 2. Diffusion of the solute concentration at different times in Example 1.

Figure 3. Comparison of numerical results with the exact solution in Example 1: (a) profiles of numerical and exact solutions; (b) errors between numerical and exact solutions.

Figure 4. Concentration profiles with different values of $\epsilon$ in Example 1.

In the second case, we take the time-independent boundary conditions $\mu _c=0.5$ at $x=0$, and $\mu _c=0.1$ at $x=1$. The parameters are taken as $\alpha =1,2$, $\beta =2,1$, $\gamma = \ln (0.1)$ and $\delta = \ln (0.2)$. At the steady state of this case, we have

(4.5)\begin{equation} -\frac{{\rm d}}{{\rm d}\kern0.7pt x} \left(c\,\frac{{\rm d} \mu_c}{{\rm d}\kern0.7pt x}\right)=0, \end{equation}

which implies that there exists a constant $\iota >0$ such that

(4.6)\begin{equation} -c\,\frac{{\rm d} \mu_c}{{\rm d}\kern0.7pt x}=\iota. \end{equation}

In the two bulk fluids, i.e. $\phi =1$ and $\phi =0$, we have

(4.7)\begin{equation} \frac{{\rm d} c}{{\rm d}\kern0.7pt x}=\left\{\begin{array}{@{}ll@{}} -\dfrac{\iota}{\alpha}, & \phi=1,\\ -\dfrac{\iota}{\beta}, & \phi=0. \end{array}\right. \end{equation}

On the other hand, the solute concentrations within two solvent fluids are expressed as

(4.8)\begin{gather} c_0=\exp\left(\frac{0.5}{\beta}+\ln(0.2))\right),\quad \phi=0, \end{gather}

(4.9)\begin{gather}c_1 =\exp\left(\frac{0.1}{\alpha}+\ln(0.1) \right),\quad \phi=1. \end{gather}

Combining (4.7)–(4.9) yields

(4.10)\begin{equation} c(x)=\left\{\begin{array}{@{}ll@{}} -\dfrac{\iota}{\alpha}\,x+c_1, & \phi=1, \\ -\dfrac{\iota}{\beta}\,(x-1)+c_0, & \phi=0, \end{array}\right. \end{equation}

where $c_0$ and $c_1$ are given in (4.8) and (4.9). Figure 5 depicts numerical results of the steady-state concentrations, which accord with the above theoretical analysis. It also indicates that the model has the ability to reflect the heterogeneity in binary solvent fluids through different choices of parameters.

Figure 5. Concentration profiles with different $\alpha$ and $\beta$.

4.2. Example 2

In this example, we simulate the single-solute problems with the Neumann boundary conditions and without velocity fields in the one-dimensional domain $\varOmega =[0,1]$. The phase variable is given by the formula

(4.11)\begin{equation} \phi = \left\{\begin{array}{@{}ll@{}} \dfrac{1}{2}\left(1+ \tanh\left(\dfrac{x-0.25}{\sqrt{2}\epsilon}\right)\right), & x\in[0, 0.5] ,\\ \dfrac{1}{2}\left(1- \tanh\left(\dfrac{x-0.75}{\sqrt{2}\epsilon}\right)\right), & x\in(0.5, 1], \end{array}\right. \end{equation}

where $\epsilon$ is the interface width. Here, we take $\epsilon =0.01$ for all tested cases in this example. The phase variable is illustrated in figure 6. The discrete solute transport equation takes the form

(4.12)\begin{equation} \frac{c^{k+1}- c^{k} }{ \tau_k} -\frac{{\rm d}}{{\rm d}\kern0.7pt x} \left(c^{k}\,\frac{{\rm d} \mu_c^{k+1}}{{\rm d}\kern0.7pt x}\right)=0, \end{equation}

associated with the boundary condition ${\textrm {d} \mu _c}/{\textrm {d}\kern0.7pt x}=0$ on the endpoints. We take $\tau _k=0.001$. The total solute mass is conserved. The fluid system will reach an equilibrium state as the total free energy within this system is dissipated with time.

Figure 6. Phase variable in Example 2.

In the first case, we take $\alpha =\beta =1$, $\gamma =\ln (0.1)$ and $\delta =\ln (0.5)$. The initial concentration is $c^0=0.3$ in the entire domain. Figure 7 illustrates the dynamics of chemical potential and concentration with time, while the total energy dissipation is depicted in figure 8. The normal diffusion makes the solute spread homogeneously in space, and it usually occurs in the single-phase fluid. The active diffusion generally stems from different chemical and physical characteristics of two solvent fluids that play a major role in dissolving a specific solute. The underlying mechanics of such scenarios is described by the free energy and chemical potential. Although the initial concentration is uniform in space, the active diffusion still takes place due to the driving force resulting from the chemical potential gradient in the space occupied by two solvent fluids, and total energy will be dissipated in this process until an equilibrium state is reached. Let $\bar {\mu }_c$ be the equilibrium-state chemical potential, which is uniform and constant in space. For a specific phase variable $\phi$, the equilibrium-state solute concentration, denoted by $c_\phi$, depends on the phase variable and satisfies

(4.13)\begin{equation} \phi\alpha \left(\ln(c_\phi)-\gamma\right) +(1-\phi)\beta \left(\ln(c_\phi)-\delta\right)=\bar{\mu}_c. \end{equation}

Solving (4.13) gives the analytical formulation of $c_\phi$ as

(4.14)\begin{equation} c_\phi=\exp\left(\frac{\bar{\mu}_c+\phi\alpha\gamma + (1-\phi)\beta\delta}{\phi\alpha+(1-\phi)\beta}\right). \end{equation}

Here, $\bar {\mu }_c$ can be calculated analytically through solving the total mass conservation equation

(4.15)\begin{equation} \int_\varOmega c_\phi\,{{\rm d}\kern0.7pt x}=\int_\varOmega c^0\,{\rm d}\kern0.7pt x. \end{equation}

In figure 9, we compare the numerical and exact solutions of the solute concentration at the equilibrium state as well as showing the numerical errors (i.e. the absolute values of the difference between numerical and exact solutions). Numerical results are in agreement with the exact solutions.

Figure 7. Dynamics of the solute chemical potential and concentration at different times in Example 2: (a) solute chemical potential; (b) solute concentration.

Figure 8. Energy dissipation curves in Example 2.

Figure 9. Comparison of numerical results with the exact solution in Example 2: (a) numerical and exact solutions; (b) numerical errors.

In the second test, we take different values of parameters $\alpha$, $\beta$, $\gamma$ and $\delta$. The initial concentration is still set to be $c^0=0.3$ in the entire domain. Figure 10 depicts the profiles of the solute concentrations at the equilibrium states with different parameters. Equation (4.14) shows that the solute concentrations at different equilibrium states can be achieved through the appropriate choices of the four parameters. This conclusion is verified clearly by figure 10. The proposed free energy and chemical potential have the capability of describing different active diffusion processes.

Figure 10. Profiles of the solute concentrations at the equilibrium states in Example 2: (a) cases with $\gamma =\delta =0$ and different values of $\alpha$ and $\beta$; (b) cases with $\alpha =\beta =1$ and different values of $\gamma$ and $\delta$.

4.3. Example 3

In this example, we simulate the dynamical process of a single solute in binary solvent fluids to show the convergence of the proposed scheme. The spatial domain is the square domain $\varOmega =[0,1]^2$, which is divided by a uniform square mesh with $50\times 50$ grid cells. The parameters in the model are chosen as $M_\phi =1$, $Re = 10^2$, $Cn = 0.03$, $We = 1$, $Pe_\phi = 10^3$, $\varrho _1/\varrho _2 = 10$ and $\eta _1/\eta _2 = 10$. The gravity effect is neglected in this example. The initial phase variable is given as

(4.16)\begin{equation} \phi^0(x,y)=\left\{\begin{array}{@{}ll@{}} \phi^*, & a\leq x\leq1-a,\ a\leq y\leq1-a, \\ \phi_*, & {\rm elsewhere}, \end{array}\right. \end{equation}

where $a=0.37$, $\phi _*=0.0735$ and $\phi ^*=0.9335$. The initial phase variable is also illustrated in figure 11. The initial solute concentration is $c^0=0.2$ in the entire domain. The reference velocity parameter is chosen as $\varsigma =1$. For the solute transport equation, we take $Pe_c=1$, $D_{1,f}=1$, $\alpha =\beta =1$, $\gamma =\ln (0.5)$ and $\delta =\ln (0.1)$. The initial velocity is zero. The total simulation time is $t_f=0.1$. Contours of the phase variable and solute concentration at $t=0.1$ are shown in figure 12. From figures 11 and 12, we can see that the interfacial tension between two bulk fluids makes the square droplet shrink into a circle, and meanwhile, a large amount of solute transfers from the surrounding into the droplet due to the active diffusion.

Figure 11. Initial phase variable in Example 3.

Figure 12. Contours of the phase variable and solute concentration at $t=0.1$ in Example 3: (a) phase variable; (b) solute concentration.

Due to strong nonlinearity and complication of the model, we cannot derive the analytical solutions, thus we use the solutions computed with a very small time step size $\tau _k=2^{-4}s$, $s=10^{-3}$, as the approximate analytical solutions, which are denoted by $\tilde {c}$, $\tilde {\phi }$, $\tilde {\boldsymbol {u}}$ and $\tilde {\mathcal {E}}$. We define the discrete norms for the errors $e_c$, $e_\phi$, $e_{\boldsymbol {u}}$ and $e_{\mathcal {E}}$ as follows:

(4.17)\begin{gather} \|e_c\|_h=\max_{k}\left(h^2\sum_{i=0}^{n_x-1}\sum_{j=0}^{n_y-1}|c^k_{i+{1}/{2},j+{1}/{2}} -\tilde{c}_{i+{1}/{2},j+{1}/{2}}(t_k)|^2\right)^{1/2}, \end{gather}

(4.18)\begin{gather}\|e_\phi\|_h=\max_{k}\left(h^2\sum_{i=0}^{n_x-1}\sum_{j=0}^{n_y-1}|\phi^k_{i+{1}/{2},j+{1}/{2}} -\tilde{\phi}_{i+{1}/{2},j+{1}/{2}}(t_k)|^2\right)^{1/2}, \end{gather}

(4.19)\begin{gather}\|e_{\boldsymbol{u}}\|_h=\max_{k}\left(h^2\sum_{i=1}^{n_x-1}\sum_{j=0}^{n_y-1}|u_{i,j+{1}/{2}}^k -\tilde{u}_{i,j+{1}/{2}}(t_k)|^2 \right.\nonumber\\ \hspace{42pt}\left.{}+h^2\sum_{i=0}^{n_x-1}\sum_{j=1}^{n_y-1}|v_{i+{1}/{2},j}^k-\tilde{v}_{i+{1}/{2},j}(t_k)|^2\right)^{1/2}, \end{gather}

(4.20)\begin{gather}\|e_{\mathcal{E}}\|_h= \max_{k}|\mathcal{E}_h^k-\tilde{\mathcal{E}}(t_k)|. \end{gather}

Numerical errors and convergence rates are displayed in table 1. It is evidently shown that the proposed scheme has first-order convergence in time.

Table 1. Temporal errors and convergence rates of the scheme.

4.4. Example 4

In this example, we simulate the dynamical process of two solutes in binary solvent fluids with a large density contrast. The spatial domain is $\varOmega =[0,1]^2$. The parameters in the model are chosen as $M_\phi =1$, $Re = 10^2$, $Cn = 0.02$, $We = 1$, $Pe_\phi = 10^3$, $\varrho _1/\varrho _2 = 100$ and $\eta _1/\eta _2 = 100$. The gravity effect is not considered in this example. The initial solute concentrations are uniform in space with $c_1^0=c_2^0=0.2$, while the initial phase variable is given as

(4.21)\begin{gather} \phi^0(x,y)=\min\left(\max\left(\hat{\phi}(x,y),\phi_*\right),\phi^*\right), \end{gather}

(4.22)\begin{gather}\hat{\phi}(x,y)=\max\left(\left( 1+\tanh\left( (0.15-r_1)/Cn \right) \right)/2, \left( 1+\tanh\left( (0.15-r_2)/Cn \right) \right)/2\right), \end{gather}

(4.23)\begin{gather}r_1 = \sqrt{ (x-0.4)^2 + (y-0.6)^2 },~~ r_2 = \sqrt{ (x-0.6)^2 + (y-0.4)^2 }, \end{gather}

where $\phi _*=0.0735$ and $\phi ^*=0.9335$. The solvent mixture is composed of two fluids with a large density contrast. There is an irregular droplet consisting mostly of the heavier fluid in the domain at the initial time. The reference velocity parameter is chosen as $\varsigma =10$. The initial velocity is zero. The solute mixture is composed of two different substances. For the solute transport equation, we take $Pe_c=1$ and $D_{1,f}=D_{2,f}=D_{1,2}=1$. The solute diffusion matrix is calculated by (2.37) and (2.38). For the solute energy parameters, we take $\alpha _1=\alpha _2=1$, $\beta _1=\beta _2=1$, $\gamma _1=\ln (0.5)$, $\gamma _2=\ln (0.1)$, $\delta _1=\ln (0.1)$ and $\delta _2=\ln (0.5)$.

The proposed scheme allows us to use non-uniform time step sizes. Here, we employ the adaptive time-stepping strategy (Qiao, Zhang & Tang 2011) to adjust the time step sizes

(4.24a,b)\begin{equation} \tau_k=\max\left(\tau_{min},\frac{\tau_{max}}{\sqrt{1+rw_k^2}}\right),\quad w_k=\frac{\mathcal{E}_h^k-\mathcal{E}_h^{k-1}}{\tau_{k-1}}, \end{equation}

where $\mathcal {E}_h^k$ is the total free energy, $\tau _{min}$ and $\tau _{max}$ are the preset lower and upper bounds of the time step sizes, and $r$ is a positive constant. In this example, we take $\tau _{min}=5\times 10^{-4}$, $\tau _{max}=5\times 10^{-3}$ and $r=10$. Figure 13 depicts the adaptive time step sizes calculated by (4.24a,b).

Figure 13. Adaptive time step sizes in Example 4.

The phase variables at different times are illustrated in figure 14, while the solute concentrations are shown in figures 15 and 16. The results show that driven by the interfacial tension between two bulk fluids, the droplets merge together and then shrink into a circle, and meanwhile, the active diffusion as well as the fluid flow makes a large amount of solute 1 transfer to the droplets from the surrounding fluid, whereas solute 2 is moving from the droplets to the surrounding fluid.

Figure 14. Phase variable contours at different times in Example 4: (a) $k=0$, (b) $k=20$, (c) $k=50$, (d) $k=150$, and (e) $k=400$.

Figure 15. Dynamics of the concentration of solute 1 in Example 4: (a) $k=20$, (b) $k=30$, (c) $k=50$, (d) $k=100$, (e) $k=150$, and ( f) $k=400$.

Figure 16. Dynamics of the concentration of solute 2 in Example 4: (a) $k=20$, (b) $k=30$, (c) $k=50$, (d) $k=100$, (e) $k=150$, and ( f) $k=400$.

The normalized total free energy is defined as

(4.25)\begin{equation} \mathcal{E}_r^k= \frac{\mathcal{E}_h^k}{|\mathcal{E}_h^0|}. \end{equation}

Figure 17 depicts the normalized total energy dissipation curve, showing that total energy is always decreasing with time. The normalized total solvent mass is defined as

(4.26)\begin{equation} {\rho}_r^k = \frac{(\rho^{k},1)_c-(\rho^{0},1)_c}{(\rho^{0},1)_c}. \end{equation}

The normalized total solute concentrations are defined as

(4.27)\begin{equation} c_{r,l}^k = \frac{(c_l^{k},1)_c-(c_l^{0},1)_c}{(c_l^{0},1)_c},\quad l=1,2. \end{equation}

Figure 18 illustrates the profiles of ${\rho }_r^k$ and $c_{r,l}^k$, and shows that the proposed method is capable of conserving total masses of solutes and solvent fluids except for roundoff errors.

Figure 17. Energy dissipation profile in Example 4.

Figure 18. Mass conservation profiles of solvents and solutes in Example 4: (a) solvent mass conservation; (b) solute mass conservation.

4.5. Example 5

This example takes into consideration the gravity effect on the solute transport in two different droplets. The spatial domain is still $\varOmega =[0,1]^2$. The initial phase variable is given as

(4.28)\begin{gather} \phi^0(x,y)=\min\left(\max\left(\hat{\phi}(x,y),\phi_*\right),\phi^*\right), \end{gather}

(4.29)\begin{gather}\hat{\phi}(x,y)=\max\left(\left( 1+\tanh\left( (0.15-r_1)/Cn \right) \right)/2, \left( 1+\tanh\left( (0.12-r_2)/Cn \right) \right)/2\right), \end{gather}

(4.30)\begin{gather}r_1 = \sqrt{ (x-0.5)^2 + (y-0.83)^2 },~~ r_2 = \sqrt{ (x-0.5)^2 + (y-0.5)^2 }, \end{gather}

where $\phi _*=0.045$ and $\phi ^*= 0.936$. The initial velocity is zero. The reference velocity parameter is chosen as $\varsigma =1$. The parameters in the fluid flow equations are chosen as $M_\phi =1$, $Re = 10^3$, $Cn = 0.02$, $We = 1$, $Pe_\phi = 10^3$, $\varrho _1/\varrho _2 = 10^3$, $\eta _1/\eta _2 = 10^3$ and $Fr=\frac {1}{5}$. In this example, we want to demonstrate the effects of solutes on the dynamics of solvents, so we simulate two cases with different single solutes. For the solute energy parameters, we take $\alpha =\beta =50$, $\gamma =\ln (0.5)$, $\delta =\ln (0.1)$ in case 1, and $\alpha =\beta =50$, $\gamma =\ln (0.1)$, $\delta =\ln (0.5)$ in case 2. The solute diffusion parameters are taken as $Pe_c=10^2$ and $D_{1,f}=1$ for both cases. The initial solute concentrations in both cases are given as

(4.31)\begin{equation} c^0(x,y)=\left\{\begin{array}{@{}ll@{}} 0.5, & y\leq\dfrac{1}{3}, \\ 0.1, & y>\dfrac{1}{3}. \end{array}\right. \end{equation}

To compare the results of the two cases, we use the uniform time step size $\tau _k=2\times 10^{-3}$ for both cases.

The phase variables of the two cases at different times are shown in figures 19 and 21. The solute concentrations of the two cases are illustrated in figures 20 and 22. In both cases, the motion of two droplets is driven by the interfacial tension, gravity and the solute chemical potentials. The droplets are moving downwards due to the gravity effect and large density contrast, and meanwhile, they are merging into a large droplet due to the interfacial tension. During the dynamical process of case 1, the solute concentration in the droplets is increasing with time due to the active diffusion, whereas in case 2, the solute behaves in the opposite pattern. It is observed clearly from figures 19 and 21 that the solutes have a significant effect on the motion of solvents; actually, in case 1, the solute dissolving in the droplets drags the downward motion of the droplets against gravity.

Figure 19. Phase variable contours of case 1 at different times in Example 5: (a) $k=0$, (b) $k=100$, (c) $k=200$, (d) $k=300$, (e) $k=400$, and ( f) $k=500$.

Figure 20. Solute concentrations of case 1 at different times in Example 5: (a) $k=0$, (b) $k=100$, (c) $k=200$, (d) $k=300$, (e) $k=400$, and ( f) $k=500$.

Figure 21. Phase variable contours of case 2 at different times in Example 5: (a) $k=0$, (b) $k=100$, (c) $k=200$, (d) $k=300$, (e) $k=400$, and ( f) $k=500$.

Figure 22. Solute concentrations of case 2 at different times in Example 5: (a) $k=0$, (b) $k=100$, (c) $k=200$, (d) $k=300$, (e) $k=400$, and ( f) $k=500$.

In this example, the gravitational potential energy is included in the total free energy. In figure 23, we plot the dissipation curves of the normalized total energy $\mathcal {E}_r^k$ defined in (4.25), which show that the total free energies of both cases are always dissipated with time. Figures 24 and 25 illustrates the profiles of ${\rho }_r^k$ and $c_{r}^k$ defined in (4.26) and (4.27). The mass loss always falls within roundoff errors. Therefore, the proposed method performs very well in not only energy stability but also mass conservation.

Figure 23. Energy dissipation profiles in Example 5: (a) case 1; (b) case 2.

Figure 24. Solvent mass conservation profiles in Example 5: (a) case 1; (b) case 2.

Figure 25. Solute mass conservation profiles in Example 5: (a) case 1; (b) case 2.

5. Conclusions

A thermodynamically consistent phase-field model has been proposed to describe the activated solute transport in binary solvent fluids. A key ingredient of the model is to introduce the mixed free-energy function for the multiple solutes, which is able to characterize the solubility difference of solutes in two solvent fluids through different solute chemical potentials in binary fluids. A general multi-component solute diffusion model is established using the Maxwell–Stefan approach, which takes into account the crossing influences between different solutes. It remains to support the free-energy function through molecular scale simulations or experimental data in the future. The hydrodynamics of two solvent fluids is described by a general phase-field model, which accounts for general average velocity and different densities. The model is derived rigorously using the second law of thermodynamics. The model system is highly nonlinear and strongly coupled. To simulate the model efficiently, we have proposed a linearized, decoupled and energy-stable numerical method. The spatial discretization is constructed using the finite difference/volume methods on staggered grids with the upwind strategy. Both semi-discrete and fully discrete schemes are proved to follow an energy dissipation law at the discrete level as well as ensuring the mass conservation law for solutes and solvents. Numerical experiments have been performed to demonstrate that the proposed model and numerical method can simulate different processes of the solute active diffusion in two-phase solvent fluids.

The proposed model has potential applications in many relevant fields, such as geological carbon sequestration with active diffusion due to the heterogeneity in underground formations, intracellular ion mass transfer, and the removal and extraction of specific substances in industrial and environmental technologies. The modified Maxwell–Stefan model can be applicable to general multi-component diffusion processes of free fluid flow and porous media flow. The proposed numerical method can provide guidance for designing efficient numerical methods for the other multiphase flow problems, for instance, the electrical ion transport in two-phase fluids coupled with electromagnetic fields.