2.1. Governing equations

We employ the sharp-interface continuum modelling approach. We start from the incompressible, two-phase Navier–Stokes equations with surface tension for Newtonian fluids under isothermal conditions (see, e.g. Slattery Reference Slattery1999; Prüss & Simonett Reference Prüss and Simonett2016). Inside the fluid phases, the governing equations are

(2.1) \begin{align} \partial _t (\rho \boldsymbol{v}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{v} \otimes \boldsymbol{v}) + \boldsymbol{\nabla }\!p &= \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} + \rho \boldsymbol{g}, \end{align}

(2.2) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v} &= 0 \end{align}

with the viscous stress tensor (we use the symbol $\boldsymbol{D}={1}/{2} (\boldsymbol{\nabla }\boldsymbol{v} + \boldsymbol{\nabla }\boldsymbol{v}^{\mathsf{T}})$ for the rate-of-deformation tensor)

(2.3) \begin{align} \boldsymbol{S} = 2 \eta \boldsymbol{D} = \eta \big(\boldsymbol{\nabla }\boldsymbol{v} + \boldsymbol{\nabla }\boldsymbol{v}^{\mathsf{T}}\big). \end{align}

These bulk equations are accompanied by jump conditions at the interface $\varSigma (t)$ . The interface is modelled as a hypersurface (i.e. it has zero thickness) and separates the domain $\varOmega$ into two bulk phases $\varOmega _{1,2}(t)$ occupied by the two fluid phases (see figure 1). Assuming that no phase change occurs in the system, the normal component of the adjacent fluid velocities $\boldsymbol{v}_{1,2}$ at the interface are coinciding and equal to the speed of normal displacement $V_\varSigma$ of the interface, resulting in the kinematic boundary condition

(2.4) \begin{align} V_\varSigma = \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{n}_\varSigma \quad \text{on} \quad \varSigma (t), \end{align}

where $\boldsymbol{n}_\varSigma$ is the interface unit normal field. Additionally, no-slip between the fluid phases is usually assumed. Assuming further that the surface tension $\sigma$ is constant, the jump conditions for mass and momentum read as

(2.5) \begin{align} \left [\![ \boldsymbol{v} \right ]\!] = 0, \quad [\![ p \boldsymbol{\mathbb{I}} - \boldsymbol{S} ]\!] \boldsymbol{n}_\varSigma = \sigma \kappa \boldsymbol{n}_\varSigma \quad \text{on} \quad \varSigma (t). \end{align}

Here, $\kappa := - \boldsymbol{\nabla}_{\!\varSigma} \boldsymbol{\cdot }\boldsymbol{n}_\varSigma$ is twice the mean curvature of the interface and

(2.6) \begin{align} \left [\![ \psi \right ]\!](t,\boldsymbol{x}) := \lim _{h \rightarrow 0^+} (\psi (t,\boldsymbol{x}+h\boldsymbol{n}_\varSigma ) - \psi (t,\boldsymbol{x}-h\boldsymbol{n}_\varSigma )) \end{align}

is the jump of a quantity $\psi$ across the interface. We assume that the solid boundary $\partial \varOmega$ is not able to store mass and we assume it to be impermeable. We consider an inertial frame of reference, where the wall is moving parallel to itself with a velocity $U_w \geqslant 0$ upwards (see figure 1). The impermeability condition in this frame of reference reads as

(2.7) \begin{align} \boldsymbol{v}_\bot = 0 \quad \text{on} \quad \partial \varOmega , \end{align}

where $\boldsymbol{v}_\bot = (\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{n}_{\partial \varOmega }) \, \boldsymbol{n}_{\partial \varOmega }$ denotes the normal part of the velocity with respect to $\partial \varOmega$ .

To obtain a closed model, the system of (2.1)–(2.7) must be complemented by (one or more) additional boundary conditions describing

1. (i) the wettability of the solid (i.e. the static and dynamic contact angle) and

2. (ii) the mobility of the contact line (i.e. the tangential velocity $\boldsymbol{v}_\parallel$ at the solid boundary).

These boundary conditions are closure relations for the continuum mechanical description and must be thermodynamically consistent, i.e. they must obey in particular the second law of thermodynamics. To arrive at a consistent closure, we consider the available energy consisting of the kinetic energy of the bulk phases and the surface energies of the liquid–gas interface as well as the wetted area $W(t) \subset \partial \varOmega$ , i.e.

(2.8) \begin{align} E(t) := \int _{\varOmega \setminus \varSigma (t)} \frac {\rho |\boldsymbol{v}|^2}{2} \, \text{d}V + \int _{\varSigma (t)} \sigma \, \text{d}A + \int _{W(t)} \sigma _w \, \text{d}A. \end{align}

Here, $\sigma = \sigma _{\textit{lg}} \gt 0$ denotes the surface tension of the liquid–gas interface and

(2.9) \begin{align} \sigma _{\textit{w}} = \sigma _{\textit{ls}} - \sigma _{\textit{sg}} \end{align}

is the specific energy density for wetting the solid surface. Note that $\sigma _{\textit{w}}$ is negative for a hydrophilic surface, as we see from Young’s equation

(2.10) \begin{align} \sigma \cos \theta _e + \sigma _{\textit{w}} = 0, \end{align}

which defines the ‘static’ or ‘equilibrium’ contact angle $\theta _e$ . It is a purely mathematical exercise (see Ren Reference Ren and W.2007, Reference Ren and W.2011b ; Fricke Reference Fricke2021 for details) to compute the rate of change $\dot {E}$ for a sufficiently regular solution of (2.1)–(2.7) (in the absence of external forces, i.e. for $\boldsymbol{g}=0$ ). The result reads as

(2.11) \begin{align} \frac {\text{d} E}{\text{d}t} = - \!\!\int _{\varOmega \setminus \varSigma (t)} \!\boldsymbol{S}:\boldsymbol{D} \, \text{d}V + \!\!\int _{\partial \varOmega } \!(\boldsymbol{v}_\parallel -\boldsymbol{U}_w) \boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel \, \text{d}A + \!\!\int _{\varGamma (t)} \!\sigma (\cos \theta _d - \cos \theta _e) \, U_{\textit{cl}} \, \text{d}l. \end{align}

In this formulation with a continuous velocity field, the scalar contact line speed (measured relative to the solid) is given as

(2.12) \begin{align} U_{\textit{cl}} = \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{n}_\varGamma -U_w. \end{align}

Closure relations are required to satisfy the second law of thermodynamics $\dot {E} \leqslant 0$ . Note that we assume an isothermal system here, so we may directly consider the change in available energy. The first contribution in (2.11) has a negative sign as we consider incompressible Newtonian fluids, i.e. $\boldsymbol{S}=2\eta \boldsymbol{D}$ . A linear closure for the second integral in (2.11) yields the well-known Navier slip condition, i.e.

(2.13) \begin{align} -\beta (\boldsymbol{v}_\parallel -\boldsymbol{U}_w) = (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel \quad \text{with a friction coefficient} \quad \beta \geqslant 0. \end{align}

Using the slip length parameter $\lambda =\eta /\beta$ , one may reformulate (2.13) as

(2.14) \begin{align} \boldsymbol{v}_\parallel + 2 \lambda (\boldsymbol{D}\boldsymbol{n}_{\partial \varOmega })_\parallel = \boldsymbol{U}_w. \end{align}

The third integral in (2.11) suggests that the dynamic contact angle $\theta _d$ , which is mathematically defined as the angle of intersection of the free surface $\varSigma$ with the solid boundary $\partial \varOmega$ , i.e.

(2.15) \begin{align} \cos \theta _d := - \boldsymbol{n}_\varSigma \boldsymbol{\cdot }\boldsymbol{n}_{\partial \varOmega } \quad \text{at} \quad \varGamma (t), \end{align}

should be linked to the contact line speed $U_{\textit{cl}}$ . Note that the definition of $\theta_{d}$ assumes that the interface has a well-defined normal field up to the solid boundary, which is the case even if the curvature has a logarithmic, hence integrable, singularity. A linear closure leads to the well-known condition

(2.16) \begin{align} -\zeta U_{\textit{cl}} = \sigma (\cos \theta _d - \cos \theta _e) \quad \text{with} \; \zeta \geqslant 0, \end{align}

where $\zeta$ can be seen as a contact line friction coefficient. Note that also more general contact angle boundary conditions are possible if a nonlinear closure relation is employed. To summarise, the ‘standard model’ based on the Navier slip condition is given by (2.1)–(2.7), (2.14) together with (2.16) or a nonlinear variant of the form

(2.17) \begin{equation} \theta _d = f(U_{\textit{cl}}) \quad \text{on} \quad \varGamma (t). \end{equation}

The mathematical model (2.1)–(2.7), (2.14), (2.16) is one of the most commonly applied models for dynamic wetting in the literature. However, there are many more modelling approaches which aim at a regularisation of the singularity and a prediction of the dynamics of wetting. For a survey of the field, we refer to de Gennes, Brochard-Wyart & Quéré (Reference de Gennes, Brochard-Wyart and Quéré2004), Blake (Reference Qian, Wang and Sheng2006), Shikhmurzaev (Reference Shikhmurzaev2008), Bonn et al. (Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), Snoeijer & Andreotti (Reference Snoeijer and Andreotti2013) and Marengo & De Coninck (Reference Marengo and De Coninck2022). To ensure thermodynamic consistency, we require that

(2.18) \begin{align} \eta \geqslant 0, \ \ \sigma \geqslant 0, \ \ \lambda \geqslant 0, \ \ U_{\textit{cl}} (f(U_{\textit{cl}})-\theta _e) \geqslant 0. \end{align}

As shown by Fricke & Bothe (Reference Fricke and Bothe2020), for regular solutions, there is a fundamental inconsistency of boundary conditions in the standard model. In fact, the evolution of the contact angle is determined by the contact angle boundary condition (say (2.16)) as well as by the flow in the vicinity of the contact line according to (1.4). As a consequence, a regular solution of the system does not exist, but a weak singularity is present at the contact line as shown already by Huh & Mason (Reference Huh and Mason1977).

2.2. Formal derivation of the GNBC

We now recall the derivation of the GNBC in the context of the sharp-interface framework starting from (2.11).

It is important to note that the GNBC was originally formulated in a diffuse interface framework (see Qian et al. Reference Qian, Wang and Sheng2003, Reference Qian, Wang and Sheng2006a ). However, the GNBC can be formally understood in the sharp interface model as a combined closure for the terms in the entropy production (2.11) that arise from the contact line motion and from slip at the solid–liquid boundary. The combined closure formally leads to a single boundary condition instead of the two independent conditions in the standard Navier slip model. Hence, the number of boundary conditions is reduced and one can show that the inconsistency at the contact line can be resolved this way (see later).

As a starting point, we consider the sum of the wall and the contact line dissipation, given as

(2.19) \begin{align} \mathcal{T} = \int _{\partial \varOmega } (\boldsymbol{v}_\parallel -\boldsymbol{U}_w)\boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel \, \text{d}A + \sigma \int _{\varGamma (t)} (\cos \theta _d - \cos \theta _e) \, U_{\textit{cl}} \, \text{d}l. \end{align}

By introducing the contact line delta distribution $\delta _\varGamma$ , one can rewrite $\mathcal{T}$ as a single integral over the entire solid boundary $\partial \varOmega$ according to

(2.20) \begin{align} \begin{aligned} \mathcal{T} = \int _{\partial \varOmega } \left ( (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma \delta _\varGamma \right ) \boldsymbol{\cdot }(\boldsymbol{v}_\parallel -\boldsymbol{U}_w) \, \text{d}A. \end{aligned} \end{align}

Note that it is possible to factor out the common co-factor $(\boldsymbol{v}_\parallel -\boldsymbol{U}_w)$ because the contact line speed can be written as $U_{\textit{cl}} = (\boldsymbol{v}_\parallel -\boldsymbol{U}_w) \boldsymbol{\cdot }\boldsymbol{n}_\varGamma$ . A linear closure relation is now provided by the GNBC (a nonlinear generalisation of the closure is discussed in § 6.1)

(2.21) \begin{align} -\beta (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) = (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma \delta _\varGamma \quad \text{on} \quad \partial \varOmega \end{align}

with a friction coefficient $\beta \gt 0$ . Notice that the boundary condition given by (2.21), called ‘delta function GNBC’ to distinguish it from other variants, should be understood in the sense of distributions.

There is an interesting link between (2.21) and the independent closure (2.14) together with (2.16): following Ren (Reference Ren and W.2007), one may assume that the friction $\beta$ in (2.21) has a singular component, i.e.

(2.22) \begin{align} \beta = \beta _\varGamma \delta _\varGamma + \beta _b \quad \text{with} \quad \beta _\varGamma , \beta _b \geqslant 0. \end{align}

In this case, (2.14) takes the form

(2.23) \begin{align} [\beta _\varGamma (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) + \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma ] \delta _\varGamma + [\beta _b (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) + (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel ] = 0 \quad \text{on} \quad \partial \varOmega . \end{align}

This distributional equation can be split into a singular and a regular part. This yields the ‘standard closure relations’ for slip and dynamic contact angle, i.e.

(2.24) \begin{align} \quad \beta _b (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) + (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel = 0, \quad \beta _\varGamma (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) \boldsymbol{\cdot }\boldsymbol{n}_\varGamma + \sigma (\cos \theta _d - \cos \theta _e) = 0. \end{align}

Hence, the sharp-interface, sharp-contact line GNBC condition (2.21), and the standard closure (2.14) and (2.16) are formally equivalent if the friction $\beta$ has a singular component at the contact line. As we will see later, this singular component $\beta _\varGamma$ may arise as a ‘lumped friction’ that results from shrinking the contact zone to a line with vanishing thickness. Let us already note at this point that we are not going to study the limit of a sharp contact line, but model it as a finite contact region.

2.3. Contact region GNBC (CR-GNBC) model

To obtain the CR-GNBC model, we replace the contact line delta distribution in (2.21) by a smooth function defined over a finite transition region with characteristic width $\varepsilon$ such that

(2.25) \begin{align} \delta _\varGamma ^\varepsilon \geqslant 0, \quad \int \delta _\varGamma ^\varepsilon (x) \, \text{d}x = 1. \end{align}

Note that this approach also requires to extend the definition of the contact angle $\theta _d$ and the contact line normal $\boldsymbol{n}_\varGamma$ away from the sharp contact line. In fact, the existence of the solid boundary, touching the interface at an angle strictly between $0$ and $\pi$ , provides a means for this extension of the contact line to a finite region. Then, the deviation of the contact angle from the equilibrium value appears in the velocity boundary condition leading to a force balance between sliding friction forces due to slip along the solid boundary, the tangential component of the viscous stress at the boundary and the uncompensated Young force:

(2.26) \begin{align} - (\boldsymbol{v}_\parallel -\boldsymbol{U}_w) = 2 \lambda (\boldsymbol{D}\boldsymbol{n}_{\partial \varOmega })_\parallel + \frac {\sigma }{\beta } \big[(\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma \delta ^\varepsilon _\varGamma \big]. \end{align}

Notably, the dynamic contact angle is not prescribed explicitly in this approach. Instead, the dynamics of the contact angle is determined by (2.26) and the kinematics of the interface transport.

2.4. Kinematics of the dynamic contact angle

We derive the evolution law for the contact angle, given a sufficiently regular solution of the CR-GNBC model (2.1)–(2.7) and (2.26). Here, we consider the limiting case of a free surface flow, where one phase is assumed to be a dynamically passive gas at a constant pressure. We assume that $\delta _\varGamma ^\varepsilon$ evaluated at the contact line yields a value of $1/\varepsilon$ . Then, the CR-GNBC condition, evaluated at the contact line, reads as

(2.27) \begin{align} \beta (\boldsymbol{v}_\parallel -\boldsymbol{U}_w) + (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \frac {1}{\varepsilon } \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma = 0 \quad \text{at} \quad \varGamma . \end{align}

By taking the inner product with the contact line normal vector $\boldsymbol{n}_\varGamma$ (normal to the contact line and tangential to the solid), we obtain the relation

(2.28) \begin{align} \frac {U_{\textit{cl}}}{\lambda } + \left \langle \boldsymbol{n}_\varGamma , (\boldsymbol{\nabla }\boldsymbol{v}) \, \boldsymbol{n}_{\partial \varOmega } \right \rangle + \left \langle (\boldsymbol{\nabla }\boldsymbol{v}) \, \boldsymbol{n}_\varGamma , \boldsymbol{n}_{\partial \varOmega } \right \rangle + \frac {\sigma }{\varepsilon \eta } (\cos \theta _d - \cos \theta _e) = 0 \quad \text{at} \quad \varGamma . \end{align}

Using the kinematic evolution equation for the contact angle derived by Fricke et al. (Reference Fricke, Köhne and Bothe2019), one can show that the rate-of-change of the contact angle ${\dot {\theta }}_{\text{d}}$ is given as

(2.29) \begin{align} 2 {\dot {\theta }}_{\textit{d}} = - \left \langle \boldsymbol{n}_\varGamma , (\boldsymbol{\nabla }\boldsymbol{v}) \, \boldsymbol{n}_{\partial \varOmega } \right \rangle . \end{align}

Moreover, it follows from the impermeability condition that the term $\left \langle \boldsymbol{\nabla }\boldsymbol{v} \, \boldsymbol{n}_\varGamma , \boldsymbol{n}_{\partial \varOmega } \right \rangle$ vanishes for a flat solid boundary. Therefore, we obtain the contact angle evolution law for a regular solution of the CR-GNBC model. It reads as

(2.30) \begin{align} {\dot {\theta }}_{\textit{d}} = \frac {U_{\textit{cl}}}{2 \lambda } + \frac {1}{\varepsilon } \, \frac {\sigma }{2\eta } (\cos \theta _d - \cos \theta _e). \end{align}

2.5. Remarks

1. (i) Compared with the standard Navier slip model (see Fricke et al. Reference Fricke, Köhne and Bothe2019 for details), the uncompensated Young stress leads to an additional term in the equation for ${\dot {\theta }}_{\textit{d}}$ , which reads as

   (2.31) \begin{align} \frac {1}{\varepsilon } \, \frac {\sigma }{2\eta } (\cos \theta _d - \cos \theta _e). \end{align}
   Obviously, the latter term is negative for $\theta _d \gt \theta _e$ (and positive for $\theta _d \lt \theta _e$ ) and, hence, drives the system towards equilibrium.

2. (ii) An important consequence of the CR-GNBC for quasi-stationary states is that it defines a functional dependence between the dynamic contact angle and the contact line speed. In fact, setting ${\dot {\theta }}_{\textit{d}} = 0$ leads to the relation

   (2.32) \begin{align} \textit {Ca}_{\textit {cl}} = \frac {\eta (-U_{\textit{cl}})}{\sigma } = \frac {\lambda }{\varepsilon } \, (\cos \theta _d - \cos \theta _e), \end{align}
   or, equivalently, to
   (2.33) \begin{align} - (\beta \varepsilon )U_{\textit{cl}} = \sigma (\cos \theta _d - \cos \theta _e). \end{align}
   By comparing (2.33) with (2.16), we see that the contact line friction parameter can be identified with the product of the ‘bulk friction’ in the Navier slip condition and the width of the contact line region, i.e.
   (2.34) \begin{align} \zeta = \beta \varepsilon . \end{align}
   The latter equation has been proposed before by Blake et al. (Reference Blake, Fernandez-Toledano, Doyen and De Coninck2015) in the context of the molecular kinetic theory. Physically, it indicates that, within the present modelling framework, there is only one friction mechanism that affects both the slip at the solid boundary and the dynamics of the microscopic contact angle. In this sense, the contact line friction $\zeta$ can be understood as the lumped wall friction of the contact region.

3. (iii) From (2.29), we conclude that the stress component $\left \langle \boldsymbol{n}_\varGamma , (\boldsymbol{\nabla }\boldsymbol{v}) \, \boldsymbol{n}_{\partial \varOmega } \right \rangle$ vanishes at the contact line for quasi-stationary states, i.e. for $\dot \theta =0$ . Hence, there appears to be ‘perfect slip’ at the contact line in that case. Actually, the concepts of the ‘apparent slip length’ $\lambda _a$ (see figure 2) and the physical slip parameter defined as $\lambda =\eta /\beta$ must be distinguished for the GNBC model. In fact, the uncompensated Young stress is able to reverse the sign of the velocity gradient at the contact line. In this case, fluid particles at the solid boundary may have a larger tangential velocity than fluid particles slightly above the boundary. This situation corresponds to a negative apparent slip length (see figure 2). It is, however, caused by the uncompensated Young stress in the velocity boundary condition. The physical slip length parameter $\lambda$ is still positive and finite in all cases.

   Figure 2.

   Different cases for the apparent slip length $\lambda _a$ : positive, perfect and negative slip (reference frame with $U_w=0$ ).

   

4. (iv) Note that (2.30) can be rephrased as a generalised mobility law of the form

   (2.35) \begin{align} U_{\textit{cl}} = f(\theta _d, {\dot {\theta }}_{\textit{d}}). \end{align}
   Therefore, the contact line speed depends on the contact angle $\theta _d$ , but also on its rate-of-change ${\dot {\theta }}_{\textit{d}}$ which, in turn, can be computed from $\boldsymbol{\nabla }\boldsymbol{v}$ (see Fricke et al. Reference Fricke, Köhne and Bothe2019). In this sense, the contact line speed in the GNBC model depends on the flow field in the vicinity of the contact line.

5. (v) Moreover, the GNBC can be understood as an inhomogeneous Robin condition for the velocity. Hence, the GNBC enforces a flow whenever $\theta _d \neq \theta _e$ . In contrast to the standard Navier slip model, the GNBC model is able to describe the relaxation process of the contact angle.

2.6. The CR-GNBC thin film equation

We now derive the thin film equation for the CR-GNBC and compare it with other known thin film equations in the context of dynamic contact lines. The coordinate system is shown in figure 1. Under the thin film assumption, we consider that the pressure remains constant along the $y$ -axis and that the Laplace pressure jump across the interface can be expressed as

(2.36) \begin{align} \Delta p = - \dfrac {\sigma H^{\prime \prime }}{(1+H^{{\prime 2 }})^{3/2}} \approx - \sigma\! H^{\prime \prime }. \end{align}

The $x$ -momentum equation

(2.37) \begin{equation} -\!\frac {1}{\rho } \frac {\partial p}{\partial x} + g + \frac {\eta }{\rho } \frac {\partial ^2 v}{\partial y^2} = 0 \end{equation}

is supplemented by a free surface condition

(2.38) \begin{align} \dfrac {\partial v}{\partial y} (x,H) = 0. \end{align}

Since we are in the reference frame of the contact line, we can write the CR-GNBC as

(2.39) \begin{equation} v(x,0) - \lambda \dfrac {\partial v}{\partial y} (x,0) = U_w \left[1 - \varepsilon \,{\textit{f}} \left (\dfrac {\textit{x}}{\varepsilon } \right )\right]\!, \end{equation}

where $\textit{f} ({\textit{x} }/{\varepsilon } )$ is the smoothed Dirac function that will be defined in (3.6). For further details on the formulation in the reference frame of the contact line, we refer to Kulkarni et al. (Reference Kulkarni, Fullana and Zaleski2023). Given the contact line boundary condition, the free surface condition and the Laplace pressure jump, the $x$ -momentum (2.37) can now be written as an ordinary differential equation in terms of $H$ ,

(2.40) \begin{equation} H^{\prime \prime \prime } + \frac {1}{\textit{l}_{c}^2} = \frac {3 \textit {Ca} \left (1- \varepsilon \,\textit{f} \left ( \dfrac {\textit{x} }{\varepsilon } \right ) \right )}{H(H+3\lambda )} - \frac {3 \eta Q}{\sigma H^2 (H+3\lambda )}, \end{equation}

where $\textit{l}_{c}$ is the capillary length and $Q = -\!\int _{0}^{H} v \, \text{d}s$ is the total flux. Assuming steady state, where $Q=0$ and $H^{\prime \prime \prime } \gg 1/\textit{l}_{c} ^2$ in the vicinity of the contact line, we obtain the CR-GNBC thin film equation

(2.41) \begin{equation} H^{\prime \prime \prime } = \dfrac {3 \textit {Ca} \left (\tanh ^2 \left (\dfrac {x}{\varepsilon }\right ) \right )}{H(H+3\lambda )}. \end{equation}

From (2.41), we can see that for $x \ll \varepsilon$ , $H^{\prime \prime }$ does not diverge and approaches a constant value at the contact line ( $x=0$ ). Hence, the equation is singularity-free. Our CR-GNBC model can therefore be viewed as Navier slip with a smoothening well of width $\varepsilon$ around the contact line where the uncompensated Young stress acts. A comparison of thin film equations from the literature and their respective smoothness is presented in table 2.

Table 2.

Thin film equations for various contact line boundary conditions.

3.1. The volume-of-fluid method

The volume-of-fluid (VOF) method for representing fluid interfaces coupled with a flow solver is well known to be suited for solving interfacial flows (see e.g. Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Popinet & Zaleski Reference Popinet and Zaleski1999; Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011; Marić et al. Reference Marić, Kothe and Bothe2020). We use the free software Basilisk, a platform for the solution of partial differential equations on adaptive Cartesian meshes (Popinet Reference Popinet2009, Reference Popinet2015, Reference Popinet2018). For a two-phase flow, the volume fraction $c(\boldsymbol{x}, t)$ is defined as the integral of the first fluid’s characteristic function in the control volume. The volume fraction $c(\boldsymbol{x}, t)$ is used to define the density and viscosity in the control volume according to

(3.1) \begin{equation} \begin{aligned} \rho ({c}) & \equiv {c} \rho _1+(1-{c}) \rho _2, \\ \mu ({c}) & \equiv {c} \mu _1+(1-{c}) \mu _2, \end{aligned} \end{equation}

with $\rho _1$ , $\rho _2$ and $\mu _1$ , $\mu _2$ the densities and viscosities of the phase 1 and 2, respectively.

The advection equation for the density is then replaced by the equation for the volume fraction

(3.2) \begin{equation} \partial _{t} c+ \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} c =0. \end{equation}

The projection method is used to solve the incompressible Navier–Stokes equations combined with a Bell–Collela–Glaz advection scheme and a VOF method for interface tracking. The resolution of the surface tension term is directly dependent on the accuracy of the curvature calculation. The height-functions method, described by Afkhami & Bussmann (Reference Afkhami and Bussmann2008, Reference Afkhami and Bussmann2009), is a VOF-based technique for calculating interface normals and curvatures. About each interface cell, fluid ‘heights’ are calculated by summing fluid volume in the grid direction closest to the normal of the interface. In two dimensions, a $7 \times 3$ stencil around an interface cell is constructed and the heights are evaluated by summing volume fractions horizontally, i.e.

(3.3) \begin{equation} h_{\!j}=\sum _{k=i-3}^{k=i+3} c_{\!j, k} \varDelta , \end{equation}

with $c_{\!j, k}$ the volume fraction and $\varDelta$ the grid spacing. The heights are then used to compute the interface normal $\boldsymbol{n}_\varSigma$ and the curvature $\kappa$ according to

(3.4) \begin{equation} \begin{array}{c}{\boldsymbol{n}_\varSigma =\left (h_{x},-1\right )}, \\ \qquad {\kappa =\dfrac {h_{x x}}{\left (1+h_{x}^{2}\right )^{3 / 2}}}, \end{array} \end{equation}

where $h_{x}$ and $h_{xx}$ are discretised using second-order central differences. The orientation of the interface, characterised by the contact angle – the angle between the normal to the interface at the contact line and the normal to the solid boundary – is imposed in the contact line cell. It is important to note that a numerical specification of the contact angle affects the overall flow calculation in two ways:

1. (i) it defines the orientation of the interface reconstruction in cells that contain the contact line;

2. (ii) it influences the calculation of the surface tension term by affecting the curvature computed in cells at and near the contact line.

We now present the numerical implementation of the GNBC as written in (2.27). The boundary condition is applied on the solid surface with a smoothing function that takes into account the relative position along the boundary with respect to the contact line, denoted $x$ ,

(3.5) \begin{align} \beta (\boldsymbol{v}_\parallel -\boldsymbol{U}_w) + (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \textit{f} \left ( \dfrac {\textit{x} }{\varepsilon } \right ) \: \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma = 0 \quad \text{on} \quad \partial \varOmega , \end{align}

with $\textit{f} ({\textit{x} }/{\varepsilon } )$ the smoothed Dirac function defined as

(3.6) \begin{equation} \textit{f} \left ( \dfrac {\textit{x} }{\varepsilon } \right ) = \dfrac {\left ( 1 - \tanh ^2 \left ( \dfrac {x}{\varepsilon } \right ) \right )}{\varepsilon }. \end{equation}

This specific function is smooth, symmetric and preserves the area for varying $\varepsilon$ , characteristic features that are necessary for the well-posedness of the discrete boundary condition. The boundary condition can be expressed as an inhomogeneous Robin boundary condition for the parallel velocity $\boldsymbol{v}_\parallel$ , as outlined previously:

(3.7) \begin{align} \boldsymbol{v}_\parallel + \dfrac {1}{\beta } (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel = \boldsymbol{U}_w + \dfrac {1}{\beta }\, \textit{f} \left ( \dfrac {\textit{x} }{\varepsilon } \right ) \: \sigma (\cos \theta _e - \cos \theta _d) \, \boldsymbol{n}_\varGamma \quad \text{on} \quad \partial \varOmega . \end{align}

We use the Navier boundary condition (Navier slip) that was implemented in the same framework by Fullana et al. (Reference Fullana, Zaleski and Popinet2020) and tested as a localised slip boundary condition by Lācis et al. (Reference Lācis, Johansson, Fullana, Hess, Amberg, Bagheri and Zaleski2020). The difference lies now in the space dependent right-hand side of (3.7). The uncompensated Young’s stress that only acts at the contact line through the discrete Dirac function needs to be computed at each grid point.

The numerical approach in this study stands out for its free contact angle method. Instead of setting the dynamic angle $\theta _d$ , we reconstruct it from the interface geometry and use it as an input parameter to calculate the right-hand side of (3.7). To reconstruct such a consistent angle from the volume fraction field, we use a Taylor expansion of the contact angle along the coordinate direction normal to the boundary (see $y$ -axis in figure 3)

(3.8) \begin{equation} \theta _d = \theta _a + \delta \, \dfrac{{\textrm{d}} \theta }{{\textrm{d}} y} + O(\delta ^2), \end{equation}

where $\theta _a$ is the above-mentioned angle located at the distance $y = \delta$ from the wall. Using that ${\textrm{d}}\theta / {\textrm{d}} s = \kappa$ and projecting onto the wall-normal direction, we obtain

(3.9) \begin{equation} \theta _d = \theta _a + \delta \, \dfrac {\kappa \, \sqrt {1 + x^\prime }}{\sin \theta _a} \, + O(\delta ^2), \end{equation}

with $x^\prime$ the local slope of the interface. Discretely, by setting $\delta = 3/2 \Delta$ , such that the above-mentioned angle is computed one layer above the wall, and using height-functions to compute the slope of the interface, such that $h_y$ is the first-order derivative of the height-functions in the $y$ direction (normal to the wall), we obtain

(3.10) \begin{equation} \theta _{d} = \theta _{a} + \dfrac {3}{2} \Delta \dfrac {\kappa \: \sqrt {1+h_y^2}}{\sin \theta _{a}} + O(\varDelta ^2). \end{equation}

Figure 3.

Extrapolation of the contact angle $\theta _d$ using the angle $\theta _a$ located $3/2 \Delta$ away from the wall. $h_0$ to $h_2$ denote the horizontal heights.

Figure 3 provides a schematic illustration of this extrapolation process. Once the extrapolated angle is computed, we enforce it through appropriate local modification of height-functions in the ghost layer. Algorithm 1 is a concise summary of the two-step procedure to apply the CR-GNBC in the VOF framework.

3.2. Kinematic transport of the contact angle

We validate the free contact angle method presented in (3.10) through an analysis of the kinematic transport of the contact angle in a simplified set-up. Leveraging kinematic considerations, Fricke et al. (Reference Fricke, Marić and Bothe2020) and Fricke (Reference Fricke2021) derived analytical solutions for the transport of the contact angle and the curvature for some specific velocity fields. To validate the present approach within the VOF framework, we conduct advection test cases for an initially circular interface in contact with the domain boundary. These advection tests are carried out for various grid sizes.

Algorithm 1:

CR-GNBC pseudo-code

The set-up involves a disk with a dimensionless diameter $D = 1$ in a $2 \times 2$ domain. The centre of the disk is shifted by $\delta _s = 0.05$ above the substrate, resulting in an initial contact angle of $\theta _0 = \arccos (\delta _s/R) \approx 84.26^\circ$ . The velocity field across the entire domain is defined as

(3.11) \begin{equation} \begin{aligned} v_x &= c_1 \cos (\pi t) \, x + c_2 \cos (\pi t)y, \\ v_y &= - c_1 \cos (\pi t) y. \end{aligned} \end{equation}

Here, $v_x$ and $v_y$ represent the $x$ and $y$ components of the velocity, while $c_1$ and $c_2$ are positive constants. We aim to validate the accuracy and reliability of the angle extrapolation method under varying grid sizes, where we only consider the advection equation of the colour function (3.2).

Figure 4.

Validation of the free angle extrapolation method for varying grid sizes. (a) Temporal evolution of the contact angle. (b) Temporal evolution of the curvature. In both panels, blue corresponds to $D/\varDelta = 32$ , orange to $D/\varDelta = 64$ , green to $D/\varDelta = 128$ , red to $D/\varDelta = 256$ and the black dashed line is the analytical solution. (c) Convergence of both the angle and curvature errors with $L_2$ and $L_\infty$ norms, showing quasi-second-order convergence for the angle and quasi-first-order convergence for the curvature.

The prescribed incompressible velocity field (3.11) will induce oscillations of the interface in both vertical and horizontal directions. The angle formed at the contact line is determined by this motion and varies in time. From the relations derived by Fricke et al. (Reference Fricke, Marić and Bothe2020), we compare the observed numerical contact angle with the analytical value $\theta _{an}$ , given by the formula

(3.12) \begin{equation} \theta _{an}(t) = \dfrac {\pi }{2} + \tan ^{-1} \left ( \dfrac {- 1}{\tan \theta _0} e^{2 c_1 S(t)} + \dfrac {c_2}{2 c_1} \big ( e^{2 c_2 S(t)} - 1 \big ) \right ) \end{equation}

with

(3.13) \begin{equation} S(t) = \dfrac {\sin (\pi t)}{\pi }. \end{equation}

Moreover, we validate the evolution of the curvature by comparison with the reference one $\kappa _{an}$ , which is given as the solution of the ordinary differential equation (see Fricke Reference Fricke2021)

(3.14) \begin{equation} \dfrac {{\textrm{d}} \kappa _{an}}{{\textrm{d}} t} = -3 \, \kappa _{an} \, \cos (\pi t) \, \big[c_1 \, \cos ^2(\theta _{an}) - c_2 \, \cos (\theta _{an}) \sin (\theta _{an}) - c_1 \, \sin ^2(\theta _{an})\big] \end{equation}

with the initial condition $\kappa _0 = 2 / D = 2$ . We use Vofi (see Bnà et al. Reference Bnà, Manservisi, Scardovelli, Yecko and Zaleski2016) to initialise the volume fraction field, which yields a more accurate initial curvature. We conduct simulations with $c_1 = 0.5$ and $c_2 = 0.2$ , using the free angle extrapolation method. The simulations run until a final dimensionless time $T = 10$ and we examine the convergence of the method with grid sizes varying from 32 to 256 points per diameter (see figure 4) with a fixed time step $\delta t = 0.2 \Delta / \max (|c_1|, |c_2|)$ . The extracted contact angles and curvatures show convergence towards the analytical solution. As shown in figure 4(c), we observe quasi-second-order convergence for the contact angle and quasi-first-order convergence for the curvature, the latter is expected since $\kappa \sim \partial \theta /\partial s$ , with $s$ the arc-length. We emphasise that this validation corresponds to a purely kinematic advection problem, in which convergence is expected to degrade by one order from interface position to contact angle and to curvature, since each quantity involves an additional spatial derivative of the interface geometry. Such behaviour is specific to purely advective test cases, while in fully coupled dynamic simulations, surface-tension-induced regularisation typically improves curvature convergence (Popinet Reference Popinet2009, Reference Popinet2018).

Figure 5.

Steady-state meniscus example for $\textit {Ca}=0.1$ using the present CR-GNBC with $\varepsilon =0.05$ . The image is in the contact line’s reference frame, where the left plate is pulled up with $\bar {U_w} = \sqrt {\textit {Ca}}$ . The inset shows a zoom around the contact line with streamlines highlighting a stagnation point in the upper phase.

6.1. A nonlinear generalisation of the GNBC

As discussed in detail in § 2.2, the CR-GNBC in the form

(6.1) \begin{align} -\beta (\boldsymbol{v}_\parallel - \boldsymbol{U}_w) = (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma \delta ^\varepsilon _\varGamma \quad \text{on} \quad \partial \varOmega \end{align}

is obtained as a linear closure relation, to render the dissipation integral

(6.2) \begin{equation} \begin{aligned} \mathcal{T} = \int _{\partial \varOmega } \left ( (\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \, \boldsymbol{n}_\varGamma \delta ^\varepsilon _\varGamma \right ) \boldsymbol{\cdot }(\boldsymbol{v}_\parallel -\boldsymbol{U}_w) \, \text{d}A \end{aligned} \end{equation}

non-positive. Since, according to kinematics, the viscous stress contribution vanishes in a quasi-stationary state (see § 2), we obtain the dynamic contact angle relation

(6.3) \begin{align} - \zeta U_{\textit{cl}} = \sigma (\cos \theta _d - \cos \theta _e) \end{align}

with the contact line friction coefficient $\zeta =\beta \varepsilon$ . Notably, (6.3) is also found in the molecular kinetic theory (MKT) in the limit of low capillary number (see, e.g. Blake et al. Reference Blake, Fernandez-Toledano, Doyen and De Coninck2015). However, for higher capillary numbers, the MKT predicts that

(6.4) \begin{align} U_{\textit{cl}}=2 \kappa ^0 \varLambda \sinh \left [\sigma \left (\cos \theta _e-\cos \theta _d\right ) / (2 n k_B T) \right ]. \end{align}

In this case, the average distance and equilibrium frequency of molecular jumps are denoted by $\varLambda$ and $\kappa ^0$ , respectively. Moreover, $n$ is the number of adsorption sites per unit area, $k_B$ is the Boltzmann constant and $T$ is the absolute temperature; see Blake et al. (Reference Blake, Fernandez-Toledano, Doyen and De Coninck2015) for more details. Therefore, it is interesting to formulate a closure relation for (6.2) that leads to the relation (6.4) in quasi-stationary states. Notice that (6.4) can be linearised for $U_{\textit{cl}} \rightarrow 0$ using $\sinh (x) = x + \mathcal{O}(x^3)$ . Hence, the contact line friction coefficient is identified as $\zeta = (n k_B T)/(\kappa ^0 \varLambda )$ .

For simplicity, let us assume that $\boldsymbol{U}_w = 0$ in the following (the generalisation to $\boldsymbol{U}_w \neq 0$ is obvious). To proceed, it is useful to decompose the integral in (6.2) into its components normal and tangential to the contact line, according to

(6.5) \begin{align} \boldsymbol{v}_\parallel = (\boldsymbol{v}_\parallel \boldsymbol{\cdot }\boldsymbol{n}_\varGamma ) \, \boldsymbol{n}_\varGamma + (\boldsymbol{v}_\parallel \boldsymbol{\cdot }\boldsymbol{t}_\varGamma )\boldsymbol{t}_\varGamma . \end{align}

Here, we denote by $\boldsymbol{t}_\varGamma$ a unit tangent vector to the contact line. We obtain the representation

(6.6) \begin{equation} \begin{aligned} \mathcal{T} = \mathcal{T}_\bot + \mathcal{T}_\parallel \end{aligned} \end{equation}

with

(6.7) \begin{align} \mathcal{T}_\parallel = \int _{\partial \varOmega } \left (\boldsymbol{t}_\varGamma \boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel \right ) (\boldsymbol{t}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ) \, \text{d}A \end{align}

and

(6.8) \begin{align} \mathcal{T}_\bot = \int _{\partial \varOmega } \left (\boldsymbol{n}_\varGamma \boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \, \delta ^\varepsilon _\varGamma \right )\! (\boldsymbol{n}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ) \, \text{d}A. \end{align}

We are now looking for closure relations to ensure that $\mathcal{T}_\bot \leqslant 0$ and $\mathcal{T}_\parallel \leqslant 0$ . A general nonlinear closure for $\mathcal{T}_\bot$ reads as

(6.9) \begin{align} {\boldsymbol{n}_\varGamma \boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel + \sigma (\cos \theta _d - \cos \theta _e) \delta ^\varepsilon _\varGamma = -f(\boldsymbol{n}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ),} \end{align}

where the scalar function $f$ satisfies the inequality

(6.10) \begin{align} x f(x) \geqslant 0 \quad \forall x \in {\mathbb{R}}. \end{align}

Such a closure is consistent with the second law of thermodynamics because it implies that

(6.11) \begin{align} \mathcal{T}_\bot = - \int _{\partial \varOmega } (\boldsymbol{n}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ) f(\boldsymbol{n}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ) \leqslant 0. \end{align}

The linear version of the GNBC is recovered as $f(x)=\beta x$ . Analogously, a closure for the component tangential to the contact line can be obtained as

(6.12) \begin{align} {\boldsymbol{t}_\varGamma \boldsymbol{\cdot }(\boldsymbol{S}\boldsymbol{n}_{\partial \varOmega })_\parallel = -\!f(\boldsymbol{t}_\varGamma \boldsymbol{\cdot }\boldsymbol{v}_\parallel ).} \end{align}

Motivated by (6.4), a special choice is

(6.13) \begin{align} f(x) = b \, \text{arcsinh}(x/a), \end{align}

with positive constants $a=2 \kappa ^0 \varLambda$ and $b=2 n k_B T$ . Clearly, (6.13) reduces by linearisation to the original GNBC (6.1) with $\beta =b/a$ if $\boldsymbol{v}_\parallel \boldsymbol{\cdot }\boldsymbol{n}_\varGamma \rightarrow 0$ . Since (6.13) corresponds to the MKT (6.4) for quasi-stationary states, it may improve the standard GNBC model for higher values of the capillary number. This shall be studied in detail in the future.

6.2. Relaxation time scale of the GNBC

We can see from figure 7(b) that after $\textit {Ca}_{\textit {loc}}$ reaches a maximum value within a short-lived impulsive slip length time scale $t_{\varepsilon }$ , there is an intermediate process of capillary relaxation before entering the slow decaying relaxation to the steady state. We refer to the time scale governing this intermediate relaxation as the GNBC relaxation time scale $t_r$ . This time scale is crucial as it controls how fast a system will relax to the quasi-stationary state locally. Hence, we now examine what determines this time scale. Our formulation of the CR-GNBC model gives a contact angle evolution law (2.30) which can be re-written in terms of capillary number as

(6.14) \begin{equation} \frac {\mu }{\sigma } {\dot {\theta }}_{\textit{d}} = - \frac {\textit {Ca}_{\textit {cl}}}{2 \lambda } + \frac {1}{2 \varepsilon } \, (\cos \theta _d - \cos \theta _e), \end{equation}

with $\textit {Ca}_{\textit {cl}}$ the contact line capillary number defined in the reference frame of the solid. In the lab frame of reference, we can use $\textit {Ca}_{\textit {loc}}$ and $\textit {Ca}$ to write the angle evolution as

(6.15) \begin{equation} \frac {\mu }{\sigma } {\dot {\theta }}_{\textit{d}} = \frac {\textit {Ca} - \textit {Ca}_{\textit {loc}}}{2 \lambda } + \frac {1}{2 \varepsilon } \, (\cos \theta _d - \cos \theta _e). \end{equation}

Assuming small deviations of $\theta _d$ from $\theta _e$ , we get

(6.16) \begin{equation} \frac {\mu }{\sigma } {\dot {\theta }}_{\textit{d}} = \frac {\textit {Ca} - \textit {Ca}_{\textit {loc}}}{2 \lambda } + \frac {\sin \theta _e}{2 \varepsilon } (\theta _e - \theta _d). \end{equation}

Note that in (6.16), the time dependent quantities as $\theta _d$ , ${\dot {\theta }}_{\textit{d}}$ and $\textit {Ca}_{\textit {loc}}$ require the solution of the Stokes flow equation subject to this boundary condition which leaves us only with numerical tools. However, motivated by the results of our DNS, we propose an asymptotic expansion of the contact line capillary number in terms of the deviation of the contact angle away from the equilibrium angle. From this expansion, we can also see that the contact line is mainly driven by the uncompensated Young stress in this relaxation time scale. In the relaxation regime, one may show that the contact-line velocity admits a regular expansion in the small angle deviation

(6.17) \begin{equation} \delta = \theta _e - \theta _d, \end{equation}

of the form

(6.18) \begin{equation} \mu U_{cl} = A_0 (\sigma \delta ) + A_1 ({\sigma \delta })^2 + A_2 ({\sigma \delta })^3 + \mathcal{O}(\delta ^4), \end{equation}

where $\mu$ is the dynamic viscosity, $\sigma$ the surface tension and the coefficients $A_{i}$ are constants. This expansion in (6.18) applies strictly on the relaxation time scale $t_r$ , over which the local capillary number $\textit {Ca}_{\textit {loc}}$ monotonically decreases from an initial value towards zero. Prior to this relaxation, an impulsive time scale of order $t_{\varepsilon }$ governs the rapid rise of $\textit {Ca}_{\textit {loc}}$ from zero to a peak value $\textit {Ca}_M$ , accompanied by a corresponding change of the dynamic contact angle from $\theta _e$ to some $\theta _{M}$ . It is important to note that (6.18) is derived in the frame in which the solid substrate is stationary. In the lab frame, where the plate translates at a constant capillary number $\textit {Ca}$ , the steady-state contact-line condition requires a finite offset in the boundary condition. The mobility law at the first order in relaxation time scale thus becomes

(6.19) \begin{equation} -\!\textit {Ca}_{\textit {loc}} = A_0 (\theta _e - \theta _d) + \varPhi (\textit {Ca}), \end{equation}

where $\varPhi (Ca)$ is the offset determined by the system’s steady-state behaviour. Substituting (6.19) in (6.16), we obtain

(6.20) \begin{equation} \frac {\mu }{\sigma } {\dot {\theta }}_{\textit{d}} = \frac { A_0 (\theta _e - \theta _d) }{2 \lambda } + \frac {\sin \theta _e}{2 \varepsilon } (\theta _e - \theta _d) + f(\textit {Ca}), \end{equation}

where the steady-state behaviour, or GNBC law (2.32), gives $f(\textit {Ca}) = - ({\varepsilon }/{\lambda }) ({\textit {Ca}}/{\sin \theta _e})$ . The solution of the above-mentioned ordinary differential equation is

(6.21) \begin{equation} (\theta _e - \theta _d) = C_{0} \, e^{-\omega t} - \frac {\varepsilon }{\lambda } \frac {\textit {Ca}}{\sin \theta _e}, \end{equation}

where

(6.22) \begin{equation} \omega = \dfrac {\sigma }{2 \mu } \left ( \dfrac {A_0}{\lambda } + \dfrac {\sin \theta _e}{\varepsilon }\right ) . \end{equation}

The solution tells us that the relaxation time constant $t_r = 1/\omega$ has a contribution from both the viscous stress $\sim \sigma\! A_0 / \mu \lambda$ and Young stress $\sim \sigma \sin \theta _e / \mu \varepsilon$ . Since $A_0$ and $\sin \theta _e$ both are of order one, the time scale is of the order of the slip length-based capillary number and a contact region-based capillary number that comes out to be a few milliseconds for our flow parameters. The constant $C_0$ , a function of $\textit {Ca}$ , needs to be obtained by asymptotic matching to the impulsive time response. We leave the solution at the impulsive time scale $t_{\varepsilon }$ and a matched asymptotic solution at the full range of scales as future scope.

We now verify the existence of this relaxation time scale $t_r$ in our DNS. Re-plotting figure 7 with shifts accounting for each $\textit {Ca}$ , we see that all the curves collapse. Hence, figure 12 verifies the region’s existence, which happens at the relaxation time scale, where the linear dependence in (6.19) is justified. We get $A_0 = 0.9$ from the best fit. At $t = 0$ for the relaxation time scale, we have $\textit {Ca}_{\textit {loc}} = \textit {Ca}_M$ and $\theta _d = \theta _M$ . Since the coordinates $(\theta _M , \textit {Ca}_M)$ are decided by the limit of $t \to \infty$ from the impulsive time scale, we directly use the numerical values of $\theta _M$ obtained from our simulations. Using the $C_0$ obtained from numerics for the highest resolution case of figure 8, the derived analytical solution (6.21) shows a good agreement with DNS results in figure 13.

Figure 12.

Collapse of shifted $\textit {Ca}_{\textit {loc}}$ as a function of the angle deviation $\theta _e - \theta _d$ .

Figure 13.

Comparison of the dynamic angle $\theta _d$ from DNS (red curve) with the analytical solution (6.21) (black dashed curve) for the CR-GNBC case with $\textit {Ca} = 0.04$ , $\varepsilon = 0.05$ and $\varepsilon /\Delta = 40$ . Inset shows the same plot in log-scale.