2.1. Phase variable and variable density

Suppose in a control volume $V$ we have two fluids (labelled by $i= 1,2$) that fill volumes $V_{i}$ separately. The volume fraction for the $i$th fluid can be given by $c_{i}= {V_{i}}/{V}$. Further we assume that the two fluids can mix within an interfacial region separating the two fluids and the volume occupied by a given amount of mass of the single fluid does not change after mixing. Therefore $c_1+c_2=1$ within the control volume $V$. Given the total mass and volume for the mixture, $M=M_{1}+M_{2}$ and $V=V_{1}+V_{2}$, we introduce the actual local mass density $\rho _{i}$ for each fluid, and the local volume-averaged mass density $\tilde {\rho }_{i}$ taken over a sufficient small volume $V$:

(2.1a,b)\begin{equation} \rho_{i}=\frac{M_{i}}{V_{i}},\quad \tilde{\rho}_{i}=\frac{M_{i}}{V}=\frac{M_{i}}{V_{i}}\frac{V_{i}}{V}=\rho_{i}c_{i}, \end{equation}

where $\rho _i$ is a positive constant. We then define the volume-averaged density of the mixture as

(2.2)\begin{equation} \rho=\tilde{\rho}_{1}+\tilde{\rho}_{2}=\rho_{1}c_{1}+\rho_{2}c_{2}. \end{equation}

Since $c_{1}+c_{2}=1$, we assume that the phase variable with respect to volume fraction is given by $c=c_{1}=1-c_{2}$. The corresponding variable density for the mixture can be rewritten as

(2.3)\begin{equation} \rho=(\rho_{1}-\rho_{2})c+\rho_{2}. \end{equation}

Suppose that the two fluids move with different velocities $\boldsymbol {u}_i(\boldsymbol {x},t)$, then mass conservation for the mixture is given by

(2.4)\begin{equation} \partial_{t}(\tilde{\rho}_{1}+\tilde{\rho}_{2})+\boldsymbol{\nabla}\boldsymbol{\cdot} (\tilde{\rho}_{1}\boldsymbol{u}_{1}+\tilde{\rho}_{2}\boldsymbol{u}_{2})=0. \end{equation}

The mass-averaged velocity $\boldsymbol {u} := {(\tilde {\rho }_{1}\boldsymbol {u}_{1}+\tilde {\rho }_{2}\boldsymbol {u}_{2})}/{\rho }$ and $\rho$, defined in (2.2), satisfy the mass conservation equation:

(2.5)\begin{equation} \partial_{t}\rho+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u})=0. \end{equation}

2.2. Balance laws

We next consider the balance laws for the two-fluid mixture. Suppose the two-phase fluid is combined in a domain $\varOmega$. We then take an arbitrary control volume $V(t)\subset \varOmega$ that moves with the mixture. Within the control volume, we define the following quantities for the two-phase fluid:

(2.6)\begin{gather} M=\int_{V(t)}\rho\,{\mathrm d}V, \end{gather}

(2.7)\begin{gather}{\boldsymbol{P}}=\int_{V(t)}\rho\boldsymbol{u}\,{\mathrm d}V, \end{gather}

(2.8)\begin{gather}E=\int_{V(t)}\left( \frac{1}{2}\rho |\boldsymbol{u}|^2+\tilde{\sigma} F(c)+ \frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}-\rho g \boldsymbol{z}\right){\mathrm d}V, \end{gather}

(2.9)\begin{gather}C=\int_{V(t)}c\,{\mathrm d}V, \end{gather}

where $M$, $\boldsymbol {P}$, $E$ and $C$ are the total mass, momentum, free energy and constituent volume of fluid 1 within the control volume $V(t)$. Further, $\rho (c)$ is the variable density of the mixture, $\boldsymbol {u}$ is the mass-averaged velocity of the mixture, $|\boldsymbol {u}|^2/2$ is the kinetic energy per unit mass, $g$ is the gravitational constant, $\boldsymbol {z}$ is the vector in the vertical direction, $\tilde {\sigma } F(c)+{\tilde {\sigma }\epsilon _c |\boldsymbol {\nabla } c|^2}/{2}$ is the Cahn–Hilliard-type interfacial free energy, $F(c)=c^2(c-1)^2/4\epsilon _c$ is double-well potential, $\epsilon _c$ is a small parameter that characterizes the thickness of the diffuse interface separating the fluid components and $\tilde {\sigma }$ is the surface tension that is related to the physical surface tension $\sigma$ through $\sigma =6\sqrt {2}\tilde {\sigma }$. The general forms of the physical balance laws associated with $M$, $\boldsymbol {P}$, $E$ and $C$ can be given as follows:

(2.10)\begin{gather} \frac{{\mathrm{d}} M}{\mathrm{d}t}=0, \end{gather}

(2.11)\begin{gather}\frac{{\mathrm{d}} \boldsymbol{P}} {\mathrm{d} t}= \int_{\partial V(t)} {\boldsymbol{\mathsf{m}}}\boldsymbol{\cdot} \boldsymbol{n}\, \mathrm{d}A -\int_{V(t)} \rho g \boldsymbol{z}\,{\mathrm{d}} V, \end{gather}

(2.12)\begin{gather}\frac{{\mathrm{d}} E}{{\mathrm{d}}t}= -\int_{\partial V(t)} {\boldsymbol{\mathsf{q}}}_{E}\boldsymbol{\cdot} \boldsymbol{n}\, {\mathrm{d}}A+\int_{V(t)} D\,{\mathrm{d}}V\quad (D\leqslant 0), \end{gather}

(2.13)\begin{gather}\frac{{\mathrm{d}} C}{{\mathrm{d}}t}= -\int_{\partial V(t)} {\boldsymbol{\mathsf{q}}}_{C}\boldsymbol{\cdot} \boldsymbol{n}\,{\mathrm{d}}A, \end{gather}

where ${\boldsymbol{\mathsf{m}}}$ is a tensor, ${\boldsymbol{\mathsf{q}}}_{E}$ and ${\boldsymbol{\mathsf{q}}}_{C}$ are fluxes and $D$ is the dissipation function, which will be specified later. The mass balance (2.10) ensures (2.5) is satisfied. The momentum balance (2.11) leads to

(2.14)\begin{equation} \rho \frac{{\mathrm{D}}\boldsymbol{u} }{{\mathrm{D}}t}={\boldsymbol{\nabla}} \boldsymbol{\cdot} {\boldsymbol{\mathsf{m}}}-\rho g \boldsymbol{z},\end{equation}

where ${{\mathrm {D}}}/{{\mathrm {D}}t}=\partial _{t}+\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla }$. The constituent volume balance equation (2.13) gives that

(2.15)\begin{equation} \frac{{\mathrm D}c}{{\mathrm D} t}+c(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} )=-{\boldsymbol{\nabla}}\boldsymbol{\cdot} {\boldsymbol{\mathsf{q}}}_{C}. \end{equation}

Further, by using the definition of $\rho$ in (2.3) and the mass conservation (2.5), we obtain

(2.16)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=- \frac{1}{\rho}\frac{{\mathrm D}\rho}{{\mathrm D}t}=- \frac{{\mathrm d}\rho}{\rho\,{\mathrm d}c}\frac{{\mathrm D}c}{{\mathrm D}t}= \frac{\rho_{2}-\rho_{1}}{\rho}\frac{{\mathrm D}c}{{\mathrm D}t}. \end{equation}

Substituting ${{\mathrm D}c}/{{\mathrm D}t}$ in (2.15), we obtain the quasi-incompressibility in the interfacial region for the two-phase flows:

(2.17)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=-\alpha\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\mathsf{q}}}_{C},\end{equation}

where $\alpha =(\rho _{2}-\rho _{1})/\rho _{2}$.

We now consider the energy balance (2.12). We first take the time derivative of the free energy (2.8) to obtain

(2.18)\begin{align} \frac{{\mathrm d} E}{{\mathrm d}t}&= \int_{V(t)} \left\{\rho \boldsymbol{u}\boldsymbol{\cdot} \frac{{\mathrm D}\boldsymbol{u}}{{\mathrm D}t}+\left(\tilde{\sigma} F'(c)- \tilde{\sigma}\epsilon_c {\rm \Delta} c\right)\frac{{\mathrm D}c}{{\mathrm D}t}+{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left(\tilde{\sigma}\epsilon_c\frac{{\mathrm D}c}{{\mathrm D}t}{ \boldsymbol{\nabla}}c\right)\right.\nonumber\\ &\quad \left.-\tilde{\sigma}\epsilon_c \left({\boldsymbol{\nabla}} c \otimes {\boldsymbol\nabla} c\right): {\boldsymbol\nabla} \boldsymbol{u}+ \left(\tilde{\sigma} F(c)+\frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}\right)\left( \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \right)-\rho g \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{z}\right\}{\mathrm d}V, \end{align}

where we have used the following identities:

(2.19)\begin{align} \frac{\mathrm D}{{\mathrm D}t}\left(F\left(c\right)+\frac{\epsilon_c|{\boldsymbol\nabla} c|^2}{2}\right)&=F'\left(c\right)\frac{{\mathrm D}c}{{\mathrm D}t}-{\rm \Delta} c\frac{{\mathrm D}c}{{\mathrm D}t}+{\boldsymbol\nabla}\boldsymbol{\cdot}\left({\boldsymbol\nabla} c\frac{{\mathrm D}c}{{\mathrm D}t}\right) \nonumber\\ &\quad -\left({\boldsymbol\nabla} c\otimes {\boldsymbol\nabla} c\right): {\boldsymbol\nabla} \boldsymbol{u}, \end{align}

(2.20)\begin{equation} {\boldsymbol\nabla} c \boldsymbol{\cdot} \frac{\mathrm D}{{\mathrm D}t}{\boldsymbol\nabla} c= {\boldsymbol\nabla} c \boldsymbol{\cdot}\left({\boldsymbol\nabla} c_{t} + \boldsymbol{u}\boldsymbol{\cdot{}\boldsymbol\nabla}\left({\boldsymbol\nabla} c\right)\right), \end{equation}

(2.21)\begin{equation}{\boldsymbol\nabla} c \boldsymbol{\cdot} {\boldsymbol\nabla} c_{t} = {\boldsymbol\nabla}\boldsymbol{\cdot} \left({\boldsymbol\nabla} c c_{t} \right)-{\rm \Delta} c c_{t}, \end{equation}

(2.22)\begin{align} {\boldsymbol\nabla} c \boldsymbol{\cdot} \left( \boldsymbol{u}\boldsymbol{\cdot}{\boldsymbol\nabla} \left({\boldsymbol\nabla} c\right)\right)&={\boldsymbol\nabla}\boldsymbol{\cdot} \left({\boldsymbol\nabla} c \left( \boldsymbol{u}\boldsymbol{\cdot}{\boldsymbol\nabla} c\right)\right)- {\rm \Delta} c \left(\boldsymbol{u}\boldsymbol{\cdot}{\boldsymbol\nabla} c\right)\nonumber\\ &\quad -\left({\boldsymbol\nabla} c \otimes {\boldsymbol\nabla} c\right): {\boldsymbol\nabla}\boldsymbol{u}. \end{align}

Substituting $\rho ({{\mathrm D}\boldsymbol {u}}/{{\mathrm D}t})$ and ${{\mathrm D}c}/{{\mathrm D}t}$ in (2.14) and (2.15), respectively, and applying the divergence theorem, we obtain

(2.23)\begin{align} \frac{{\mathrm D}E}{{\mathrm D}t}&= \int_{V(t)} \left\{{\boldsymbol\nabla} \boldsymbol{\cdot}\left( {\boldsymbol{\mathsf{m}}} \cdot \boldsymbol{u} +\tilde{\sigma}\epsilon_c\frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)- \left({\boldsymbol{\mathsf{m}}}_0 +\tilde{\sigma}\epsilon_c ({\boldsymbol\nabla} c \otimes {\boldsymbol\nabla} c)+{\tau}\right): {\boldsymbol \nabla}\boldsymbol{u}\right. \nonumber\\ &\quad - \left(-\tilde{\sigma} F(c)-\frac{\tilde{\sigma}\epsilon_c | \boldsymbol{\nabla} c|^2}{2}+\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c {\rm \Delta} c\right)c\right) {\boldsymbol\nabla}\boldsymbol{\cdot} \boldsymbol{u}\nonumber\\ &\quad \left.\vphantom{\left( {\boldsymbol{\mathsf{m}}} \cdot \boldsymbol{u} +\tilde{\sigma}\epsilon_c\frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)}- {\boldsymbol\nabla} \boldsymbol{\cdot}\left( \tilde{\sigma}\left(F'(c)- \epsilon_c{\rm \Delta} c\right) {\boldsymbol{\mathsf{q}}}_{C}\right)+{\boldsymbol{\mathsf{q}}}_{C} \boldsymbol{\cdot} {\boldsymbol\nabla} \left(\tilde{\sigma}\left(F'(c)-\epsilon_c{\rm \Delta} c\right) \right)\right\}{\mathrm d}V\nonumber\\ &=\int_{V(t)} \left\{{\boldsymbol\nabla} \boldsymbol{\cdot}\left( {\boldsymbol{\mathsf{m}}} \boldsymbol{\cdot} \boldsymbol{u} +\tilde{\sigma}\epsilon_c \frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)- \left({\boldsymbol{\mathsf{m}}}_0 +\tilde{\sigma}\epsilon_c ({\boldsymbol\nabla} c \otimes {\boldsymbol\nabla} c) \right): {\boldsymbol D}\boldsymbol{u}-{\tau}:{\boldsymbol\nabla} \boldsymbol{u} \right. \nonumber\\ &\quad - \left({\boldsymbol{\mathsf{m}}}_0 +\tilde{\sigma} \epsilon_c ({\boldsymbol\nabla} c \otimes {\boldsymbol\nabla} c) -\tilde{\sigma} F(c) {\boldsymbol{\mathsf{I}}}-\frac{\tilde{\sigma}\epsilon_c | \boldsymbol{\nabla} c|^2}{2}{{\boldsymbol{\mathsf{I}}}}\!+\! \left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c {\boldsymbol{\mathsf{I}}}\right): \frac{{\boldsymbol{\mathsf{I}}}}{3}{\boldsymbol\nabla} \boldsymbol{\cdot} \boldsymbol{u} \nonumber\\ &\quad\left. \vphantom{\left( {\boldsymbol{\mathsf{m}}} \boldsymbol{\cdot} \boldsymbol{u} +\tilde{\sigma}\epsilon_c \frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)}-{\boldsymbol\nabla} \boldsymbol{\cdot}\left(\tilde{\sigma}\left(F'(c)- \epsilon_c{\rm \Delta} c\right) {\boldsymbol{\mathsf{q}}}_{C}\right)+ {\boldsymbol{\mathsf{q}}}_{C}\boldsymbol{\cdot}{\boldsymbol\nabla} \left(\tilde{\sigma}\left(F'(c)-\epsilon_c{\rm \Delta} c\right)\right)\right\}{\mathrm d}V, \end{align}

where we denote the general stress tensor by ${\boldsymbol{\mathsf{m}}} ={\boldsymbol{\mathsf{m}}}_0 + {\tau }$, in which ${\tau }$ is the deviatoric stress tensor with zero trace, ${\boldsymbol{\mathsf{m}}}_0$ is the unknown part to be defined later and ${\boldsymbol{\mathsf{D}}} \boldsymbol {u}=\boldsymbol {\nabla } \boldsymbol {u}-({\boldsymbol \nabla } \boldsymbol {\cdot } \boldsymbol {u}){\boldsymbol{\mathsf{I}}} /3$ is the deviatoric part of ${\boldsymbol \nabla } \boldsymbol {u}$. Moreover, using the following identity and (2.17):

(2.24)\begin{align} &-\left({\boldsymbol{\mathsf{m}}}_0 +\tilde{\sigma}\epsilon_c ({\boldsymbol\nabla} c \otimes {\boldsymbol\nabla} c) -\tilde{\sigma} F(c){\boldsymbol{\mathsf{I}}}-\frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}{\boldsymbol{\mathsf{I}}}+\left(\tilde{\sigma} F'(c)\epsilon_c-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c {\boldsymbol{\mathsf{I}}}\right): \frac{{\boldsymbol{\mathsf{I}}}}{3}{\boldsymbol\nabla}\cdot \boldsymbol{u} \nonumber\\ &\quad =-\frac{1}{3}\left({\mathrm{tr}}\,{\boldsymbol{\mathsf{m}}}_0 - \frac{\tilde{\sigma} \epsilon_c|{\boldsymbol\nabla} c|^2}{2}- 3\tilde{\sigma} F(c)+3\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c\right) ({\boldsymbol\nabla}\boldsymbol{\cdot} \boldsymbol{u}) \nonumber\\ &\quad = \frac{1}{3}\left({\mathrm{tr}}\,{\boldsymbol{\mathsf{m}}}_0 - \frac{\tilde{\sigma} \epsilon_c|{\boldsymbol\nabla} c|^2}{2}- 3\tilde{\sigma} F(c)+3\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c \right) \alpha {\boldsymbol\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathsf{q}}}_{C} \nonumber\\ &\quad = \frac{1}{3}{\boldsymbol\nabla} \boldsymbol{\cdot} \left(\left({\mathrm{tr}}\, {\boldsymbol{\mathsf{m}}}_0 - \frac{\tilde{\sigma} \epsilon_c|{\boldsymbol\nabla} c|^2}{2}- 3\tilde{\sigma} F(c)+3\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c \right) \alpha{\boldsymbol{\mathsf{q}}}_{C}\right)\nonumber\\ &\qquad - \frac{1}{3}{\boldsymbol{\mathsf{q}}}_{C}\boldsymbol{\cdot} {\boldsymbol\nabla}\left(\left({\mathrm{tr}}\,{\boldsymbol{\mathsf{m}}}_0 - \frac{\tilde{\sigma} \epsilon_c|{\boldsymbol\nabla} c|^2}{2}- 3\tilde{\sigma} F(c)+3\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c\right) \alpha \right), \end{align}

we define the pressure $p$ as

(2.25)\begin{equation} -p=\frac{1}{3}{\mathrm{tr}}\,{\boldsymbol{\mathsf{m}}}=\frac{1}{3}\left({\mathrm{tr}}\, {\boldsymbol{\mathsf{m}}}_0 - \frac{\tilde{\sigma} \epsilon_c|{\boldsymbol\nabla} c|^2}{2}- 3\tilde{\sigma} F(c) +3\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c\right) , \end{equation}

such that

(2.26)\begin{equation} {\boldsymbol{\mathsf{m}}}_0=-\tilde{\sigma}\epsilon_c({\boldsymbol\nabla} c \otimes{\boldsymbol\nabla}c)+\left(-p+\tilde{\sigma} F(c)+\frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}-\left(\tilde{\sigma} F'(c)-\tilde{\sigma}\epsilon_c{\rm \Delta} c\right)c\right){\boldsymbol{\mathsf{I}}}. \end{equation}

Moreover, the term $({\boldsymbol{\mathsf{m}}}_0 +\tilde {\sigma }\epsilon _c ({\boldsymbol \nabla } c \otimes {\boldsymbol \nabla } c) ): {\boldsymbol {D}}\boldsymbol {u}=0$ due to the above identity. By using (2.24), we rewrite (2.23) as

(2.27)\begin{align} \frac{{\mathrm D}E}{{\mathrm D}t}&= \int_{V(t)} \left\{{\boldsymbol\nabla} \boldsymbol{\cdot}\left( {\boldsymbol{\mathsf{m}}} \boldsymbol{\cdot} \boldsymbol{u} +\tilde{\sigma}\epsilon_c \frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)-{\tau}:{\boldsymbol\nabla} \boldsymbol{u}- {\boldsymbol\nabla} \boldsymbol{\cdot}\left( \left(\tilde{\sigma} F'(c)-\epsilon_c\tilde{\sigma}{\rm \Delta} c+\alpha p\right) {\boldsymbol{\mathsf{q}}}_{C}\right)\right.\nonumber\\ &\quad \left.\vphantom{\left( {\boldsymbol{\mathsf{m}}} \boldsymbol{\cdot} \boldsymbol{u} +\tilde{\sigma}\epsilon_c \frac{{\mathrm D} c}{{\mathrm D}t}{\boldsymbol\nabla} c\right)}+ {\boldsymbol{\mathsf{q}}}_{C}\boldsymbol{\cdot} {\boldsymbol\nabla} \left(\tilde{\sigma} F'(c)-\epsilon_c\tilde{\sigma}{\rm \Delta} c+\alpha p\right)\right\}{\mathrm d}V, \end{align}

(2.28)\begin{equation} \boldsymbol{q}_{E}= {\boldsymbol{\mathsf{m}}} \boldsymbol{\cdot} \boldsymbol{u} + \tilde{\sigma}\epsilon_c \frac{{\mathrm D} \phi}{{\mathrm D}t}{\boldsymbol\nabla} c- \left(\tilde{\sigma} F'(c)-\epsilon_c\tilde{\sigma}{\rm \Delta} c+\alpha p\right) {\boldsymbol{\mathsf{q}}}_{C}, \end{equation}

(2.29)\begin{equation}D=-{\tau}:{\boldsymbol\nabla} \boldsymbol{u}+{\boldsymbol{\mathsf{q}}}_{C} \boldsymbol{\cdot} {\boldsymbol\nabla} \left(\tilde{\sigma} F'(c)-\epsilon_c \tilde{\sigma}{\rm \Delta} c+\alpha p\right). \end{equation}

Therefore, to ensure the non-positivity of the dissipation term ($D \leqslant 0$, the second law of thermodynamics), we specify the unknown terms in the form

(2.30)\begin{gather} {\tau}= \nu(c) (\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T}) - \tfrac{2}{3}\nu(c)(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}) {\boldsymbol{\mathsf{I}}}, \end{gather}

(2.31)\begin{gather}{\boldsymbol{\mathsf{q}}}_{C}=-M(c){\boldsymbol\nabla}\left(\mu+\alpha p\right), \end{gather}

where $M(c)\geqslant 0$ is the mobility function, $\nu (c)=\nu _{1}c+\nu _{2}(1-c)$ is variable viscosity for the two-phase fluid and $\nu _{i}$ is the viscosity for the $i$th fluid, $\mu =\tilde {\sigma } F'(c)-\epsilon _c\tilde {\sigma }{\rm \Delta} c$. Combining (2.29) and (2.30) we obtain

(2.32)\begin{equation} {\boldsymbol{\mathsf{m}}}={\tau}-\tilde{\sigma}\epsilon_c({\boldsymbol\nabla} c \otimes{\boldsymbol\nabla} c)+\left(-p+\tilde{\sigma} F(c)+\frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}-c\mu\right){\boldsymbol{\mathsf{I}}} .\end{equation}

2.3. Governing equations

By substituting (2.30)–(2.32) in (2.14) and (2.15), and using the identity

(2.33)\begin{equation} {\boldsymbol\nabla} \left(\tilde{\sigma} F(c)+\frac{\tilde{\sigma}\epsilon_c |\boldsymbol{\nabla} c|^2}{2}-c\mu \right)=\tilde{\sigma}\epsilon_c{\boldsymbol\nabla}\boldsymbol{\cdot} \left({\boldsymbol\nabla} c \otimes{\boldsymbol\nabla}c\right)-c{\boldsymbol\nabla}\mu, \end{equation}

we obtain the q-NSCH system equations governing two-phase flows:

(2.34)\begin{gather} \rho \boldsymbol{u}_{t}+\rho\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} = -{\boldsymbol\nabla} p+{\boldsymbol\nabla} \boldsymbol{\cdot} \left(\nu\left(c\right) {\boldsymbol\nabla}\boldsymbol{u}\right)+\tfrac{1}{3}{\boldsymbol\nabla}\left(\nu\left( c\right){\boldsymbol\nabla}\boldsymbol{\cdot} \boldsymbol{u}\right) +\mu{\boldsymbol\nabla}c-\rho g \boldsymbol{z}, \end{gather}

(2.35)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=\alpha{\boldsymbol\nabla} \boldsymbol{\cdot}\left(M(c){\boldsymbol\nabla}\mu\right)+\alpha^2{\boldsymbol\nabla}\boldsymbol{\cdot}\left( M(c){\boldsymbol\nabla} p\right), \end{gather}

(2.36)\begin{gather}c_{t}+\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{u}c)={\boldsymbol\nabla}\boldsymbol{\cdot}\left( M(c) {\boldsymbol\nabla} \mu\right)+\alpha{\boldsymbol\nabla}\boldsymbol{\cdot}\left( M(c){\boldsymbol\nabla}p\right), \end{gather}

(2.37)\begin{gather}\mu=\tilde{\sigma} F'(c)-\epsilon_c \tilde{\sigma}{\rm \Delta} c. \end{gather}

Here $M(c)=\sqrt {c^2(1-c)^2}$ is the degenerate mobility, which becomes zero in the bulk regions ensuring that the velocity field is divergence free for single-phase flow. Moreover, the following boundary conditions on $\partial \varOmega$ are imposed:

(2.38)\begin{gather} \boldsymbol{u}= \boldsymbol{u}_{g}\quad (\textrm{no slip}), \end{gather}

(2.39)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}(\mu+\alpha p)=0\quad (\textrm{no flux}), \end{gather}

(2.40)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}c= \frac{1}{\epsilon_c\sqrt{2}}\cos (\theta) (1-c)c\quad (\textrm{contact angle}). \end{gather}

Here $\boldsymbol {u}_{g}=(g_{u},g_{v})$ is the velocity on the boundary $\partial \varOmega$ and $\boldsymbol {n}$ is the outward normal vector to $\varOmega$. Assuming $\partial \varOmega$ is closed, the boundary velocity needs to satisfy $\int _{\partial \varOmega }\boldsymbol {u}_{g}\boldsymbol {\cdot }\boldsymbol {n}\,{\mathrm d}S=0$. Further, the contact angle $\theta$ is the static contact angle between the fluid–fluid interface and the physical boundary. This form of the boundary condition (2.40) was derived in Jacqmin (Reference Jacqmin1999) (see also Gránásy et al. (Reference Gránásy, Pusztai, Saylor and Warren2007), Do-Quang & Amberg (Reference Do-Quang and Amberg2009) and Aland et al. (Reference Aland, Lowengrub and Voigt2010) in the context of the DD method). In particular, for an interface in equilibrium, a contact angle condition (Young's law) may be posed: $\cos (\theta )=(({\sigma _{2S}-\sigma _{1S}})/{\sigma })$, where $\theta$ is the (equilibrium) contact angle and $\sigma _{iS}$ is the interfacial energy between the solid and component $i$. The contact angle characterizes the wettability of the surface with $\theta =0$ denoting complete wetting. However, in the dynamical setting with the no-slip condition, the force needed to move a contact line is infinite making it necessary to regularize the system. Numerous methods have been implemented to regularize the problem including introducing a slip velocity, microscopic interactions, a precursor film or a diffuse interface.

Recently, two other quasi-incompressible two-fluid models have been presented in Gong et al. (Reference Gong, Zhao and Wang2017) and Shokrpour Roudbari et al. (Reference Shokrpour Roudbari, Şimşek, van Brummelen and van der Zee2018), where a different phase variable $\phi$ is introduced leading to a different Cahn–Hilliard free energy and different systems of equations. In the paper by Shokrpour Roudbari et al. (Reference Shokrpour Roudbari, Şimşek, van Brummelen and van der Zee2018), in order to ensure the energy conservation for both the model equations and the corresponding numerical method, they introduced a dual-density form (two formulations of the variable density of the mixture) to the modelling framework. These two forms of the density appear in the model simultaneously and cannot be used interchangeably, leading to an extra mass conservation equation that needs to be solved along with the model equations. Our model equations are more compact and satisfy the energy conservation naturally. There is only one form of the variable density, and no extra mass conservation equation. In the paper by Gong et al. (Reference Gong, Zhao and Wang2017), the modelling framework they follow is based on the free energy dissipation which means that the framework is mainly focused on the isothermal case. Our modelling framework is more general, as the total energy, free energy and entropy can be introduced to the system in a thermodynamically consistent way, so that we can derive more general models to simulate the non-isothermal case or two-phase flows with non-uniform surface tension (surface tension gradient induced by the gradient of temperature or concentration of surfactants). Moreover, in their modelling framework, they did not provide the specific formulation for the Cahn–Hilliard free energy, which shortens the model derivation. However, to ensure the mass conservation of the fluid mixture and each fluid component, a parameter ‘$B$’ is introduced to the Cahn–Hilliard equation, serving as the coefficient of the mobility of the chemical potential. In their work, the calculation of ‘$B$’ seems quite long and complicated, and the resulting ‘$B$’ also makes the model equations more complicated compared to our model.

4.1. Cavity flow

We first consider a two-phase fluid flow in a driven cavity (Boyer Reference Boyer2002). The numerical results from the q-NSCH system (2.34)–(2.37) on the physical domain $\varOmega = [0,1]^2$ are compared to those from the q-NSCH-DD system (3.1)–(3.4) on an extended domain $\tilde \varOmega = [-0.25,1.25]^2$. For the original q-NSCH system, the boundary conditions on $\partial \varOmega$ are imposed as

(4.1)\begin{gather} \boldsymbol{u}= \boldsymbol{u}_{g}, \end{gather}

(4.2)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}(\mu+\alpha p)=0, \end{gather}

(4.3)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}c= 0, \end{gather}

where $\boldsymbol {u}_{g}=(g_{u},g_{v})=(g_{u},0)$, and

(4.4)\begin{equation} g_{u} = \begin{cases} {\max} \left(0, 1-4 \times (x-0.5)^2\right) , & y=1 \\ 0, & y<1. \end{cases} \end{equation}

For the q-NSCH-DD system, the function $\phi$, given by

(4.5)\begin{align} \phi&=0.5\left(1+\tanh {\left( \frac{3x}{\epsilon}\right)}\right) \times 0.5\left(1+\tanh {\left( \frac{3(1-x)}{\epsilon}\right)}\right)\nonumber\\ &\quad \times0.5\left(1+\tanh {\left( \frac{3y}{\epsilon}\right)}\right) \times 0.5\left(1+\tanh {\left( \frac{3(1-y)}{\epsilon}\right)}\right), \end{align}

is used to approximate the characteristic function of $\varOmega$, and $g_{u}$ is extended to $\tilde \varOmega$ as follows:

(4.6)\begin{equation} g_{u} = \begin{cases} {\max} \left(0, 1-4 \times (x-0.5)^2\right) , & y\geqslant 1-\epsilon \\ 0, & y<1-\epsilon. \end{cases} \end{equation}

For both systems, the initial $c$ is

(4.7)\begin{equation} c=\frac{1}{2}\left( 1- \tanh \left(\frac{0.5-y}{2\sqrt{2}\epsilon_c}\right)\right), \end{equation}

such that $c=0$ for the fluid at the bottom and $c=1$ for the top fluid. In addition, we set

(4.8a–f)\begin{align} \rho_{1}=\rho_{2}=1,\quad \nu_{1}=\nu_{2}=0.002,\quad g=0,\quad \tilde{\sigma} = 0,\quad \theta = 90^{\circ}\quad {\mathrm{and}}\quad \epsilon_c= 0.05. \end{align}

We present the evolution of the phase-field variable $c$ obtained from the q-NSCH-DD system in figure 2 and compare $c$ at $t=15$ with the one from the q-NSCH system at the same time in figure 3.

Figure 2.

Evolution of the phase-field variable $c$ in a cavity flow. The results are obtained using the q-NSCH-DD system with a grid size $h=1.5/256$, $\textrm {d}t = 2.5\times 10^{-3}$ and $\epsilon =0.025$ in $\tilde \varOmega$. The white box represents the physical domain $\varOmega$.

Figure 3.

Comparisons between the phase-field variable obtained from the q-NSCH-DD system (a) and that from the q-NSCH system (b) at $t=15$. The white box in (a) marks the physical domain $\varOmega$.

Moreover, we compare the results between the q-NSCH and q-NSCH-DD solutions to test the convergence of the q-NSCH-DD system as $\epsilon \rightarrow 0$. In particular, we use the results, including the velocities $u_{exc}$, $v_{exc}$, and phase variable $c_{exc}$, from the q-NSCH system with a grid size $h=1/1024$ and a time step $\textrm {d}t=3.90625\times 10^{-5}$ as the exact solutions. For the q-NSCH-DD system, we obtain the corresponding numerical solutions $u_\epsilon$, $v_\epsilon$ and $c_\epsilon$ using different values of $\epsilon$, namely $0.1$, $0.05$, $0.025$, $0.0125$, $0.00625$. We start with a grid size $h = 1.5/64$ and time step $\textrm {d}t=10^{-2}$ and $\epsilon =0.1$, and refine them together. For each $\epsilon$, we restrict the solutions from q-NSCH-DD to the original domain $\varOmega$ (see the smaller region bounded by the white solid line in figure 2) and compute the $L^2$ ($E^{(2)}$) and $L^\infty$ ($E^{(\infty )}$) errors as

(4.9a,b)\begin{equation} E_{y}^{(2)}=||y_\epsilon-y_{exc}||_{L^{2}(\varOmega)}\quad \textrm{and}\quad E_{y}^{(\infty)}=||y_\epsilon-y_{exc}||_{L^{\infty}(\varOmega)}, \end{equation}

where $y_\epsilon$ represents $u_\epsilon$, $v_\epsilon$, $c_\epsilon$ and $y_{exc}$ represents $u_{exc}$, $v_{exc}$, $c_{exc}$. Note that $E_{y}^{(\infty )}$ is computed where $\phi \geqslant 0.5$. The errors at $t=0.1$ are presented in tables 1 and 2. We observe that our DD approximation converges to the q-NSCH model with $\mathcal {O}(\epsilon )$.

Remark Another choice of $BC_u = -\epsilon ^{-3}(1-\phi )(\boldsymbol {u}-\boldsymbol {u_g})$ was presented by Li et al. (Reference Li, Lowengrub, Rätz and Voigt2009) and later the scheme was shown to be $\mathcal {O}(\epsilon \ln (\epsilon ))$ accurate by Yu et al. (Reference Yu, Guo and Lowengrub2020). We here provide the numerical results using $BC_u = -\epsilon ^{-3}(1-\phi )(\boldsymbol {u}-\boldsymbol {u}_{\boldsymbol {g}})$ for the cavity flow example in tables 3 and 4. This form of the velocity boundary condition was also used by Aland et al. (Reference Aland, Lowengrub and Voigt2010). The numerical set-up is identical to that considered above. We find that the errors are of a similar order of magnitude but are larger than those using $BC_u = -\epsilon ^{-2}(1-\phi )(\boldsymbol {u}-\boldsymbol {u}_{\boldsymbol {g}})$ and are sub-first-order accurate, which is consistent with the results of Yu et al. (Reference Yu, Guo and Lowengrub2020). In the following, we use $BC_u = -\epsilon ^{-2}(1-\phi )(\boldsymbol {u}-\boldsymbol {u}_{\boldsymbol {g}})$ exclusively. Moreover, our numerical results agree with the convergence results shown in Franz et al. (Reference Franz, Roos, Gärtner and Voigt2012).

Table 1.

The $L^2$ errors and convergence rate ($k$) of the velocities $u$, $v$ and phase-field variable $c$ with different values of $\epsilon$ in a cavity flow in § 4.1. Here $k$ stands for the convergence rate. See text for details.

Table 2.

The $L^\infty$ errors and convergence rate ($k$) of the velocities $u$, $v$ and phase-field variable $c$ in a cavity flow in § 4.1. Here $k$ stands for the convergence rate. See text for details.

Table 3.

The $L^2$ errors and convergence rate ($k$) of the velocities $u$, $v$ and phase-field variable $c$ in a cavity flow using another choice of $BC_u = -\epsilon ^{-3}(1-\phi )(\boldsymbol {u}-\boldsymbol {u_g})$. See the remark in § 4.1 for details.

Table 4.

The $L^\infty$ errors of the velocities $u$, $v$ and phase-field variable $c$ in a cavity flow using another choice of $BC_u = -\epsilon ^{-3}(1-\phi )(\boldsymbol {u}-\boldsymbol {u_g})$. See the remark in § 4.1 for details.

4.2. Droplet spreading

In order to test the accuracy of the q-NSCH-DD system with the moving contact line, we consider an example of a water droplet spreading on a flat substrate. In particular, we assume that there is a water droplet in the air, and is placed above a flat substrate. Due to the effect of gravity, the droplet falls down and spreads on the substrate. For the original q-NSCH system, the computational domain is $\varOmega =[0,1]^2$ with the following boundary conditions on $\partial \varOmega$:

(4.10)\begin{gather} \boldsymbol{u}= 0\quad (\textrm{no slip}), \end{gather}

(4.11)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}(\mu+\alpha p)=0\quad (\textrm{no flux}), \end{gather}

(4.12)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}c= \frac{1}{\epsilon_c\sqrt{2}}\cos (\theta) (1-c)c\quad (\textrm{contact angle}). \end{gather}

For the q-NSCH-DD system, we extend the computational domain to $\tilde {\varOmega }=[-0.25,1.25]^2$ with the following boundary conditions on $\partial \tilde {\varOmega }$:

(4.13)\begin{gather} \boldsymbol{u}= 0\quad (\textrm{no slip}), \end{gather}

(4.14)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}(\mu+\alpha p)=0\quad (\textrm{no flux}), \end{gather}

(4.15)\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\nabla}c= 0\quad (\textrm{no flux}). \end{gather}

Moreover, we use the function $\phi$, given by

(4.16)\begin{align} \phi&=0.5\left(1+\tanh {\left( \frac{3x}{\epsilon}\right)}\right) \times 0.5\left(1+\tanh {\left( \frac{3(1-x)}{\epsilon}\right)}\right)\nonumber\\ &\quad \times 0.5\left(1+\tanh {\left( \frac{3y}{\epsilon}\right)}\right)\times 0.5\left(1+\tanh {\left( \frac{3(1-y)}{\epsilon}\right)}\right), \end{align}

to approximate the characteristic function of $\varOmega$. For both systems, the initial condition for the phase variable $c$ is given by

(4.17)\begin{equation} c=\frac{1}{2}\left( 1- \tanh \left(\frac{r-0.2}{2\sqrt{2}\epsilon_c}\right)\right), \end{equation}

with $r=\sqrt {(x-0.5)^2+(y-0.21)^2}$, such that $c=0$ for the water droplet and $c=1$ for the air. In addition, we set

(4.18)\begin{align} \left.\begin{gathered} \rho_{1}=1,\quad \rho_{2}=1000,\quad \nu_{1}=0.1,\quad \nu_{2}=10,\quad g=0.98, \quad \tilde{\sigma} = 10,\quad \theta = 60^{\circ},\\ \epsilon= 0.0025\quad {\mathrm{and}}\quad \epsilon_c= 0.005. \end{gathered}\right\} \end{align}

In figure 4, we show the domains $\tilde {\varOmega }$ and $\varOmega$ together with the initial $c$ for both systems. In figure 5, we show the dynamics of the contact line of the moving droplet, where the solid (dotted) lines represent the contour line ($c=0.5$) of the numerical results obtained from the q-NSCH (the q-NSCH-DD) system. A good agreement can be observed between the two results. Again, we compare the results obtained from the two systems, and test the convergence rate of the q-NSCH-DD system as $\epsilon \rightarrow 0$. As in the previous section, we use the results of phase variable $c_{exc}$ from the q-NSCH system with agrid size $h=1/1024$ and a time step $\textrm {d}t=3.90625\times 10^{-5}$ as the exact solutions. For the q-NSCH-DD system, we obtain the corresponding numerical solutions $c_\epsilon$ using different values of $\epsilon$, namely $0.1$, $0.05$, $0.025$, $0.0125$, $0.00625$. We start with a grid size $h = 1.5/64$ and time step $\textrm {d}t=10^{-2}$ and $\epsilon =0.1$, and refine them together. For each $\epsilon$, we restrict the solutions from q-NSCH-DD to the original domain $\varOmega$ and compute the $L^2$ ($E^{(2)}$) and $L^\infty$ ($E^{(\infty )}$) errors as shown in (4.9a,b). The errors at $t=0.1$ are presented in table 5, where we observe that our DD approximation converges to the q-NSCH system with $\mathcal {O}(\epsilon )$.

Figure 4.

The computational domain for the q-NSCH-DD system (a) and the q-NSCH system (b), and the initial condition for phase variable $c$. The white box in (a) marks the physical domain $\varOmega$.

Figure 5.

The dynamics of moving contact line of a water droplet on a solid substrate, where the solid (dotted) lines are the contour lines ($c=0.5$) of the numerical results obtained from the q-NSCH (q-NSCH-DD) system.

Table 5.

The $L^2$ and $L^{\infty }$ errors and convergence rate ($k$) of the phase-field variable restricted on $\varOmega$, $\phi c$, at time $t=0.1$ with contact angle $\theta = 60^{\circ }$. See § 4.2 for details.

4.3. Droplet on cylinder

We next consider a circular droplet of diameter $d$ wetting and dewetting on a solid cylinder. Initially, the droplet is placed on the cylinder of the same diameter $d$ with contact angle $\theta _{0}=90^{\circ }$ (see figure 6). We mention that careful studies for wetting and dewetting phenomena using a phase-field model have been presented in Gránásy et al. (Reference Gránásy, Pusztai, Saylor and Warren2007), Dziwnik, Münch & Wagner (Reference Dziwnik, Münch and Wagner2017) and Rainer et al. (Reference Rainer, Steven, Marco and Axel2019). With $\theta$ being the static contact angle, the contact line of the droplet recedes if $\theta <\theta _{0}$ or spreads if $\theta >\theta _{0}$, and eventually reaches an equilibrium with the interface intersecting the solid cylinder at an angle of $\theta$. In the absence of gravity, the equilibrium state of the droplet can be obtained analytically by carving the cylinder out of a circle (Liu & Ding Reference Liu and Ding2015; O'Brien & Bussmann Reference O'Brien and Bussmann2018); see figure 7. The distance between the centres of the droplet and the cylinder, $d_{c}$, and the diameter of the circle, $d_{l}$, can be computed by solving

(4.19)\begin{gather} \pi d^2/4=(\alpha + \theta )(d_{l}/2)-\alpha(d/2)^2+(d_{l}d/4)\sin{\theta}, \end{gather}

(4.20)\begin{gather}d_{l}=d\cos{\theta}+\sqrt{d^{2}\cos^2{\theta}-d^2+4d_{c}^{2}}, \end{gather}

where

(4.21)\begin{equation} \alpha=\arccos{\left(\frac{d^{2}+4d_{c}^{2}-d^{2}_{l}}{4d d_{c}}\right)}.\end{equation}

The exact solution $c_{exc}$ is given by

(4.22)\begin{equation} c_{exc}=\frac{1}{2}\left( 1+ \tanh \left(\frac{r_{exc}-0.5d_{l}}{2\sqrt{2}\epsilon_c}\right)\right) \end{equation}

and $r_{exc}=\sqrt {(x-1.25d)^2+(y-0.6d-d_{c})^2}$. To apply the q-NCSH-DD system, we extend $\varOmega$ to a larger domain $\tilde \varOmega$ = $[0,2.5d]\times [0,2.5d]$. The initial $\phi$ and $c$ are given by

(4.23)\begin{equation} \left.\begin{gathered} \phi=\frac{1}{2}\left( 1+\tanh \left(\frac{3(r_{\phi}-0.5d)}{\epsilon}\right)\right),\\ c=\frac{1}{2}\left( 1+ \tanh \left(\frac{r_c-0.5d}{2\sqrt{2}\epsilon_c}\right)\right),\end{gathered}\right\} \end{equation}

where $r_{\phi }\!=\!\sqrt {(x\!-\!1.25d)^2\!+\!(y\!-\!0.6d)^2}$ and $r_c\!=\!\sqrt {(x\!-\!1.25d)^2\!+\!(y\!-\!0.6d\!-\!0.5\sqrt {2}d)^2}$, such that $c=1$ for the droplet (indicated by 1) and $c=0$ for the surrounding fluid medium (indicated by 2). All the other parameters are set as follows:

(4.24a–h)\begin{equation} \left.\begin{gathered} \rho_{1}=1000\ (\mathrm{droplet}),\quad \rho_{2}=1\ (\mathrm{air}), \quad\nu_{1}=0.1,\quad \nu_{2}=0.001,\quad \tilde{\sigma}=1, \\ g=0,\quad \epsilon_c=\epsilon,\quad d=1. \end{gathered}\right\} \end{equation}

Figure 6.

Set-up for droplet on cylinder.

Figure 7.

Comparisons between the numerical results at equilibrium (solid lines) and the analytic solutions (dashed lines) with different static contact angles, where the angles are indicated inside the blue cylinders. See text for details.

We solve the q-NSCH-DD system with different values of $\epsilon$ ($=\epsilon _{c}$) and various $\theta$ ($=30^\circ , 90^\circ ,120^\circ$) to obtain the numerical solutions $c$ at equilibrium and then compare with the corresponding analytical solutions $c_{exc}$ in both $L^2$ and $L^{\infty }$ norm:

(4.25a,b)\begin{equation} E_{c}^{(2)}=||c-c_{exc}||_{L^{2}(\varOmega)}\quad \textrm{and} \quad E_{c}^{(\infty)}=||c-c_{exc}||_{L^{\infty}(\varOmega)}. \end{equation}

We start with a grid size $h=2.5/128$, a time step $\textrm {d}t=0.01$ and $\epsilon =\epsilon _c=0.02$ and conduct the error analysis by refining them together. In tables 6 and 7, we present the $L^2$ and $L^\infty$ errors and find that our q-NSCH-DD system provides a first-order approximation to the exact solution (sharp interface).

Table 6.

The $L^2$ errors and convergence rate ($k$) of the phase-field variable restricted on $\varOmega$, $\phi c$, at equilibrium with different static contact angles $30^{\circ }$, $90^{\circ }$ and $120^{\circ }$.

Table 7.

The $L^{\infty }$ errors of the phase-field variable restricted on $\varOmega$, $\phi c$, at equilibrium with different static contact angles $30^{\circ }$, $90^{\circ }$ and $120^{\circ }$. See text for details.

In figure 7, we also present the numerical solutions $\phi c$ with analytical solutions at equilibrium from the q-NSCH-DD system with various $\theta$ ($=30^{\circ }, 90^\circ ,120^\circ$). The analytical solutions are the dashed curves. Here the grid size is $n=2048$ and $\epsilon =\epsilon _{c}=0.00125$. There is excellent agreement between the numerical and analytical solutions.

4.4. Penetration of a liquid droplet into porous media

As a final test, we consider a two-dimensional droplet falling through a porous substrate consisting of three cylinders ($C_i$, $i=1,2,3$) with different wettabilities, as shown in figure 8. We also fix the contact angle of the fluid–fluid interface at the walls to be $90^{\circ }$. The extended domain $\tilde \varOmega = [0,1]\times [0,1]$ is discretized on a uniform mesh of $[1024\times 1024]$. The original physical domain is $\varOmega = \tilde \varOmega /(\cup C_i)$. There are three cylinders in the domain with radii $r_1=0.1$, $r_2=0.12$, $r_3=0.15$, and centres $(0.3,0.6)$, $(0.4,0.25)$, $(0.7,0.5)$, respectively. Initially, a heavy liquid drop with radius $r_d=0.15$ is placed above the cylinders with the centre $(0.5,0.8)$. Note that the initial $\phi$ and $c$ can be given by a slight modification of (4.23). Again, we let $c=1$ stand for the heavy fluid droplet and $c=0$ stands for the light fluid medium. All other parameters are set as follows:

(4.26a–h)\begin{equation} \left.\begin{gathered} \rho_{1}=1000,\quad \rho_{2}=1,\quad\nu_{1}=10,\quad\nu_{2}=0.1,\quad g=1,\quad\tilde{\sigma}=1,\\ \epsilon_c=0.001,\quad\epsilon=0.002. \end{gathered}\right\} \end{equation}

The dynamics of droplet impact is shown in figure 8 via snapshots at different times. Initially, the droplet starts moving downwards due to gravitational forces. Dramatic deformations are observed when the droplet comes into contact with the left and right cylinders. After impacting on the cylinder at the bottom, the droplet splits into several parts. In the end, it is interesting to see that there is no liquid on the cylinders with contact angle of $\theta =90^{\circ }$ and $\theta =120^{\circ }$. All the liquid rests either on the cylinder with $\theta =60^{\circ }$ or on the walls.

Figure 8.

Snapshots of the penetration of a liquid droplet into porous media with different static contact angles, where the angles are indicated inside the white cylinders.

4.5. Moving domains

In the next example, we simulate solid hydrophobic and hydrophilic balls impacting an air–water interface in two dimensions. The extended domain is $\tilde {\varOmega }=[0,1]\times [0,1]$, and the air is placed on top of the water, with a flat interface at $y=0.5$. Figure 9 shows the initial $\phi$ (the complex domain) and $c$ (the air–water interface), which are given by

(4.27)\begin{equation} \left.\begin{gathered} \phi=\frac{1}{2}\left( 1+\tanh \left(\frac{3(r_{\phi}-0.1)}{\epsilon}\right)\right),\\ c=\frac{1}{2}\left( 1- \tanh \left(\frac{y-0.5}{2\sqrt{2}\epsilon_c}\right)\right), \end{gathered}\right\} \end{equation}

where $r_{\phi }=\sqrt {(x-0.5d)^2+(y-0.65d)^2}$ and such that $\phi =0$ for the solid ball (in white) and $\phi =0$ for the bulk regions with water (in red) or air (in blue). Here we assume that the ball moves with a constant velocity $\boldsymbol {v}_{ball}=(0,-1)$. All the other parameters are set as follows:

(4.28)\begin{equation} \left.\begin{gathered} \rho_{1}=1000\ (\mathrm{water}, \mathrm{lower\ fluid}),\quad \rho_{2}=1\quad (\mathrm{air}, \mathrm{upper\ fluid}),\quad \nu_{1}=0.1, \\ \nu_{2}=0.001,\quad \tilde{\sigma}=4,\quad g=0.98,\quad \epsilon_c=\epsilon=0.003. \end{gathered}\right\} \end{equation}

Two contact angle boundary conditions, namely $\theta _{0}=60^{\circ }$ (hydrophilic) and $\theta _{0}=120^{\circ }$ (hydrophobic), are imposed at the boundary of the solid ball. Figures 9 and 10 show the results of the simulations. For the hydrophilic ball, the contact line moves continuously on the surface of the ball until reaching the ball peak, at which the liquid merges and detaches from the surface of the ball. In contrast, for the hydrophobic ball, the contact line initially moves along the ball surface, but eventually rests on the ball. As a consequence of the pinned contact line, the length of the air cavity keeps growing. A neck then appears at the cavity and breaks up, thereby generating a bubble attached at the top of the ball.

Figure 9.

Snapshots of water entry of a hydrophilic ball with contact angle $60^{\circ }$. Here the white region represents the solid falling ball, the blue region represents the air and the red region represents the water.

Figure 10.

Snapshots of water entry of a hydrophobic ball with contact angle $120^{\circ }$. Here the white region represents the solid falling ball, the blue region represents the air and the red region represents the water.