2.1. Governing equations

We consider an incompressible Newtonian fluid of dynamic viscosity $\mu$, density $\rho$ and surface tension $\gamma$. The kinematic viscosity is denoted by $\nu$. We study the propagation of a droplet upward on an inclined plane under the action of a constant shear stress of intensity $\tau _e$. The inclination angle is denoted by $\alpha$. The $Ox$-axis is oriented upward with an angle $\beta$ with the line of the greatest slope $OX$ and the $Oz$-axis is normal to the plane. Finally, the $Oy$-axis is chosen to form a direct orthonormal basis (see figure 1). The liquid depth (in the $Oz$ direction) is denoted by $h$. The equations are written in the reference frame of the plane, which is supposed to be Galilean. The gas is supposed to impose a constant pressure and a constant shear stress. The fluid velocity $\boldsymbol {v}$ satisfies the continuity equation

(2.1)\begin{equation} \mathrm{div}\, \boldsymbol{v}=0 . \end{equation}

Denoting by $\boldsymbol {\tau }$ the viscous stress tensor and by $p$ the pressure, the Navier–Stokes equation can be written as

(2.2)\begin{equation} \rho\left[\dfrac{\partial \boldsymbol{v}}{\partial t} +{\boldsymbol{div}} \left(\boldsymbol{v} \otimes \boldsymbol{v}\right)\right] = \rho \, \boldsymbol{g}-\boldsymbol{grad}\, p+{\boldsymbol{div}} \boldsymbol{\tau}, \end{equation}

where $\boldsymbol {g}$ is the weight acceleration and $\otimes$ the tensor product. The constitutive relation is $\boldsymbol {\tau }=2 \mu \boldsymbol{\mathsf{D}}$ where $\boldsymbol{\mathsf{D}}$ is the strain rate tensor defined by

(2.3)\begin{equation} \boldsymbol{\mathsf{D}}= (1/2) [ \boldsymbol{{grad}}\, \boldsymbol{v}+(\boldsymbol{grad} \, \boldsymbol{v})^{\rm T} ]. \end{equation}

The components of the velocity are denoted by $u$, $v$ and $w$ in the $Ox$, $Oy$ and $Oz$ directions, respectively, and the components of the tensor $\boldsymbol {\tau }$ are denoted by

(2.4)\begin{equation} \boldsymbol{ \tau}= \begin{pmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{xy} & \tau_{yy} & \tau_{yz} \\ \tau_{xz} & \tau_{yz} & \tau_{zz} \end{pmatrix}. \end{equation}

Note that the continuity equation implies that $\tau _{zz}=-\tau _{xx}-\tau _{yy}$.

Figure 1. Definition sketch.

At the bottom, the no-slip and no-penetration boundary conditions hold

(2.5)\begin{equation} \boldsymbol{v}(0)=0. \end{equation}

At the free surface, the kinematic boundary condition can be written

(2.6)\begin{equation} \dfrac{\partial h}{\partial t} +u(h)\dfrac{\partial h}{\partial x} +v(h) \dfrac{\partial h}{\partial y}=w(h). \end{equation}

The dynamic boundary condition at the free surface gives three scalar equations, which are

(2.7)\begin{gather} \tau_{xz}(h)+[ p(h)-\tau_{xx}(h)]\dfrac{\partial h}{\partial x} -\tau_{xy}(h) \dfrac{\partial h}{\partial y} +\gamma \dfrac{\partial h}{\partial x} K = \tau_{ex} , \end{gather}

(2.8)\begin{gather} \tau_{yz}(h)+[p(h)-\tau_{yy}(h)]\dfrac{\partial h}{\partial y} -\tau_{xy}(h) \dfrac{\partial h}{\partial x} +\gamma \dfrac{\partial h}{\partial y} K = \tau_{ey} , \end{gather}

(2.9)\begin{gather} p(h)+\tau_{xz}(h)\dfrac{\partial h}{\partial x} + \tau_{yz}(h) \dfrac{\partial h}{\partial y} -\tau_{zz}(h)+\gamma K= \tau_{ex} \dfrac{\partial h}{\partial x} + \tau_{ey} \dfrac{\partial h}{\partial y}. \end{gather}

In these equations and in the following, $K$ denotes the total curvature

(2.10)\begin{equation} K= \mathrm{div} \frac{\boldsymbol{grad} \, h}{ \sqrt{ 1 + \boldsymbol{grad} \, h \boldsymbol{\cdot} \boldsymbol{grad} \, h}} \end{equation}

and $\tau _{ex}$ and $\tau _{ey}$ are the components in the $Ox$ and $Oy$ directions, respectively, of the shear stress $\boldsymbol {\tau _e}$ imposed on the free surface. In the following, this shear stress is supposed to be a constant.

2.2. Scaling

The equations are derived with several assumptions concerning the order of magnitude of the dimensionless parameters of the problem. The shallow-water parameter is supposed to be a small parameter so that

(2.11)\begin{equation} \varepsilon = \frac{h_0}{L} \ll 1, \end{equation}

where $h_0$ is the characteristic depth and $L$ the characteristic length in the $Ox$ and $Oy$ directions. The order of magnitude of all other dimensionless parameters will be defined by comparison with $\varepsilon$. The imposed shear stress at the free surface is large compared with the hydrostatic pressure. More precisely, we suppose that

(2.12)\begin{equation} \delta= \dfrac{\rho g h_0}{\tau_e}= O(\varepsilon^2). \end{equation}

The characteristic velocity is defined from the imposed shear stress at the surface as

(2.13)\begin{equation} u_0=\dfrac{h_0 \tau_e}{\mu}. \end{equation}

The Reynolds number defined with this velocity and with the characteristic fluid depth is supposed to be of $O(1)$, i.e.

(2.14)\begin{equation} Re=\dfrac{h_0 u_0}{\nu}=\dfrac{\rho h_0^2 \tau_e}{\mu^2}= O(1). \end{equation}

The Weber number is defined by

(2.15)\begin{equation} We =\dfrac{\rho h_0 u_0^2}{\gamma}=\dfrac{\rho h_0^3 \tau_e^2}{\gamma \mu^2}=O(\varepsilon^2). \end{equation}

It will be convenient to use the number

(2.16)\begin{equation} \kappa=\dfrac{\varepsilon^2}{We}=O(1), \end{equation}

which is of $O(1)$. The angle $\alpha$ is such that $\sin {\alpha }=O(1)$ and the same assumption is made for $\beta$, although the equations will be used in practice most of the time with $\beta = 0$. We will keep a non-zero value for $\beta$ in all the derivation process, so that the final equations do not depend on the particular orientation of the axes. We define the following dimensionless numbers:

(2.17a–c)\begin{align} \lambda_x = \frac{\delta \sin{\alpha} \cos{\beta}}{Re}=\varepsilon^2 \lambda_1 , \quad \lambda_y = \frac{\delta \sin{\alpha} \sin{\beta}}{Re}=\varepsilon^2 \lambda_2 , \quad \lambda_z = \frac{\delta \cos{\alpha} }{Re}=\varepsilon^2 \lambda_3 , \end{align}

where $\lambda _1$, $\lambda _2$ and $\lambda _3$ are of $O(1)$. The dimensionless quantities (denoted with a tilde) are defined with the following scaling:

(2.18)\begin{equation} \left. \begin{gathered} \tilde{u}=\dfrac{u}{u_0}=\dfrac{\mu u}{h_0 \tau_e}, \quad \tilde{v}=\dfrac{v}{u_0}=\dfrac{\mu v}{h_0 \tau_e}, \quad \tilde{w}=\dfrac{u_0}{\varepsilon u_0}=\dfrac{\mu w}{\varepsilon h_0 \tau_e},\\ \tilde{x}=\dfrac{x}{L}, \quad \tilde{y}=\dfrac{y}{L}, \quad \tilde{z}=\dfrac{z}{h_0}, \quad \tilde{h}=\dfrac{h}{h_0}, \quad \tilde{t}=\dfrac{u_0}{L}t =\dfrac{h_0 \tau_e}{\mu L} t, \quad \tilde K=\frac{L}{\varepsilon}K,\\ \tilde{\tau}_{xz}=\dfrac{\tau_{xz}}{\tau_e}, \quad \tilde{\tau}_{yz}=\dfrac{\tau_{yz}}{\tau_e}, \quad \tilde{\tau}_{ex} = \frac{\tau_{ex}}{\tau_e}, \quad \tilde{\tau}_{ey} = \frac{\tau_{ey}}{\tau_e},\\ \tilde{\tau}_{xy}=\dfrac{\tau_{xz}}{\varepsilon \tau_e}, \quad \tilde{\tau}_{xx}=\dfrac{\tau_{xx}}{\varepsilon \tau_e}, \quad \tilde{\tau}_{yy}=\dfrac{\tau_{yy}}{\varepsilon \tau_e}, \quad \tilde{\tau}_{zz}=\dfrac{\tau_{zz}}{\varepsilon \tau_e}. \end{gathered} \right\} \end{equation}

The pressure is dominated by the Laplace pressure and scaled accordingly

(2.19)\begin{equation} \tilde{p}=\dfrac{L^2}{\gamma h_0} p. \end{equation}

To lighten the notation, the tildes are now dropped. In dimensionless form, we can write the mass conservation equation

(2.20)\begin{equation} \dfrac{\partial u}{\partial x}+ \dfrac{\partial v}{\partial y}+\dfrac{\partial w}{\partial z} = 0 \end{equation}

and the momentum equation, in the $Ox$, $Oy$ and $Oz$ directions, respectively, as

(2.21)\begin{gather} \dfrac{\partial u}{\partial t}+\dfrac{\partial u^2}{ \partial x} + \dfrac{\partial uv}{\partial y}+\dfrac{\partial uw}{\partial z} ={-}\kappa \dfrac{\partial p}{\partial x} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{xx}}{\partial x} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{xy}}{\partial y} + \dfrac{1}{\varepsilon R_e} \dfrac{\partial \tau_{xz}}{\partial z} - \varepsilon \lambda_1 , \end{gather}

(2.22)\begin{gather} \dfrac{\partial v}{\partial t}+\dfrac{\partial uv}{ \partial x} + \dfrac{\partial v^2}{\partial y}+\dfrac{\partial vw}{\partial z} ={-}\kappa \dfrac{\partial p}{\partial y} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{xy}}{\partial x} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{yy}}{\partial y} + \dfrac{1}{\varepsilon R_e} \dfrac{\partial \tau_{yz}}{\partial z} - \varepsilon \lambda_2 , \end{gather}

(2.23)\begin{gather} \varepsilon^2\left[\dfrac{\partial w}{\partial t}+\dfrac{\partial uw}{ \partial x} + \dfrac{\partial vw}{\partial y}+\dfrac{\partial w^2}{\partial z} \right] ={-}\kappa \dfrac{\partial p}{\partial z} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{xz}}{\partial x} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{yz}}{\partial y} + \dfrac{\varepsilon}{R_e} \dfrac{\partial \tau_{zz}}{\partial z}\nonumber\\ - \varepsilon^2 \lambda_3. \end{gather}

The boundary conditions at the bottom plane become $u(0) = v (0) = w(0) =0$. The boundary conditions at the free surface can be written

(2.24)\begin{equation} w (h) =\dfrac{ \partial h }{\partial t }+ u (h)\dfrac{ \partial h }{ \partial x }+ v (h) \dfrac{\partial h }{ \partial y} \end{equation}

and

(2.25)\begin{gather} \tau_{xz}(h)+[\varepsilon \kappa Re\, p(h)-\varepsilon^2 \tau_{xx}(h)] \dfrac{\partial h}{\partial x}-\varepsilon^2 \tau_{xy}(h) \dfrac{\partial h}{\partial y}+\varepsilon \kappa Re \dfrac{\partial h}{\partial x} K = \tau_{ex} , \end{gather}

(2.26)\begin{gather} \tau_{yz}(h)+[\varepsilon \kappa Re\, p(h)-\varepsilon^2 \tau_{yy}(h)] \dfrac{\partial h}{\partial y}-\varepsilon^2 \tau_{xy}(h) \dfrac{\partial h}{\partial x}+\varepsilon \kappa Re \dfrac{\partial h}{\partial y} K = \tau_{ey}, \end{gather}

(2.27)\begin{gather} p(h)+ \dfrac{\varepsilon }{\kappa Re} \tau_{xz}(h) \dfrac{\partial h}{\partial x}+\dfrac{\varepsilon }{ \kappa Re} \tau_{yz}(h)\dfrac{\partial h}{\partial y} -\dfrac{\varepsilon }{\kappa Re} \tau_{zz}(h)+ K\nonumber\\ = \dfrac{\varepsilon }{\kappa Re} \tau_{ex}\dfrac{\partial h}{\partial x}+\dfrac{\varepsilon }{\kappa Re} \tau_{ey}\dfrac{\partial h}{\partial y} . \end{gather}

Finally, the constitutive relation leads to

(2.28)\begin{equation} \left. \begin{gathered} \tau_{xy} = \dfrac{\partial u}{\partial y}+\dfrac{\partial v}{ \partial x} ; \quad \tau_{xz} = \dfrac{\partial u}{\partial z}+ \varepsilon^2 \dfrac{\partial w}{ \partial x} ; \quad \tau_{yz} = \dfrac{\partial v}{\partial z}+ \varepsilon^2 \dfrac{\partial w}{ \partial y} ; \\ \tau_{xx} = 2 \dfrac{\partial u}{\partial x} ; \quad \tau_{yy} = 2 \dfrac{\partial v}{\partial x}; \quad \tau_{zz} ={-} \tau_{xx}- \tau_{yy} . \end{gathered} \right\} \end{equation}

2.3. Asymptotic expansions

To derive a consistent first-order model, accurate to within $O(\varepsilon ^2)$, an asymptotic method is used. The fields (velocity, pressure, viscous stress) are expanded as $X=X^{(0)}+ \varepsilon X^{(1)} +O(\varepsilon ^2)$ where $X$ refers to other $u, v, w, p, \tau _{xx} , \tau _{yy}, \tau _{zz}, \tau _{xz}, \tau _{yz}, \tau _{xy}.$

These expansions are inserted into the dimensionless equations of the flow to calculate the fields at orders $0$ and $1$.

At order $0$, the momentum equations (2.21) and (2.22) lead to

(2.29a,b)\begin{equation} \dfrac{\partial \tau_{xz}^{(0)}}{\partial z}=0 ; \quad \dfrac{\partial \tau_{yz}^{(0)}}{\partial z}=0. \end{equation}

This gives $\tau _{xz}^{(0)}=\tau _{xz}^{(0)}(h)$ and $\tau _{yz}^{(0)}=\tau _{yz}^{(0)}(h)$. These expressions can be found from the dynamic boundary conditions (2.25) and (2.26). At order $0$, we have simply

(2.30a,b)\begin{equation} \tau_{xz}^{(0)}= \tau_{ex} ; \quad \tau_{yz}^{(0)}= \tau_{ey}. \end{equation}

Since we suppose that the imposed shear stress at the free surface is a constant, these two components of the viscous stress tensor are uniform in the droplet at order $0$. From the constitutive relation (2.28) and the no-slip condition, we obtain the components $u^{(0)}$ and $v^{(0)}$

(2.31a,b)\begin{equation} u^{(0)} = \tau_{ex} z ; \quad v^{(0)} = \tau_{ey} z. \end{equation}

The linear profile of the velocity is characteristic of a planar Couette flow. The mass conservation equation (2.20) enables us to calculate $w^{(0)}$. Since $u^{(0)}$ and $v^{(0)}$ do not depend on $x$ or $y$, the integration is straightforward. Taking into account the no-penetration boundary condition, we find $w^{(0)} = 0$. At order $0$, the momentum equation (2.23) reduces to $-K \partial p^{(0)} / \partial z = 0$. This implies that $p^{(0)} = p^{(0)}(h)$. At this order, the dynamic boundary condition (2.27) is simply $p^{(0)}(h) = -K$, which gives

(2.32)\begin{equation} p^{(0)} ={-}K. \end{equation}

The pressure in the droplet is, at order $0$, entirely determined by the Laplace pressure. The remaining components of the viscous stress tensor are calculated from $u^{(0)}$, $v^{(0)}$ and the constitutive relation (2.28). The result is simply

(2.33a–d)\begin{equation} \tau_{xx}^{(0)} = 0 ; \quad \tau_{yy}^{(0)} = 0 ; \quad \tau_{zz}^{(0)} = 0 ; \quad \tau_{xy}^{(0)} = 0 . \end{equation}

At order $1$, the momentum equation (2.21) gives

(2.34)\begin{equation} \frac{\partial \tau_{xz}^{(1)}}{\partial z}=\kappa Re \frac{\partial p^{(0)}}{\partial x}. \end{equation}

This equation can be integrated together with the expression (2.32) of the pressure at order 0 and the dynamic boundary condition (2.25), which at order $1$ may be written

(2.35)\begin{equation} \tau_{xz}^{(1)} (h) ={-} \kappa Re \, p^{(0)} (h) \frac{\partial h}{\partial x }-\kappa Re \frac{\partial h}{\partial x}K = 0. \end{equation}

This leads to

(2.36a,b)\begin{equation} \tau_{xz}^{(1)} = \kappa Re \frac{\partial K}{\partial x}\left(h-z\right), \quad \tau_{yz}^{(1)} = \kappa Re \frac{\partial K}{\partial y}\left(h-z\right)\!. \end{equation}

The constitutive relation (2.28) gives at order $1$

(2.37a,b)\begin{equation} \frac{\partial u^{(1)}}{\partial z}=\tau_{xz}^{(1)} ; \quad \frac{\partial v^{(1)}}{\partial z}=\tau_{yz}^{(1)} . \end{equation}

These equations can be integrated with the no-slip condition to obtain

(2.38a,b)\begin{equation} u^{(1)}=\kappa \, Re \dfrac{\partial K}{\partial x} z \left( h-\dfrac{z}{2} \right)\!; \quad v^{(1)}=\kappa \, Re \dfrac{\partial K}{\partial y} z \left( h-\dfrac{z}{2} \right)\!. \end{equation}

It is not necessary to calculate $w^{(1)}$, $p^{(1)}$ or the other components of the viscous stress tensor at order $1$.

2.4. Depth-averaging procedure

2.4.1. Average velocity

The model is obtained by averaging over the depth the equations of the flow. For any quantity $A$, its depth-averaged value is defined by

(2.39)\begin{equation} \left\langle A \right\rangle = \dfrac{1}{h} \, \int_0^h A\, \mathrm{d}z. \end{equation}

Furthermore, we use the notation $U = \langle u \rangle \; ;\; V = \langle v \rangle$. It is necessary to expand also $U$ and $V$ as $U=U^{(0)} + \varepsilon U^{(1)} + O(\varepsilon ^2), \; V=V^{(0)} + \varepsilon V^{(1)} + O(\varepsilon ^2)$. The expressions (2.29a,b) and (2.36a,b) enable us to calculate

(2.40a–d)\begin{equation} U^{(0)} = \dfrac{1}{2} \tau_{ex}h ; \quad V^{(0)}=\dfrac{1}{2} \tau_{ey}h ; \quad U^{(1)}= \kappa \, Re \dfrac{\partial K}{\partial x} \dfrac{h^2}{3} ; \quad V^{(1)} = \kappa \, Re \dfrac{\partial K}{\partial y} \dfrac{h^2}{3}. \end{equation}

2.4.2. Mass and momentum equations

Integrating the mass conservation equation (2.20) together with the no-penetration condition at the bottom and the kinematic boundary condition at the free surface leads to the equation

(2.41)\begin{equation} \frac{\partial h}{\partial t} + \mathrm{div}\left(h\boldsymbol U \right)=0, \end{equation}

where $\boldsymbol U = (U,V)^{\mathrm {T}}$ is the depth-averaged velocity vector. The momentum equation (2.21) in the $Ox$ direction is integrated over the depth, with the no-penetration condition and the kinematic boundary condition, to obtain

(2.42)\begin{equation} \frac{\partial hU}{\partial t} + \frac{\partial h\left\langle u^2 \right\rangle}{\partial x}+\frac{\partial h \left\langle uv \right\rangle }{\partial y}+\kappa \int_0^h \frac{\partial p}{ \partial x}\, \mathrm{d}z=\frac{1}{\varepsilon Re}\left(\tau_{xz}(h)-\tau_{xz}(0)\right)+O(\varepsilon). \end{equation}

Note that the terms with the derivatives of the components $\tau _{xx}$, $\tau _{xy}$ and $\tau _{yy}$ are negligible at this order of accuracy. This implies that there is no diffusive term due to the viscosity in the model, which, apart from the capillary terms, is hyperbolic. The effect of viscosity is represented by the terms with $\tau _{xz}$, which give, by integration over the depth, a viscous friction and a driving force due to the shear at the free surface. In this equation, the pressure can be evaluated at order zero to calculate the integral

(2.43)\begin{equation} \kappa \int_0^h \frac{\partial p}{\partial x}\, \mathrm{d}z ={-}\kappa h\frac{\partial K}{\partial x}+O(\varepsilon). \end{equation}

The right-hand side of (2.42) is evaluated with the expressions at order $0$ (2.29a,b) and order $1$ of $\tau _{xz}$ and with the expression of $U^{(1)}$:

(2.44a,b)\begin{equation} \tau_{xz}(h) = \tau_{ex} + O(\varepsilon^2) ; \quad \tau_{xz}(0) = \tau_{ex} + \varepsilon \frac{3U^{(1)}}{h}+O(\varepsilon^2). \end{equation}

Since $U^{(1)}$ can be written

(2.45)\begin{equation} U^{(1)} = \frac{U-U^{(0)}}{\varepsilon}+O(\varepsilon). \end{equation}

Equation (2.42) can be written with a relaxation term as

(2.46)\begin{equation} \frac{\partial hU}{\partial t} + \frac{\partial h\left\langle u^2 \right\rangle}{\partial x}+\frac{\partial h \left\langle uv \right\rangle }{\partial y} = \frac{1}{\varepsilon Re}\left(\frac{3}{2}\tau_{ex}-\frac{3U}{h} \right)+\kappa h\frac{\partial K}{\partial x}+O(\varepsilon). \end{equation}

To calculate $\langle u^2 \rangle$, $\langle uv \rangle$ and thereafter $\langle v^2 \rangle$, $u$ and $v$ are expanded as

(2.47a,b)\begin{equation} u=U+u' ; \quad v=V+v', \end{equation}

where $u'$ and $v'$ are the deviations of $u$ and $v$, respectively, with respect to their depth-averaged values $U$ and $V$. Then $\langle u^2 \rangle = U^2 + \langle u'^2 \rangle$ since, by definition, $\langle u' \rangle = 0$. In vector form, the velocity $\boldsymbol u = (u,v)^{\mathrm {T}}$ is written $\boldsymbol u = \boldsymbol U + \boldsymbol {u}'$ where $\boldsymbol {u}'=(u',v')^{\mathrm {T}}$ is the deviation to the average velocity. Then we define the tensor

(2.48)\begin{equation} \boldsymbol{\varPhi}=\dfrac{1}{h^3} \int_0^h \boldsymbol{u}' \otimes \boldsymbol{u}' \, \mathrm{d}z, \end{equation}

which will be called thereafter enstrophy because it is homogeneous to the square of a vorticity. The components of this two-dimensional symmetrical and anisotropic tensor are defined by $\boldsymbol \varPhi = \varphi _{11} \boldsymbol {e_x}\otimes \boldsymbol {e_x} + \varphi _{12} \boldsymbol {e_x} \otimes \boldsymbol {e_y} + \varphi _{12} \boldsymbol {e_y} \otimes \boldsymbol {e_x} + \varphi _{22} \boldsymbol {e_y} \otimes \boldsymbol {e_y}.$ We can write

(2.49a–c)\begin{equation} \langle u^2 \rangle = U^2 + h^2 \varphi_{11} ; \quad \langle uv \rangle = UV + h^2 \varphi_{12} ; \quad \langle v^2 \rangle = V^2 + h^2 \varphi_{22}. \end{equation}

The enstrophy terms are not negligible because the velocity $\boldsymbol {u}$ is not constant in the depth. This implies that $\langle \boldsymbol {u}\otimes \boldsymbol {u}\rangle \neq \boldsymbol {U}\otimes \boldsymbol {U}$. In contrast, at order $0$, as shown previously, the variations of the velocity with the depth is linear as in a planar Couette flow. Furthermore, at order 1, the velocity profile can be different from a linear profile. With the enstrophy tensor, the nonlinear term is written $\langle \boldsymbol {u}\otimes \boldsymbol {u}\rangle = \boldsymbol {U}\otimes \boldsymbol {U} +h^2 \boldsymbol {\varPhi }$. The introduction of the enstrophy as an additional variable of the model guarantees a well-posed model, i.e. with an energy conservation equation (see the discussion of the two-equation model (2.56) below), in the case of a non-constant velocity profile.

In tensor form, the depth-averaged momentum equation can be written

(2.50)\begin{equation} \dfrac{\partial h \boldsymbol{U} }{\partial t} +\boldsymbol{div}\left(h\boldsymbol{U}\otimes \boldsymbol{U} +h^3 \boldsymbol{\varPhi}\right) = \dfrac{3}{ \varepsilon R_e} \left( \dfrac{\boldsymbol{\tau_{e}}}{2} -\dfrac{\boldsymbol{U}}{h} \right)+\kappa h\, \boldsymbol{grad}\,{K}+O(\varepsilon). \end{equation}

2.4.3. Enstrophy equation

The derivation of the conservative part of the equations for the tensors $h^3 \boldsymbol \varPhi$ and $h^2 \boldsymbol \varPhi$ can be found in Teshukov (Reference Teshukov2007) and for the tensor $\boldsymbol \varPhi$ in Richard, Duran & Fabrèges (Reference Richard, Duran and Fabrèges2019a), in both cases under the approximation of a weakly sheared flow, which means that the tensor $\langle \boldsymbol {u}' \otimes \boldsymbol {u'} \otimes \boldsymbol {u'} \rangle$ is negligible. In the present case, the third-order tensor $\langle \boldsymbol {u}' \otimes \boldsymbol {u'} \otimes \boldsymbol {u'} \rangle$ can be consistently evaluated at order 0. At this order, the flow is a plane Couette flow with a linear velocity profile, which implies that this third-order tensor is equal to zero at order zero. The structure of the conservative part of the equations of Teshukov can thus be found consistently even if this flow is not weakly sheared.

The enstrophy tensor is expanded as $\boldsymbol {\varPhi } = \boldsymbol {\varPhi ^{(0)}} + \varepsilon \boldsymbol {\varPhi ^{(1)}}+O(\varepsilon ^2)$ with at order $0$,

(2.51)\begin{equation} \boldsymbol{\varPhi^{(0)}}=\dfrac{1}{12} \boldsymbol{\tau_e}\otimes \boldsymbol{\tau_e}, \end{equation}

and at order $1$,

(2.52)\begin{equation} \boldsymbol{\varPhi^{(1)}}= \kappa \, Re \, \frac{h}{24} \left(\boldsymbol{\tau_e}\otimes \boldsymbol{grad}\,K+\boldsymbol{grad}\, K \otimes \boldsymbol{\tau_e} \right). \end{equation}

It follows that

(2.53)\begin{equation} \boldsymbol{\varPhi} - \frac{\boldsymbol U \otimes \boldsymbol U}{3h^2}+\frac{1}{12h^2}\left(\boldsymbol U \otimes \boldsymbol U -\frac{h^2\boldsymbol{\tau_e}\otimes \boldsymbol{\tau_e}}{4}\right)=O(\varepsilon^2). \end{equation}

The enstrophy equation can be consistently written at first order as

(2.54)$$\begin{gather} \dfrac{\partial h \boldsymbol{ \varPhi}}{\partial t} +\boldsymbol{div}( h \boldsymbol{ \varPhi} \otimes \boldsymbol{U})-2~h \left(\mathrm{div}\, \boldsymbol{U}\right)\boldsymbol{ \varPhi}+\boldsymbol{grad}\, \boldsymbol{U} \boldsymbol{\cdot} h \boldsymbol{ \varPhi} + h \boldsymbol{ \varPhi} \boldsymbol{\cdot} \left(\boldsymbol{grad}\, \boldsymbol{U}\right)^{\mathrm{T}}\nonumber\\ ={-} \dfrac{1}{ \varepsilon Re} \dfrac{B}{h}\left[ \boldsymbol{ \varPhi} -\frac{ \boldsymbol{U} \otimes \boldsymbol{U}}{3 h^2} +\dfrac{1}{12 h^2} \left(\boldsymbol{U}\otimes \boldsymbol{U} -\dfrac{h^2}{4} \boldsymbol{\tau_e} \otimes \boldsymbol{ \tau_e} \right)\right]+O(\varepsilon), \end{gather}$$

where $B$ is an arbitrary dimensionless constant. Details on this derivation are given in Appendix A.

Physically, $B$ controls the relaxation of the enstrophy $\boldsymbol { \varPhi }$ on the tensor $\boldsymbol {U}\otimes \boldsymbol {U}$. Considering a large value of $B$ yields

(2.55)\begin{equation} \boldsymbol{\varPhi} \approx \frac{1}{4h^2}\boldsymbol U \otimes \boldsymbol U + \frac{1}{48}\boldsymbol{\tau_e}\otimes \boldsymbol{\tau_e}, \end{equation}

and leads back to a two-equation system of equation for $h$ and $\boldsymbol {U}$ where the averaged momentum balance reads

(2.56)\begin{equation} \dfrac{\partial h \boldsymbol{U} }{\partial t} +\boldsymbol{div}\left(\frac{5}{4} h\boldsymbol{U}\otimes \boldsymbol{U} +\frac{h^3}{48}\boldsymbol{\tau_e}\otimes \boldsymbol{\tau_e}\right) = \dfrac{3}{ \varepsilon R_e} \left( \dfrac{\boldsymbol{\tau_{e}}}{2} -\dfrac{\boldsymbol{U}}{h} \right)+\kappa h\, \boldsymbol{grad}\,{K} . \end{equation}

Unfortunately, this two-equation system does not admit an energy conservation equation because of the factor $5/4$ instead of $1$ in the momentum flux (the justification is similar to the case studied in Richard et al. Reference Richard, Gisclon, Ruyer-Quil and Vila2019b). The only consistent two-equation system with a factor 1 in the momentum flux in front of $h\boldsymbol {U}\otimes \boldsymbol {U}$ has a flux equal to $h\boldsymbol {U}\otimes \boldsymbol {U}+(1/12)h^3\boldsymbol {\tau _e}\otimes \boldsymbol {\tau _e}$. Due to the anisotropy of the tensor $\boldsymbol {\tau _e}\otimes \boldsymbol {\tau _e}$, the second term of this flux does not behave as a pressure as, for example, the term $gh^2 \boldsymbol{\mathsf{I}}\cos \theta /2$ in the usual nonlinear shallow-water equations ($\boldsymbol{\mathsf{I}}$ being the identity tensor). As a result, even this two-equation system does not admit an energy conservation equation in conservative form.

In contrast, the three-equation system (2.41), (2.50) and (2.54) does admit an energy conservation equation, which can be written

(2.57)$$\begin{gather} \frac{\partial he}{\partial t}+\mathrm{div}(he\boldsymbol{U}+h^3 \boldsymbol{\varPhi}\boldsymbol{\cdot}\boldsymbol{U})=\left[\frac{3}{\varepsilon Re}\left(\frac{\boldsymbol{\tau_e}}{2}-\frac{\boldsymbol{U}}{h}\right)+\kappa h \boldsymbol{grad}\,K \right]\boldsymbol{\cdot} \boldsymbol{U}\nonumber\\ -\frac{1}{\varepsilon Re}\frac{B h}{2}\left(\mathrm{tr}\,\boldsymbol{\varPhi}-\frac{\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{U}}{4~h^2}-\frac{\boldsymbol{\tau_e}\boldsymbol{\cdot} \boldsymbol{\tau_e}}{48}\right) +O(\varepsilon), \end{gather}$$

where the energy $e$ is $e=\boldsymbol {U}\boldsymbol{\cdot} \boldsymbol {U}/2+h^2\mathrm {tr}\,\boldsymbol {\varPhi }/2$. Because of the existence of this energy conservation equation in conservative form, the three-equation system is clearly preferable to all consistent two-equation systems.

2.5. Inertialess limit

Before completing our flow description with a model for partial wetting conditions, let us underline the link between the shallow-water three-equation model (2.41), (2.50) and (2.54), which we have derived so far, and the surface equations that are usually employed based on lubrication theory. A more complete justification of the link between shallow-water and lubrication equations can be found in Bresch & Noble (Reference Bresch and Noble2007).

In the limit of a vanishing Reynolds number, the momentum balance (2.50) reduces to

(2.58)\begin{equation} \boldsymbol{U} = \frac{h}{2}\boldsymbol{\tau_e} + \epsilon^3 C\!a \frac{h^2}{3} \boldsymbol{grad}\,{K}, \end{equation}

where $C\!a ={\gamma }/{\mu u_0}$ is a capillary number. Following the usual approximation of the lubrication theory, the curvature $K$ of the free surface can be approached by $\Delta h$. As a consequence, the mass balance (2.41) yields (Oron et al. Reference Oron, Davis and Bankoff1997)

(2.59)\begin{equation} \partial_t h + \boldsymbol{div} \left(\frac{h^2}{2}\boldsymbol{\tau_e} + \epsilon^3 C\!a\, \frac{h^3}{3} \boldsymbol{grad}\, \Delta h \right) = 0, \end{equation}

which corresponds to the surface equation used by Espín & Kumar (Reference Espín and Kumar2017) in the appropriate limit (absence of gravity, absorption and disjoining pressure).

Interestingly, shallow-water systems of equations, equivalent to lubrication equations of the form (2.59) in the limit of a vanishing Reynolds number, have been used to prove global existence of non-negative weak solutions of lubrication equations (Bresch et al. Reference Bresch, Colin, Msheik, Noble and Song2019).

2.6. Disjoining pressure

In order to account for dewetting phenomena and the displacement of contact lines, we introduce a regularisation of the jump of surface energy from the liquid–gas interface ($\gamma _{lg} = \gamma$) to the solid–gas interface ($\gamma _{sg}= \gamma \cos (\theta _s) < \gamma _{sl}$) by means of a disjoining energy $e_d(h)$ function of the free surface elevation. For convenience as will be discussed later in § 2.7, we set the reference of the dimensionless surface energy to zero at a flat solid–gas interface, so that the dimensionless surface energy at the wall is equal to $\kappa (\cos (\theta _s) -1)$.

We adopt the classical formulation proposed by Derjarguin (Churaev & Sobolev Reference Churaev and Sobolev1995), where the disjoining energy density reads

(2.60)\begin{equation} e_d(h)=\frac{(n-1)(m-1)}{n-m}\, \kappa [\cos(\theta_s)-1] \Bigg[\dfrac{1}{1-n} \left( \dfrac{h^*}{h}\right)^{n-1} - \dfrac{1}{1-m}\left( \dfrac{h^*}{h}\right)^{m-1}\Bigg], \end{equation}

where $n> m$, and varies accordingly from $0$ in the bulk of the liquid ($h \gg h^*$) to $\kappa [\cos (\theta _s) -1]$ at the precursor film thickness $h^*$. The precursor film $h=h^*$ thus plays the role of the solid–gas interface. The introduction of a precursor film is an elegant way to deal with the singularity of the viscous stress at the contact line. Another approach is the introduction of a slip at the wall. Comparisons of the slip and precursor-film approaches show that they are more or less equivalent (Sibley et al. Reference Sibley, Nold, Savva and Kalliadasis2015). Yet Diez et al. (Reference Diez, Kondic and Bertozzi2000) found the former less numerically demanding.

Associated to the disjoining energy is a disjoining pressure

(2.61)\begin{equation} \varPi(h) ={-}\frac{{\rm d} e_d}{{\rm d}h}= \frac{(n-1)(m-1)}{n-m}\, \frac{\kappa [1- \cos(\theta_s)]}{h^*} \left[\left( \dfrac{h^*}{h}\right)^{n} -\left( \dfrac{h^*}{h}\right)^{m}\right] \end{equation}

so that the pressure in the liquid reads $p= -K - \varPi (h)$.

The disjoining pressure is negative for $h> h^*$ and positive for $h < h^*$ which guarantees the stability of the precursor film. At a contact line, the disjoining pressure promotes a gradient of pressure which drives the liquid out of the precursor film.

In line with the observations of Diez et al. (Reference Diez, Kondic and Bertozzi2000), we observe that the Frumkin–Derjaguin disjoining energy (2.60) enables to take quite large values $h^*$ of the precursor film while modelling correctly the apparent static contact angle. In addition, the precursor-film model acts as a low-pass filter: drops with a height which is less than $2 h^*$ are absorbed by the precursor film. As a consequence, with this model, the minimal size of the droplets that are represented is controlled by $h^*$.

2.7. Augmented formulation

In order to mimic the dynamics of shear driven sliding droplets, a mathematical framework has been developed. The derived formulation needs to overcome the limitations of previous attempts, especially the correct inclusion of surface tension, which is required in order to capture correctly the shape of the droplets. A second requirement is to obtain a conservative hyperbolic formulation and to guarantee the consistency with the long-wave expansion in the appropriate limit. These properties would enable to develop efficient numerical schemes which preserve the energy of the flow.

Lastly, partial wetting has to be accounted for and must enable to capture the hysteresis of the static contact angle. This last requirement will allow for a correct description of droplets shedding and water accumulation. Bresch et al. (Reference Bresch, Cellier, Couderc, Gisclon, Noble, Richard, Ruyer-Quil and Vila2020) developed an augmented formulation for surface tension, accounting for the full curvature of the free surface, thus improving over the initial formulation of Noble & Vila (Reference Noble and Vila2014) and Bresch et al. (Reference Bresch, Couderc, Noble and Vila2016) which was limited to linearised curvature of the free surface in the long wave limit.

In short, the idea is to introduce an additional variable with the dimension of a velocity whose kinetic energy is equal to the surface energy of the film. By writing a transport equation for this additional variable, we are able to recast shallow-water equations with full surface tension terms into an ‘augmented’ system of equations with a skew-symmetric structure with respect to the $L^{2}$ scalar product which makes the proof of energy estimates. In addition, the design of compatible numerical scheme is made easier as surface tension terms are then recast as generalised diffusion terms.

Following Bresch et al. (Reference Bresch, Cellier, Couderc, Gisclon, Noble, Richard, Ruyer-Quil and Vila2020), we introduce a vector variable $\boldsymbol {W}$ which is colinear to the gradient of the free-surface location, $\boldsymbol {p} = \boldsymbol {grad} h$ and verifies

(2.62a)\begin{gather} \boldsymbol{W} = \dfrac{\sqrt{\kappa}}{\sqrt{h}} \, \alpha(q^2) \boldsymbol{p}, \end{gather}

(2.62b)\begin{gather}\text{with}\quad \alpha(q^2)=\dfrac{\sqrt{2}}{q}\left(\sqrt{1+q^2}-1\right)^{1/2}=\sqrt{2}\left(\sqrt{1+q^2}+1\right)^{{-}1/2}, \end{gather}

(2.62c)\begin{gather}q=\|\boldsymbol{grad} h\| = \|\boldsymbol{p}\|. \end{gather}

Note that, in this context, $\boldsymbol {W}$ has the dimension of a velocity and transforms the capillary energy density into a virtual kinetic energy as

(2.63)\begin{equation} \tfrac{1}{2}h\left\Vert \boldsymbol{W}\right\Vert ^{2}= \kappa\left(\sqrt{1+q^{2}}-1\right). \end{equation}

Again, in order to make reversible the relation between the additional velocity $\boldsymbol {W}$ and the gradient $\boldsymbol {grad}\, h$, the constant $\kappa =\gamma /\rho$ has been subtracted to the usual definition of the surface capillary energy, which does not affect the dynamics of the flow.

The system of dimensionless equations (2.41), (2.50) and (2.54) is then modified into

(2.64a)\begin{gather} \frac{\partial h}{\partial t} + \mathrm{div}\left(h\boldsymbol U \right)=0, \end{gather}

(2.64b)\begin{gather} \dfrac{\partial h \boldsymbol{U} }{\partial t} +\boldsymbol{div}(h\boldsymbol{U}\otimes \boldsymbol{U} +h^3 \boldsymbol{\varPhi}) = \dfrac{3}{R_e} \left( \dfrac{\boldsymbol{\tau_{e}}}{2} -\dfrac{\boldsymbol{U}}{h} \right) + h \boldsymbol{grad} \varPi(h) \nonumber\\ +\,\boldsymbol{div}[h \left( \boldsymbol{grad}( \boldsymbol{\mathsf{f}}_{\boldsymbol{\mathsf{1}}} \boldsymbol{\cdot} \boldsymbol{W})\right)^{\mathrm{T}}]-\boldsymbol{grad}(\,\boldsymbol{f}_{\boldsymbol{2}}\boldsymbol{\cdot} \boldsymbol{W}), \end{gather}

(2.64c)\begin{gather} \dfrac{\partial h \boldsymbol{W}}{\partial t} +\boldsymbol{div}\left(h\boldsymbol{W}\otimes \boldsymbol{U} \right)={-}\boldsymbol{\mathsf{f}}_{\boldsymbol{\mathsf{1}}} \boldsymbol{\cdot} \boldsymbol{div} [h \left( \boldsymbol{grad}\, \boldsymbol{U}\right)^{\rm T}]-\boldsymbol{f_2}\, \mbox{div} \,\boldsymbol{U}\,, \end{gather}

(2.64d)\begin{gather} \dfrac{\partial h \boldsymbol{ \varPhi}}{\partial t} +\boldsymbol{div}( h \boldsymbol{ \varPhi} \otimes \boldsymbol{U})-2 h \left(\mathrm{div}\, \boldsymbol{U}\right)\boldsymbol{ \varPhi}+\boldsymbol{grad}\, \boldsymbol{U} \boldsymbol{\cdot} h \boldsymbol{ \varPhi} + h \boldsymbol{ \varPhi} \boldsymbol{\cdot} \left(\boldsymbol{grad}\, \boldsymbol{U}\right)^{\mathrm{T}} \nonumber\\ ={-} \dfrac{1}{Re} \dfrac{B}{h}\left[ \boldsymbol{ \varPhi} -\frac{ \boldsymbol{U} \otimes \boldsymbol{U}}{3~h^2} +\dfrac{1}{12~h^2} \left(\boldsymbol{U}\otimes \boldsymbol{U} -\dfrac{h^2}{4} \boldsymbol{\tau_e} \otimes \boldsymbol{ \tau_e} \right)\right], \end{gather}

where the second-order tensor $\boldsymbol{\mathsf{f}}_{\boldsymbol{\mathsf{1}}}$ is given by

(2.64e)\begin{equation} \boldsymbol{\mathsf{f}}_{\boldsymbol{\mathsf{1}}}(h, \boldsymbol{W})= \sqrt{\kappa} \sqrt{h}\left(1+\dfrac{h}{4 \kappa}\| \boldsymbol{W}\|^2 \right) ^{{-}1/2} \left( \boldsymbol{\mathsf{I}} - \dfrac{h}{4 \kappa} ( 1+\dfrac{h}{2 \kappa} \|\boldsymbol{W}\|^2)^{{-}1} \boldsymbol{W} \otimes \boldsymbol{W} \right) \end{equation}

and the vector $\boldsymbol {f_2}$ by

(2.64f)\begin{equation} \boldsymbol{f_2} (h, \boldsymbol{W})= \dfrac{h \boldsymbol{W}}{2} \left( 1+\dfrac{ h}{2 \,\kappa} \|\boldsymbol{W}\|^2 \right)^{{-}1}. \end{equation}

An advantage of the above formulation is to commute the third-order capillary term $\boldsymbol {grad}\,{K}$ in the averaged momentum balance into second-order generalised diffusion terms in (2.64b), a simplification that turns out to be useful in numerical simulations, especially on unstructured meshes. The principal advantage of (2.64) is the skew-symmetry of these second-order terms with respect to the $L^2$ scalar product and allows for the construction of numerical schemes which conserve the energy (see the following).

Note that the definitions (2.62a), (2.64e) and (2.64f) of the additional velocity $\boldsymbol {W}$ and functions $\boldsymbol{\mathsf{f}}_{\boldsymbol{\mathsf{1}}}$ and $\boldsymbol {f_2}$ remain unchanged when the kinematic surface tension $\gamma /\rho$ is substituted for $\kappa$.

Dimensional equations are given in Appendix B.

3.1. Contact angle hysteresis

The adherence of the sessile drops is the hallmark of static contact angle hysteresis. Small drops have the ability to deform and resist the onset of motion due to the shear of the gas flow since a pressure gradient may be sustained between the front and the rear of the drop as a result of the difference between the contact angle at its front and back.

Numerical implementation of contact angle hysteresis is generally based on an evaluation of the triple line orientation with respect to the orientation of the flow (Ding & Spelt Reference Ding and Spelt2008). For instance, using a precursor film formulation and a Frumkin–Derjaguin disjoining pressure, Ahmed et al. (Reference Ahmed, Sellier, Jeremy and Taylor2014) tracked the location of the contact line point at which the contact line was orthogonal to the flow direction, thus defining the back and front of the drop. However, Ahmed et al. computed only one drop at a time and this approach seems difficult to apply for several drops when coalescence or splitting may occur.

As a first approach, the orientation of the contact line with respect to the flow has been evaluated by computing ${\boldsymbol U}\boldsymbol{\cdot} \boldsymbol {grad}\, h$. The front of the drop is identified as ${\boldsymbol U} \boldsymbol{\cdot} \boldsymbol {grad}\, h<0$ and the static contact angle $\theta _s$ is then set to it advancing value $\theta _a$. Conversely, if ${\boldsymbol U} \boldsymbol{\cdot} \boldsymbol {grad}\, h>0$, $\theta _s$ is adjusted to the receding contact angle $\theta _r$. However, this method leads to numerical difficulties. The first is the sharp jump of surface energies at the edges of the drop, i.e. for ${\boldsymbol U}\perp \boldsymbol {grad}\, h$. Small fluctuations of the orientation of the contact line at these edges thus promote sharp surface forces, which generate oscillations of the contact line. These numerical spurious waves propagate, grow and may lead to failures of the numerical simulations.

In order to avoid this problem, the front and back of the moving drops are identified based on $\textrm {div}(h {\boldsymbol U}) = -\partial h / \partial t$. Thus, the front of the drop is identified if mass is gained, i.e. $\textrm {div}(h {\boldsymbol U}) <0$. The limit of the front and back regions of the drop does not depend on the orientation of the contact line with this definition. A direct consequence is that the local contact angle may not switch sharply from $\theta _a$ to $\theta _r$ by changing the orientation of the contact line. This proves to be sufficient to eliminate the occurrence of spurious oscillations at the contact line.

In practice, the discontinuity of the value of the static contact angle has been regularised with a hyperbolic tangent:

(3.1)\begin{equation} \theta_s = \frac{\theta_a+\theta_r}{2} + \frac{\theta_r-\theta_a}{2}\tanh\left[\frac{{\rm div}(h {\boldsymbol U})}{\varepsilon'}\right], \end{equation}

where $\varepsilon '$ has been taken sufficiently low.

3.2. Simulations of single drops

Simulations of the dynamics of sessile drops under the action of a shear stress have been conducted initially without taking into account the hysteresis of the contact angle.

Convergence of our numerical scheme with respect to the grid mesh has been checked by conducting simulations in one dimension (i.e. assuming all derivatives $\partial _y=0$ with only one cell in the $y$ direction) with periodic boundary conditions. A spherical cap drop is placed at initial time in the numerical domain. The characteristic height of the drop is $h_0=0.1$ mm. Physical properties correspond to water ($\mu = 10^{-3}$ Pa s, $\rho =10^3\,\textrm {kg}\,\textrm {m}^{-3}$ and $\gamma = 75\,\textrm {mN}\,\textrm {m}^{-1}$) and the wall shear stress $\tau _e$ is adjusted to $\tau _e=8$ Pa s. The Reynolds and Weber numbers are thus $Re=80$ and $We=0.85$. These values are unchanged throughout this section.

The numerical domain extension is $7.2$ mm. The precursor film thickness is set to $h^* = 0.05$, which corresponds to $5\,\mathrm {\mu }\textrm {m}$. Figure 2 presents the result of our convergence test with respect to the mesh size $\varDelta$, or number of mesh nodes $N$. The final shape of the droplet is reported at the end of the simulation when a constant shape and speed are achieved. However, low values of $N$ are associated to spurious oscillations of the drop. A number $N=720$ of nodes is sufficient to approach satisfactorily the shape of the sliding drop, as simulations with $N=720$ and $N=1440$ provides close final shapes of the drop. Numerical tests with different precursor film thicknesses (not shown) indicate that the mesh size $\varDelta$ must be of the order of $h^*$ in order to prevent the onset of spurious oscillations which emanate from an inadequate representation of the contact line region of the drop. We thus conclude that a ratio $\varDelta /h^*=2$ is sufficient to capture the transition region which replaces a sharp contact line within the disjoining pressure description of partial wetting. A value of $\varDelta /h^* = 2$ has thus been used hereinafter for every simulations presented in this work. However, a lower resolution, i.e. $\varDelta /h^* =4$ is still acceptable.

Figure 2. Convergence test in one dimension. Here $L=7.2$ mm, contact angle $\theta _s=30^\circ$ and $\varDelta /h^*= 2$ for $N=720$.

Figure 3 presents the distribution of the velocity $U$ and enstrophy $\varphi$ within the droplet at final time. The velocity $U$ is almost constant in the bulk of the droplet and varies only at the front and back of the drop as the film thickness reduces to precursor film thickness. Indeed, for travelling-wave solutions, i.e. solutions whose shape and speed do not vary with time, the mass balance (2.41) can be rewritten in the moving frame $\xi = x - u_{drop} t$ and integrated to yield

(3.2)\begin{equation} h (U- u_{drop}) = h^* \left( \frac{\tau_e h^*}{2} - u_{drop} \right), \end{equation}

where the constant flow rate in the moving frame of reference is expressed on the right-hand side of (3.2) from the Nusselt solution $h=h^*$. Therefore, $h^* \ll 1$ implies that $U\approx u_{drop}$ whenever $h \gg h^*$.

Figure 3. Distribution of averaged velocity $U$, enstrophy weighted by the square of the thickness, $h^2 \varphi$ and shape factor $\alpha = 1 + h^2 \varphi /U^2$ (the red horizontal line shows the value $\alpha =4/3$) for the 1-D droplet of figure 2 ($N=1440$, $h^*=0.05$). The dashed black curves show the depth profile of the droplet.

Instead of plotting $\varphi$, we choose to show $h^2 \varphi$ which corresponds to the contribution of the deviations of the velocity profile $u-U$ to the kinetic energy. We compare $h^2 \varphi$ to $U^2$ by computing the shape factor

(3.3)\begin{equation} \alpha = \frac{\langle u^2 \rangle}{\langle u\rangle^2} = 1 + \frac{h^2\varphi}{U^2}. \end{equation}

The shape factor $\alpha$ does not deviate much from the value $4/3$ corresponding to the Couette flow solution (2.29a,b). As a consequence, $\varphi$ adopts large values at the edges of the droplet, so that displaying $\varphi$ instead of $h^2 \varphi$ gives the false impression that the enstrophy is negligible in the bulk of the droplet. A value of the shape factor equal to $4/3$ indicates a linear velocity profile in the depth. A shape factor smaller than $4/3$ indicates a convex velocity profile, i.e. with a larger shear near the bottom than close to the free surface, while the opposite is true for a shape factor larger than $4/3$, which denotes a concave profile, with a stronger shear near the free surface than near the bottom. The distribution of the shape factor shows that there is a strong shear effect close to the free surface in the centre of the droplet, resulting in a shape factor larger than $4/3$. At the edges of the droplet, the curvature of the velocity profile has the opposite sign.

Figure 4 illustrates the influence of the precursor film thickness on the shape of one-dimensional (1-D) drops which have reached a steady state. For a large value of $h^*$, the extension of the drop length is larger, and the front is more rounded. The front shape of the drop is more impacted than its back shape by the precursor film thickness. The maximum elevation of the drop also decreases as $h^*$ is raised.

Figure 4. Convergence test in one dimension. Here $L=7.2$ mm, contact angle $\theta _s=30^\circ$ and $\varDelta /h^*= 2$. Comparison of the drop shapes for different precursor film thicknesses $h^*$. All simulations have reached steady states. The drop profiles have been aligned to allow for an easy comparison. The red dashed line is the shape of the drop for $h^*=0.01$.

Figure 5 discusses the influence of the precursor film thickness $h^*$ on the drop velocity $u_{drop}$ and its extension length $L$. Both of them decrease monotonically with the precursor film thickness. In agreement with Schwartz & Eley (Reference Schwartz and Eley1998) and Sellier et al. (Reference Sellier, Grayson, Renbaum-Wolff, Song and Bertram2015), we observe that $1/u_{drop}$ varies linearly with the logarithm of $h^*$ as expected from (C6) as a consequence of the dependency of the wall shear dissipative work with respect to $h^*$ (see Appendix C).

Figure 5. Convergence test in one dimension. Here $L=7.2$ mm, contact angle $\theta _s=30^\circ$ and $\varDelta /h^*= 2$: (a) inverse of drop velocity $u_{drop}$ as a function of the precursor film thickness $h^*$; (b) spread length versus $h^*$.

We next turn to the simulation of two-dimensional drops. The precursor film thickness is set to $h^*=0.05$, a reasonable trade-off between accuracy and numerical cost as the mesh size $\varDelta$ must be chosen of the order of $h^*$. Two-dimensional simulations have been performed in a rectangular domain with periodic boundary conditions. At initial time, a single drop is placed in the numerical domain. Figure 6 illustrates the evolution of the drop as time is increased. The drop has initially a spherical cap shape, corresponding to the equilibrium solution in the absence of a shear, with an apparent radius $R_\theta =0.8$ mm. A constant static contact angle $\theta _s=30^\circ$ is again imposed.

Figure 6. Snapshots of the free surface elevation at different times in a domain of size $7.2\,\textrm {mm}\times 2.4\,\textrm {mm}$ with $3N\times N\approx 2\times 10^6$ nodes for a contact angle $\theta _s=30^\circ$.

The numerical domain is $7.2\,\textrm {mm}\times 2.4\,\textrm {mm}$, or $72$ by $24$ units. The mesh size to precursor film thickness is fixed to $\varDelta /h^*=2$ so that the numerical domain is discretised with $N=2160$ nodes in the $x$ direction and $720$ nodes in the $y$ direction. Mesh cells are thus squares of dimension $10\,\mathrm {\mu }\textrm {m}$. The sessile drop is initially elongated by the shear, accelerates and, after some oscillations, achieves a constant shape and speed as can be observed from figure 6. The final shape of the drop is reminiscent of the teardrop shape of sessile drops entrained by gravity at moderate inclination angles (Ahmed et al. Reference Ahmed, Sellier, Jeremy and Taylor2014).

Figure 7 presents the components of the velocity field $\boldsymbol {U}$ as well as $h^2\boldsymbol {\varPhi }$. The velocity field is nearly constant everywhere and equal to the velocity of the droplet, except in the vicinity of the contact line. Similarly $h^2\varphi _{12}$ and $h^2\varphi _{22}$ are much smaller than the component $h^2\varphi _{11}$. As for the 1-D droplets discussed in the previous section, the contribution of the enstrophy to kinetic energy $h^2\mathrm {tr}\,\varPhi /2$ is significant throughout the droplet.

Figure 7. Velocity components, $U$ and $V$, and enstrophy components $h^2\varphi _{11}$, $h^2\varphi _{22}$ and $h^2\varphi _{12}$ for the droplet in steady state presented in figure 6.

3.3. Hysteresis effect on the drop shape

We next turn to the simulations of single sliding drops in the presence of a hysteresis of the static contact angle. We replicate the simulation presented in figure 6 with a hysteresis range $2\delta _s = 10^\circ$ so that the advancing and receding contact angles are $\theta _a = 35^\circ$ and $\theta _r = 25^\circ$. We again observe a strong initial elongation of the drop as it starts moving under the action of the gas shear stress, followed by oscillations leading finally to a drop propagating with a constant shape and speed (see figure 8). In contrast to the no hysteresis case, we observe a very different shape of the drop at the final stage of the simulation, which we may labelled as inverted-teardrop shape or bullet shape, with a rounded back and an ogival front. Similar observations have been made by Ding & Spelt (Reference Ding and Spelt2008) in their simulations of the onset of motion of shear-driven droplets with contact-angle hysteresis.

Figure 8. Snapshots of the free surface elevation at different times in a domain of size $7.2\,\textrm {mm}\times 2.4\,\textrm {mm}$ with $3N\times N\approx 2\times 10^6$ nodes for a hysteresis $2\delta \theta _s=10^\circ$ ($\theta _r=25^\circ$ and $\theta _a=35^\circ$).

As discussed in § 3.1, the simulation of sliding drops in the presence of a hysteresis has been conducted considering the dependence of the static contact angle on $\textrm {div}\, h {\boldsymbol U}$ given by (3.1) with a regularisation parameter $\varepsilon =10^{-3}$. We checked that the front and rear of the drop at the end of the simulation were identified similarly by spotting the accumulation of mass in time (sign of $\textrm {div}\, h {\boldsymbol U}$) or spotting the orientation of the slope of the free surface (sign of ${\boldsymbol U}\boldsymbol{\cdot}\boldsymbol {grad}\, h$). Figure 9 presents the drop elevation $h$, and the distributions of ${\boldsymbol U}\boldsymbol{\cdot}\boldsymbol {grad}\, h$ and $\textrm {div}\, h {\boldsymbol U}$ within the drop, demonstrating that the locations of the contact line at which mass is accumulated ($\textrm {div}\, h {\boldsymbol U}<0$) coincides with the front of the drop (${\boldsymbol U}\boldsymbol{\cdot}\boldsymbol {grad}\, h<0$). In particular, the region at which the static contact angle varies continuously from $\theta _r$ to $\theta _a$ is restricted to the side edges of the drop (see figure 9b).

Figure 9. Stationary drop at the end of a simulation for a hysteresis $2\delta \theta _s=10^\circ$ ($\theta _r=25^\circ$ and $\theta _a=35^\circ$): (a) elevation $h$; (b) locations where $1.2 h^* < h<2 h^*$ and $|\textrm {div}(hU)|<\epsilon$; (c) $\textrm {div}(hU)$; (d) $U\boldsymbol{\cdot} \textrm {grad} h$.

As the hysteresis of the static contact angle is increased, the sliding drop continues to modify its shape, the widest section of the drop moving further to the back. The front adopts a more and more ogival shape. The width and maximum elevation of the drop increase, whereas its length diminishes (see figure 10). The variation of the drop length $L$, width $l$ and wetting area $S$ is reported in table 1.

Figure 10. Snapshots of the free surface elevation at $t=350$ in a domain of size $2.4\,\textrm {mm}\times 7.2\,\textrm {mm}$ with $N\times 3N=3\times 10^4$ nodes for different hysteresis of the static contact angle $2\delta \theta _s=0^\circ$, $10\,^\circ$ and $18\,^\circ$.

Table 1. Drop shape parameters for different hysteresis of the static contact angle $2\delta \theta _s=0^\circ$, $10^\circ$ and $18^\circ$.

3.4. Dynamics of coarsening

We next investigate the coarsening dynamics of a cloud of small droplets deposited at initial stage on a planar substrate. The parameter set is again chosen to correspond to water droplets, $\mu =10^{-3}$ Pa s, $\rho =1000\,\textrm {kg}\,\textrm {m}^{-3}$, $\gamma =75\,\textrm {mN}\,\textrm {m}^{-1}$, undergoing a constant shear stress $\tau _e = 8$ Pa. The numerical domain has a size of $8\,\textrm {mm} \times 24\,\textrm {mm}$. The initial condition consists in a set of $1000$ droplets for a total mass $0.22\,\textrm {mm}^3$, so that the mean apparent radius of a sessile droplet is $80\,\mathrm {\mu }\textrm {m}$. The mesh size is $\varDelta =10\,\mathrm {\mu }\textrm {m}$ and the precursor film thickness is again $h^*=5\,\mathrm {\mu }\textrm {m}$. These simulations are costly as the number of nodes reach $N \times 3 N=2\times 10^6$.

Figure 11 illustrates the observed dynamics in the absence of contact angle hysteresis. After an initial absorption of the droplets into the precursor film, droplets of minimal thickness roughly equal to $15\,\mathrm {\mu }\textrm {m}$, that is $3\times h^*$ start to emerge. This behaviour results from the disjoining pressure modelling of partial wetting, which acts as a low-pass filter. The attractive part of the disjoining pressure leads to the absorption of droplets whose thickness is below two to three times the precursor film thickness $h^*$. However, the mass of the film is conserved, and the absorbed droplets raise the mean thickness of the film above $h^*$ leading to its instability and the formation of larger sliding drops. These drops have a typical maximal height ranging from $5$ to $10$ times the thickness $h^*$ of the precursor film. The later evolution of the film is characterised by the merging of drops as large drops tend to move faster than smaller ones, catch them and coalesce. As a result, fewer and fewer drops are observed.

Figure 11. Snapshots of the free surface elevation at different times in a domain of size $24$ mm $\times\, 8$ mm with $3N\times N=2\times 10^6$ nodes for a constant contact angle $\theta _s=30^\circ$.

In the presence of a contact angle hysteresis, a similar evolution of the coarsening dynamics is observed. However, this scenario seems to be delayed in that case. For instance, comparing figures 11 and 12, one may note that the evolution of the merging process observe without contact angle hysteresis at time $t=250$ (seven drops) is reached only at time $t=750$ for $2\delta \theta _s=14^\circ$. This slowing down is a direct consequence of the lowering of the speed of the drops in the presence of a contact angle hysteresis, as one can anticipate from the drop speed estimate (C6).

Figure 12. Snapshots of the free surface elevation at different times in a domain of size $24\,\textrm {mm}\times 8\,\textrm {mm}$ with $3N\times N=2\times 10^6$ nodes for a hysteresis $2\delta \theta _s=14^\circ$ ($\theta _r=23^\circ$ and $\theta _a=37^\circ$).

To further understand the effect of the contact angle hysteresis on the drop merging process, we have performed a parametric study of the speed of a single sliding drop for different values of $2\delta \theta _s$ as the mass of the drop is varied. Figure 13 reports the result of this parametric study and shows a dramatic effect of the contact angle hysteresis on small drops, whose speed is lowered significantly. However, for large drops, this effect is mitigated by their size. As a consequence, the slowing down of the coarsening dynamics by the contact angle hysteresis is more efficient in the early stages of the coalescence process.

Figure 13. Droplet terminal velocity in a domain of size $24\,\textrm {mm}\times 8\,\textrm {mm}$ with $3N\times N=2\times 10^6$ nodes for hysteresis $2\delta \theta _s=0^\circ$, $2\delta \theta _s=5^\circ$ and $2\delta \theta _s=9^\circ$ according to the droplet initial mass.