---
abstract: 'In this work, an analytical Volume Of Fluid (VOF) implementation of the Generalized Navier Boundary Condition is presented based on the Brackbill surface tension model. The model is validated by simulations of droplets on a smooth surface in a planar geometry. Looking at the static behavior of the droplets, it is found that there is a good match between the droplet shape resolved in the simulations and the theoretically predicted shape for various values of the Young’s angle. Evaluating the spreading of a droplet on a completely wetting surface, the Voinov-Tanner-Cox law ($\theta \propto {\text{Ca}}^{1/3}$) can be observed. At later times scaling follows $r \propto t^{1/2}$, suggesting spreading is limited by inertia. These observations are made without any fitting parameters except the slip length.'
author:
- 'A.M.P. Boelens'
- 'J.J. de Pablo'
title: 'Generalized Navier Boundary Condition for a Volume Of Fluid approach using a Finite-Volume method.'
---

Introduction
============

For any multi-phase flow over a wall where dissipation at the contact line becomes of the same order as the bulk dissipation, a good understanding of contact line behavior is essential [@squires2005]. However, despite many investigations [@dussan1979; @degennes1985; @deconinck2008; @eral2013], for applications including adhesion of liquids to a solid surface [@furmidge1962], transport of liquid water in fuel cells [@zhu2008], liquid infused surfaces [@wexler2015], and coating [@vandre2013], contact line behavior is still not well understood.

The physics of a static contact angle between a liquid and a gas on a smooth solid surface is well established [@young1805; @gauss1830]. However, real surfaces are not completely smooth. They are not chemically homogeneous, and/or have roughness. This causes static contact angle hysteresis and contact line pinning, both of which are difficult to model. When looking at a moving contact line instead of a static contact angle things get even more complicated. Dynamic contact angle behavior is not well understood, even on a completely smooth solid surface. The origin of this poor understanding of the moving contact line is twofold: there is the contact line singularity problem, and the question of how the contact angle depends on contact line velocity.

The contact line singularity problem was first identified by @huh1971. While normally it is a good approximation to use the no-slip condition as boundary condition on the wall, for corner flow this assumption causes the viscous stress and pressure to scale as $r^{-1}$ and thus to diverge as $r \to
0$ at the contact line. Numerous methods have been proposed to resolve this discontinuity or work around it. @hocking1977 showed that, using domain perturbation method in cylindrical coordinates, any slip-velocity model [@navier1823; @thompson1997] resolves the velocity singularity. Another method to circumvent the contact line singularity has been to use the Cahn-Hilliard-van der Waals model [@cahn1977] as the basis for either diffusive interface models [@seppecher1996; @jacqmin2000], or for precursor film models [@hervet1984; @degennes1985]. Precursor film models use a disjoining pressure [@derjaguin1955; @teletzke1987; @herring2010] to model the van der Waals forces that cause the formation of a precursor film ahead of the interface, removing the singularity. More exotic models have suggested local shear-thinning [@cox1998] and non-constant surface tension [@shikhmurzaev1993; @shikhmurzaev1997] as possible solutions.

For the question of how the contact angle depends on the contact line velocity there are also various models. Typically a distinction is made between the local microscopic contact angle at the contact line and the macroscopic apparent contact angle, which is observed in experiments. Due to experimental limitations to access sufficiently small length scales, the apparent contact angle is measured away from the contact line, and curvature of the interface causes this angle to be different from the microscopic contact angle [@chen1989]. Arguably, the easiest method to define the dynamic microscopic contact angle is to simply assume it is fixed and the same as the static contact angle [@eggers2004a]. Instead of a fixed dynamic contact angle, Molecular Kinetic Theory (MKT) [@blake1969; @blake2006] predicts a microscopic contact angle which changes with the contact line velocity. The Voinov-Tanner-Cox [@voinov1976; @tanner1979; @cox1986] law describes the relation between the apparent contact angle and the microscopic contact angle. This law is based on the assumption that it is possible to choose a length scale arbitrary close to the contact line. This makes it impossible to identify a characteristic length scale of the contact line geometry, and reduces its physics to a balance between capillary and viscous stresses [@kafka1979]. Using the lubrication approximation one can now derive the @voinov1976 equation for some specific asymptotic limits, and matching solutions of the Voinov equation with the mesoscopic hydrodynamic solution further away from the wall then gives the Cox-Tanner-Voinov law [@voinov1976; @tanner1979; @cox1986].

Although there are a lot of different models to describe contact lines, there is no consensus on what is the correct description of the physics. Many models have multiple fitting parameters which can be tuned to give the same results [@sibley2014]. Because impurities and surface heterogeneities have a large effect on measurements, experiments are very difficult to reproduce. On top of this, one needs access to microscopic length scales to get to the details and the outcome of experiments on a macroscopic level only depends very weakly on these small length scales [@bonn2009]. With the advent of Molecular Dynamics (MD) [@koplik1988; @koplik1989; @thompson1989; @thompson1990] simulations, and new experimental techniques, such as Atomic Force Microscopy (AFM) [@checco2006], it is now possible to access length and time scales at the contact line that were not accessible before. While MD simulations still can only probe small systems for short times, a couple of fundamental discoveries have been made. Although it had been argued that continuum models break down all together at the contact line [@dussan1979], it was found that both the Navier-Stokes equations and Young’s equation hold up even down to the nanometer scale [@bocquet2010; @seveno2013]. Furthermore, support was found for contact line slip [@thompson1989], precursor films [@he2003], and non-static dynamic microscopic contact angles [@shen1998]. Apart from simulations, precursor films have also been found experimentally [@hardy1919; @kavehpour2003].

In addition to the experimental difficulties and often contradicting findings, another reason that there is no consensus is that all of the above mentioned models have some fundamental short comings and/or work in different regimes. Further analysis of contact line slip models by @huh1977 showed that, even though stresses are not diverging anymore, the pressure still shows a weak singularity i.e. the pressure diverges, but becomes finite when integrated. While using a slip length model is sufficient from a modeling perspective, a divergent pressure cannot be right from a physics perspective. The precursor film on the other hand does successfully work around the contact line singularity for both viscous stresses and pressure. However, they are typically not seen under partial wetting conditions [@snoeijer2013]. The different models for microscopic and macroscopic contact angles have some limitations too. While a model which describes the contact angle as a function of the capillary number might describe dynamic contact angle hysteresis correctly, there is also static contact angle hysteresis, which causes contact line pinning of non-moving contact lines. Any model that lets the contact angle only depend on the local capillary number will not be able to capture static contact angle hysteresis. Assuming a static angle ,on the other hand, does not capture static contact angle hysteresis, is not able to properly describe the flow of a liquid over a chemically patterned surface with different wetting properties, and does not accurately predict contact line velocity [@yamamoto2013].

In this work, a validation study of a Volume Of Fluid (VOF) implementation of a contact line model called the Generalized Navier Boundary Condition (GNBC) is presented. While the non-sharp interface of the VOF method implicitly resolves the contact line singularity problem, even with a no-slip condition at the wall [@seppecher1996; @chen2000; @sibley2013], the question of what is the right contact angle is still valid for this method. In addition to applying the Navier slip condition at the wall, the GNBC uses the reduced Young’s stress as a restoring force when the contact angle deviates from its equilibrium value, and is informed by Molecular Dynamics [@qian2003; @gentner2004]. While previous implementations used a friction factor to link the reduced Young’s stress to a contact line velocity [@gerbeau2009], this approach does not work for the VOF model, because of the implicit slip of the interface [@ashish2009]. Instead, in this work, the reduced Young’s stress is used directly as a source term in the navier stokes equations [@mahady2015]. A consequence of using the GNBC is that the contact angle no longer is a constraint imposed onto the system, but that the contact angle is a self-selecting variable. Using this approach has a couple of advantages over existing models. Because there is no enforcement of a model that relates the contact angle to contact line speed, this model can reproduce static contact angle hysteresis, without artificially fixing the position of the contact line [@fang2008; @dupont2010]. Additionally, a variable contact angle model can describe flow over a chemically patterned surface [@dupuis2004; @wang2008]. While the Generalized Navier Boundary Condition has been implemented for continuum simulations using a diffuse interface Cahn Hillard method [@qian2006; @ren2007; @wang2008; @cai2015], using a Volume Of Fluid approach has the advantage that many less grid points are needed to resolve the interface [@ding2007]. In contrast to the body force term derived by @mahady2015, the model presented in this work describes wetting of a dry surface. Opposed to @deganello2011, the model presented here is a analytical derivation of the line tension force.

In our work, we uncover the following findings: when looking at the static behavior of droplets, it is found there is a good match between the droplet shape resolved in the simulations and the theoretically predicted shape for various values of the Young’s angle. Investigating the spreading of a droplet on a completely wetting surface, the Voinov-Tanner-Cox law [@voinov1976; @tanner1979; @cox1986] ($\theta \propto {\text{Ca}}^{1/3}$) can be observed. Late time scaling follows $r \propto t^{1/2}$, suggesting spreading is limited by inertia. These observations are made without any fitting parameters except the slip length.

Theory
======

Body force reformulation of uncompensated Young’s stress
--------------------------------------------------------

The traditional form of the Generalized Navier Boundary Condition (GNBC) used in, for example, Arbitrary Lagrangian-Eulerian simulations [@gerbeau2009], looks like: $$\beta {\left( \vec{u} - \vec{u}_{w} \right)} \cdot \tau
+ \mu {\left( \nabla \vec{u} + \nabla \vec{u}^{T} \right)} \hat{n}_{w} \cdot \tau
+ \sigma {\left( \cos{\theta} - \cos{\theta_{0}} \right)} \hat{t}_{w} \cdot \tau \delta_{{\text{cl}}}
=
  0
\label{eqn:gnbc}$$ Here the first two terms on the left hand side represent the Navier slip boundary condition [@navier1823], and the third term represents the unbalanced Young’s stress. $\beta$ is the slip coefficient, $\vec{u}$ is velocity, $\vec{u}_{w}$ is the velocity of the wall, $\tau$ is any vector tangent to the wall, $\mu$ is the dynamic viscosity, $\hat{n}_{w}$ is the normal of the wall pointing outward, $\sigma$ is the surface tension coefficient, $\theta$ is the dynamic contact angle, $\theta_{0}$ is the equilibrium Young’s contact angle, $\hat{t}_{w}$ is the vector tangent to the wall and normal to the contact line, and the distribution $\delta_{{\text{cl}}}$ is defined as [@gerbeau2009]: $${\left<{\delta_{{\text{cl}}}, \psi}\right>}
=
  \int_{{\text{cl}}} \psi d \lambda_{{\text{cl}}}$$ where $\psi$ is any smooth function, and $\lambda_{{\text{cl}}}$ denotes the Lebesgue measure (i.e. the length measure) on the contact line. The slip length in equation \[eqn:gnbc\] is equal to $l_{s} = \mu / \beta$.

Analogous to the methods of @brackbill1992 a Volume Of Fluid expression is derived for the above equation for the Generalized Navier Boundary Condition. However, instead of relating the uncompensated Young’s stress to a velocity using a slip coefficient, the uncompensated Young’s stress is modeled as an extra body force acting at the contact line in the Navier-Stokes equation. The reason for this approach is that the intrinsic contact line slip in a VOF code is large enough that converting the contact line tension directly to a velocity on the wall does not move the contact line. The Navier-slip boundary condition component of the above equation is left unchanged in our implementation of the Generalized Navier Boundary Condition.

The goal of this derivation is to find an expression that rewrites the contact line force as a body or volume force, so it can be treated as a momentum source term in the Navier-Stokes equations. The first step is to find a relation, that rewrites the line force $\vec{F}_{\tau L}$ as surface force $\vec{F}_{\tau A}$ acting on the wall: $$\lim_{h \to 0} \int_{\Delta A_{w}} \vec{F}_{\tau A} {\left( \vec{x} \right)} d^{2} \vec{x}
=
  \int_{\Delta L} \vec{F}_{\tau L} {\left( \vec{x}_{L} \right)} d \vec{x}
\label{lineForce}$$ The points, $\vec{x}_{L}$, form the contact line or triple point, $\Delta L$ is the length of the line segment in a the small volume of integration $\Delta V$, and $\Delta A_{w}$ is the side of $\Delta V$ that is part of the wall, and in which the points $\vec{x}_{L}$ lay. An additional constraint for $\vec{F}_{\tau
A} {\left( x \right)}$ is that it is zero outside of the interface region: $$\vec{F}_{\tau A} {\left( x \right)}
=
  0
{\text{ for }}
  {\left|{
    \hat{n}_{2D} {\left( \vec{x}_{L} \right)}
    \cdot
    {\left( \vec{x} - \vec{x}_{L} \right)}
  }\right|}
\geq
  h$$ where $\hat{n}_{2D} {\left( \vec{x}_{L} \right)}$ is the normal to the contact line in the plane of the wall, and $\delta {\left[ \hat{n}_{2D} {\left( \vec{x}_{L} \right)} \cdot
{\left( \vec{x} - \vec{x}_{L} \right)} \right]}$ describes the plane $\Delta A_{w}$.

Consider a system of two fluids, fluid $1$, and fluid $2$, separated by an interface, and define a discontinuous function, $c {\left( \vec{x} \right)}$, to distinguish between the two phases: $$c {\left( \vec{x} \right)}
=
  \left\{
    \begin{array}{ll}
      c_{1}                           & {\text{in liquid 1}}  \\
      {\left<{c}\right>} = {\left( c_{1} + c_{2} \right)}/2 & {\text{on interface}} \\
      c_{2}                           & {\text{in liquid 2}}
    \end{array}
  \right.$$ An example of such a function would be the density of two different incompressible liquids, $\rho_{1}$, and $\rho_{2}$. In this case the position of the contact line can be found with: $$\rho {\left( \vec{x}_{L} \right)}
=
  {\left<{\rho}\right>}$$ In order to change this problem from a boundary value problem on the contact line to an approximate continuous model, one can define a continuous function $\tilde{c}$, which varies smoothly over thickness $h$ going from $c_{1}$ to $c_{2}$ , $c_{1} \le \tilde{c} \le c_{2}$. $h$ is a length scale of the order of the grid size $\Delta x$, and defines the width of the transition region from $c_{1}$ to $c_{2}$. The two functions $c$ and $\tilde{c}$ are related through the interpolation function $\mathcal{P} {\left( \vec{x} \right)}$: $$\tilde{c} {\left( \vec{x} \right)}
=
  \frac{1}{h^{2}} \int_{A} c {\left( \vec{x}' \right)} \mathcal{P} {\left( \vec{x}'-\vec{x} \right)} d^{2} \vec{x}$$ which is normalized as: $$\int_{A} \mathcal{P} {\left( \vec{x} \right)} d^{2} \vec{x}
=
  h^{2}$$ is bounded as: $$\mathcal{P} {\left( \vec{x} \right)}
=
  0 
  {\text{ for }}
  {\left|{\vec{x}}\right|}
\ge
  h/2$$ is differentiable, and decreases monotonically with ${\left|{\vec{x}}\right|}$. The continuous function is defined such that: $$\lim_{h \to 0} \tilde{c} {\left( \vec{x} \right)}
= 
  c {\left( \vec{x} \right)}$$ i.e. the function $\tilde{c} {\left( \vec{x} \right)}$ approaches $c {\left( \vec{x} \right)}$ as the interface thickness goes to zero. $\tilde{c}$ is differentiable because $\mathcal{P}$ is, and: $$\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}
=
  \frac{1}{h^{2}} \int_{A} c {\left( \vec{x}' \right)} \nabla_{2D} \mathcal{P} {\left( \vec{x}'-\vec{x} \right)} d^{2} \vec{x}$$ where $\nabla_{2D}$ is the two-dimensional gradient in the plane of the wall. Using Gauss’ theorem and the realization that $c {\left( \vec{x} \right)}$ is constant within each fluid, the above integral can be written as: $$\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}
=
  \frac{{\left[ c \right]}}{h^{2}} \int_{L} \hat{n}_{2D} {\left( \vec{x}_{L} \right)} \mathcal{P} {\left( \vec{x}-\vec{x}_{L} \right)} d L
\label{gauss}$$ where ${\left[ c \right]} = c_{2} - c_{1}$, thus converting the surface integral to a line integral. To pull the normal out of the integral, its weighted mean is calculated. Since $\mathcal{P}$ is bounded, its maximum contribution to the line integral is $\mathcal{O} {\left( h \right)}$. Integral \[gauss\] can thus be approximated as: $$\frac{1}{h^{2}} \int_{L} \hat{n}_{2D} {\left( \vec{x}_{L} \right)} \mathcal{P} {\left( \vec{x}-\vec{x}_{L} \right)} d L
\approx
  \frac{1}{h^{2}} \hat{n}_{2D} {\left( \vec{x}_{L0} \right)} \int_{L} \mathcal{P} {\left( \vec{x}-\vec{x}_{L} \right)} d L
+
  \mathcal{O} {\left( \frac{h}{R} \right)}
\label{approximate}$$ where $\vec{x}_{L0}$ is the point on $L$ closest to $x$, and $R$ is the radius of the contact line.

The integral in equation \[approximate\] can be bounded by: $$\frac{1}{h} \int_{L} \mathcal{P} {\left( \vec{x}-\vec{x}_{L} \right)} d L
\le
  \mathcal{P} {\left( \vec{x}-\vec{x}_{L0} \right)}$$ where in the limit $h \to 0$, $\mathcal{P} {\left( \vec{x}-\vec{x}_{L0} \right)}$ is zero everywhere except for $\vec{x} = \vec{x}_{L0}$. Taking the corresponding limit of $\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}$ over the interface gives: $$\lim_{h \to 0} \int \hat{n}_{2D} {\left( \vec{x}_{L0} \right)} \cdot \nabla_{2D} \tilde{c} {\left( \vec{x} \right)} dx
= 
  {\left[ c \right]}$$ As $h$ goes to $0$, $\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}$ is thus equivalent to: $$\lim_{h \to 0} \nabla_{2D} \tilde{c} {\left( \vec{x} \right)}
=
  \hat{n}_{2D} {\left[ c \right]} \delta_{2D} {\left[ \hat{n}_{2D} \cdot {\left( \vec{x} - \vec{x}_{L} \right)} \right]} 
=
  \nabla_{2D} c {\left( \vec{x} \right)}
\label{delta}$$

Because the Brackbill surface tension already takes care of the $\cos(\theta)$ component of the reduced Young’s stress, only the $\cos{\theta_{0}}$ component needs to be modeled. Using the delta function, the contact line force can now be written as a surface force as follows: $$\begin{aligned}
  \int_{L} \vec{F}_{\tau L} {\left( \vec{x}_{L} \right)} d L 
&= 
  \int_{A} \vec{F}_{\tau L} {\left( \vec{x} \right)} \delta_{2D} {\left[ \hat{n}_{2D} {\left( \vec{x}_{L} \right)} \cdot {\left( \vec{x} - \vec{x}_{L} \right)} \right]} d^{2} \vec{x} \nonumber \\
&= 
  \int_{A} \sigma \cos{\theta_{0}} \hat{n}_{2D} {\left( \vec{x} \right)} \delta_{2D} {\left[ \hat{n}_{2D} {\left( \vec{x}_{L} \right)} \cdot {\left( \vec{x} - \vec{x}_{L} \right)} \right]} d^{2} \vec{x}\end{aligned}$$ converting the line integral over the contact line into an integral over the wall surface. Equation \[delta\] can be used as an approximation for the delta function when the interface has a finite thickness. Substitution gives: $$\int_{\Delta L} \vec{F}_{\tau L} {\left( \vec{x}_{L} \right)} d L 
=
  \lim_{h \to 0} \int_{\Delta A_{w}} \sigma \cos{\theta_{0}} \frac{\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}}{{\left[ c \right]}} d^{2} \vec{x}
\label{surfaceForce}$$ and comparing equation \[lineForce\] with \[surfaceForce\], the surface force $\vec{F}_{\tau A}$ can be identified as: $$\vec{F}_{\tau A}
=
  \sigma \cos{\theta_{0}} \frac{\nabla_{2D} \tilde{c} {\left( \vec{x} \right)}}{{\left[ c \right]}}$$ As a last step, the surface integral needs to be converted to a volume integral. Since $\nabla_{2D} \tilde{c}$ is independent of the distance away from the wall, integrating over $\Delta V$ is the same as multiplying $\vec{F}_{\tau A}$ with mesh size, i.e. $\Delta A_{w} / \Delta V$: $$\vec{F}_{\tau V}
=
  \vec{F}_{\tau A} \frac{\Delta A_{w}}{\Delta V}
=
  \sigma \cos{\theta_{0}} \frac{\nabla_{2D} \tilde{c}
  {\left( \vec{x} \right)}}{{\left[ c \right]}} \frac{\Delta A_{w}}{\Delta V}
\label{eqn:volumeForce}$$ where $\Delta A_{w}$ is the surface area of the wall in volume $\Delta V$.

While the above equation has the correct limiting behavior, it was found that, due to the diffuse nature of the interface, at too low resolution the interface gets pulled apart. To counter this phenomenon, a modification is proposed to localize the contact line force more to the interface: $$g {\left( \tilde{c} \right)}
=
  H {\left( \pi/2-\theta_{0} \right)} \frac{5}{4} H {\left( \tilde{c}-\frac{1}{5} \right)}
+ H {\left( \theta_{0}-\pi/2 \right)} \frac{5}{4} H {\left( \frac{4}{5}-\tilde{c} \right)}$$ Where $H(x)$ is the Heaviside step function. Depending on the value of the Young’s angle, this function truncates the contact line force. For $\theta_{0} <
\pi/2$ this means that the contact line does not get pulled apart, while for $\theta_{0} > \pi/2$ this prevents the gas phase from being pulled into the droplet at the contact line. The value $5/4$ is to normalize the function. The new function for the contact line restoring force at low resolution thus becomes: $$\vec{F}_{\tau V}
=
  \vec{F}_{\tau A} \frac{\Delta A_{w}}{\Delta V} g {\left( \tilde{c} \right)}
\label{eqn:volumeForceCor}$$ where the same definition is used for $\vec{F}_{\tau A}$ as in equation \[eqn:volumeForce\].

Numerical implementation
------------------------

This section focuses on the numerical implementation of the above derived equation in the Volume Of Fluid (VOF) solver that comes with OpenFOAM [@jasak1996; @openfoam]. This involves implementing the uncompensated Young’s stress and the Navier slip boundary condition. In this code the general phase parameter $\tilde{c} {\left( \vec{x} \right)}$ is called $\alpha$, and has the following properties: $$\alpha
=
  \left\{
  \begin{array}{ll}
  0      & {\text{in phase 1}}    \\
  (0, 1) & {\text{on interface}}    \\
  1      & {\text{in phase 2}}
  \end{array}
  \right.$$ Phase parameter $\alpha$ is stored as a separate field, just like velocity and pressure, and its evolution is calculated using the following transport equation: $$\frac{\partial \alpha}{\partial t}
+ \nabla \cdot {\left( \alpha \vec{v} \right)}
+ \nabla \cdot {\left( \alpha {\left( 1 - \alpha \right)} \vec{v}_{lg} \right)}
=
  0$$ where $\vec{v} = \alpha \vec{v}_{l} + {\left( 1 - \alpha \right)} \vec{v}_{g}$ is the phase averaged velocity, and $\vec{v}_{lg} = \vec{v}_{l} - \vec{v}_{g}$ is the velocity difference between the liquid and the gas phase. This equation is equivalent to a material derivative, but rewritten to minimize numerical diffusion [@rusche2002].

The volume fraction is used to calculate phase-averaged densities, velocities, and viscosities, which are used in the momentum balance $$\frac{\partial \rho \vec{v}}{\partial t}
+ \nabla \cdot {\left( \rho \vec{v} \otimes \vec{v} \right)}
=
- \nabla p
+ \nabla \cdot {\left( \mu \nabla \vec{v} \right)}
+ \rho \vec{g}
+ \vec{f}_{{\text{st}}}
- \vec{f}_{{\text{cl}}}$$ and the continuity equation: $$\nabla \cdot \vec{v}
=
  0$$ In the above equations $\rho$ is the density, $\vec{v}$ is the velocity, $t$ time, $p$ is pressure $\mu$ is the viscosity, $g$ is gravity, $\vec{f}_{{\text{st}}}$ is the surface tension force, $\vec{f}_{{\text{cl}}}$ is the contact line tension force, and $\otimes$ is the dyadic product. The density $\rho$, velocity $\vec{v}$, and viscosity $\mu$, are all phase averaged using $\alpha$.

The surface tension force is calculated using the expression: $$\vec{f}_{{\text{st}}}
=
  \sigma_{{\text{st}}} \kappa \nabla{\alpha}$$ where $\sigma_{{\text{st}}}$ is the surface tension coefficient, $$\kappa
=
- {\left( \nabla \cdot \vec{n} \right)}$$ is the curvature of the interface, and $$\vec{n}
=
  \frac{\nabla \alpha}{{\left|{\nabla \alpha}\right|}}$$ is the normal of the interface [@brackbill1992].

Using $\alpha$, the line tension force defined in equation \[eqn:volumeForceCor\] is rewritten as: $$\vec{f}_{{\text{cl}}}
=
  \sigma_{{\text{st}}} \cos{\theta_{0}} \;
  g {\left( \alpha \right)}
  \frac{\Delta A_{w}}{\Delta V}$$ and $\nabla_{2D} \alpha$ is implemented as: $$\nabla_{2D} \alpha 
=
  \nabla \alpha - {\left( \hat{n}_{w} \cdot \nabla \alpha \right)} \hat{n}_{w}$$ where $\hat{n}_{w}$ again is the normal of the wall pointing outward. The equilibrium Young’s angle $\theta_{0}$ can be defined uniquely for any grid cell along the wall, so an arbitrary wettability pattern can be created. $\Delta A_{w}$ and $\Delta V$ are properties of the mesh that can be accessed directly in OpenFOAM. The contact line tension source term is solved explicitly along with the surface tension source term.

The Navier slip condition is implemented as: $$\vec{v}
=
  l_{s} {\left( {I} - \hat{n}_{w}^{2} \right)} \nabla \vec{v}$$ where $l_{s}$ is the slip length ${I}$ is the identity matrix, and $\hat{n}_{w}^{2}$ is the dyadic product of the wall normal, $\hat{n}_{w}$, with itself. Using this formulation, the velocity perpendicular to the wall is always set to zero.

Results
=======

The first section of the results focuses on the static behavior of the implementation of the Generalized Navier Boundary Condition presented above. The second section covers the dynamic contact angle behavior. To speed up simulation times all results are for 2D droplets in a planar geometry, and using the symmetry of the system only half of the droplet is simulated. All simulations are performed without gravity and represent water droplets in air at room temperature.

Static
------

![The interface of various droplets at time $0.1 \si{\s}$ \[figure:contourEquil\]](fig1){width="60.00000%"}

Figure \[figure:contourEquil\] shows the interface ($\alpha = 0.5$) of various droplets with different Young’s angles at a resolution of $512 \times 256$ in a box of $1.5 \si{\milli\m} \times 0.75 \si{\milli\m}$. The half droplets in the figure have a surface area of about $0.2 \si{\milli\m^{2}}$ (i.e. a radius of $0.5 \si{\milli\m}$ when the Young’s angle is $\pi/2$). As initial condition these droplet where given their equilibrium shape, and it can be seen that after a simulation time of $0.1 \si{\s}$ the droplets have maintained their shape.

![Radius of the base of the droplet and droplet height as a function of different contact angles.\[figure:angleRadiusEquil\]](fig2){width="60.00000%"}

To further quantify these droplets, figure \[figure:angleRadiusEquil\] shows the height and radius of the base of the droplets as a function of the Young’s angle. The radius of the base of the droplet $r_{B}$ is calculated as: $$r_{b} =  R \sin{\theta_{0}},$$ and the height of the droplet is equal to: $$h_{d} = R {\left( 1 - \cos{\theta_{0}} \right)}.$$ In the above equations: $$R = \sqrt{\frac{2 A}{\theta_{0} - \cos{\theta_{0}} \sin{\theta_{0}}}}$$ is the radius of the droplet, $\theta_{0}$ is the Young’s angle, and $A$ is the surface area of the droplet. While there are small deviations between the theoretically predicted droplet shapes and the simulations, both match well.

![Error in the radius of the base of the droplet as a function of resolution for a contact angle of $\theta_{0} = \pi/4$. \[figure:resolution\]](fig3){width="60.00000%"}

The convergence of the error as a function of resolution is shown in figure \[figure:resolution\]. The x axis shows the resolution. The values of the time steps at these resolutions are: $0.1 \si{\mu \s}$, $0.05 \si{\mu \s}$, $0.025 \si{\mu \s}$, and $0.01 \si{\mu
\s}$, to keep the Courant number of the same order between simulations. The y axis shows the absolute error, $\epsilon_{{\text{abs}}}$, between the theoretical value for the radius of the base of the droplet and the simulated value for a Young’s angle of $\theta_{0} = \pi/4$. The graph shows that the code is close to 2nd order accurate.

![\[figure:pressureEquil\] Pressure at the wall at $t = 50 \si{\milli
\s}$ for different resolutions and a contact angle of $\theta_{0} = \pi/4$](fig4){width="60.00000%"}

Because the interface gets thinner as resolution increases the Volume Of Fluid method is mesh dependent. As can be seen in figure \[figure:pressureEquil\], this results in the pressure peak on the surface getting sharper with increasing resolution. The increased pressure inside the drop, left of the pressure peak, is the Laplace pressure.

![Maximum of pressure peak at the contact line as a function of different contact angles.\[figure:anglePresEquil\]](fig5){width="60.00000%"}

The reason there is a pressure peak at the contact line for a stationary droplet in the first place can be seen in figure \[figure:anglePresEquil\]. This figure shows the value of the pressure peak as a function of the Young’s angle. From the figure it is clear that there is a strong dependence of the pressure peak on the Young’s angle. Since the line tension force only acts parallel to the surface, this suggest that the pressure peak is the result of the surface tension force calculated by the @brackbill1992 model. If the vertical component of the reduced Young’s stress also were implemented, an extra term, proportional to $\sigma \sin{\theta_{0}}$, would have been present in the plot as an additional contribution to the pressure peak. Since we are only concerned about solid surfaces, this term was not incorporated in the model.

![\[figure:velocityEquil\] Velocity at the wall at $t = 50 \si{\milli \s}$ for different resolutions](fig6){width="60.00000%"}

Due to the large pressure peak there is also a small residual slip velocity at the contact line, which can be observed in figure \[figure:velocityEquil\]. As was the case with the pressure, this velocity peak gets sharper at higher resolutions. As can be seen in figure \[figure:angleRadiusEquil\] the residual slip velocity does not negatively affect the shape of the drop. However, spurious currents are a well known issue of Volume Of Fluid solvers [@deshpande2012] and, if needed, can be controlled by making the time step sufficiently small.

![Integral of pressure peak at the contact line as a function of resolution.\[figure:resolutionPressure\]](fig7){width="60.00000%"}

![Integral of slip velocity at the contact line as a function of resolution.\[figure:resolutionVelocity\]](fig8){width="60.00000%"}

To investigate to what extent the pressure and velocity peaks affect the solution they are both integrated over the wall for various resolutions. Figure \[figure:resolutionPressure\] shows this integral for the pressure, and figure \[figure:resolutionVelocity\] for the slip velocity. The integration limits in both plots are from $\alpha = 0.01$ to $\alpha = 0.99$. For the integral over the pressure peak this approach makes sure that the integral is not affected by the Laplace pressure inside the droplet, and the velocity integral uses the same limits to be consistent with the pressure. It can be appreciated how the pressure integral converges to a constant value and the velocity integral approaches zero as resolution increases, showing a convergent solution for both the pressure and velocity.

Dynamic
-------

For the validation of the dynamic case, the starting point is again a 2D droplet in planar geometry. However, in this case, the equilibrium Young’s contact angle is set to $\theta_{0} = 0$, and the spreading behavior as a function of time is investigated.

![\[figure:contourTanner\] Interface of a droplet at different times at a resolution of $512 \times 256$.](fig9){width="60.00000%"}

Figure \[figure:contourTanner\] shows the interface ($\alpha = 0.5$) at various times for a box of $1.5 \si{\milli\m} \times 0.75 \si{\milli\m}$ and with a resolution of $512 \times 256$. As expected the droplet keeps spreading until the edge of the simulation box is reached and the simulation is stopped.

![\[figure:diameterTanner\] Radius of a droplet as a function of time.](fig10){width="60.00000%"}

Figure \[figure:diameterTanner\] shows how the droplet spreads by showing the radius of the droplet as a function of time for two different resolutions of $256
\times 128$ and $512 \times 256$. For lower resolutions the contact line was pulled apart and the results are not shown in the graph. For reference, the power law $r \propto t^{1/2}$ is also shown. The $r \propto t^{1/2}$ power law describes the late time spreading behavior for low viscosity axisymmetric droplets [@biance2004; @ding2007b]. At a resolution of $256 \times 128$ initially the contact line hardly moves, but then the spreading radius as a function of time becomes proportional to: $r \propto t^{1/2}$. At the larger resolution of $512 \times
256$ it can be seen that the scaling also converges to $r \propto t^{1/2}$.

![\[figure:pressureTanner\] Pressure at the wall for various times at a resolution of $512 \times 256$](fig11){width="60.00000%"}

.

![\[figure:velocityTanner\] Velocity at the wall for various times at a resolution of $512 \times 256$](fig12){width="60.00000%"}

.

Figures \[figure:pressureTanner\] and \[figure:velocityTanner\] show the pressure profile and velocity at the wall for various times. The pressure shows a sharp peak already, but especially at $t=0.5 \si{\milli\s}$ there is quite some noise in both the curves for pressure and velocity. The peaks are expected to become even sharper with increasing resolution [@afkhami2009]. While it is known that for the Navier-slip boundary condition pressure is divergent [@devauchelle2007], whether the reduced Young’s stress provides a cut-off mechanism for the pressure has not yet been analytically determined.

Because of the difference in the spreading curves in figure \[figure:diameterTanner\] an additional simulation was performed at a much larger resolution. The simulation domain for this simulation is $0.4\si{\milli\m} \times 0.2\si{\milli\m}$ with a resolution of $256 \times 128$, but the mesh is refined at the wall to better capture the curvature at the contact line. The grid cells at the wall are about $10 \si{\nano\m}$ cubed. The droplet has an initial radius of $0.15\si{\milli\m}$ and makes an angle to $\theta = \pi/4$ with the surface. The function $g$ which was used in the above simulations to keep the contact line tension localized was omitted in this simulation because it is not needed at larger resolutions.

![\[figure:contourTannerRes\] Interface of a droplet at different times at high resolution.](fig13){width="60.00000%"}

![\[figure:diameterTannerRes\] Radius of a droplet as a function of time at high resolution.](fig14){width="60.00000%"}

Figure \[figure:contourTannerRes\] shows the interface ($\alpha = 0.5$) of the droplet for various times. Figure \[figure:diameterTannerRes\] shows the corresponding radial position of the contact line as a function of time. As was observed in Figure \[figure:diameterTanner\], this smaller droplet also shows inertia dominated spreading.

![\[figure:angleTannerRes\] Angle of the inflection point (i.e. the apparent contact angle) as a function of the capillary number at high resolution.](fig15){width="60.00000%"}

Because of the larger resolution in these simulations the apparent contact angle can accurately be determined. Figure \[figure:angleTannerRes\] shows the apparent contact angle as a function of the capillary number. It can be appreciated that the simulation recovers the Voinov-Tanner-Cox [@voinov1976; @tanner1979; @cox1986] law: ($\theta \propto {\text{Ca}}^{1/3}$).

![\[figure:pressureTannerRes\] Pressure at the wall for various times at high resolution.](fig16){width="60.00000%"}

![\[figure:velocityTannerRes\] Velocity at the wall for various times at high resolution.](fig17){width="60.00000%"}

Figures \[figure:pressureTannerRes\] and \[figure:velocityTannerRes\] show the pressure peak and slip velocity at the wall at various times. Because of the high resolution at the wall both are very sharp peaks confined to the interface. The fluctuations of the pressure and velocity that was observed in figures \[figure:pressureTanner\] and \[figure:velocityTanner\] is no longer there, suggesting the simulations have fully converged.

Conclusions & Discussion
========================

An implementation of the Generalized Navier boundary condition for the Volume Of Fluid method is presented in this work. In analogy with the Brackbill surface tension model, a body force representation is developed for the contact line tension, while the Navier slip condition is applied on the wall. A validation of the code is presented for both a static case of a droplet maintaining its equilibrium shape, and a dynamic case of a spreading droplet.

It is shown how, on a completely smooth solid surface, in a system without gravity, the shape of the simulated droplets matches with their theoretically predicted shape for various Young’s angles. In addition, it is shown how the pressure peak and corresponding velocity peak at the interface converge with increasing resolution.

For the dynamic case it is found that the spreading of the droplet scales as $r
\propto t^{1/2}$. This suggests that spreading is limited by inertia. Also the Voinov-Tanner-Cox law is observed ($\theta \propto {\text{Ca}}^{1/3}$). This behavior is observed without using any fitting parameters.

An interesting future application of this model is flow over pattered surfaces [@dupuis2004; @wang2008]. Another topic of interest is the possibility to investigate both static and dynamic contact angle pinning. One can apply any pattern of equilibrium contact angles on a surface, and study how different patterns pin the contact line, either in a static or in a dynamic system.

Many thanks to Sidney Nagel and the Nagel group at the Department of Physics of the University of Chicago for stimulating discussions and feedback.

[74]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, , , , , , ().

, , , ****, ().

, pp. ().

, ** (, ).

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , pp. ().

, , , , , ****, ().

, , , ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , ****, ().

, , , , ****, ().

, , , ****, ().

, , , , ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, , , , , ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, , , ****, ().

, ****, ().

, , , , ****, ().

, , , ****, ().

, , , , , , ****, ().

, , , ****, ().

, Ph.D. thesis, ().

, ** (), <http://www.openfoam.org>.

, Ph.D. thesis, ().

, , , ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, , , ****, ().
