# Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized Navier boundary condition

Yash Kulkarni, Tomas Fullana, Stephane Zaleski

► To cite this version:

Yash Kulkarni, Tomas Fullana, Stephane Zaleski. Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized Navier boundary condition. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2023, 479 (2278), pp.20230141. 10.1098/rspa.2023.0141. hal-05417267

HAL Id: hal-05417267

https://hal.science/hal-05417267v1

Submitted on 15 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PROCEEDINGS A

royalsocietypublishing.org/journal/rspa

Research

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

**Cite this article:** Kulkarni Y, Fullana T, Zaleski S. 2023 Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized Navier boundary condition. *Proc. R. Soc. A 479:* 20230141.

https://doi.org/10.1098/rspa.2023.0141

Received: 22 February 2023

Accepted: 14 September 2023

**Subject Areas:**

applied mathematics, fluid mechanics

**Keywords:**

Stokes flow, stream function, three-phase
contact line, bi-harmonic equation

**Author for correspondence:**

Yash Kulkarni
e-mail: yash.kulkarni@dalembert.upmc.fr

Stream function solutions for
some contact line boundary
conditions: Navier slip, super
slip and the generalized Navier
boundary condition

Yash Kulkarni<fnref n="1" />, Tomas Fullana<fnref n="1" /> and Stephane
Zaleski<fnref n="1" />,<fnref n="2" />

<figure><img src="image_4.png" /></figure>

YK, 0000-0001-6228-2865

The stream function solution for the inner region Stokes flow, for a locally plane moving fluid interface near the triple point, is derived considering three different boundary conditions: the Navier slip boundary condition (NBC), the super-slip boundary condition and the generalized Navier boundary condition (GNBC). The NBC, incorporating a slip length parameter λ, is a well-known method for regularization in the context of the three-phase dynamic contact line problem. It is demonstrated that the velocity field solution under this boundary condition maintains a C⁰ continuity at the contact line, resulting in a logarithmic divergence of the pressure at the contact line. By contrast, the super-slip boundary condition establishes a proportional relationship between the wall velocity and the normal derivative of the shear stress, leading to a C<fnref n="1" /> velocity field. Furthermore, the GNBC, which introduces an uncompensated Young stress to drive the contact line, yields a C<fnref n="2" /> velocity field. The dominant terms are explicitly derived, and the analytical approach presented here can be extended to other bi-harmonic problems as well.

Electronic supplementary material is available
online at https://doi.org/10.6084/m9 figshare.
c.6858438.

© 2023 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

<footnote n="1">¹Sorbonne Université and CNRS, UMR 7190, Institut Jean le Rond
∂'Alembert, 75005 Paris, France</footnote>

<footnote n="2">²Institut Universitaire de France, UMR 7190, Institut Jean le Rond
∂'Alembert, 75005 Paris, France</footnote>
# 1. Introduction

A three-phase dynamic contact line emerges at the intersection of a fluid-fluid interface and a solid surface. It is widely recognized that the traditional no-slip boundary condition exhibits a force singularity at the contact line [1] and results in a discontinuous velocity field [2]. The proof of velocity field discontinuity was demonstrated by Dussan & Davis [2] solely based on the continuum hypothesis and mass conservation, without invoking momentum balance. To address this critical flaw in the no-slip model, one proposed solution is to allow for slip velocity at the fluid-solid interface. Various mechanisms for fluid slip exist, but the most popular one, suggested by Henry Navier (even before the discovery of the no-slip paradox), assumes that the slip velocity is proportional to the shear itself. This proportionality constant is referred to as the slip length, typically of the order of a few nanometres. In the vicinity of a three-phase contact line at this scale, gravity and inertial effects can be neglected, and the interface can be assumed to be locally planar. According to Hocking [3], the assumption of a locally planar interface is equivalent to the presence of surface tension between the two fluids, which is sufficiently large to compensate for the normal stress imbalance across the interface, resulting in negligible curvature. This enables the formulation of a Stokes flow problem, and by defining a stream function, the governing partial differential equation reduces to the bi-harmonic equation for the stream function, subject to appropriate boundary conditions. Hocking [3] solved such a flow problem for a three-phase moving contact line and demonstrated the absence of force singularity in the no-slip model. However, the general case considered by Hocking [3] does not yield an explicit stream function solution since, at that time, the primary focus for fluid mathematicians regarding contact line motion was resolving the force singularity. In parallel to Hocking's work [3], Huh & Mason [4] obtained an axisymmetric stream function solution for a capillary rise problem using the slip boundary condition, demonstrating the elimination of the force singularity for the Navier slip model (Model 1 in their paper [4]). Unlike Hocking [3], Huh & Mason [4] provided an explicit expression for the stream function but only for a specific case involving a 90° contact angle with a free surface. Recently, in a kinematic analysis of a sharp interface by Fricke et al. [5], it was shown that assuming a $C^1$ velocity field up to the contact line (as stated in definition 7 of their paper) implies that the velocity gradient in the direction normal to the solid boundary is equal to the rate of rotation of the normal vector at the contact line. This implies that the shear stress must vanish at the triple point. For the Navier slip model, this implies that the contact angle must change with time if the contact line moves, rendering steady-state motion of the contact line impossible. This motivates us to re-examine Hocking's analysis [3] and calculate the complete stream function to obtain the velocity field, investigating how the observed steady-state motion reported by Hocking is possible, despite the apparent contradiction raised by Fricke et al. [5]. In the subsequent section, we demonstrate that the Navier slip boundary condition (NBC) solution presented by Hocking [3] is non-differentiable at the contact line, thus violating the sufficiently smooth ($\ge C^1$) assumption made by Fricke et al. [5]. This non-differentiability issue also arises in the inner region solution of Huh & Mason [4] for their Model 1, which corresponds to the Navier slip model and exhibits a logarithmic singularity. As a non-differentiable velocity field implies an ill-posed nature of the boundary condition, we explore the possibility of introducing an additional level of smoothness by incorporating the slip velocity's dependence on the shear derivative. Hocking [6] and Devauchelle et al. [7] investigated such a boundary condition in the lubrication approximation. We demonstrate that this modified boundary condition provides an additional degree of smoothness to the Navier slip model, although the resulting velocity field remains non-$C^2$. We then present the expression for the generalized Navier boundary condition (GNBC) for a locally planar interface, and solving the Stokes flow with this boundary condition yields a $C^2$ velocity field. The outline of this paper is as follows: §2 defines the domain and governing equation and provides an overview of Hocking's [3] solution method. In §3, we solve for the stream function under the NBC. Subsequently, §4 presents the stream function solution for the super-slip boundary condition, followed by §5, where we obtain the stream function for the GNBC.
## 2. The Hocking solution for a moving fluid interface

For length scales below the slip length $\lambda$, the inertial and gravity effects are negligible and hence Hocking [3] solves the Stokes equation. The interface is assumed to be flat and a steady-state motion implies that the velocity of the solid as well as the contact angle are also constant. A polar coordinate system, as shown in figure 1a is used as per the details written in the caption. The flow is solved in the reference frame of the contact line (origin O).

Defining the stream functions $\Psi_1$ and $\Psi_2$ for each fluid, we get the bi-harmonic equation representing Stokes flow.<fnref n="1" />

$$ \Delta^2 \Psi_i = 0, \quad i=1,2. \qquad (2.1) $$

The boundary conditions include the NBC at the solid surface $\theta = \alpha_i$ ($i=1,2$). The slip length $\lambda$ is assumed to be the same for both liquid and gas.

$$ \left. \begin{array}{l} \frac{1}{r} \frac{\partial \Psi_1}{\partial \theta}(r, \alpha_1) + \frac{\lambda}{r^2} \frac{\partial^2 \Psi_1}{\partial \theta^2}(r, \alpha_1) = U \\ \\ \frac{1}{r} \frac{\partial \Psi_2}{\partial \theta}(r, \alpha_2) - \frac{\lambda}{r^2} \frac{\partial^2 \Psi_2}{\partial \theta^2}(r, \alpha_2) = -U \end{array} \right\} , \qquad (2.2) $$

and

and the streamline boundary condition at the interface and the solid boundary

$$ \Psi_i(r, 0) = \Psi_i(r, \alpha_i) = 0, \quad i=1,2. \qquad (2.3) $$

The continuity of the velocity and tangential stresses across the interface can be written as

$$ \partial_{\theta} \Psi_{1}(r, 0) = \partial_{\theta} \Psi_{2}(r, 0) \qquad (2.4) $$

and

$$ \partial_{\theta\theta} \Psi_1(r, 0) = \chi \partial_{\theta\theta} \Psi_2(r, 0), \qquad (2.5) $$

equations (2.4) and (2.5), respectively.

Solution of equation (2.1) can be written in the following form:

$$ \Psi_i = U r \Phi_i(\rho, \theta), \qquad (2.6) $$

where $\rho = \log(r/\lambda)$.

Under this transformation equation (2.1) becomes

$$ [(\partial_{\rho} + 1)^2 + \partial_{\theta\theta}] [(\partial_{\rho} - 1)^2 + \partial_{\theta\theta}] \Phi_i = 0, \qquad (2.7) $$

and the streamline boundary conditions (2.3) become

$$ \Phi_i(\rho, 0) = \Phi(\rho, \alpha_i) = 0. \qquad (2.8) $$

The slip boundary condition (2.2) becomes

$$ \left. \begin{array}{l} \partial_{\theta} \Phi_{1}(\rho, \alpha_{1}) + e^{-\rho} \partial_{\theta\theta} \Phi_{1}(\rho, \alpha_{1}) = 1 \\ \\ \partial_{\theta} \Phi_{2}(\rho, \alpha_{2}) - e^{-\rho} \partial_{\theta\theta} \Phi_{2}(\rho, \alpha_{2}) = -1 \end{array} \right\} . \qquad (2.9) $$

and

A two-sided Laplace transform or Bi-Laplace transform of a function $f(r)$ with respect to the variable $r$ is defined by $B(f(r)) = \tilde{f}(s) = \int_{-\infty}^{\infty} e^{-sr} f(r) dr$. Now Hocking [3] takes the two-sided Laplace transform w.r.t. variable $\rho$ and equation (2.7) subject to the streamline boundary condition (2.8) has the following solution (2.10) in the Laplace domain:

$$ \tilde{\Phi}_i(s, \theta) = A_i(s) \sin(s\theta) \sin(\theta - \alpha_i) + B_i(s) \sin(\theta) \sin(s(\theta - \alpha_i)). \qquad (2.10) $$

The coefficients $A_i(s)$ and $B_i(s)$ need to be determined. In the paper, from here on, Hocking [3] solves for the wall shear stress, which can be expressed in terms of the boundary values of $\Phi(\rho, \alpha)$ and its derivatives. These are later substituted to the transform of the slip boundary condition along with the no-slip matching condition and the full expression for the Wall-shear stress can

<footnote n="1">¹Details on the stream function formulation can be found in ch. 4 of [8].</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Illustrative figures: (a) polar coordinate system. The origin is the contact line and the interface is the plane $\theta = 0$. The sector $0 < \theta < \alpha_1$ is the fluid with viscosity $\mu_1$ and vice versa. $\alpha_2 = \alpha_1 - \pi$ (b) The illustration of various length scales for the GNBC.</figcaption></figure>

be obtained in the Laplace domain. Analytical inversion for the wall shear stress is possible for $\alpha_1 = \pi/2$ and for asymptotic values close to $\pi$. The conclusion is that the Wall shear stress is finite and the force singularity is removed. We call equation (2.10), the Hocking Solution.

### 3. The Navier slip model solution in the whole domain: complete form of the Hocking solution

The full solution of the flow field in the above solution requires the knowledge of the stream function in the full domain $\theta \in [\alpha_2, \alpha_1]$. This requires direct calculation of the coefficients $A_i(s)$ and $B_i(s)$ in equation (2.10). This was not done in [3] and we obtain them in this section. The continuity of velocity along the interface gives

$$s \sin \alpha_1 (A_1(s) + A_2(s)) = -B_1(s) \sin s\alpha_1 - B_2(s) \sin s(\pi - \alpha_1), \quad (3.1)$$

and the continuity of the shear stress across the interface gives

$$\cos \alpha_1 (A_1(s) - A_2(s)) = -B_1(s) \cos s\alpha_1 + \chi B_2(s) \cos s(\pi - \alpha_1). \quad (3.2)$$

Note that the solution (2.10) already satisfies the Stokes equation and the streamline boundary condition. We now must use the slip boundary condition to obtain the remaining two functional equations (one for each phase) for $A_i(s)$ and $B_i(s)$. To do this, we take the Bi-Laplace transform of the slip boundary condition (2.9) for each fluid.

The transform of (2.9) is written as

$$\left. \begin{array}{l} \partial_{\theta} \tilde{\Phi}_{1}(s, \alpha_{1}) = \mathcal{B}[1 - e^{-\rho} \partial_{\theta\theta} \tilde{\Phi}_{1}(\rho, \alpha_{1})] \\ \\ \partial_{\theta} \tilde{\Phi}_{1}(s, \alpha_{2}) = -\mathcal{B}[1 - e^{-\rho} \partial_{\theta\theta} \tilde{\Phi}_{1}(\rho, \alpha_{2})] \end{array} \right\} . \quad (3.3)$$

For the integration on RHS, one can write $1 \sim e^{-\epsilon|\rho|}$ for some positive $\epsilon \ll 1$. The integration then gives the following functional equation after substituting the boundary values of the derivatives of $\Phi$ w.r.t. $\theta$:

$$\left. \begin{array}{l} A_1(s) \sin s\alpha_1 + B_1(s) s \sin \alpha_1 \\ \\ + 2(s+1)[A_1(s+1) \cos(s+1)\alpha_1 + B_1(s+1) \cos \alpha_1] = \left( \frac{1}{s+\epsilon} - \frac{1}{s-\epsilon} \right) \\ \\ A_2(s) \sin s\alpha_2 + B_2(s) s \sin \alpha_2 \\ \\ - 2(s+1)[A_2(s+1) \cos(s+1)\alpha_2 + B_2(s+1) \cos \alpha_2] = - \left( \frac{1}{s+\epsilon} - \frac{1}{s-\epsilon} \right) \end{array} \right\} . \quad (3.4)$$

The solution of the inverse transform of the above four coupled equations (3.1), (3.2) and (3.4) would provide the needful expressions of the stream function for each fluid.
(a) Solution for the liquid–gas interface

The condition (3.2) for continuity of the shear stress along the interface decouples for fluid 1 if $\chi = 0$, which could be used as an approximation for the liquid–gas interface. In such a case, we can solve the decoupled equation obtained for fluid 1.

Put $\chi = 0$ in equation (3.2) to get the relation between the coefficients for fluid 1:

$$A_1(s) = C_0(s)B_1(s), \quad (3.5)$$

where

$$C_0(s) = \frac{-\cos(s\alpha)}{\cos\alpha}.$$

The slip boundary condition (3.4), now becomes a functional equation for $B_1(s)$

$$\begin{aligned} B_1(s+1) + B_1(s) \frac{(C_0(s) \sin s\alpha + s \sin \alpha)}{2(s+1)\left(C_0(s+1) \cos(s+1)\alpha + \cos \alpha\right)} \\ = \frac{1}{2(s+1)\left(C_0(s+1) \cos(s+1)\alpha + \cos \alpha\right)} \left(\frac{1}{s+\epsilon} - \frac{1}{s-\epsilon}\right). \end{aligned} \quad (3.6)$$

Knowing the function $B_1(s)$, the full solution for the streamfunction $\Phi_1(\rho, \theta)$ is obtained by doing the inverse transform of equation (2.10), which in this case becomes

$$\Phi_1(\rho, \theta) = \frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} e^{s\rho} B_1(s) \left( C_0(s) \sin(s\theta) \sin(\theta-\alpha) + \sin(\theta) \sin(s\theta-\alpha) \right) ds. \quad (3.7)$$

We now use the residue theorem to get the inverse transform. We examine the poles for a particular solution of the functional equation for $\Re(s) > 0$, which would correspond to the solution in the inner region ($r < \lambda$). There must not be a pole at $s = \epsilon$ as it is prohibited by the existence condition of the two-sided Laplace transform. Hence, there is a simple pole at $s = 1 + \epsilon$ and by induction, there are simple poles at $s = n + \epsilon$ for $n = 1, 2, 3, \dots$ From the functional equation, the inversion at these poles gives the particular solution

$$\Phi_{1_p}(r, \theta) = \left(\frac{r}{\lambda}\right) \mathcal{K}_1(\theta) + \left(\frac{r}{\lambda}\right)^2 \mathcal{K}_2(\theta) + \left(\frac{r}{\lambda}\right)^3 \mathcal{K}_3(\theta) + \dots, \quad (3.8)$$

where the coefficient functions $\mathcal{K}_i(\theta)$ for $i = 1, 2, 3, \dots$ are determined recursively using the functional equation (3.6). We write $\mathcal{K}_1(\theta)$ explicitly, due to its dominance and other coefficient functions can then be determined recursively. Note that the following notation is used in general:

$$\mathcal{T}_n(\theta, \alpha) = C_0(n) \sin(n\theta) \sin(\theta - \alpha) + \sin(\theta) \sin(n\theta - \alpha). \quad (3.9)$$

The coefficient functions are

$$\left. \begin{aligned} \mathcal{K}_1(\theta) &= \left(-\frac{1}{2\alpha \sin \alpha}\right) \left( \sin(\theta - \alpha)(\alpha \sin \theta \tan \alpha - \theta \cos \theta) + (\theta - \alpha) \sin \theta \cos(\theta - \alpha) \right) \\ \mathcal{K}_2(\theta) &= \mathcal{F}_2(\alpha)\mathcal{T}_2(\theta, \alpha) \quad \text{where } \mathcal{F}_2(\alpha) = -\frac{1}{8 \cos^2 \alpha - \cos^2 2\alpha} \left( \frac{\tan \alpha}{2} - \frac{1}{2 \tan \alpha} + \frac{1}{2\alpha} \right). \end{aligned} \right\} \quad (3.10)$$

For $n \ge 3$,

$$\mathcal{K}_n(\theta) = \mathcal{F}_n(\alpha)\mathcal{T}_n(\theta, \alpha). \quad (3.11)$$

The constant $\mathcal{F}_n(\alpha)$ is obtained recursively as follows (for $m \ge 2$):

$$\mathcal{F}_{m+1}(\alpha) = \mathcal{F}_m(\alpha) \frac{(C_0(m) \sin m\alpha + m \sin \alpha)}{2(m+1)\left(C_0(m+1) \cos(m+1)\alpha + \cos \alpha\right)}. \quad (3.12)$$

Note that the above solution exhibits a singularity at $\alpha = \pi/2$. Such spurious singularities often arise while solving Stokes equation [9] and this implies the presence of an inner region
eigensolution that is asymptotically dominant for angles greater than $\pi/2$. To deal with it in the present case, one needs to consider the double poles of the functional equation (3.6). We do not attempt that here as the qualitative description is already contained in our solution for $\theta < \pi/2$.

The solution for gas phase $\Phi_2(r, \theta)$ can now be obtained from the functional equations arising from the continuity of velocity along the interface and the slip boundary condition. From repeated substitution, the continuity of velocity along the interface condition becomes

$$A_2(s) = -(C_0 B_1(s) + C_1 B_1(s) + C_2 B_2(s)), \quad (3.13)$$

where

$$C_1(s) = \frac{\sin s\alpha_1}{s \sin \alpha_1} \quad \text{and} \quad C_2(s) = \frac{\sin s(\pi - \alpha_1)}{s \sin \alpha_1}.$$

The slip boundary condition takes the form

$$\left.
\begin{aligned}
& -B_1(s)(C_0(s) + C_1(s)) \sin s\alpha_2 + B_2(s)(-C_2(s) \sin s\alpha_2 + s \sin \alpha_2), \\
& -2(s+1)\left(-B_1(s+1)(C_0(s+1) + C_1(s+1)) \cos(s+1)\alpha_2\right), \\
& +B_2(s+1)(-C_2(s+1) \cos(s+1)\alpha_2 + \cos \alpha_2)
\end{aligned}
\right\} = -\left(\frac{1}{s+\epsilon} - \frac{1}{s-\epsilon}\right). \quad (3.14)$$

and

The function $B_1(s)$ is already known from the liquid solution. Using the functional equation for $B_2(s)$, we can invert the full flow field in the gas region

$$\Phi_2(\rho, \theta) = \frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} e^{s\rho} \left( -B_1(s)(C_0(s) + C_1(s)) \sin s\theta \sin(\theta - \alpha_2) + B_2(s)(-C_2(s) \sin s\theta \sin(\theta - \alpha_2) + \sin\theta \sin s(\theta - \alpha_2)) \right) ds. \quad (3.15)$$

A qualitative description of the stream function can be obtained without calculating the exact coefficients. Expanding for $K(1+\epsilon) = -(B_1C_0 + B_1C_1 + B_2C_2)/B_2$, the pole at $s=1+\epsilon$ it can be shown that the leading order is $\epsilon^3$. The solution thus obtained will be of the following form:

$$T_{\text{gas}}(\alpha, \theta) = -\frac{\theta^2}{2!} \alpha \tan \alpha \sin \theta \sin \phi + \frac{\theta^3}{3!} \sin \phi \cos \theta - \frac{\phi^3}{3!} \sin \theta \cos \phi, \quad (3.16)$$

where $\phi = (\theta - \alpha)$.

Now we are in a position to plot the flow fields in both the liquid and the gas regions. The plots are shown in figure 2 (liquid corresponds to $K_1(\alpha, \theta)$ and gas to $T_{\text{gas}}(\alpha, \theta)$) and are seen to depict the practical branching features of the flow solution.

One can now check the leading-order term for the liquid region from equations (3.8) and (3.10). There is a presence of $r\theta$ term, which makes the solution non-differentiable at the contact line $r=0$.

## 4. The super slip boundary condition

Further regularization of the contact line singularity can be done by introducing a second-order slip law, which we call the super slip boundary condition. In reference with figure 1a this boundary condition,<fnref n="2" /> for the liquid region, is written as

$$\frac{1}{r} \frac{\partial \Psi_1}{\partial \theta}(r, \alpha_1) + \frac{\lambda}{r^2} \frac{\partial^2 \Psi_1}{\partial \theta^2}(r, \alpha_1) + \frac{S^2}{r^3} \frac{\partial^3 \Psi_1}{\partial \theta^3}(r, \alpha_1) = U, \quad (4.1)$$

where $S$ is the super slip length.

From the molecular kinetics [7,10], we know that $S \sim O(\lambda)$. Setting $S = \lambda$ allows us to solve for the above boundary condition. Using the previous transformation of $\Psi(r, \theta) = Ur\Phi(\rho, \theta)$, the

<footnote n="2">²We continue to solve in the liquid-gas approximation of zero viscosity ratio ($\chi=0$).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Flow field corresponding to $K_1(\alpha, \theta)$, which represents the viscous liquid region. $\theta = 0$ represents the liquid-gas interface and $\theta = 45^\circ$ represents the solid-liquid wall. (b) Flow field visualization for $T_{gas}(\alpha, \theta)$, which corresponds to the gas region. $\theta = 0$ represents the liquid-gas interface and $\theta = 45^\circ$ represents the solid-air wall.</figcaption></figure>

above boundary condition becomes

$$ \partial_{\theta} \Phi_{1}(\rho, \alpha_{1}) + e^{-\rho} \partial_{\theta\theta} \Phi_{1}(\rho, \alpha_{1}) + e^{-2\rho} \partial_{\theta^{3}} \Phi_{1}(\rho, \alpha_{1}) = 1. \quad (4.2) $$

Note that in the Laplace domain, the form of the Stokes flow solution is still the same as (2.10), and the constants $A_1(s)$ and $B_1(s)$ will now be calculated as per the super slip boundary condition (4.2). Taking the two-sided Laplace transform of (4.2) and substituting the boundary values of $\Phi$ and its derivatives, we get a functional equation. With the super slip length dominating at the smallest scale, we can simplify the leading order solution by neglecting the first slip length and considering only the super slip length. The functional equation then becomes

$$ B_1(s + 2, \alpha_1) + \frac{\mathcal{N}}{D} B_1(s, \alpha_1) = \frac{1}{D} \left( \frac{1}{s + \epsilon} - \frac{1}{s - \epsilon} \right), \quad (4.3) $$

where

$$ \mathcal{N} = C_0(s) \sin s\theta + s \sin \theta $$

and

$$ \mathcal{D} = - \left( (3s^2 + 12s + 13) \frac{-\sin(2s + 4)\alpha_1}{2 \cos \alpha_1} + (s^3 + 6s^2 + 11s + 6) \sin \alpha_1 \right). $$

The coefficient function $B(s)$ is obtained from (4.3) and the full solution, which is the inverse transform of (2.10), becomes (3.7). Using the Cauchy residue theorem, it is possible to write the inner solution corresponding to the inverse (3.7), knowing the poles of (4.3). There must not be a pole at $s = \epsilon$ due to Hocking's existence condition. Now there are poles at $s = 2 + \epsilon, 4 + \epsilon, 6 + \epsilon, \dots$ giving the particular solution in the inner region. The leading order solution is the following:

$$ \Phi(r, \theta) = \frac{r^2}{S^2} \left( \frac{2 \cos \alpha_1}{2 \sin 2\alpha_1 (-13 \cos 2\alpha_1 + 3)} \right) \left( \sin \theta \sin 2(\theta - \alpha_1) - \frac{\cos 2\alpha_1}{\cos \alpha_1} \sin 2\theta \sin(\theta - \alpha_1) \right) + \frac{r^4}{S^4} K_2(s, \theta) + \dots \quad (4.4) $$

The term $r^2 \sin \theta \sin 2\theta$ increases the degree of smoothness of the solution as compared with the Navier slip model. This makes the solution differentiable at the contact line, however, it is still not $C^2$.

## 5. The generalized Navier boundary condition

The GNBC is a contact line boundary condition that incorporates the uncompensated Young's stress, in addition to the typical viscous stresses. This was suggested in [11] considering the
irreversible thermodynamics and for a brief derivation from entropy principles in the sharp interface limit, one can refer to [12]. The popular way of writing the GNBC in the reference frame of solid, also mentioned in [12] is

$$u - \lambda \frac{\partial u}{\partial y} = \frac{\sigma \lambda}{\mu} [\cos \theta - \cos \theta_{eq}] f_{\epsilon}(r). \quad (5.1)$$

Here, $f_{\epsilon}(x)$ is an approximation of the delta function concentrated near the contact line. As per the kinematic analysis of [5], for a well-defined $C^2$ velocity field, $\partial u/\partial y \propto d\theta/dt$, so the steady-state contact line velocity becomes

$$u_{cl} = \frac{\sigma \lambda}{\mu} [\cos \theta - \cos \theta_{eq}] \frac{1}{\epsilon'}. \quad (5.2)$$

Here, $\epsilon'$ is the length scale over which we smear out the delta function. This also means that the definition of the contact angle and the contact line itself is extended on this length scale. To formulate the problem in a flat interface approximation, we substitute the result (5.2) into (5.1) and we now approximate the local function as $f_{\epsilon}(r) = 1 - \tanh^2 r/\epsilon$. This is illustrated in figure 1b. In the reference frame of the contact line, this gives us the following form of the GNBC at the leading order valid for $r > \epsilon \gg \epsilon'$:

$$u + \epsilon' \frac{\partial u}{\partial y} = U_s \frac{r^2}{\epsilon^2}. \quad (5.3)$$

Using the polar coordinates and after defining the streamfunction, this is written as

$$\frac{1}{r} \frac{\partial \Psi}{\partial \theta} (r, \alpha) + \frac{\epsilon' \partial^2 \Psi}{r^2 \partial \theta^2} (r, \alpha) = U_s \frac{r^2}{\epsilon^2}. \quad (5.4)$$

Substitute $\rho = \log(r/\epsilon')$ and the form of solution $\Psi = Ur\Phi(\rho, \theta)$, we get the following equation at the boundary:

$$\Phi_{\theta}(\rho, \alpha) + e^{-\rho} \Phi_{\theta\theta}(\rho, \alpha) = e^{2\rho} \delta^2, \quad (5.5)$$

where $\delta = \epsilon'/\epsilon$.

Taking the Laplace transform of the above w.r.t. variable $\rho$, we get the functional equation in the Laplace domain<fnref n="3" />

$$\tilde{\Phi}_{\theta}(s, \alpha) + \tilde{\Phi}_{\theta\theta}(s + 1, \alpha) = \delta^2 \left( \frac{1}{s - 2} \right). \quad (5.6)$$

Using Cauchy's residue theorem, the inverse transform is written as equation (3.7). In the limit $\chi \to 0$, we have for fluid 1, the functional equation

$$B_1(s + 1) + B_1(s) \frac{(C_0(s) \sin s\alpha + s \sin \alpha)}{2(s + 1)(C_0(s + 1) \cos(s + 1)\alpha + \cos \alpha)} \\ = \frac{\delta^2}{(s - 2)} \frac{1}{2(s + 1)(C_0(s + 1) \cos(s + 1)\alpha + \cos \alpha)}. \quad (5.7)$$

The above functional equation gives a particular solution in the inner region corresponding to the poles at $s = 3, 4, 5, \dots$. The solution thus becomes

$$\Phi_p(r, \theta) = \left(\frac{r}{\epsilon'}\right)^3 \delta^2 \mathcal{G}_3(\theta) + \left(\frac{r}{\epsilon'}\right)^4 \delta^2 \mathcal{G}_4(\theta) + \dots, \quad (5.8)$$

where $\mathcal{G}_3(\theta)$ is given as

$$\mathcal{G}_3(\theta) = -\frac{1}{6} \frac{\cos \alpha}{\cos^2 3\alpha (1 + \cos \alpha)} \left( \frac{-\cos 3\alpha}{\cos \alpha} \sin 3\theta \sin(\theta - \alpha) + \sin \theta \sin 3(\theta - \alpha) \right). \quad (5.9)$$

Note that the above solution is $C^2$.

<footnote n="3">³RHS converges iff $\Re(s) > 2$.</footnote>
## 6. A note on asymptotic matching to the far-field solution

For a full-scale description of the flow field, we need to match the inner-scale stream function solutions we obtained to the outer solution. A reasonable choice for the outer solution is an intermediate length scale $r \gg \lambda$ that is large enough to assume the continuum hypothesis and approximate no-slip boundary condition and small enough to have viscous-capillary balance and neglect the inertial effects. This can be realized at small Ca. In that case, the Stokes flow assumption remains valid. Hence, one can observe this particular region at low capillary numbers, for example in receding contact lines, where the interface bending, as discussed by Cox [13] yields the curvature $\kappa \sim Ca/r$.

The outer solution for $r \gg \lambda$, that is a Stokes flow subject approximately to the no-slip boundary condition could be obtained by substituting $\lambda/r^2 = 0$ in the boundary condition 2.2. For the free surface case, the stream function solution can be written as

$$ \Psi_o = U r \Phi_o(\theta), \qquad (6.1) $$

where subscript 'o' indicates the outer region, and

$$ \Phi_o(\theta) = \frac{2}{2\alpha - \sin 2\alpha} \theta \sin(\theta - \alpha) + \frac{-2 \cos \alpha}{2\alpha - \sin 2\alpha} (\theta - \alpha) \sin \theta. \qquad (6.2) $$

The equation above corresponds to a corner flow solution and is well-known in the literature [1,3,9]. To take into account the effect of the bending, one can use this corner flow stream function and compute the latter from the Laplace pressure jump. The bending of the interface thus obtained in terms of the local slope $\theta$ is given by the famous Cox–Voinov log law [13,14]

$$ g(\theta) = g(\theta_m) - Ca \ln \left( \frac{r}{r_m} \right), \qquad (6.3) $$

where

$$ g(\theta) = \int_{0}^{\theta} \frac{t - \cos t \sin t}{2 \sin t} dt, \qquad (6.4) $$

$\theta_m$ is the microscopic contact angle and $r_m$ is the microscopic length scale. Note that $\theta_m$ is not necessarily the equilibrium contact angle obtained from the Young's law. For example, for the GNBC, this angle is the $\theta$ obtained in equation (5.2). The matching would therefore be between the inner region solutions we obtained with the Cox–Voinov log law (6.3). For the GNBC, this was done in [15] using a generalized lubrication approach. The microscopic length scale, in this case, becomes $r_m = c\lambda$, where coefficient $c$ was found by matching. Note that equation (6.3) represents a family of curves on the $g(\theta)$ versus $\ln(r)$ plane and the problem of matching is the problem to calculating the scale $r_m$. This is illustrated in figure 3.

This problem of matching is a non-trivial problem and one needs to use a numerical solution to get the parameter $c$. In our new solutions, this parameter $c$ is a function of super slip length $S$ for the super-slip boundary condition and $\epsilon$ for the GNBC. This is left as a future scope of the paper. To summarize, there is nothing in the new expressions for a boundary condition near the contact line that prevents the matching of the inner solution found in the present work with the Cox outer solution. This matching would provide the functional relation between the Cox–Voinov microscopic length and the microscopic lengths present in the models discussed here.

However, one must be careful with the steady state advancing contact lines for dense fluids such as water. In those cases, the contact line capillary number might be of order one and inertial effects would play a role therefore altering the stream function solution even at the scales of a few hundred nanometres as shown in [16]. In that case, the Cox region is absent and the inner region of the Stokes flow solutions with slip models we obtained should be matched to the first-order inertially corrected Stokes flow with the no-slip boundary condition.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Example of a matching procedure. $\epsilon$ is the microscopic length scale corresponding to the slip length $\lambda$ in the NBC and to $\epsilon$ in the GNBC.</figcaption></figure>

# 7. Conclusion

We have obtained the leading order Stokes flow solutions for the steady-state contact line motion subject to the NBC, the super slip boundary condition and the GNBC. The velocity flow field thus obtained for the NBC is $C^0$, for the super slip boundary condition is $C^1$ and for the GNBC is $C^2$. As shown in appendix A, our stream function solution asymptotically makes the pressure $p \sim \log r$ for the NBC, $p \sim r$ for the super slip and $p \sim r^2$ for the GNBC, where $r$ is the distance from the contact line. This means that the Navier slip model is singular at the contact line because the pressure diverges, therefore making the model thermodynamically ill-posed. The logarithmic singularity is however integrable, and hence weak. This weakly singular behaviour is regularized by the super slip and the GNBC models.

We would like to emphasize that our goal is not to design the boundary condition that would regularize the flow, instead, it is to understand the flow subject to already hypothesized boundary conditions derived from physics. The super slip boundary condition can be obtained using fundamental principles, similar to how we derive the Navier slip condition from the Molecular Kinetic Theory, but by considering an additional order. Our analysis of the super slip boundary condition reveals a substantial modification in the pressure behaviour at the contact line, effectively smoothing it out. However, it is important to note that the velocity field remains non-$C^2$. The GNBC is derived from entropy principles in continuum sharp interface limit [12], and has a clear meaning in the phase-field model [11]. This boundary condition not only regularizes the pressure but also makes the velocity flow field $C^2$. Hence, the increased smoothness that we obtain is an outcome of our analysis and in no way an a priori statement.

Another advantage these improved boundary conditions may provide benefits is in the numerical modelling of the contact line in steady-state wetting where obtaining convergence is quite a challenging task [17,18]. A preliminary attempt to use these improved models numerically can be seen in the appendix of the doctoral thesis [19], where a much faster convergence is seen.

We thus show direct calculations that increase our understanding of the smoothness of these boundary conditions. The analytical procedure presented here might also be used for other Bi-Laplacian problems. To justify the correctness of our analytical procedure, we present a corner flow between the slip and no-slip wall in appendix B. Another insight is that the boundary conditions we have, result in a *slip* like behaviour for the fluid at the contact line, rather than a *rolling* like behaviour, which is observed for the no-slip boundary condition. This means that, unlike the no-slip model, the contact line is a material line because the fluid particle approach
with the velocity $u \sim r^m$ where $m \ge 1$ and hence would take infinite time to reach the contact line. This is a well-known behaviour and was commented upon in [20]. In physical reality, the slip lengths are quite small and of the order of nanometres, so the slip models would still preserve the rolling type of kinematics at the experimental resolutions, which are often of the order of a few microns. As we zoom in more at the contact line, we would eventually enter a length scale smaller than the mean free path where the continuum mechanics will no longer hold. What we present are mathematical models that are useful for mathematicians studying the contact line in the continuum sharp interface limit.

**Data accessibility.** The data are provided in the electronic supplementary material [21].

**Declaration of AI use.** Yes, we have used AI-assisted technologies in creating this article.

**Authors' contributions.** Y.K.: conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing—original draft, writing—review and editing; T.F.: conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing—review and editing; S.Z.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 883849).

**Acknowledgements.** We would like to thank Mathis Fricke for the discussions. Yash Kulkarni would also like to thank the TU Darmstadt staff for their hospitality during his visit to Germany where several ideas on the contact line flowed freely. The slip kinematics and rolling kinematics discussed in the Conclusion were pointed out by an anonymous reviewer.

## Appendix A. Obtaining velocity fields and pressure fields

From the full stream function solutions, one can obtain the velocity and pressure fields.

The velocity field in the polar coordinates is obtained as

$$u_r = \frac{1}{r} \frac{\partial \Psi(r, \theta)}{\partial \theta} \quad \text{and} \quad u_\theta = -\frac{\partial \Psi(r, \theta)}{\partial r} \qquad (A 1)$$

Here, $u_r$ and $u_\theta$ represent the radial and polar components of the velocity field. For the Navier slip boundary condition, one can see from the solution (3.8) that $u \sim r/\lambda$ at the leading order. This is multiplied with the appropriate function of $\theta$ and plotted in figure 2. For the super slip boundary condition, one can see that $u \sim r^2$ and hence increasing the smoothness. For the GNBC one can see that $u \sim r^3$ increases the smoothness further and makes the solution $C^2$.

Once we have the velocity field, the pressure field is obtained as

$$\frac{\partial p}{\partial r} = \Delta u_r - \frac{u_r}{r^2} - \frac{2}{r^2} \frac{\partial u_\theta}{\partial \theta} \quad \text{and} \quad \frac{1}{r} \frac{\partial p}{\partial \theta} = \Delta u_\theta - \frac{u_\theta}{r^2} + \frac{2}{r^2} \frac{\partial u_r}{\partial u_\theta}, \qquad (A 2)$$

where $\Delta$ is the Laplacian operator. With the leading order stream function solutions obtained by us, one can easily conclude the asymptotic behaviour of the pressure scalar, that is $p \sim \log(r)$ for the Navier slip, $p \sim r$ for the super slip and $p \sim r^2$ for the GNBC and the conclusions follow.

## Appendix B. Stokes corner flow between perpendicular slip and no-slip walls

As proof that the analytical procedure works, we provide an example of a Stokes flow in a perpendicular corner formed by a slip wall and a no-slip wall. We compare the analytical result at the leading order with the full Navier-Stokes simulations done using the Basilisk flow solver.

The flow domain is shown in figure 4a and a formulation can be done in polar coordinates, such that $\theta = 0$ is the no-slip wall and $\theta = \alpha = \pi/2$ is the slip wall where streamline boundary conditions are applied. One needs to solve the bi-Laplacian of the stream function equation in the
<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) The flow domain schematic showing perpendicular slip wall and a no-slip wall. (b) Comparison for the velocity along the slip wall as obtained numerically and analytically. $\lambda$ is the slip length, the blue points represent the simulations and the orange line represents the leading-order term, asymptotically dominant for $y \ll \lambda$.</figcaption></figure>

polar coordinates written as,

$$ \Delta^2 \Psi(r, \theta) = 0. \qquad (B1) $$

Proceeding on the same lines as in §2, an ansatz is introduced $\Psi(r, \theta) = U_s r \Phi(\rho, \theta)$ where $\rho = \log(r/\lambda)$. We again take two-sided Laplace transform with respect to variable $\rho$ and the solution in the Laplace domain is written as

$$ \tilde{\Phi}(s, \theta) = A(s) \sin(s\theta) \sin(\theta - \alpha) + B(s) \sin(\theta) \sin(s(\theta - \alpha)). \qquad (B2) $$

No slip condition at $\theta = 0$ gives $\tilde{\Phi}_{\theta}(r, 0) = 0^4$. This gives a relation between the coefficient function, which is,

$$ A(s) = \mathcal{M}_0(s)B(s), \qquad (B3) $$

where

$$ \mathcal{M}_0(s) = - \frac{\sin s\alpha}{s \sin \alpha}. \qquad (B4) $$

We now must satisfy the slip boundary condition at $\theta = \alpha$, where $\alpha = \pi/2$. In the $\rho$ domain, this gives,

$$ \Phi(\rho, \alpha) + e^{-\rho} \Phi_{\theta\theta}(\rho, \alpha) = 1. \qquad (B5) $$

Again, we write $1 \approx e^{-\epsilon|\rho|}$ and this allows us to take the two-sided Laplace transform and after substituting appropriate derivatives of $\Phi$ and using coefficient relation (B3), results in the following functional equation,

$$ B(s + 1) - \left( \frac{\sin^2(s\pi/2) - s^2}{s \sin s\pi} \right) B(s) = \frac{1}{\sin s\pi} \left( \frac{1}{s + \epsilon} - \frac{1}{s - \epsilon} \right). \qquad (B6) $$

The solution of $\Phi(\rho, \theta)$ is obtained by doing an inverse transform of $\tilde{\Phi}(s, \theta)$, which is written as

$$ \Phi_1(\rho, \theta) = \frac{1}{2\pi i} \int_{\sigma - i\infty}^{\sigma + i\infty} e^{s\rho} B(s) \left( \mathcal{M}_0(s) \sin(s\theta) \sin(\theta - \alpha) + \sin(\theta) \sin(s(\theta - \alpha)) \right) ds. \qquad (B7) $$

Note that we are now only attempting a special case of $\alpha = \pi/2$. A solution like (3.8) can be constructed by doing term-wise inversion at poles at $s = 1 + \epsilon, 2 + \epsilon, 3 + \epsilon, \dots$. For the special case

<footnote>$^4Y_x$ represent $\partial Y/\partial x$.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 5. The Stokes flow field obtained by (a) Basilisk simulation and (b) leading order term of the solution (B 8).</figcaption></figure>

of $\alpha = \pi/2$, the leading order term, corresponding to the pole at $s = 1 + \epsilon$ is given as

$$ \Phi(r, \theta) = \frac{r}{\lambda} \frac{1}{\pi} \left( \frac{\sin 2\theta}{2} + \frac{\pi}{2} \sin^2 \theta - \theta \right) + \mathcal{O}\left(\frac{r}{\lambda}\right)^2. \quad (B8) $$

The advantage of having perpendicular walls is that we can perform numerical simulation of a single-phase Stokes flow using the open-sourced CFD solver Basilisk (basilisk.fr) and compare the result against our calculations. The flow fields obtained by numerics and that given by the first term of the solution (B8), are shown in figure 5, where a qualitative agreement can be seen. A quantitative comparison is possible for the velocity along any given angular coordinate. We show the variation of velocity along the slip wall ($\theta = \pi/2$) in figure 4b where one may see excellent agreement in the region $r \ll \lambda$. Not only the linear behaviour is reproduced, but the pre-factor of $2/\pi$ as predicted by solution (B8) is also coherent with the simulation result.

## References

1. Huh C, Scriven L. 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. *J. Colloid Interface Sci.* **35**, 85–101. (doi:10.1016/0021-9797(71)90188-3)

2. Dussan VEB, Davis SH. 1974 On the motion of a fluid-fluid interface along a solid surface. *J. Fluid Mech.* **65**, 71–95. (doi:10.1017/S0022112074001261)

3. Hocking LM. 1977 A moving fluid interface. Part 2. The removal of the force singularity by a slip flow. *J. Fluid Mech.* **79**, 209–229. (doi:10.1017/S0022112077000123)

4. Huh C, Mason SG. 1977 The steady movement of a liquid meniscus in a capillary tube. *J. Fluid Mech.* **81**, 401–419. (doi:10.1017/S0022112077002134)

5. Fricke M, Köhne M, Bothe D. 2019 A kinematic evolution equation for the dynamic contact angle and some consequences. *Physica D* **394**, 26–43. (doi:10.1016/j.physd.2019.01.008)

6. Hocking LM. 2001 Meniscus draw-up and draining. *Eur. J. Appl. Math.* **12**, 195–208. (doi:10.1017/S0956792501004247)

7. Devauchelle O, Josserand C, Zaleski S. 2007 Forced dewetting on porous media. *J. Fluid Mech.* **574**, 343–364. (doi:10.1017/S0022112006004125)

8. Batchelor GK. 2000 Flow of a uniform incompressible viscous fluid. In *An introduction to fluid dynamics* (ed. GK Batchelor), Cambridge Mathematical Library, pp. 174–263. Cambridge, UK: Cambridge University Press.

9. Moffatt HK, Duffy BR. 1980 Local similarity solutions and their limitations. *J. Fluid Mech.* **96**, 299–313. (doi:10.1017/S0022112080002133)

10. Hadjiconstantinou NG, Al-Mohssen HA. 2005 A linearized kinetic formulation including a second-order slip model for an impulsive start problem at arbitrary Knudsen numbers. *J. Fluid Mech.* **533**, 47–56. (doi:10.1017/S0022112005003897)
11. Qian T, Wang XP, Sheng P. 2006 A variational approach to moving contact line hydrodynamics. *J. Fluid Mech.* **564**, 333–360. (doi:10.1017/S0022112006001935)

12. Fricke M. 2021 *Mathematical modeling and Volume-of-Fluid based simulation of dynamic wetting*. PhD thesis, Technische Universität Darmstadt, Germany.

13. Cox RG. 1986 The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. *J. Fluid Mech.* **168**, 169–194. (doi:10.1017/S0022112086000332)

14. Voinov OV. 1976 Hydrodynamics of wetting. *Fluid Dyn.* **11**, 714–721. (doi:10.1007/BF01012963)

15. Chan TS, Kamal C, Snoeijer JH, Sprittles JE, Eggers J. 2020 Cox–Voinov theory with slip. *J. Fluid Mech.* **900**, A8. (doi:10.1017/jfm.2020.499)

16. Varma A, Roy A, Puthenveettil BA. 2021 Inertial effects on the flow near a moving contact line. *J. Fluid Mech.* **924**, A36. (doi:10.1017/jfm.2021.582)

17. Fullana T, Zaleski S, Popinet S. 2020 Dynamic wetting failure in curtain coating by the volume-of-fluid method volume-of-fluid simulations on quadtree meshes. *Eur. Phys. J. - Spec. Top.* **229**, 1923–1934. (doi:10.1140/epjst/e2020-000004-0)

18. Afkhami S, Zaleski S, Bussmann M. 2009 A mesh-dependent model for applying dynamic contact angles to VOF simulations. *J. Comput. Phys.* **228**, 5370–5389. (doi:10.1016/j.jcp.2009.04.027)

19. Fullana T. 2022 *Simulation and optimization of complex phenomena in multiphase flows*. PhD thesis, Sorbonne University, France.

20. Dussan VEB. 1976 The moving contact line: the slip boundary condition. *J. Fluid Mech.* **77**, 665–684. (doi:10.1017/S0022112076002838)

21. Kulkarni Y, Fullana T, Zaleski S. 2023 Stream function solutions for some contact line boundary conditions: Navier slip, super slip and the generalized Navier boundary condition. Figshare. (doi:10.6084/m9 figshare.c.6858438)