A level-set method for moving contact lines with contact angle hysteresis

Highlights

• Mesh-independent results can be achieved by properly choosing the friction coefficients.

• Microscopic slip length can be achieved by tuning the contact line friction.

• The contact angle hysteresis is automatically captured with a simple modification.

• Contact line is pinned when the dynamic angle is within the hysteresis window.

• Contact angle condition is easily imposed.

Abstract

We develop a level-set method in the finite-element framework. The contact line singularity is removed by the slip boundary condition proposed by Ren and E (2007) [6], which has two friction coefficients: ${\beta }_N$ that controls the slip between the bulk fluids and the solid wall and ${\beta }_{CL}$ that controls the deviation of the microscopic dynamic contact angle from the static one. The predicted contact line dynamics from our method matches the Cox theory very well. We further find that the same slip length in the Cox theory can be reproduced by different combinations of $({\beta }_N,{\beta }_{CL})$, based on which we come up with a computational strategy for mesh-independent results that can match the experiments. There is no need to impose the contact angle condition geometrically, and the dynamic contact angle automatically emerges as part of the numerical solution. With a little modification, our method can also be used to compute contact angle hysteresis, where the tendency of contact line motion is readily available from the level-set function. Different test cases, including code validation and mesh-convergence study, are provided to demonstrate the efficiency and capability of our method.

Introduction

The moving contact line problem has attracted intensive research in the past few decades due to its importance in many natural processes and industrial applications. This problem is difficult due to the stress singularity at the contact line caused by the discrepancy between the no-slip boundary condition and the moving interface. In continuum numerical simulations, different models have been adopted to relax the stress singularity, e.g., Navier slip [1], diffusion (such as those in the Cahn-Hilliard model [2], conservative level set method [3], and lattice Boltzmann method [4]), and the generalized Navier slip [5], [6]. The readers are referred to [7] for a comprehensive review on this topic and we will focus on the generalized Navier slip in this work.

From molecular dynamics (MD) simulations, Qian et al. found that the slip velocity at the wall was proportional to the sum of the tangential viscous stress and the uncompensated Young's stress (a.k.a. the unbalanced Young's stress), based on which they developed the generalized Navier boundary condition (GNBC) in the phase-field framework [5]. The velocity profiles in the vicinity of the contact line from their phase-field simulations agreed very well with the MD results. Ren and E later developed a similar slip condition [6], which is no longer restricted to the phase-field method. Their continuum modeling based on the immersed boundary method compared favorably with the MD results. Theoretically, these two slip conditions can also be derived from Onsager's minimum energy dissipation rate principle [8] or simply from the second law of thermodynamics [9]. It should be noted that Ren and E's slip condition is different from the sharp-interface limit of Qian et al.'s GNBC [10]. We follow Ren and E's slip condition in this work.

With the support from MD simulations and thermodynamic principles, the generalized Navier slip condition has gained popularity in recent years. In the phase-field community, the GNBC has been frequently adopted for contact line problems, e.g., [11], [12], [13], [14], [15]. Meanwhile, the GNBC has also been adopted in many other numerical methods for interfacial flows. For example, Gerbeau and Lelièvre incorporated the GNBC into a variational arbitrary Lagrangian-Eulerian (ALE) formulation which is well suited for energy stability analysis [16]. Li et al. developed an efficient augmented immersed interface method to implement Ren and E's slip condition [17]. Ren and E applied their slip condition to the level-set method and investigated contact line dynamics on heterogeneous surfaces [18]. This level-set work was later extended to moving contact lines with insoluble surfactants [19]. Recently, Zhang and Ren also investigated the influence of viscoelasticity on contact line dynamics using an immersed boundary method combined with the generalized slip condition [20]. The implementation of the GNBC in the front-tracking method can be found in [21]. Most recently, the GNBC was also extended to the volume-of-fluid method [22], [23].

Although the generalized Navier slip has been widely used, it is still challenging to obtain mesh-independent results, because the physical slip length is usually at the nanoscale and cannot be resolved by the computational mesh. It has been shown that Ren and E's slip condition cannot remove the weak singularity at the contact line [24], [25]. Furthermore, Ren and E's slip condition introduces two friction coefficients (${\beta }_{CL}$ and ${\beta }_N$ in Section 2.1) rather than a single slip length, and it is unclear how to choose them for predictive simulations that can match the experiments. In this work, we aim to address these issues based on a level-set method. In the literature, a standard treatment to remove mesh dependency, as proposed in [26], [27], [28], [29], is to determine a numerical contact angle at the grid scale based on macroscale models such as the Cox-Voinov model [30], [31]; this numerical angle is then applied at the contact line in place of the static contact angle. A drawback of this method is that it requires the input of contact line velocity, which may be difficult to obtain, especially in three dimensions. The similar idea was also used in the GNBC, however, in a different flavor [21], [32]: the grid-scale contact angle from the simulation is used to determine a microscopic dynamic contact angle, which is then fed to the GNBC to compute the slip velocity. In this work, we propose a different approach which does not rely on hydrodynamic models and is thus much easier to implement. Meanwhile, by properly choosing the friction coefficients, we will show that Ren and E's slip condition itself is sufficient to reproduce the well-established Cox theory [31] with realistic slip lengths.

Another challenging issue is the contact angle hysteresis, since most solid surfaces are intrinsically rough or chemically heterogenous. In this case, the contact line stays pinned when the microscopic dynamic contact angle is between a receding contact angle ${\theta }_R$ and an advancing contact angle ${\theta }_A$. The most popular approach for contact angle hysteresis was developed by Spelt for a level-set method [33]. An intermediate contact angle is obtained such that the contact line is pinned. If this angle is within the hysteresis window, the solution is accepted; otherwise, the solution is abandoned and the contact line is moved with prescribed contact angles. This idea was later extended to different methods, e.g., the phase-field method [34], the volume-of-fluid method [35], the Lattice Boltzmann method [36], and the front-tracking method [37]. However, this approach relies on ghost cells outside the boundary to pin the contact line or to impose the contact angle condition, which can be challenging on curved boundaries and unstructured meshes. Recently, we developed a thermodynamically consistent phase-field model for contact angle hysteresis [38]. Since the dynamic contact angle is part of the solution instead of being imposed, this method is easy to implement and automatically captures the pinning, advancing, and receding of the contact line. Motivated by [38], we will show that Ren and E's slip condition can also be easily modified to capture the contact angle hysteresis.

The rest of this paper is organized as follows. We first introduce the governing equations and numerical methods in Section 2. We then explain how to incorporate contact angle hysteresis in Section 3. The numerical results, including code validation and mesh convergence studies, are given in Section 4.