An improved threshold dynamics method for wetting dynamics

Highlights

• Total surface energy is approximated using different Gaussian kernel convolutions.

• The method is derived based on the relaxation and linearization of the energy.

• Asymptotic analysis is performed to show the consistency of the method.

• Gamma convergence of the approximate energy is rigorously proved.

• Numerical experiments show significant improvements in the accuracy.

Abstract

We propose a modified threshold dynamics method for wetting dynamics, which significantly improves the behavior near the contact line compared to the previous method (J. Comput. Phys. 330 (2017) 510–528). The new method is also based on minimizing the functional consisting of weighted interface areas over an extended domain including the solid phase. However, each interface area is approximated by the Lyapunov functional with a different Gaussian kernel. We show that a correct contact angle (Young's angle) is obtained in the leading order by choosing correct Gaussian kernel variances. We also show the Gamma convergence of the functional to the total surface energy. The method is simple, unconditionally stable with $O(NlogN)$ computational complexity per time step and is not sensitive to the inhomogeneity or roughness of the solid surface. It is also shown that the dynamics of the contact point is consistent with the dynamics of the interface away from the contact point. Numerical examples have shown significant improvements in the accuracy of the contact angle and the hysteresis behavior of the contact angle.

Introduction

Wetting describes how a liquid drop spreads on a solid surface. The study of wetting is of critical importance for many applications and has attracted much interest in the physics and applied mathematics communities [2], [11], [16], [35], [50]. The equilibrium configuration of the liquid drop can be obtained by minimizing the total interface energy:

$$\mathcal{E}={\gamma }_{LV}|{\mathrm{\Sigma }}_{LV}|+{\gamma }_{SL}|{\mathrm{\Sigma }}_{SL}|+{\gamma }_{SV}|{\mathrm{\Sigma }}_{SV}|$$

where ${\gamma }_{SV}$, ${\gamma }_{SL}$ and ${\gamma }_{LV}$ are the solid-vapor, solid-liquid and liquid-vapor surface energy densities, respectively and $|{\mathrm{\Sigma }}_{SV}|,|{\mathrm{\Sigma }}_{SL}|$ and $|{\mathrm{\Sigma }}_{LV}|$ are the corresponding interface areas. When the solid surface is homogeneous, the contact angle for a static drop is given by the famous Young's equation:

$$\cos {\theta }_Y=\frac{{\gamma }_{SV}-{\gamma }_{SL}}{{\gamma }_{LV}},$$

where ${\theta }_Y$ is the so-called Young's angle [51]. Analytic solution of the minimization problem of (1) is difficult and the numerical solution is also challenging. There have been many numerical methods proposed for simulating the free interface problem using front-tracking [25], [48], level set method [52] or the phase-field method [9], [17].

The threshold dynamics method developed by Merriman, Bence, and Osher (MBO) [29] is an efficient numerical method for the motion of the interface driven by the mean curvature. The method alternately diffuses and sharpens characteristic functions of regions and is easy to implement and highly efficient. The MBO method has been shown to converge to the continuous motion by mean curvature [3], [5], [15], [42] when the interface is away from the solid boundary. Esedoglu and Otto [13] generalized this type of method to multiphase flow with arbitrary surface tensions. The method has attracted much attention and becomes very popular due to its simplicity and unconditional stability. It has been subsequently extended to deal with many other applications including the problem of area or volume preserving interface motion [19], [21], [41], [48], image processing [12], [28], [45], problems of anisotropic interface motions [4], [10], [31], [39], generating quad mesh [43], and foam bubble problems [44]. Various algorithms and rigorous error analysis have been carried out to refine and extend the original MBO method and related methods for the aforementioned problems (see, for example, [14], [18], [26], [30], [37], [38], [40]). Some mesh free methods are also considered to accelerate this type of method [20] based on non-uniform fast Fourier transform (NUFFT) [8], [24]. Laux et al. [22], [23] rigorously proved the convergence of the method proposed in [13]. Recently, a generalized target-valued diffusion generated method was studied in [33], [34], [46], [47].

In [49], we proposed an efficient threshold dynamics method for the wetting and interface motion on the rough solid surface. The domain is extended to include the solid phase as the third phase and the method is based on the minimization of the approximate energy to (1) (as $h\to 0$)

$${\mathcal{E}}^h({\chi }_{D_1},{\chi }_{D_2})=\frac{{\gamma }_{LV}\sqrt{\pi }}{\sqrt{h}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_1}{G}_h⁎{\chi }_{D_2}\mathrm{d}\mathbf{x}+\frac{{\gamma }_{SL}\sqrt{\pi }}{\sqrt{h}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_1}{G}_h⁎{\chi }_{D_3}\mathrm{d}\mathbf{x}+\frac{{\gamma }_{SV}\sqrt{\pi }}{\sqrt{h}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_2}{G}_h⁎{\chi }_{D_3}\mathrm{d}\mathbf{x},$$

where

$$G_h(\mathbf{x})=\frac{1}{{(4\pi h)}^{n/2}}\exp (-\frac{|\mathbf{x}{|}^2}{4h})$$

is the Gaussian kernel and ${\chi }_{D_1}$ and ${\chi }_{D_2}$ are characteristic functions of domain $D_1,D_2$ in Fig. 1. An efficient iterative algorithm is then designed to find the minimizer of (3) (with volume constraints on $D_1$ and $D_2$). The method is simple, efficient, unconditionally stable and insensitive to the inhomogeneity of the solid surface. However, numerical experiments in [49] have shown that, although the apparent (macroscopic) contact angle satisfies the Young's equation, the microscopic contact angle at the contact point deviates from the correct Young's angle. There seems to be a boundary layer on the solid surface around the contact points.

In this paper, we show that the method can be improved by using heat kernel with different variances for different surface energy terms in (3), i.e.,

$${\mathcal{E}}^{h_1,h_2}({\chi }_{D_1},{\chi }_{D_2})=\frac{{\gamma }_{LV}\sqrt{\pi }}{\sqrt{h_1}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_1}{G}_{h_1}⁎{\chi }_{D_2}\mathrm{d}\mathbf{x}+\frac{{\gamma }_{SL}\sqrt{\pi }}{\sqrt{h_2}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_1}{G}_{h_2}⁎{\chi }_{D_3}\mathrm{d}\mathbf{x}+\frac{{\gamma }_{SV}\sqrt{\pi }}{\sqrt{h_2}}\underset{\tilde{\mathrm{\Omega }}}{\int }{\chi }_{D_2}{G}_{h_2}⁎{\chi }_{D_3}\mathrm{d}\mathbf{x},$$

where we use $h_1$ for approximating the liquid-vapor interface energy and $h_2$ for approximating the solid-liquid and the solid-vapor interface energy. We perform asymptotic analysis and show that, to remove the boundary layer near the contact point and obtain the correct Young's angle ${\theta }_Y$, we need to have $h_2={\lambda }^2{h}_1$ with $\lambda =\frac{\pi \cos {\theta }_Y}{\pi -2{\theta }_Y}$. We then derive the dynamic of the contact point which is consistent with the dynamic of the interface away from the contact point. We show that the improved threshold dynamics method still enjoys the energy-decaying property and is unconditionally stable. Furthermore, we also prove the Γ-convergence of the weighted functional (5) with $h_2={\lambda }^2{h}_1$ to the functional (1). This extends the analysis in [13].

This paper is organized as follows. In Section 2, we derive the modified threshold dynamics method and prove that the modified method has energy-decaying property which implies the unconditional stability. In Section 3, we use asymptotic analysis to derive the dynamic law of the contact point. In Section 4, we prove the Γ-convergence result. We present several numerical examples to verify the improvement of our modified method in Section 5. We then draw a conclusion and make some discussions in Section 6.