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 (Reference Patankar2016) established the equations for the solute equilibrium and the phase equilibrium. The logarithmic chemical potentials were introduced by Patankar (Reference Patankar2016) 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 Reference Patankar2016). 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 Reference Wells, Kuhl and Garikipati2006; Wang, Kou & Cai Reference Wang, Kou and Cai2020b). 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 (Reference Brenner2009a,Reference Brennerb, Reference Brenner2013), 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 Reference Brenner2009a,Reference Brennerb, Reference Brenner2013). 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 Reference Kou and Sun2018a,Reference Kou and Sunb) 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. Reference Kou, Wang, Zeng and Cai2020b). 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 Reference Wesselingh and Krishna2000), 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. Reference Qin, Huang, Zhu, Liu and Xu2022) 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.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 Reference Kou, Sun and Wang2020a; Wang et al. Reference Wang, Kou and Cai2020b; Wang, Kou & Gao Reference Wang, Kou and Gao2021) 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. Reference Wang, Kou and Cai2020b) 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. (Reference Wang, Kou and Gao2021). 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 Reference Tryggvason, Scardovelli and Zaleski2011; Kou, Sun & Wang Reference Kou, Sun and Wang2018) 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. Reference Kou, Sun and Wang2018):

(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.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 Reference Garcke, Lam and Stinner2014; Qin et al. Reference Qin, Huang, Zhu, Liu and Xu2022), 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 Reference Qiao, Zhang and Tang2011) 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.