Phase field modeling and simulation of three-phase flow on solid surfaces

Abstract

Phase field models are widely used to describe the two-phase system. The evolution of the phase field variables is usually driven by the gradient flow of a total free energy functional. The generalization of the approach to an N phase ($N\ge 3$) system requires some extra consistency conditions on the free energy functional in order for the model to give physically relevant results. A projection approach is proposed for the derivation of a consistent free energy functional for the three-phase Cahn–Hilliard equations. The system is then coupled with the Navier–Stokes equations to describe the three-phase flow on solid surfaces with moving contact line. An energy stable scheme is developed for the three-phase flow system. The discrete energy law of the numerical scheme is proved which ensures the stability of the scheme. We also show some numerical results for the dynamics of triple junctions and four phase contact lines.

Introduction

Phase field method has been widely used for numerical simulations of two-phase systems. The basic idea of the phase field method is to introduce an order parameter ϕ which takes two distinct values (+1 and −1, for instance) in each of the phases, with a smooth transition between both values in the zone around the interface, which is then diffuse with a finite width. The key advantage of the method is that the need to explicitly track the interface between different phases is removed. The phase field concept can be generalized to systems with three phases. In such three-phase system, the scalar variable ϕ is replaced by a vector c where the i-th element $c_i$ represents the volume fraction of i-th phase. These phases are linked through the constraint

$$c_1+c_2+c_3=1.$$

The evolution of the three-phase system is then driven by the gradient of a total free energy of the system subject to the constraint Eq. (1.1). If the free energy of the solid boundary is taken into account, the total free energy is then a sum of three terms: a bulk free energy which is usually taken as a multi-well potential describing the free energy density of the bulk of each phase, an interface energy term depending on the gradient of c and a surface energy on the solid boundary. Considering a three-phase system in domain Ω with a solid boundary ∂Ω, the total free energy functional can be written as

$${\mathfrak{F}}_3(\mathbf{c})=\underset{\mathrm{\Omega }}{\int }F(\mathbf{c})d\mathbf{x}+\underset{\partial \mathrm{\Omega }}{\int }{\gamma }_3(\mathbf{c})d\mathbf{s},\text{where}F(\mathbf{c})=F_3(\mathbf{c})+G_3(\mathbf{c},\mathrm{\nabla }\mathbf{c}).$$

Here, $G_3(\mathbf{c},\mathrm{\nabla }\mathbf{c})$ is the interfacial free energy density, and $F_3(\mathbf{c})$ and ${\gamma }_3(\mathbf{c})$ are the free energy densities in the bulk and on the solid boundary, respectively. The bulk free energy $F_3(\mathbf{c})$ should be a “triple-well” potential with minima at $c_i=1$ for all $i\in \{1,2,3\}$. It is natural to require that the three-phase free energy functional is consistent with that of two-phase flow, i.e., given $i\in \{1,2,3\}$, if $c_i=0$ and $j\ne i,k\ne i$, then

$$F_3(c_j,c_k,0)=F_2(c_j,c_k),G_3(c_j,c_k,0)=G_2(c_j,c_k),{\gamma }_3(c_j,c_k,0)={\gamma }_2(c_j,c_k),$$

where $F_2,G_2$ and ${\gamma }_2$ are the corresponding free energy densities in two-phase flow on the solid surfaces. If consider the time evolution of the three-phase system driven by the minimization of the free energy functional, the dynamic equations are

$$\frac{\partial c_i}{\partial t}=-{\mu }_i$$

with a boundary condition on ∂Ω

$$\frac{\partial c_i}{\partial t}=-L_i$$

for $i=1,2,3$, where ${\mu }_i$ and $L_i$ are defined by the variation of the total free energy functional ${\mathfrak{F}}_3(\mathbf{c})$, that is

$$\delta {\mathfrak{F}}_3(\mathbf{c})=\underset{\mathrm{\Omega }}{\int }d\mathbf{x}[{\mu }_i\delta c_i]+\underset{\partial \mathrm{\Omega }}{\int }d\mathbf{s}[L_i\delta c_i].$$

For (1.2) and (1.3), the following extra consistency conditions (P.1) and (P.2) are also needed in order for the system to give physically relevant results, see also [3], [9],

• When the phase i does not present in the mixture at the initial time, the phase should not appear during the time evolution of the system, i.e. given $i\in [1,3]$,

  $$c_i(0)=0⟹c_i(t)=0,\forall t\ge 0.$$

• When there are only two phases in the model, the model of three-phase problem should always degenerate into the corresponding two-phase model.

Boyer et al. [3] shows that, for the three-phase problem, one must carefully choose the form of the bulk free energy so that the model can be well-posed and satisfy the algebraically and dynamically consistency conditions. However, their method is not easy to be generalized to general N-phase problems with $N>3$ and to problems with general boundary conditions such as (1.3).

When the multiphase fluid dynamics is considered, the multiphase system is then coupled with the Navier–Stokes equations. For two-phase flow on solid boundary, the generalized Navier boundary condition is proposed for the moving contact line problem [12]. The two-phase flow model has been studied extensively [13], [14] with many efficient numerical methods developed [6], [8]. There are, however, very few results on the three-phase (or components) flow with solid boundary modeling and simulations [3], [11], [16], [17].

In this paper, a new approach is developed to derive the consistent free energy functional for three-phase flow with moving contact line problem. The idea of the approach is to introduce a projection operator to enforce constraints Eq. (1.1), (P.1) and (P.2). We then employ the operator to derive the three-phase Cahn–Hilliard Navier–Stokes equations with generalized Navier boundary condition. In addition, we develop an energy stable scheme for the coupled Cahn–Hilliard Navier Stokes equations with boundary conditions on solid surface. We also show that the scheme has the total energy decaying property which ensures the stability of the scheme. Numerical examples are carried out to verify capability of our model and numerical scheme.

The rest of this paper is organized as follows. In Section 2, the consistent bulk free energy functional and surface energy functional is derived from the projection approach. The energy functional is then incorporated into the Cahn–Hilliard Navier Stokes system for three-phase flow with the generalized Navier boundary condition for the motion of the contact line. In Section 3, the energy law for the model of three-phase flows on solid surfaces is proved. The numerical scheme is developed and discrete energy laws of the numerical scheme is discussed in Section 4. In Section 5, we give some numerical examples for dynamics of triple junctions and four phase contact lines (points). The discrete energy law and mass conservation law is verified through numerical examples. The paper concludes in Section 6.