Sliding and merging of strongly sheared droplets

Abstract

A mathematical and numerical framework is proposed to compute the displacement and merging dynamics of sliding droplets under the action of a constant shear exerted by a gas flow. An augmented formulation is implemented to model surface tension including the full curvature of the free surface. A set of shallow-water evolution equations is obtained for the film thickness, the averaged velocity, an additional quantity (with dimension of a velocity) taking into account the capillary effects and a tensor called enstrophy. The enstrophy accounts for the deviation of the velocity profile from a constant velocity distribution. The formulation is consistent with the long-wave expansion of the basic equations with a conservative part and source terms including the effect of viscosity, in the form of a viscous friction and the effect of the shear stress. The model is hyperbolic with generalised diffusion terms due to capillarity. Finally, our model is completed with a disjoining pressure formulation that is able to account for the hysteresis of the static contact angle. In this formulation, the advancing or receding nature of the contact line is assessed by the accumulation or reduction of mass of the droplet at the contact line. Simulations of sliding water droplets are performed with periodic boundary conditions in a domain of limited size. Hysteresis of the static contact angle causes a slowdown of the drops and a delay in the sequence of coalescence of the drops.

1. Introduction

Transport of water droplets or films at the surface on a rigid wall under the action of a shear exerted by the surrounding atmosphere (shear flow) is drawing an increasing attention due to its importance in aeronautics in the context of icing phenomena.

Some experimental studies describe the phenomenon of shedding, which is the detachment of sessile droplet induced by the gas flow, and considered either the effect of the presence of other droplets (Razzaghi & Amirfazli 2019) or the influence of the surface wettability (Milne & Amirfazli 2009; Fan, Wilson & Kapur 2011). Fan et al. (2011) observed that the droplet may retain a footprint similar to that at the point of motion or exhibit a tail. In some cases, a trail remains behind the droplet (that can shed smaller droplets). In the context of aerodynamic applications, Moghtadernejad et al. (2016) further considered the coalescence of identical water droplets, on an aluminium plate. Their study reveals the formation of a rivulet as the result of the coalescence process. These authors further considered the evolution of rivulets on substrates of different wettabilities sheared by a high-velocity air stream (Moghtadernejad et al. 2014).

In the context of the wind-driven water run-back process on an airfoil, Zhang, Wei & Hu (2015); Zhang, Rothmayer & Hu (2016) considered the formation of water rivulets on a NACA0012 airfoil. The rivulets width and rivulets distributions along the wing span were found to depend on the airflow velocity and to result from the destabilisation of the liquid film which forms at the leading edge of the air foil. The instability mechanism is related to the dynamics of the advancing contact line during the water film run-back.

In most cases, numerical studies focus on the averaged properties of the film (liquid hold-up) and did not resolve the wavy nature of the liquid interface (see, for instance, Lan et al. 2008). However, some recent studies have been proposed based on crude low-dimensional modelling of the film flow yielding conservative, hyperbolic and two-dimensional equations (Gosset 2017; Lallement et al. 2018). Gosset (2017) used the numerical framework proposed by Meredith et al. (2011) which is based on a VOF formulation and the continuum surface force (CSF) model introduced by Brackbill, Kothe & Zemach (1992). Using Openfoam software, her simulations reproduced quite satisfactorily the onset of rivulets observed by Zhang et al. (2015, 2016). However, reasonable comparisons are impaired by the dependency of the results on the parameters of the partial-wetting model. In addition, the curvature of the interface is linearised using the long-wave approximation in this study.

Lallement et al. (2018) proposed instead to model surface tension with an augmented formulation. They introduced a transport equation for the gradient of the film thickness, ${\boldsymbol p} = \boldsymbol {{grad}}\, h$. This formulation enables one to lower the order of the derivatives in the averaged momentum balance from third order to second order. Such an augmented formulation has been initially proposed in Bresch et al. (2016) but was limited to the linearised long-wave approximation of the curvature of the free surface. They also introduced an innovative approach of partial wetting within the framework of the disjoining pressure model proposed by Derjaguin (1940). A literature review of the different formulations of disjoining pressure to model long-range surface forces within the framework of long-wave thin-film equations is available in Oron, Davis & Bankoff (1997). Lallement et al. thus introduced a disjoining force which accounts for partial wetting and enables to regularise the discontinuity of surface energy at the contact line at a scale $h_\star$, which corresponds physically to the range of microscopic forces, but is taken much larger at the order of the mesh size. The evolution equations derived by Lallement et al. is compatible with a conservation law of the total energy of the equation at the macroscopic level. However, a discretisation which preserves this property at the discrete level has not yet been achieved (Lallement 2019).

Within the framework of long-wave approximation, the Derjaguin disjoining pressure model has been used extensively to simulate the spreading and sliding of droplets (Bertozzi & Pugh 1994; Brandon, Wachs & Mamur 1997; Schwartz & Eley 1998; Ahmed et al. 2014; Espín & Kumar 2017) as well as the fingering instability of a liquid film front (Zhao & Marshall 2006). One advantage of the Derjaguin formulation is the introduction of a precursor film which alleviates the divergence of the viscous stresses at the contact line, also known as ‘the contact-line paradox’ (Dussan 1979). Precursor films are known to be observable for static and spreading non-volatile droplets and emerge from the absorption of the liquid at the solid surface by long-range forces such as van der Vaals forces. Typical measured thicknesses of precursor films are $0 \, (100\,{\unicode{x00C5}}$) (Popescu et al. 2012).

An alternative to this approach is the assumption of a slip at the wall boundary (see, e.g., Haley & Miksis 1991; Savva & Kalliadasis 2009). This approach is justified by molecular dynamics simulations (Ren & E 2007). The two approaches have been compared by Diez, Kondic & Bertozzi (2000) who reported that much lower values of the slip length than the precursor film thickness are required to capture satisfactorily the moving contact line characteristics, either speed or shape. In addition, the precursor film formulation does not require the contact line location to be tracked. The Derjaguin approach has been employed by Ahmed et al. (2014) to study the sliding of drops. Their formulation adjusted the Hamacker constant of the disjoining pressure at the front of the rear of the drops, to prescribe an advancing and a receding contact angles. Indeed, a hysteresis of the static contact angles at the advancing and receding fronts may be observed as a result of the surface inhomogeneities (Schwartz & Eley 1998; Zhao & Marshall 2006), either roughness (Savva & Kalliadasis 2009) or chemical heterogeneities of the surface (Brandon et al. 1997), or the presence of solutes (polymers or surfactants) in the liquid which may contaminate the surface and form a film (de Gennes 1985). Typical contact angle hysteresis is $10^\circ$ (de Gennes 1985) but may be higher. This phenomenon generates the adhesion force of a sessile droplet and determines the onset of shedding (Mahé et al. 1988; Milne & Amirfazli 2009; Moghtadernejad et al. 2014), as well as the speed and shape of sliding droplets on an inclined plane (Ahmed et al. 2014).

In this paper, we wish to formulate a set of averaged equations which models the evolution of a gas liquid film sheared by a gas flow within the framework of the shallow-water equations, which is well adapted for large-Reynolds-number flows. This approach enables us to depart from the usual creeping flow assumption used in lubrication theory, as in Espín & Kumar (2017) for instance, which is not valid for typical experimental conditions of shear-driven droplets.

As an example, Fan et al. (2011) conducted experiments of water, water–glycerin and glycerin droplets of typical thickness which can be estimated to be $h \approx 1$ mm, sheared by an air flow in a wind tunnel of effective parameter $L_{eff} = 0.032$ m. They reported typical gas velocities of $V_a \approx 10\,{\rm m}\,{\rm s}^{-1}$. With the help of the correlation

(1.1)\begin{equation} \frac{2\rho_a\tau_e}{V_a^2} = 0.0792 Re_G^{{-}1/4}, \end{equation}

where the gas Reynolds number defined as $Re_G = V_a L_{eff}/\mu _a$. This gives $\tau _e \approx 0.4$ Pa and consequently a liquid Reynolds number $Re \approx 400$ for water, in which case the creeping flow assumption is no more admissible.

In addition, our formulation is based on an augmented formulation proposed by Bresch et al. (2020) which enables one to avoid the long-wave approximation of the free-surface curvature, a questionable assumption for large contact angles.

Our formulation alleviates the limitations of previous attempts by conserving the surface and disjoining energy and preserving the consistency with the long-wave asymptotics. This proposed derivation will account for the displacement of contact lines with the Derjaguin formula. The evolution of droplets under the action of a flow is investigated. In particular, we are interested in the effect of the contact angle hysteresis on the droplets speed and shape. In this work, we limit ourselves to consider only a constant external shear stress. For a boundless gas flow, this assumption is admissible as long as the wake of the droplet does not generate a boundary-layer separation in the gas flow. Indeed, the shear stress exerted by the gas flow is determined by the boundary layer at the gas–liquid interface. Since, boundary-layer equations obey the Prandtl transformation, the response of the shear stress to the wall geometry is negligible as long as the potential flow far from the wall remains unaffected. For open flows, wall roughness thus only affect the wall shear stress through second-order modifications of the pressure in the outer potential region (Luchini & Charru 2019).

Section 2 is devoted to the mathematical formulation that is able to describe the sliding of the drops under a constant shear. This mathematical framework is based on the Saint-Venant approach, i.e. an in-depth averaging of the basic equations which is made possible by the long-wave nature of the flow as the typical thickness of the water layer is much smaller than the extension of the drop on the wall. Section 3 presents a numerical investigation of the characteristics of single droplets and the coalescence dynamics of a cloud of small droplets in a periodic domain. Concluding remarks are given in § 4.

2. Mathematical developments

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 (2007) and for the tensor $\boldsymbol \varPhi$ in Richard, Duran & Fabrèges (2019a), 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. 2019b). 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 (2007).

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. 1997)

(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 (2017) 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. 2019).

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 1995), 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. 2015). Yet Diez et al. (2000) 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. (2000), 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. (2020) 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 (2014) and Bresch et al. (2016) 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. (2020), 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. Simulations

An in-house two-dimensional numerical solver implementing the proposed model (2.64) has been developed. The code is written in Julia, and uses the method of line. This method adopts a semidiscrete form of the model: an ordinary differential equation (ODE) system is obtained after the spatial discretisation is implemented. The spatial discretisation of the model uses both the finite-volume method for the hyperbolic terms at the left-hand side of the equations and the finite-difference method for the remaining terms on the right-hand side. The MUSCL method is employed along with a slope limiter (optional), ensuring a numerical scheme of second order in space for the hyperbolic part of the equations. The ‘minmod’, ‘superbee’ and ‘Van-Leer’ limiters have been implemented. Their objective is to ensure the connection between the zones of strong and weak gradients as finite-volume schemes of order greater than one are known to be unstable at strong gradients (Godunov & Bohachevsky 1959). The choice of limiters enables to switch from a scheme of second order in the regions where the solution is smooth to a scheme of first order where the gradients are important to ensure total variation diminishing property of the global scheme. The choice of the limiter has a major impact on the numerical diffusion: this diffusion artificially smoothes the solution, which is to be avoided, but makes the diagram more robust to sudden oscillations near the discontinuities. In implemented limiters, ‘minmod’ is the most diffusive (and stable) and ‘superbee’ the most precise (but can generate digital oscillations). ‘Van-Leer’ is a compromise between these two extremes. The remaining terms are discretised by centered finite differences of second order. Integration in time of the obtained system of ODE requires the selection of a time solver which satisfies the following constraints.

1. (i) The presence of a steep fronts at the contact lines implies a ‘stiff’ problem, which requires the use of an implicit solver. This means that a large-dimensional linear system $(N_x \times N_y \times 8) ^ 2$ must be solved.

2. (ii) Solving this linear system requires the computation of the associated Jacobian. Its finite-difference approximation is expensive and imprecise. This requires us to make use of sparse matrices in order to limit its memory storage.

In comparison to the scheme developed by Bresch et al. (2020), we have avoided a semi-implicit scheme for the surface tension terms. Instead of using an implicit solver, we have chosen to use an explicit strong stability-preserving Runge–Kutta (SSPRK) method (Gottlieb, Shu & Tadmor 2001). This method is a Runge–Kutta method that preserves the stability of first-order methods. They have proved helpful in solving hyperbolic partial differential equations. Being explicit, this method does not require the computation of the Jacobian, neither the solution of a linear system. This leads to a significant reduction of the computational time as well as easy-to-write parallelisation of the code. A CUDA compatible version of the code has been developed, which allows to perform simulations on a GPU. The SSPRK method enables to control the time step according to a local truncation error, which can be used to control the accuracy of the solution. The time step is controlled by a proportional–integral algorithm that adjusts the time step according to the local truncation error. Absolute local threshold is set to $10^{-6}$ and relative local threshold is set to $10^{-3}$.

A reprojection routine is used to ensure that the augmented variable $\boldsymbol {W}$ stay close to its definition (see (2.62a)): this is mandatory when the variable become discontinuous, which will be the case near the triple point. The reprojection is done at each timestep.

Simulations have been carried out on an AMD Ryzen 9 3950X 16-core processor and a GeForce RTX 3080 Ti (10’240 CUDA core). The performance of the numerical code has been evaluated for a benchmark consisting of the simulation of a two-dimensional droplet for the flow conditions corresponding to figure 10 with no hysteresis and $N\times 3N=2 \times 10^6$ nodes, aborting the simulations at time $t=20$. Single-thread, mutlithread and GPU runs have been compared. A single-thread run lasts $8018$ s, an execution time, which is lowered to $2771$ for 4 threads, but only reduces to $2001$ s for 8 threads and actually increase to $2050$ s for 16 threads. In contrast, the GPU simulation lasted only $189$ seconds. Tests of convergence show a quadratic dependence of the time step $\delta t$ on the mesh size $\varDelta$ in the form $\delta t = O( \varDelta ^2)$ as expected from the second-order generalised-diffusion nature of the capillary terms.

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 2008). For instance, using a precursor film formulation and a Frumkin–Derjaguin disjoining pressure, Ahmed et al. (2014) 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 (1998) and Sellier et al. (2015), 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. 2014).

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 (2008) 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.

4. Conclusion and perspectives

In this work, a mathematical and numerical framework has been developed to study the displacement and merging dynamics of sliding droplets under the action of a constant shear exerted by a gas flow. The mathematical development has been developed in § 2. An augmented formulation has been implemented to model surface tension including correctly the full curvature of the free surface. Details can be found in Bresch et al. (2020). A set of shallow-water evolution equations for the film thickness $h$, the averaged velocity ${\boldsymbol U}$, the additional velocity ${\boldsymbol W}$ and an enstrophy tensor $\varPhi$, have been obtained. The enstrophy $\varPhi$ accounts for the deviation of the velocity profile from the constant velocity distribution that is usually assumed in the Saint-Venant equations and is not verified for viscous flows at low Reynolds numbers. The formulation is consistent with the long-wave expansion of the basic equations, remains hyperbolic and conservative. Finally, our model has been completed with a disjoining pressure formulation that is able to account for the hysteresis of the static contact angle. In this formulation, the advancing or receding nature of the contact line is assessed by the accumulation or reduction of mass of the droplet at the contact line (sign of $\textrm {div}(h{\boldsymbol U})$). Section 3 is devoted to an in-house two-dimensional numerical solver implementing the proposed model. This code is written in Julia and uses the method of line. We discuss the contact angle hysteresis and provide simulation of single droplet as well as dynamics of coarsening. The derivation of the enstrophy equation and a discussion on drop speed is given in Appendix A.

Simulations of sliding water droplets have been performed with periodic boundary conditions in a domain of limited size. Our simulations show a transition of the shape of the sliding droplets from a tear-drop shape to a bullet shape with a round back and an ogival front as the hysteresis of the contact angle is raised. Simulations of the coarsening dynamics of the drop demonstrate a slowdown of the drops and a delay in the sequence of coalescence of the drops due to the contact angle hysteresis.

We have limited ourselves to consider only a constant shear stress at the free surface. However, several studies suggest a retroaction of the droplet geometry on the gas flow (Pozrikidis 1997; Ding & Spelt 2008; Sellier et al. 2019). Razzaghi & Amirfazli (2019) thus showed that the shedding of a sessile droplet by an airflow is affected by the presence of other droplets nearby. Hooshanginejad & Lee (2017) revealed that this interaction results from the modification of the pressure field induced by the wake of the first droplet on the second one. It is therefore important to extend the present liquid-side formulation to a two-phase flow one for which the gas motion is resolved (Lavalle et al. 2017). This will be investigated in a subsequent work.