Efficient Second Order Unconditionally Stable Schemes for a Phase-field Moving Contact Line Model Using Invariant Energy Quadratization Approach

In this paper, we consider numerical approximations for a hydrodynamically-coupled phase field model [34, 35] with moving contact line (MCL) boundary conditions. Phase field method is a popular approach that is widely used to simulate the interfacial dynamics of multiple material components due to its versatility in modeling as well as numerical simulations. The fluid-fluid interface in this method is considered as a continuous but steep change of some physical property of two fluid components, e.g., density or viscosity, etc. An order parameter (phase field variable) is introduced to label the two fluid components, the interface is then represented by a thin but smooth transition layer that can remove the singularities in practice. The standard phase field model for incompressible immiscible fluid mixture is a nonlinear system that couples the Cahn-Hilliard equation and the Navier-Stokes equations via convective and stress terms. Once the fluid-fluid interface touches a solid wall, a MCL problem is induced. This phenomenon exists in many physical and engineering processes such as wetting, coating, or painting, etc. In this situation, the no-slip boundary condition for the Navier-Stokes equations is no longer applicable (cf. [32, 7, 6]). Simulations in [22, 23, 49] using molecular dynamics simulations showed that nearly complete slip happens near the MCL. In the context of phase field method, a set of accurate boundary conditions for the MCL problem was derived by Qian et. al. in [34, 35], resulting in a standard two phase model supplemented with a generalized Navier boundary condition (GNBC) and a dynamic contact line condition (DCLC), where the nonlinear couplings also show up in the boundary conditions.

Numerically, although the phase field variable is continuous and smooth, the model is still very stiff since a small dimensionless parameter related to the thickness of the interface layer, is involved, for which certain numerical methods like fully implicit or explicit methods for nonlinear terms, need a very small time marching step size (see the analysis in [11, 54, 42]). Hence a challenging issue for solving the model is to develop numerically stable methods, namely, to establish efficient numerical schemes that can verify the so-called energy stable property at the discrete level, irrespective of the coarseness of the discretization. Such a kind of algorithms is usually called unconditionally energy stable or thermodynamically consistent. Schemes with this property are specially preferred since it is critical for the numerical schemes to use larger time steps to capture the correct long time dynamics that can reduce the time cost in computing. Nonetheless, we need mention a basic fact that a larger time step will definitely induce larger numerical errors. In other words, schemes with unconditional energy stability can allow arbitrarily large time step only for the sake of the stability concern. To measure whether a scheme is reliable or not, the controllable accuracy is another important factor in addition to the stability. Therefore, if one attempts to use a time step as large as possible while maintaining the desirable accuracy, the only possible choice is to develop more accurate schemes, e.g., second order energy stable schemes, that is the main focus of this paper.

It is remarkable that, unlike the enormous numerical scheme developments on the standard phase field model for two phase fluid flows system with easy boundary conditions (without MCLs), e.g., see [27, 10, 21, 30, 48, 20, 26], almost all developed time marching schemes for solving the model with MCLs are first order, e.g., see [19, 5, 4, 12, 36, 1, 13, 46, 29, 62]. More precisely, to the best of the authors’ knowledge, no schemes can be claimed to posses the following three properties, namely, linear, unconditionally energy stable and second order accuracy. This is because two additional numerical difficulties, the discretization of the GNBC and DCLC conditions on the boundary, emerge besides the regular stiffness issue induced by the nonlinear double well potential in the Cahn-Hilliard equation. At the very least, even for the Cahn-Hilliard equation, the algorithm design is still challenging. To overcome the stiffness issue, many efforts had been implemented to remove the time step constraint, including, e.g., the nonlinear convex splitting approach [9, 52, 39, 47], and the linear stabilization approach [56, 42, 44, 45, 46, 3, 28, 33, 43]. About the pros and cons of these two methods, we give some detailed discussions in Section 3.

Therefore, the aim of this paper is to develop some more efficient and accurate schemes for solving the phase-field MCL model proposed in [34, 35]. We shall construct second order time stepping schemes which satisfy a discrete energy law by combining several successful approaches including the projection method for the Navier-Stokes equations to decouple the velocity and pressure, the invariant energy quadratization (IEQ) method (cf. [57, 17, 58, 55, 60]) for nonlinear gradient terms that appear in the bulk as well as the boundary of the phase field equation, and a subtle implicit-explicit treatment for the stress and convective terms. At each time step, one can solve a linear elliptic system for the phase variable and the velocity field, and a Poisson equation for the pressure. We shall give a rigorous proof of the well-posedness of the linear system together with numerical results to verify the second order accuracy in time and the efficiency.

The rest of the paper is organized as follows. In Section 2, we briefly describe the phase field model with MCL boundary conditions and its associated energy dissipation law. In Section 3, we present the numerical schemes, and prove the well-posedness of the semi-discretized linear system and their discrete energy dissipation law rigorously. In Section 4, we describe the implementation based on a Fourier-Legendre spectral-Galerkin spatial discretization. In section 5, we present various numerical examples to illustrate the accuracy and efficiency of the proposed schemes. Some concluding remarks are given in Section 6.

Consider a mixture of two immiscible, incompressible fluids in a confined domain $\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) with matched density and viscosity. We introduce a phase field variable $\phi({\boldsymbol{x}},t)$ such that

We consider the following Ginzburg-Landau free energy for the mixture

where $\lambda$ denotes rescaled characteristic strength of phase mixing energy. The first gradient term in $E_{mix}$ contributes to the hydro-philic type (tendency of mixing) of interactions between the materials and the second part, the double well bulk energy $F(\phi)=\frac{1}{4\varepsilon}(\phi^{2}-1)^{2}$, represents the hydro-phobic type (tendency of separation) of interactions. As the consequence of the competition between the two types of interactions, the equilibrium configuration will include a diffusive interface with a thickness proportional to the parameter $\varepsilon$.

The total bulk energy of the hydrodynamic system is a sum of the kinetic energy $E_{k}$ together with the mixing energy $E_{mix}$:

Here ${\boldsymbol{u}}$ is the fluid velocity field, and we assume the fluid density is $1$.

The evolution of the phase function is governed by the Cahn-Hilliard equation

where $\mu$ is the chemical potential, $M$ is a mobility parameter, and $f(\phi)=F^{\prime}(\phi)=\frac{1}{\varepsilon}\phi(\phi^{2}-1)$. The momentum equation for the hydrodynamics takes the usual form of the Navier-Stokes equation as follows (cf. e.g. [27, 21, 35, 46])

where $p$ is the pressure, $\nu$ is the kinetic viscosity of the mixture.

On the boundary $\Gamma$, we use the GNBC for the velocity as follows [34, 35],

and together with the DCLC for the phase field variable on the boundary $\Gamma$,

where $\dot{\phi}=\phi_{t}+{\boldsymbol{u}}_{\tau}\cdot\nabla_{\tau}\phi$, $\ell(\phi)\geq 0$ is a given coefficient function that is the ratio of domain length to the slip length, $\gamma$ is a boundary relaxation coefficient, ${\boldsymbol{u}}_{w}$ is the wall velocity, ${\boldsymbol{u}}_{\tau}$ is the tangential velocity along the boundary tangential direction $\tau$, $\nabla_{\tau}=\nabla-({\boldsymbol{n}}\cdot\nabla){\boldsymbol{n}}$ is the gradient along $\tau$, $g(\phi)=G^{\prime}(\phi)$ and $G(\phi)$ is the interfacial energy that is defined as

where $\theta_{s}$ is the static contact angle. From (2.8), we have ${\boldsymbol{u}}={\boldsymbol{u}}_{\tau}$ on boundary $\Gamma$.

We now describe the energy dissipation law for the above model. Here and after, for any function $f,g\in L^{2}({\Omega})$, we use $(f,g)=\int_{\Omega}f({\boldsymbol{x}})g({\boldsymbol{x}})\mathrm{d}{\boldsymbol{x}}$ to denote the $L^{2}$ inner product between functions $f({\boldsymbol{x}})$ and $g({\boldsymbol{x}})$, $(f,g)_{\Gamma}$ to denote $\int_{\Gamma}f(s)g(s)\mathrm{d}s$, and $\|f\|^{2}=(f,f)$ and $\|f\|^{2}_{\Gamma}=(f,f)_{\Gamma}$.

We now construct time marching schemes to solve the NSCH system (2.4)-(2.5)-(2.6)-(2.7) with boundary conditions of GNBC (2.8)-(2.9) and DCLC (2.10)-(2.11). With the aim of constructing schemes that are linear, second order, and unconditionally energy stable, we notice that there are several numerical challenges, including (i) how to decouple the computations of velocity and pressure; (ii) how to discretize $f(\phi)$; (iii) how to discretize $g(\phi)$; and (iv) how to develop proper discretizations for convective and stress terms.

The first difficulty (i) actually has been well studied during the last forty years, e.g., the projection type methods are one of the best ways to solve it (cf. the review in [14] and the references therein). The difficulty (ii) is also well studied recently by two class of methods where one is the nonlinear convex splitting method (cf. e.g. [9, 52, 39, 47]), and the other is the linear stabilization approach(cf. e.g. [56, 42, 44, 45, 46, 3, 28, 33, 43]). The convex splitting approach is energy stable, however, it produces nonlinear schemes at most cases, thus the implementations are often complicated and the computational costs are high. The linear stabilization approach introduces purely linear schemes, thus it is easy to implement. But, its stability usually requests a special property (generalized maximum principle)[42, 53] satisfied by the classical PDE solution and the numerical solution, that is very hard to prove in general. Recently, some theoretical progress has been made to overcome this barrier for first order [25] and second order [24] stabilization methods, using subtle Fourier analysis. However, to prove the unconditional stability, large stabilization constants are required. In this paper, we use a newly developed IEQ approach, which has been successfully applied to solve several gradient flow type models (cf. [55, 57, 17, 58, 61, 60, 59]). Its idea is to make the free energy quadratic in terms of new variables via the change of variables. Then the free energy and the PDE system are transformed into equivalent forms and thus the nonlinear terms can be treated semi-explicitly.

More precisely, we define two new variables

where $C$ is a constant to ensure $G(\phi)+C$ be positive, for instance, $C=\frac{\sqrt{2}}{3}+\eta$ with any $\eta>0$. Hence we can rewrite the bulk and total free energy as

Thus, we have an equivalent new PDE system as follows

with the GNBC on $\Gamma$ as

and DCLC as

where $Z(\phi)={g(\phi)}/{\sqrt{G(\phi)+C}}$. The initial conditions read as

We derive the energy dissipative law for this new system (3.18)-(3.27) as follows.

We emphasize that the new transformed system (3.18)-(3.27) is exactly equivalent to the original system (2.4)-(2.5)-(2.6)-(2.7), (2.8)-(2.9)-(2.10)-(2.11) that can be easily obtained by integrating (3.22) and (3.27) with respect to the time. Therefore, the energy law (3.29) for the transformed system is exactly the same as the energy law (2.14) for the original system for the time-continuous case. We will develop time-marching schemes for the new transformed system (3.18)-(3.27) that in turn follows the new energy dissipation law (3.29) instead of the energy law (2.13) for the original system.

We fix some notations here. We define two Sobolev space $H_{c}^{1}(\Omega)=\{\phi\in H^{1}(\Omega):\int_{\Omega}\phi=0\}$ and $H_{{\boldsymbol{u}}}(\Omega)=\big{\{}{\boldsymbol{u}}\in[H^{1}(\Omega)]^{d}:{\boldsymbol{u}}\cdot{\boldsymbol{n}}|_{\Gamma}=0\big{\}}$. Let $\delta t>0$ be a time step size and set $t^{n}=n\delta t$. For any function $S({\boldsymbol{x}},t)$, let $S^{n}$ denotes the numerical approximation to $S(\cdot,t)|_{t=t^{n}}$, and $S^{n+\frac{1}{2}}:=\frac{S^{n+1}+S^{n}}{2}$, $S_{\star}^{n+\frac{1}{2}}:=\frac{3}{2}S^{n}-\frac{1}{2}S^{n-1}$, $S_{\star}^{n+1}:=2S^{n}-S^{n-1}$, $\delta S^{n+1}:=S^{n+1}-S^{n}$, $\delta^{2}S^{n+1}:=S^{n+1}-2S^{n}+S^{n-1}$.