Droplet settling on solids coated with a soft layer

Abstract

Gravitational settling of a droplet in air onto a soft substrate is a ubiquitous event relevant to many natural processes and engineering applications. We study this phenomenon by developing a three-phase lubrication model of droplet settling onto a solid substrate coated by a thin soft layer represented by a viscous film, an elastic compressible layer and an elastic sheet supported by a viscous film. By combining scaling analysis, analytical methods and numerical simulations we elucidate how the resulting droplet dynamics is affected by the nature of the soft layer. We show that these soft layers can significantly affect the droplet shape during gravitational settling. When there is a linear response of the deformations of the soft layer, the air layer takes longer to drain as compared with the case of a droplet settling onto a rigid substrate. Our results provide new insight into the coupled interactions between droplets and solids coated by a thin film of a soft material.

1. Introduction

Impacts of drops onto solid substrates or liquids are ubiquitous in many natural and industrial processes. They include inkjet printing, spray coating, forensic analysis, air–sea transfer and epidemiology of foliar diseases to name but a few examples. The broad relevance of droplet impact has made it a widely studied topic, for which the effects of the interfacial and bulk properties of the drop as well as the substrate properties have been characterized (Rein 1993; Neitzel & Dell'Aversana 2002; Yarin 2006; Kavehpour 2015; Josserand & Thoroddsen 2016; Ajaev & Kabov 2021). The drop interface dynamics has been elucidated through asymptotic analysis, scaling laws and simulations when the impact speed is very slow (Yiantsios & Davis 1990; Duchemin & Josserand 2020), a situation we refer to as settling. When the drop approaches the solid it must drain the air layer separating the two interfaces, which leads to a build-up of pressure in the air film. As a consequence, the drop interface deforms and takes the shape of a dimple before making direct contact with the solid. The dimple-shaped interface has a maximum thickness at the axis of symmetry of the drop and a minimum near its outer edge where a neck radius can be defined. The formation of a dimple-shaped interface is not only observed for settling droplets but also when a bubble slowly approaches a rigid surface (Chan, Klaseboer & Manica 2011), as well as for inertial drop impacts where it significantly affects the dynamics (Thoroddsen et al. 2005; Xu, Zhang & Nagel 2005; Mani, Mandre & Brenner 2010; Hendrix et al. 2016).

Drop impacts onto compliant soft substrates (see figure 1) are also common. This soft substrate can for instance be a viscous liquid film, a soft elastic layer or an elastic sheet supported by a viscous film. These are the three cases we describe herein. Over recent years there has been an emerging interest in how soft materials influence capillary flows (Bico, Reyssat & Roman 2018; Andreotti & Snoeijer 2020). Problems involving elastohydrodynamic lubrication, also known as soft lubrication (Skotheim & Mahadevan 2005; Essink et al. 2021), are also critical for a wide variety of systems ranging from stereolithography to biological adhesion (Wang et al. 2017; Chan & Carlson 2019; Wang, Feng & Frechette 2020). In this context, experiments have demonstrated that the dynamics of drop impacts can be controlled by the softness of the solid (Pepper, Courbin & Stone 2008; Chen & Li 2010; Chen et al. 2016; Howland et al. 2016; Langley, Castrejón-Pita & Thoroddsen 2020). The deformations of the substrate can affect the dynamics after contact, either by absorbing some of its energy or through the contact line motion (Andreotti & Snoeijer 2020; Dervaux, Roché & Limat 2020). The high lubrication pressure can also deform the substrate before contact occurs, as recently observed in experiments (Langley et al. 2020) as well as in numerical and theoretical work (Pegg, Purvis & Korobkin 2018; Henman, Smith & Tiwari 2021). How the compliance of the substrate affects the air drainage and dimple formation during settling of droplets has so far been overlooked, but could produce significant effects for the resulting interfacial flow. Recent experimental advances (Lo, Liu & Xu 2017; Pack et al. 2017; Lakshman et al. 2021; Zhang et al. 2021) now also allow us to probe the influence of these surface deformations for drop impacts on thin liquid films down to the nanoscale.

Figure 1. (a) Schematic of a droplet of radius $a^\star$, density $\rho _1$, viscosity $\mu_1$, surface tension coefficient $\sigma _1$, suspended in a fluid at rest (typically air) of density $\rho _2$, viscosity $\mu_2$ and settling on a soft material under the influence of the gravity field $-g\boldsymbol {e}_z$. Deformations are not to scale. The inset shows the definition of the droplet profile $h_1^\star (r^\star,t^\star )$, the profile of the soft material $h_2^\star (r^\star,t^\star )$ and the thickness of the air layer $H^\star (r^\star,t^\star )=h_1^\star (r^\star,t^\star )-h_2^\star (r^\star,t^\star )$ between the two. We consider three different soft substrates: (b) a compressible Hookean solid characterized by its Lamé coefficients $G$ and $\lambda$; (c) a thin viscous film with surface tension coefficient $\sigma _3$; (d) an elastic sheet with thickness $d^\star$, Young's modulus $E$, and Poisson's ratio $\nu$, giving a bending stiffness $B=Ed^{\star 3}/12(1-\nu ^2)$. In all cases the height of the undeformed soft layer is $h_s^\star$, and for (c,d) the liquid film has a viscosity $\mu_3$. The dimensional heights denoted by stars ()$^\star$ are non-dimensionalized with the initial air layer thickness at $r^\star =0, H_0^\star$, whilst the radial coordinate $r^\star$ is non-dimensionalized with $(H_0^\star a^\star )^{1/2}$. The dimensional time $t^\star$ is non-dimensionalized with $t_0^\star =\mu_2/{\rm \Delta} \rho g a^\star$.

To describe the settling of a droplet onto a soft surface, represented by a thin compliant layer, we study a minimal model considering a very viscous flow (very small Reynolds numbers) and droplets small enough for the interface deformations due to the flow to be localized near the substrate and to not affect the overall drop shape; this requires a small Bond number, i.e. capillary effects to dominate over gravity. These assumptions allow us to describe the settling dynamics based on the lubrication approximation, and we also consider a regime where deformations of the drop and the interface only appear once the lubrication assumptions hold. Our analysis builds on the previous works from Yiantsios & Davis (1990) and Duchemin & Josserand (2020). Yiantsios & Davis (1990) described the settling of a droplet onto a bath of liquid and derived the long-term asymptotics of the quantities defining the dimple. They also considered the effects of slip at the droplet interface by coupling the lubrication equations with boundary integral equations to account for the flow inside the droplet. More recently, Duchemin & Josserand (2020) presented numerical simulations and scaling analysis using the lubrication approximation to rationalize the settling of a large drop onto a thin liquid film, showing how slip at the film interface can accelerate the settling process. However, the effects of the deformations of a soft substrate on the settling dynamics of droplets has so far been overlooked. In this article we consider how the settling dynamics is affected by the deformations of either a compressible elastic layer, a thin viscous liquid film or an elastic sheet supported by a viscous film, as outlined in figure 1.

2. Problem set-up and droplet settling onto a rigid substrate

2.1. Lubrication equations

Figure 1(a) illustrates the settling of a droplet towards a solid substrate coated with a thin soft layer. The initially spherical droplet of radius $a^\star$, viscosity $\mu_1$, density $\rho _1$, is suspended in a fluid (typically air) with density $\rho _2<\rho _1$, viscosity $\mu_2$ and surface tension coefficient $\sigma _1$. The droplet settles towards the soft surface by gravity, characterized by the gravitational acceleration $g$. To describe this interface dynamics, we assume an axisymmetric flow and start with the lubrication theory derived by Yiantsios & Davis (1990) for a rigid substrate, before moving on to a description of how the dynamics is altered when the substrate is coated with a soft layer.

The governing equations describing the flow in the air layer are made dimensionless; we denote dimensional lengths, times, pressures and velocities with a star $(\ )^\star$. Since the droplet is initially spherical and the substrate is undeformed, the air layer thickness is well approximated by a parabolic profile $H^\star (r^\star,0)=H_0^\star + r^{\star 2}/2a^\star$ near the axis of symmetry ($r^\star =0$) at time $t^\star =0$. Therefore the characteristic vertical length scale is the initial thickness of the film at $r^\star =0, H_0^\star$, whilst the characteristic radial length scale is $r_0^\star =(H_0^\star a^\star )^{1/2}$. The droplet motion is driven by its weight, giving a characteristic pressure in the air film $p_0^\star ={\rm \Delta} \rho g a^{\star 2}/H_0^\star$ where ${\rm \Delta} \rho =\rho _1-\rho _2$. The time scale of the process is $t_0^\star =\mu_2/{\rm \Delta} \rho ga^\star$. We assume a small aspect ratio, $\varepsilon = H_0^\star /r_0^\star =(H_0^\star /a^\star )^{1/2}\ll 1$, as well as small Reynolds numbers both in the droplet and in the air film. This allows us to take advantage of the lubrication approximation, i.e. we neglect any inertial effects and consider the equations at leading order in $\varepsilon$. There are also additional assumptions we make in order to simplify the interfacial boundary conditions. First, we only consider small interfacial slopes. Second, we consider the no-slip condition for the air film at the substrate and at the droplet interface. The vertical length scale in the air layer is $H_0^\star$ whilst it is the same as the horizontal length scale, $r_0^\star$, in the droplet. If $U_2^\star$ is the scale of the radial velocity in the air layer, the continuity of the shear stress at the interface shows that the radial velocity in the droplet scales as $U_1^\star = \varepsilon \mu_2 U_2^\star /\mu_1$. The continuity of radial velocity at the interface gives, in dimensionless form, $U_1^\star u_1 = U_2^\star u_2$, which simplifies to $\mu_2 u_2=\mu_1 u_1 / \varepsilon$. When the droplet is very viscous, i.e. when $\mu_2/\mu_1 \ll \varepsilon$, we obtain at the interface $u_2\rightarrow 0$. This approximation becomes inaccurate when the droplet is very close to the substrate, but since slip at one of the interfaces only changes the prefactor of the governing equations, the dynamics is expected to be similar regardless. Finally, we assume that the Bond number is small enough, namely $\delta ={Bo}/\varepsilon ^2 \ll 1$ with ${Bo}={\rm \Delta} \rho g a^{\star 2}/\sigma _1$, to considerably simplify the normal stress balance. Considering a typical fluid with a capillary length of approximately 3 mm, and setting $\varepsilon =0.1$, this last condition amounts to considering a droplet of radius $a^\star \ll 300\ {\rm \mu}{\rm m}$.

By combining these assumptions we obtain a set of three equations. The link between the pressure $p_2^\star (r^\star,t^\star )$ and the thickness $H^\star (r^\star,t^\star )$ of the air layer is given by Reynolds’ thin film equation, derived from the Navier–Stokes equations by considering the classical assumptions of lubrication theory (Batchelor 1967), $\partial H^\star /\partial t^\star = (1/12\mu ) \nabla ^\star (H^{\star 3} \nabla ^\star p_2^\star )$, where $\nabla ^\star$ is the gradient operator. The linearized normal stress balance simplifies to the Young–Laplace equation and links the pressure $p_2^{\star }(r^\star,t^\star )$, measured relative to the undeformed drop, to the drop height $h_1^\star (r^\star,t^\star )$: $p_2^{\star }=\sigma _1 (2/a^\star - \nabla ^2 h_1^{\star } )$. Finally, a constraint on the pressure field is given by a quasi-steady force balance on the droplet obtained by neglecting its inertia: the force from the pressure field balances the weight of the drop, which gives the integral condition $2{\rm \pi} \int _{0}^{+\infty }p_2^\star r^\star \, {\rm d}r^\star =4{\rm \pi} a^{\star 3} g{\rm \Delta} \rho /3.$ In dimensionless form, this leads to the following set of equations (Yiantsios & Davis 1990):

(2.1a)$$\begin{gather} {\frac{\partial H}{\partial t}}(r,t)=\frac{1}{12}\frac{1}{r} {\frac{\partial }{\partial r}}\left(rH^3(r,t){\frac{\partial p_2}{\partial r}}(r,t)\right), \end{gather}$$

(2.1b)$$\begin{gather}\delta p_2(r,t) = 2-\frac1r{\frac{\partial }{\partial r}}\left(r{\frac{\partial h_1}{\partial r}}(r,t)\right), \end{gather}$$

(2.1c)$$\begin{gather}\int_0^{+\infty} p_2(r,t) r\, {\rm d}r = \frac{2}{3}, \end{gather}$$

where $H(r,t)=h_1(r,t)-h_2(r,t)$ is the air film thickness and $h_1(r,t)$ represents the droplet interface (figure 1a). As mentioned above, the force balance (2.1c) neglects the inertia of the droplet. This model is valid for small approach speeds $V^\star$, when $V^\star < H_0^\star /t_0^\star$. This translates to a condition on the Capillary number ${Ca}=\mu_2 V^\star /\sigma _1$: ${Ca}\ll \delta \varepsilon ^2$, which also ensures that the droplet is initially spherical.

The system (2.1) must be coupled with a governing equation for the height of the soft substrate $h_2(r,t)$. In § 2.3 we consider a rigid surface with $h_2(r,t)=0$; in § 3 we consider a compressible elastic layer with (3.1); in § 4 we consider a thin liquid film with (4.2) and (4.3); and in § 5 we consider an elastic sheet supported by a thin viscous film with (5.1). We recall that the heights $H(r,t), h_1(r,t), h_2(r,t)$ are scaled with $H_0^\star$, the radial coordinate $r$ with $r_0^\star =(H_0^\star a^\star )^{1/2}$, the time $t$ with $t_0^\star =\mu_2/{\rm \Delta} \rho g a^\star$ and the pressure $p_2(r,t)$ with $p_0^\star = {\rm \Delta} \rho g a^{\star 2}/H_0^\star$. Table 1 summarizes the definitions of the dimensionless numbers characterizing the settling dynamics, along with those relevant for describing settling onto soft layers.

Table 1. Dimensionless numbers characterizing the settling of a droplet on soft layers. We consider a compressible elastic layer, a viscous liquid film (viscous film, capillary) and an elastic sheet supported by a viscous film (viscous film, elastic). The last column indicates the range of values used in our numerical simulations.

2.2. Numerical procedure

The systems of equations we study, (2.1) supplemented with equations for $h_2(r,t)$, are solved numerically using the finite element method implemented in the FEniCS code (Alnæs et al. 2015) with quartic polynomial elements and an implicit time integration procedure. For all cases studied, the initial condition for the air layer thickness is fixed, presenting an initially small deviation to the spherical droplet following Yiantsios & Davis (1990): $H(r,0)=1+r^2/2 + \delta (5/18 - \ln (1+r^2/2)/3)$. The integral condition (2.1c) is implemented as a boundary condition using (2.1b): $(\partial h_1/\partial r) (r,t) \overset {r\rightarrow \infty }{\sim } r- 2\delta /3r$.

2.3. Droplet settling onto a rigid substrate

For a droplet settling onto a rigid substrate the system (2.1) is closed since $h_2(r,t)$ is a constant, taken here as $0$ without loss of generality, and $h_1(r,t)=H(r,t)$. A remarkable feature of the air film dynamics is that the droplet interface evolves into a dimple, a small region centred around the axis of symmetry where its deformations are localized. Assuming that the droplet becomes almost flat in this dimple region, the pressure there is then almost uniform and given by the Laplace pressure of the droplet: $p_2(r< r_d,t) \simeq 2/\delta$. The force balance (2.1c) then gives the dimple radius as $r_d\simeq (2\delta /3)^{1/2}$ (Derjaguin & Kussakov 1939; Frankel & Mysels 1962). This radius corresponds to the neck of the dimple where the air layer thickness $H(r,t)$ is minimum. Outside the dimple, for $r>r_d$, the droplet is nearly undeformed and spherical. The dimple geometry is illustrated by numerical simulations of (2.1) in figure 2(a).

Figure 2. Results from a simulation of (2.1) for droplet settling on a rigid substrate, with $\delta =0.05$. (a) Droplet profile at $t=50$, 150, 400, 1000 and 20 000 in the dimple. The expected dimple radius is $r_d=(2\delta /3)^{1/2}\simeq 0.18$. (b) Evolution of the height at the axis of symmetry, $H(0,t)$, and of the minimum height, $H_{min}(t)$. (c) Profiles of the pressure in the air layer at $t=10$, 100, 1000 and 10000. The dotted profile represents a uniform pressure $2/\delta$ over a region of size $r_d$. The inset shows the time evolution of the pressure at $r=0$, with the dotted lined representing the asymptotic value $2/\delta$.

Further insight into the interfacial dynamics may be gained through scaling analysis (Duchemin & Josserand 2020). Following the above discussion we assume a dimple profile to be formed with $r_d\sim \delta ^{1/2}$ and a uniform pressure $p_2(r< r_d,t)\sim \delta ^{-1}$. Integrating (2.1a) then gives a mass balance as

(2.2)\begin{equation} r_d^2 {\frac{\partial H}{\partial t}}(0,t) \sim r_d H_{min}^3(t) {\frac{\partial p_2}{\partial r}}(r_d,t),\end{equation}

where the gradient of pressure $(\partial p_2/\partial r)(r_d,t) \sim p_2/\ell (t) \sim 1/\delta \ell (t)$ is localized within the radial extent of the dimple neck $\ell (t)$, and the minimum thickness is $H_{min}(t)=H(r_d,t)$. Matching the curvature of the neck with the curvature of the undeformed droplet yields $H_{min}(t)/\ell (t)^2 \sim 1$, whilst the slope is matched between the dimple and the neck: $H_{min}(t)/\ell (t) \sim H(0,t)/r_d$. These two matching conditions simplify (2.2) to $({\rm d}H/{\rm d}t)(0,t) \sim \delta ^{-4}H(0,t)^5$, which finally allows us to find the following scaling laws: $H(0,t) \sim \delta t^{-1/4}, H_{min}(t) \sim \delta t^{-1/2}, \ell (t) \sim \delta ^{1/2} t^{-1/4}$. A more rigorous mathematical analysis has also been conducted by Yiantsios & Davis (1990), who derived and validated numerically the following long-time behaviours from an asymptotic expansion of the governing equations (2.1)

(2.3a)$$\begin{gather} H(0,t)=(0.3273)\delta t^{{-}1/4} , \end{gather}$$

(2.3b)$$\begin{gather}H_{min}(t)=(0.4897)\delta t^{{-}1/2}, \end{gather}$$

where $r_d$ approaches a constant value as $(2\delta /3)^{1/2} (1-0.1652 t^{-1/4})$. These asymptotic results are in very close agreement with the numerical simulations of (2.1) as shown in figure 2. In the next sections we describe how this physical picture changes when considering the settling of a droplet onto soft surfaces.

3. Solid substrate coated with a thin compressible elastic layer

We move on to describe the case when the rigid substrate is coated with a thin and compressible elastic layer (figure 1$b$). At rest, the layer has uniform thickness $h_s^\star$, such that $h_2(r,0)=h_2(r\rightarrow \infty,t)=h_s^\star /H_0^\star =h_s$. The response of the layer is assumed to follow Winkler's model (Dillard et al. 2018; Chandler & Vella 2020), i.e. $h_2^\star (r^\star,t^\star )=h_s^\star -\eta ^\star p_2^\star (r^\star,t^\star )$, with $\eta ^\star$ a proportionality constant relating the displacement to the external pressure acting on the layer. Winkler's model can be formally derived for compressible Hookean layers using lubrication theory, in which case $\eta ^\star =h_s^\star /(2G+\lambda )$ with $G$ and $\lambda$ the Lamé coefficients of the elastic material, assumed to be of the same order of magnitude. Such a simple relationship is obtained in the limit $\varepsilon \ll 1, \varepsilon h_s \ll 1$, as shear stresses are negligible and only the pressure contributes to deformations (Skotheim & Mahadevan 2005; Chandler & Vella 2020). Winkler's model can also be used to describe thin liquid-infused poroelastic layers (Skotheim & Mahadevan 2005) as well as soft polymer brushes (Gopinath & Mahadevan 2011; Davies et al. 2018). In dimensionless form, the relation between the pressure and displacement becomes

(3.1a)$$\begin{gather} h_2(r,t)=h_s-\eta p_2(r,t), \end{gather}$$

(3.1b)$$\begin{gather}\eta= \frac{{\rm \Delta} \rho g a^{\star2} h_s^\star}{H_0^{\star2} (2G+\lambda)}. \end{gather}$$

The dimensionless number $\eta$ is a softness parameter measuring the relative importance of the pressure exerted by the drop on the compressible layer compared with its Lamé coefficients, and which controls the importance of the layer deformations compared with its thickness (Skotheim & Mahadevan 2005). As $\eta \rightarrow 0$, i.e. as the layer becomes very stiff, we recover $h_2(r,t)=h_s$, a constant, and the results for a rigid solid presented in § 2.3 apply.

To illustrate the effect of a soft layer on the droplet settling dynamics, we perform numerical simulations of (2.1) and (3.1). When the layer is stiff compared with the droplet, i.e. when $\eta$ is sufficiently small, we expect the droplet settling dynamics to approach the analytical solution for a rigid wall (2.3). Figure 3 also shows that when $\eta$ is large, the droplet profile $h_1(r,t)$ becomes nearly undeformed while the elastic layer $h_2(r,t)$ also adopts the parabolic shape of the droplet. Accordingly, we start the analysis in the limit of a very soft compressible layer by assuming the following ansatz:

(3.2)\begin{equation} h_s-h_2(r,t)=\frac12(r_e^2-r^2),\end{equation}

valid up to a radius $r_e$, and an undeformed elastic layer $h_2(r,t) = h_s$ for $r>r_e$. From (3.1), the corresponding pressure at the axis of symmetry is $\eta p_2(0,t) = r_e^2/2$. The force balance (2.1c), limiting the integration to $0\leq r\leq r_e$, gives: $r_e^4-4\eta p_2(0,t) r_e^2 + 16\eta /3=0$. These two relations yield

(3.3a)$$\begin{gather} r_e=\left(\frac{16}{3}\eta\right)^{1/4}, \end{gather}$$

(3.3b)$$\begin{gather}p_2(r=0,t)=\left(\frac{3}{4}\eta\right)^{{-}1/2}. \end{gather}$$

To arrive at these results we have assumed that the droplet maintains its spherical shape and that the elastic layer follows a similar profile with the same curvature as the droplet. This means that the thickness of the air film is uniform: $H(r,t)\simeq H(t)$ for $r< r_e$. The evolution equation (2.1a), with the parabolic pressure profile given by (3.1), (3.2) and (3.3) as $p_2(r,t)=2/\sqrt {3\eta }-r^2/2\eta$, then yields for $t \gg \eta ^{-1}$ and $r< r_e$

(3.4)\begin{equation} H(r,t)=\left(\frac{t}{3\eta}\right)^{{-}1/2}.\end{equation}

Figure 3. Profiles of the droplet, $h_1(r,t)$, and of the compressible elastic layer $h_2(r,t)$, at $t=1000$ and with $\delta =0.05$ and (a) $\eta =10^{-5}$, (b) $\eta =10^{-3}$. The initial height of the elastic layer, $h_s$, is only significant in the definition of $\eta$ (3.1b) and $h_2(r,t)$ is translated in the figures so that it is 0 in its undeformed state.

To find the range of applicability of the above results we can integrate (2.1b), which we have not yet used, with the pressure profile derived herein

(3.5)\begin{equation} h_1(r,t)=h_1(0,t)+\frac{r^2}{2}\left(1-\frac{1}{\sqrt{3}} \frac{\delta}{\sqrt{\eta}}\right)+\frac{1}{32}\frac{\delta}{\eta}r^4. \end{equation}

The analysis is consistent only if the corrections to the ansatz (3.2) are small, i.e. when $h_1(r,t)\simeq h_1(0,t) + r^2/2$, which requires $\eta \gg \delta ^2$. Cancelling the quartic term up to $r=r_e=(16\eta /3)^{1/4}$ reduces to $\delta \ll 1$, consistent with our assumptions discussed in § 2.1.

In order to verify these predictions, we show in figure 4 the time evolution of the air film thickness at $r=0$, of the minimum air film thickness, of the pressure in the air film and of the radial position of the neck. These numerical results confirm the asymptotic behaviour derived herein as well as their application range, $\eta \gg \delta ^2$, while the effects of elasticity are negligible when $\eta \ll \delta ^2$. In particular, the minimum height of the air film scales as $H_{min}(t)\sim t^{1/2}$ with a prefactor always larger than for a rigid wall. This suggests that a soft enough material would delay direct contact between the droplet and the solid.

Figure 4. Air layer dynamics during a droplet settling onto a thin compressible elastic layer. (a,b) Time evolution of the air layer thickness $H$ at $r=0$ (a) and at its minimum value (b), for $\delta =0.05$ and $\eta =10^{-7}, 10^{-5}, 10^{-3}, 10^{-1}$. (c) Radius of the neck, $r_{neck}(t)={\rm argmin}_r[H(r,t)]$ at $t=10^3$. (d) Pressure profile in the air layer at $t=10^3$. Solid lines represent the numerical results, the dotted line are the profiles $p_2(r,t)=2/\sqrt {3\eta }-r^2/2\eta$ derived in the text assuming the dynamics is dominated by the elasticity of the compressible layer. The inset shows the pressure at $r=0$ and at the same time for various values of $\eta$.

To place these results in the context of relevant material parameters, we consider a rescaled Bond number $\delta ={O}(10^{-1})$, an aspect ratio $\varepsilon ={O}(10^{-1})$, a liquid–gas density difference ${\rm \Delta} \rho ={O}(10^3\ {\rm kg\ m}^{-3})$, a surface tension coefficient $\sigma _1={O}(50\ {\rm mN\ m}^{-1})$ and a compressible layer thickness $h_s^\star \sim H_0^\star$. The condition $\eta \gg \delta ^2$ translates into the following condition for the material properties of the soft layer for compressibility effects to be dominant: $(2G+\lambda ) \ll 10^4\ {\rm Pa}$. Consequently, the analysis is primarily reserved for very soft materials, for instance substrates grafted by a layer of polymer brushes which can verify this criterion (Davies et al. 2018) and for which the Winkler mode applies (Gopinath & Mahadevan 2011).

4. Solid substrate coated with a thin viscous liquid film

Next we consider a droplet settling on a rigid substrate coated with a thin viscous liquid film (figure 1$c$) of viscosity $\mu_3=\mu_2/\lambda$ and surface tension coefficient $\sigma _3$. We continue to consider no slip of air at the droplet interface, as discussed in § 2. At the interface between the air layer and the liquid film we consider the continuity of radial velocity, $u_2^\star =u_3^\star,$ as well as the shear stress balance, $\mu_2 \partial u^\star _2/\partial z^\star = \mu_3 \partial u^\star _3/\partial z^\star$, where $u_2^\star$ and $u_3^\star$ are the radial velocities in the air and liquid layer, respectively. There is also no slip at the base of the liquid film. With these boundary conditions and using the lubrication approximations, the Navier–Stokes and mass conservation equations written in both layers yield the following coupled equations for the air film thickness $H(r,t)$ and the liquid film height $h_2(r,t)$

(4.1a)\begin{align} {\frac{\partial H}{\partial t}}(r,t)&= \frac{1}{12}\frac{1}{r}{\frac{\partial }{\partial r}} \left[\frac{r}{H(r,t)+\lambda h_2(r,t)}\left( H^4(r,t){\frac{\partial p_2}{\partial r}}(r,t) \right.\right. \nonumber\\ &\quad + \left.\left.4\lambda H^3(r,t) h_2(r,t) {\frac{\partial p_2}{\partial r}}(r,t) + 3\lambda H^ 2(r,t)h_2^2(r,t) {\frac{\partial p_3}{\partial r}}(r,t) \right)\right], \end{align}

(4.1b)\begin{align} {\frac{\partial h_2}{\partial t}}(r,t)&= \frac{\lambda}{3}\frac{1}{r}{\frac{\partial }{\partial r}}\left[\frac{r}{H(r,t)+\lambda h_2(r,t)}\left(H(r,t)h_2^3(r,t) {\frac{\partial p_3}{\partial r}}(r,t) \right.\right. \nonumber\\ &\quad + \left.\left.\frac\lambda4 h_2^4(r,t) {\frac{\partial p_3}{\partial r}}(r,t) + \frac34H^2(r,t)h_2^2(r,t){\frac{\partial p_2}{\partial r}}(r,t)\right)\right]. \end{align}

Equation (4.1a) replaces (2.1a), which was derived assuming no slip of air at the interface with the soft substrate, whilst the force balance (2.1c) and normal stress balance at the air–drop interface (2.1b) remain unaltered. The pressure $p_3(r,t)$ in the viscous film is given by the normal stress balance at the interface between the air layer and the viscous film

(4.2)\begin{equation} \xi (p_3(r,t)-p_2(r,t))={-}\frac1r {\frac{\partial }{\partial r}}\left(r{\frac{\partial h_2}{\partial r}}(r,t)\right),\end{equation}

where $\xi = \varepsilon ^{-2} {\rm \Delta} \rho ga^{\star 2}/\sigma _3$ is a rescaled Bond number defined similarly to $\delta, \xi =\delta \sigma _1/\sigma _3$. At $t=0$ and far from the droplet, the height of the layer is $h_2(r,0)=h_2(r\rightarrow \infty,t)=h_s$. For the most common liquids, e.g. aqueous substances and common oils, we expect that surface tension coefficients remain in the same range, and hence that $\xi$ and $\delta$ have the same order of magnitude. We set $\delta =0.05$ for the numerical results presented below and vary $\xi$ from 0.01 to 0.2.

We focus on cases where the viscosity ratio $\lambda =\mu_2/\mu_3 \ll 1$, i.e. we consider the liquid film to be very viscous as compared with the air layer. In the limit $\lambda \rightarrow 0$ the no-slip boundary condition effectively applies in the air layer at the interface with the liquid film, and (4.1a) simplifies to (2.1a). In the liquid layer, we expect the terms in (4.1b) proportional to $\lambda ^2$ and to $(\partial p_2 / \partial r)(r,t)$ to be subdominant: they are the terms appearing from the shear stress imposed by the air layer, which should be negligible for a small viscosity ratio $\lambda$. Neglecting them is equivalent to considering free slip in the liquid layer at the air–liquid interface. By making these assumptions, (2.1) continues to hold and the evolution equation of the viscous film height (4.1b) simplifies to

(4.3)\begin{equation} {\frac{\partial h_2}{\partial t}}(r,t) = \frac{\lambda}{3} \frac{1}{r} {\frac{\partial }{\partial r}}\left(r h_2^3(r,t) {\frac{\partial p_3}{\partial r}}(r,t)\right).\end{equation}

We discuss next what can be learned from analysing these simplified equations at small $\lambda$ for the air layer profile and deformation of the liquid film. We will also verify in § 4.4 and Appendix B that neglecting the terms that are ${O}(\lambda )$ in (4.1a) and those that are ${O}(\lambda ^2 ; \partial p_2/\partial r)$ in (4.1b), i.e. considering no slip of air and free slip of liquid at the air-film interface, is justified for $\lambda \ll 1$.

4.1. Behaviour of the air layer

In the limit of very viscous films $\lambda \ll 1$, insight into the behaviour of the air layer thickness $H(r,t)$ can be gained since the expression of the pressures $p_2(r,t)$ in the air layer (2.1b) and $p_3(r,t)$ in the liquid film (4.2) are similar. Indeed, subtracting (4.2) from (2.1b) yields

(4.4)\begin{equation} (\delta+\xi) p_2(r,t) - \xi p_3(r,t) = 2-\frac1r{\frac{\partial }{\partial r}}\left(r{\frac{\partial H}{\partial r}}(r,t)\right).\end{equation}

Assuming $\xi p_3(r,t) \ll (\delta +\xi )p_2(r,t)$, the system (2.1a), (4.4) and (2.1c) for the air layer thickness $H(r,t)$ is closed and can be solved similarly to (2.1). The discussion of § 2 for the deposition on a rigid wall can then also be applied here, except that $\delta$ now becomes $\delta +\xi$ to account for the liquid film. In particular the dimple radius approaches $(2(\delta +\xi )/3)^{1/2}$, the pressure in the dimple approaches $2/(\delta +\xi )$, and the results given by (2.3) for the height at the axis of symmetry and minimum height apply as well (with $\delta \rightarrow \delta +\xi$). This was recognized by Yiantsios & Davis (1990) for the deposition of a droplet on a bath of its own fluid ($\xi =\delta$), neglecting any influence of the pressure in the bath ($p_3$ implicitly assumed to be zero), and was also used to study the approach of a droplet towards another droplet (Yiantsios & Davis 1991).

4.2. Response of a thin viscous film under constant load

Before describing the full problem we consider a simplified, generic, situation where the thin fluid film is exposed to an external load $p_e$ and with no external shear stresses. To account for this situation we change the non-dimensional units to scale (4.2) and (4.3) using $t=\hat {t} t_c, h_2=\hat {h}_2(\hat {\boldsymbol {x}},t) h_s, \boldsymbol {\nabla }=\hat {\boldsymbol {\nabla }}/x_c, p_e=\hat {p}_e(\hat {\boldsymbol {x}},t) p_c$; with $p_c$ the characteristic magnitude of the load which is distributed over a characteristic length $x_c, h_s$ the initial height of the liquid layer, and $t_c=3\xi x_c^{4}/\lambda h_s^{3}$ the characteristic time. We temporarily use Cartesian coordinates and denote by $\hat {\boldsymbol {x}}$ the spatial variable. Using these scalings, we obtain the following dimensionless equation for the liquid film height:

(4.5)\begin{equation} {\frac{\partial \hat{h}_2}{\partial \hat{t}}}(\hat{\boldsymbol{x}},\hat{t}) = \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} ( \hat{h}_2^3(\hat{\boldsymbol{x}},\hat{t}) \hat{\boldsymbol{\nabla}} (-\hat{\nabla}^2 \hat{h}_2(\hat{\boldsymbol{x}},\hat{t}) + \beta_{cap} \hat{p}_e(\hat{\boldsymbol{x}},\hat{t}) ) ).\end{equation}

The parameter $\beta _{cap}=p_c x_c^{2} \xi /h_s$ represents the ratio of the characteristic pressure force to the characteristic capillary force. We consider a long and initially flat viscous film, $\hat {h}_2(\hat {\boldsymbol {x}},0)=1$, exposed to a load uniform in time: $\hat {p}_e(\hat {\boldsymbol {x}},\hat {t})=\hat {p}_e(\hat {\boldsymbol {x}})\mathcal {H}(\hat {t})$ where $\mathcal {H}$ is the Heaviside function. By assuming small deformations, we can derive the film height profile analytically by adapting a procedure recently used to study flows induced in glassy polymer films (Pedersen et al. 2021a). Indeed, as long as the film deformations are small, (4.5) can be linearized defining $\hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t})$ such that $\hat {h}_2(\hat {\boldsymbol {x}},\hat {t})=1-\beta _{cap}\hat {\epsilon } (\hat {\boldsymbol {x}},\hat {t})$. By assuming $|\beta _{cap}\hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t})|\ll 1$ and neglecting terms that are ${O}((\beta _{cap}\hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t}))^2)$, we obtain from (4.5) the following linear equation:

(4.6)\begin{equation} \left({\frac{\partial }{\partial \hat{t}}}+\hat{\nabla}^{4}\right)\hat{\epsilon}(\hat{\boldsymbol{x}},\hat{t}) ={-}\hat{\nabla}^2 \hat{p}_e(\hat{\boldsymbol{x}},\hat{t}).\end{equation}

The associated Green's function is

(4.7)\begin{equation} G_{cap}(\hat{\boldsymbol{x}},\hat{t})=\frac{\mathcal{H}(\hat{t})}{(2{\rm \pi})^{{2}}} \int_{\mathbb{R}^{{2}}} \, {\rm e}^{-|\boldsymbol{k}|^{4} \hat{t}} \, {\rm e}^{{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \hat{\boldsymbol{x}}} \, {\rm d}\boldsymbol{k},\end{equation}

and the solution of (4.6) for an arbitrary load $\hat {p}_e(\hat {\boldsymbol {x}})$ is given by

(4.8)\begin{equation} \hat{\epsilon}(\hat{\boldsymbol{x}},\hat{t})={-}\int_{\mathbb{R}^{{2}}} \int_0^{\hat{t}} G_{cap}(\hat{\boldsymbol{x}}-\boldsymbol{x}^\prime; \hat{t}-t^\prime)\hat{\nabla}^2\hat{p}_e(\boldsymbol{x}^\prime)\, {\rm d} t^\prime\,{\rm d} \boldsymbol{x}^\prime.\end{equation}

When the load is a Dirac delta function, $\hat {p}_e(\hat {\boldsymbol {x}})=\delta _{Dirac}(\hat {\boldsymbol {x}})$, this simplifies to

(4.9)\begin{equation} \hat{\epsilon}(\hat{\boldsymbol{x}},\hat{t})={-}\int_0^{\hat{t}} \hat{\nabla}^2 G_{cap}(\hat{\boldsymbol{x}};t^\prime) \,{\rm d} t^\prime.\end{equation}

For a confined uniform load, $\hat {p}_e(\hat {\boldsymbol {x}})=1/{\rm \pi}$ if $\|\hat {\boldsymbol {x}}\|< 1$, the solution (4.8) becomes

(4.10)\begin{equation} \hat{\epsilon}(\hat{\boldsymbol{x}},\hat{t})={-}\frac{1}{\rm \pi} \int_0^{\hat{t}} \int_{\mathcal{C}}\hat{\boldsymbol{\nabla}} G_{cap}(\hat{\boldsymbol{x}}-\boldsymbol{x^\prime};t^\prime)\boldsymbol{\cdot}\boldsymbol{n} \,{\rm d}\ell^\prime\,{\rm d} t^\prime, \end{equation}

where $\boldsymbol {n}$ is the outward normal to $\mathcal {C}=\{\hat {\boldsymbol {x}}, \|\hat {\boldsymbol {x}}\|=1\}$.

It is tempting to try to compute the solution for a Dirac load (4.9): this generally gives the self-similar intermediate asymptotic solution (Barenblatt 1996) towards which solutions to any confined load would converge, similarly to the levelling scenario (Benzaquen, Salez & Raphaël 2013; Benzaquen et al. 2014) where an initial deformation is allowed to relax. However, the calculations show a singularity in time that cannot be integrated, which suggests that (4.6) does not possess a self-similar universal attractor. Instead we compute the response to a confined uniform load, which is relevant for the droplet settling case. By using (4.7) and writing spatial variables in polar coordinates, in particular letting $\hat {\boldsymbol {x}}=(\hat {r}\cos \theta,\hat {r}\sin \theta )$, and after using the identities $\int _{0}^{2{\rm \pi} }\exp ({\rm i}k\cos \theta )\,{\rm d}\theta =2{\rm \pi} J_0(k), \int _{0}^{2{\rm \pi} } \cos (\theta )\exp (-{\rm i}k\cos \theta )\,{\rm d}\theta =-2\,{\rm i}{\rm \pi} J_1(k)$ where $J_n$ is the Bessel function of the first kind of order $n$, (4.10) simplifies to the following Hankel transform:

(4.11a,b)\begin{equation} \hat{\epsilon}(\hat{r},\hat{t}) = \psi\left(\frac{\hat{r}}{\hat{t}^{1/4}},\hat{t}\right), \quad \psi(y,\hat{t})=\frac{\hat{t}^{1/4}}{\rm \pi} \int_0^{+\infty} J_0(uy)J_1 \left(\frac{u}{\hat{t}^{1/4}}\right) \frac{1- {\rm e}^{{-}u^4}}{u^2}\, {\rm d}u.\end{equation}

As expected, this does not converge towards a self-similar solution as illustrated in figure 5(a,b). This is a peculiar property of this problem and is at odds with the levelling scenario mentioned above and with the elastohydrodynamic case discussed in the next section. In particular, this means that the long-term behaviour of the interface not only depends on the total weight applied on it but also on its specific distribution.

Figure 5. (a) Deformation profiles of a viscous film under an external load for $\hat {t}=10^{-1}, 10, 10^{3}$ and $10^{7}$ in the linear limit, using an initially flat profile and a confined uniform load defined numerically as $\hat {p}_e(\hat {r})=p_0(1+\text{erf}[c(1-\hat {r})])$ with $c=50, \text{erf}$ the error function, $p_0$ chosen such that $\int _{\mathbb {R}^+} 2{\rm \pi} p_e(\hat {r}) \hat {r}\, {\rm d}\hat {r}=1$. (b,c) Evolution at the axis of symmetry ($\hat {r}=0$) of the (b) height and (c) pressure from simulations of the nonlinear evolution equation (4.5).

In the droplet settling dynamics, we are particularly interested in what happens at the axis of symmetry. The integral in (4.11a,b) admits a closed form expression at $\hat {r}=0$ in terms of hypergeometric functions, shown in Appendix A, from which the following asymptotic expansion as $\hat {t}\rightarrow +\infty$ can be found:

(4.12a)$$\begin{gather} \hat{\epsilon}(0,\hat{t}) = \frac{\ln(\hat{t})}{8{\rm \pi}} + k + \mathcal{O}(\hat{t}^{{-}1/2}), \quad k = \frac{2-3\gamma_{em}+4\ln(2)}{8{\rm \pi}}\simeq0.121, \end{gather}$$

(4.12b)$$\begin{gather}\hat{h}_2(0,\hat{t}) \simeq 1-\beta_{cap} k - \frac{\beta_{cap}}{8{\rm \pi}}\ln(\hat{t}), \end{gather}$$

with $\gamma _{em}\simeq 0.577$ the Euler–Mascheroni constant. It is also interesting to look at the capillary pressure at the axis of symmetry, -$\hat {\nabla }^2 \hat {h}(0,\hat {t})$, which in the linear approximation is given by $\hat {p}(0,\hat {t})=-\beta _{cap} \hat {\nabla }^2 \hat {\epsilon }(0,\hat {t})$. By using (4.11a,b), we can find that the asymptotic expansion of this quantity as $\hat {t}\rightarrow \infty$ is

(4.13)\begin{equation} \hat{p}(0,\hat{t})= \frac{\beta_{cap}}{8\sqrt{\rm \pi}}\hat{t}^{{-}1/2} + \mathcal{O}(\hat{t}^{{-}1}).\end{equation}

In this linear and asymptotic approximation, the time to reach $\hat {h}_2(0,t)=0$ is $\hat {\tau }= {\rm e}^{8{\rm \pi} (1/\beta _{cap}-k)}\simeq 0.048\, {\rm e}^{8{\rm \pi} /\beta _{cap}}$. We therefore expect (4.12) to be valid up to $\hat {t}\approx \hat {\tau }$, after which nonlinear effects become important and prevent this singularity. We also expect that for this asymptotic linear solution to have enough time to develop before nonlinear effects appear we need an upper bound on $\beta _{cap}$, i.e. a weak enough load.

We verify the results from this minimal model by comparing them against simulations of the nonlinear evolution equation (4.5) in figure 5(b,c). We indeed observe that the linear approximation (4.6) and expansions (4.12) and (4.13) are accurate up to $\hat {t}\simeq \hat {\tau }$ when $\beta _{cap} \lesssim 1.5$, while for larger values of $\beta _{cap}$ the asymptotic regime is not reached before nonlinearities appear. We note that the fact that $\hat {\epsilon }(0,t)$ does not evolve as a power law is consistent with the lack of self-similarity of the solution. Yet, the asymptotic pressure evolution (4.13) is nevertheless scale invariant and universal. We hypothesize that the logarithmic evolution given by (4.12), $\hat {\epsilon }(0,\hat {t}) \simeq \ln (\hat {t})/8{\rm \pi} + k$, is also universal but that the prefactor $k$ depends on the functional form of the load. We verified this numerically for a few other loads, but these are not presented here.

4.3. Response of a viscous film to a settling droplet

We now come back to the complete droplet settling case, where we anticipated from the discussion of § 4.1 that the dimple radius evolves towards a constant, $(2(\delta +\xi )/3)^{1/2}$, and that the pressure in the dimple is also constant, equal to $2/(\delta +\xi )$. Assuming that this is indeed the case, we expect the results derived in § 4.2 to apply in the case of droplet settling as well. Accordingly, we define $r_c=(2(\delta +\xi )/3)^{1/2}$ and $p_c=2{\rm \pi} /(\delta +\xi )$; this gives $t_c=4\xi (\delta +\xi )^2/3h_s^3\lambda, \beta _{cap}=4{\rm \pi} \xi /3h_s$ and the rescaled height $\hat {h}_2(\hat {r},\hat {t})$ is governed by (4.5) where the external pressure $\hat {p}_e$ is now that in the air layer, $\hat {p}_2(\hat {r},\hat {t})$.

We first need to verify that (4.4) applies, i.e. that $\xi p_3(r,t) \ll (\delta +\xi )p_2(r,t)$. Assuming $\beta _{cap} \lesssim 1$ for an asymptotic regime to be reached and $\hat {t}\lesssim \hat {\tau } \simeq 0.05 \, {\rm e}^{8{\rm \pi} /\beta _{cap}}$, i.e. ${t\lesssim 0.1 \xi (\delta +\xi )^2 h_s^{-3} \lambda ^{-1} \, {\rm e}^{6h_s/\xi }}$ for the deformations of the liquid film interface to remain small, the full film profile can be derived as shown in § 4.2. In particular, (4.13) gives the pressure at $r=0$ in the liquid film as

(4.14)\begin{equation} p_3(0,t) = \frac{\beta_{cap}}{8\sqrt{\rm \pi}} \frac{h_s}{\xi r_c^2} \left(\frac{t}{t_c}\right)^{{-}1/2} = \frac{\sqrt{\rm \pi}}{2\sqrt{3}} \left(\frac{\xi}{h_s^3 \lambda t}\right)^{1/2} . \end{equation}

In this regime ($\beta _{cap} \lesssim 1, \hat {t}\lesssim \hat {\tau }$), we therefore expect the pressure in the liquid to continuously decrease and the condition $\xi p_3(r,t) \ll (\delta +\xi ) p_2(r,t)$ will eventually be satisfied: this justifies a posteriori the discussion of § 4.1 and the scaling we have chosen for $r_c$ and $p_c$. We can now expect from (4.11a,b) and (4.12) that the sheet profile $h_2(r,t)$ evolves according to

(4.15a)$$\begin{gather} 1-\frac{h_2(r,t)}{h_s}=\beta_{cap}\psi\left(\frac{\hat{r}}{\hat{t}^{1/4}},\hat{t}\right) =\frac{4{\rm \pi}\xi}{3h_s}\psi\left(\frac{r}{(h_s^3\lambda t)^{1/4}},\frac{3h_s^3\lambda}{4\xi(\delta+\xi)^2} t\right), \end{gather}$$

(4.15b)$$\begin{gather}1-\frac{h_2(0,t)}{h_s} = \beta_{cap} \left(\frac{\ln(t/t_c)}{8{\rm \pi}}+k\right) \simeq \frac{\xi}{6 h_s}\left[ \ln\left(\frac{h_s^3\lambda}{\xi(\delta+\xi)^2}t\right) + 2.75\right]. \end{gather}$$

4.4. Numerical results

In figure 6 we show the time evolution of the pressure in the liquid film and its height from numerical solutions of the complete nonlinear equations, i.e. considering the normal stress balances (2.1b) and (4.13), the force balance (2.1c) and the governing equations (4.1) accounting for the full shear stress balance. We present numerical results for $\lambda =10^{-5}$. The results are in close agreement with the theoretical expectations (4.14) and (4.15) for $\beta \lesssim 1$, a condition required for the asymptotic regime to be reached. We also show in Appendix B that $\lambda =10^{-5}$ is small enough for (2.1a) and (4.3) to apply, i.e. to consider free slip in the liquid film and no slip in the air layer at the interface between the two. Indeed, figures 12 and 13 show that considering the full viscous stress balance through (4.1) is equivalent to solving the simplified system considering the no-slip and free-slip conditions.

Figure 6. Height evolution of the coated liquid film for a droplet settling with $\delta =0.05$. Pressure (a) and deformation (b) of the liquid film at the axis of symmetry ($r=0$) for $\delta =0.05, \lambda =10^{-5}$ with $h_s=0.1, \xi =0.2, 0.05, 0.01$ and $h_s=1, \xi =0.2$, corresponding to $\beta =8.4, 2.1, 0.4$ and 0.8, respectively. (c) Normalized deformation of coated the liquid film shown here for $\lambda =10^{-5}, h_s=0.1, \xi =0.2$, corresponding to $\beta =0.8$, and at times $t=10^2, 10^3, 10^4$ and $10^5$.

The results presented in figure 7 for the air layer dynamics show that the discussion of § 4.1 indeed applies. The long-term asymptotics of the air film thickness $H(r,t)$ is similar to that of the solid case, with $\delta \rightarrow \delta +\xi$. We note that the air layer thickness profile has a dimple structure while both the droplet and the liquid film profiles are monotonically increasing near the axis of symmetry and that the neck of the dimple is located at the inflection point of the film profile. As discussed by Duchemin & Josserand (2020) when a droplet settles on a solid surface, we expect that this neck approaches a self-similar form, $H(r,t)=H_{min}(t) F(({r-r_d(t)})/{\ell (t)})$, where $F$ is linear for $r< r_d$ and quadratic for $r>r_d$, and where $\ell (t)$ is the radial extent of the neck region. This is verified in figure 7(b), where we used $\ell (t)=(\delta +\xi )^{1/2} t^{-1/4}$ following the scaling analysis presented in § 2.3.

Figure 7. Evolution of the air layer for the settling of a droplet with $\delta =0.05$ on a viscous liquid film. (a) Profiles for $\xi =0.2, h_s=1, \lambda =10^{-5}$ of the droplet $h_1(r,t)$ and the liquid film $h_2(r,t)$ at $t=5000$. (b) Air film profile $H(r,t)=h_1(r,t)-h_2(r,t)$ for $\xi =0.05, h_s=1, \lambda =10^{-5}$ at $t=50, 10^2, 10^3, 10^4$ and $10^5$. The neck is predicted to be located at $r=(2(\delta +\xi )/3)^{1/2}\simeq 0.26$. The inset shows the rescaled neck structure near $r=r_d(t)$ with $\ell (t)=(\delta +\xi )^{1/2} t^{-1/4}$. $(c,d)$ Thickness of the air film at $r=0$ (c) and at its minimum ($d$) for $\delta =0.05, \lambda =10^{-5}$ with $h_s=0.1, \xi =0.2, 0.05, 0.01$ and $h_s=1, \xi =0.2$, corresponding to $\beta _{cap}=8.4, 2.1, 0.4$ and 0.8, respectively.

When the deformations of the liquid film become large, the pressure $p_3(r,t)$ increases and can no longer be neglected in (4.4). This leads the air layer into another regime. For $\beta _{cap} \lesssim 1$, the time at which deformations become large can be estimated from § 4.2 as $\tau \approx 0.05 t_0 \, {\rm e}^{8{\rm \pi} /\beta _{cap}}\simeq 0.1 \xi (\delta +\xi )^2 h_s^{-3}\lambda ^{-1} \, {\rm e}^{6h_s/\xi }$. The effect of large deformations on the liquid film evolution can be seen in figure 6(a,b) for the case $\beta _{cap}=8.3$ at time $t\gtrsim 10^4$. As shown in figure 7(c), the resulting effect on the air film is characterized by a sudden decrease of its thickness at the axis of symmetry. In fact, figure 8(a) shows that when nonlinearities of the liquid film are present the profile evolves towards a more complicated shape presenting a local maximum not located at the axis of symmetry and a minimum not located at the neck. The pressure profile in the air, shown in figure 8(b), also significantly deviates from the near-uniform pressure characteristic of the classical dimple (figure 2$c$), and the assumption $(\delta +\xi )p_2(r,t) \gg \xi p_3(r,t)$ is no longer valid. Such a shape of the film is not uncommon in thin film drainage and has been observed, both experimentally and numerically, for droplets and bubbles approaching a rigid substrate (Clasohm et al. 2005; Ajaev, Tsekov & Vinogradova 2007, 2008), and has been referred to as rippled deformation, or as a wimple (Chan et al. 2011). This wimple occurs due to a competition of two effects in the lubrication pressure: the capillary-driven deformations of the drop/bubble, along with an additional physical effect. In prior work this additional effect originated from surface forces, such as van der Waals interaction or electrostatic effects. Here, we do not consider such effects, but it is the nonlinear deformations of the liquid film towards which the droplet settles that lead to a wimple. Such a rippled shape also seems to develop in the numerical simulations of Duchemin & Josserand (2020) during the late stage of a large drop settling on a liquid film. The results we presented above are valid before the appearance of a wimple; an understanding of this shape and the associated drainage dynamic would be an important step towards an understanding of droplet settling when the parameter $\beta _{cap}=4{\rm \pi} \xi /3h_s$ becomes large, i.e. on very thin films.

Figure 8. (a) Interface profile of the air layer thickness for a droplet settling on a viscous liquid film at $t=500,$ $10^4$ and $10^5$, for $\xi =0.2, h_s=0.1, \lambda =10^{-3}$. (b) Corresponding pressure profiles in the air layer ($p_2$) and in the liquid film ($p_3$) at $t=10^5$. (c) Interface profile of the air layer thickness for a droplet settling on an elastic sheet supported by a thin viscous layer at $t=10\ 000, 20\ 000$ and $40\ 000$, for $\alpha =10, h_s=1, \lambda =10^{-3}$. (d) Corresponding pressure profiles in the air layer ($p_2$) and in the viscous film ($p_3$) at $t=40\ 000$.

5. Solid substrate coated with a thin viscous film and an elastic sheet on top

We now consider the case where the interface of the thin viscous film is substituted by an elastic sheet (figure 1$d$). We neglect any tension in the sheet and only consider its bending, an assumption valid when the deformation of the sheet remains small compared with its thickness (Landau & Lifshitz 1959; Tulchinsky & Gat 2016). This means that the effect of the viscous shear stress by the air layer on the elastic sheet is very small compared with the normal stress pressure exerts on it. However, when the deformation becomes of the same order as the thickness of the sheet, tension in the sheet should be taken into account. This can be done using a nonlinear plate model such as the Föppl–von Kármán equations (Landau & Lifshitz 1959). Here, we consider that we are always in the case when deformations are small compared with the thickness of the elastic sheet. This leads to a thin film equation for $h_2(r,t)$ similar to (2.1a) and (4.3), where the pressure enters the equation as $p_3^\star (r^\star,t^\star )=p_2^\star + B \nabla ^{\star 4}h_2^\star (r^\star,t^\star )$. The operator $\nabla^4(\cdot ) = \nabla^2(\nabla^2(\cdot ))$ is the bi-Laplacian and $B=Ed^{\star 3}/12(1-\nu ^2)$ is the bending stiffness of the sheet, with $d^\star, E$ and $\nu$ its thickness, Young's modulus and Poisson ratio, respectively. Using the scaling defined in § 2, the non-dimensional governing equation for the height of the film reads

(5.1a)$$\begin{gather} {\frac{\partial h_2}{\partial t}}(r,t) = \frac{\lambda}{12} \frac{1}{r} {\frac{\partial }{\partial r}}\left(r h_2^3(r,t) {\frac{\partial p_3}{\partial r}}(r,t)\right), \end{gather}$$

(5.1b)$$\begin{gather}\alpha (p_3(r,t)-p_2(r,t)) =\nabla^4 h_2 (r,t), \end{gather}$$

with $\lambda =\mu_2/\mu_3, \mu_3$ the viscosity of the fluid supporting the elastic sheet. The effect of bending is quantified by the parameter $\alpha ={\rm \Delta} \rho g a^{\star 4}/B=(a^\star /\ell _g^\star )^4$, analogous to an elastic Bond number, which compares the droplet radius with the elastogravity length $\ell _g^\star =(B/ g {\rm \Delta} \rho )^{1/4}$. At the initial time $t=0$ and far from the droplet, the height of the sheet is $h_2(r,0)=h_2(r\rightarrow \infty,t)=h_s$. We set $\delta =0.05$ for the numerical results presented below.

We note that for a typical density ratio between gas and liquids ${\rm \Delta} \rho =10^3\ {\rm kg\ m}^{-3}$, and considering very soft elastic sheets, e.g. $\nu \simeq 0.5, E={O}(10^9\ {\rm Pa}), d={O}(100\ {\rm nm}), \ell _g$ is approximately 50 ${\rm \mu} $m, which is also the typical upper bound we expect for droplet size in order for the lubrication scalings presented in § 2.1 to be valid. We are therefore mostly interested in the cases where $\alpha \leq 1.$

5.1. Response of an elastic sheet supported by a viscous film under constant load

We start by looking at the general case of an external load $p_e(\boldsymbol {x})$ exerted over an elastic sheet supported by a thin viscous film, before coming back to the complete droplet settling situation. The rescaled governing equation of the sheet height for a load $\hat{p}_e(\hat {\boldsymbol {x}})$ reads

(5.2)\begin{equation} {\frac{\partial \hat{h}_2}{\partial \hat{t}}}(\hat{\boldsymbol{x}},\hat{t}) = \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} ( \hat{h}_2^3(\hat{\boldsymbol{x}},\hat{t})\hat{\boldsymbol{\nabla}} (\hat{\nabla}^4 \hat{h}_2(\hat{\boldsymbol{x}},\hat{t}) + \beta_{el} \hat{p}_e(\hat{\boldsymbol{x}}) ) ),\end{equation}

with $t=\hat {t} t_c, h_2=\hat {h}_2(\hat {\boldsymbol {x}},t) h_s, \boldsymbol {\nabla }=\hat {\boldsymbol {\nabla }}/x_c, p_e=\hat {p}_e(\hat {\boldsymbol {x}}) p_c$; with $p_c$ the characteristic magnitude of the load which is distributed over a characteristic length $x_c, h_s$ the initial height of the elastic sheet, and $t_c=12\alpha x_c^6/h_s^3 \lambda$ the characteristic time. The parameter $\beta _{el}=\alpha x_c^4 p_c / h_s$ characterizes the ratio of the force due to the applied pressure to the characteristic bending force. For small deformations of the elastic sheet, (5.2) can be linearized by defining $\hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t})$ such that $\hat {h}_2(\hat {\boldsymbol {x}},\hat {t})=1-\beta _{el}\hat {\epsilon } (\hat {\boldsymbol {x}},\hat {t})$ and considering $\vert \beta _{el} \hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t})\vert \ll 1$. This leads to the following equation when neglecting terms that are ${O}((\beta _{el} \hat {\epsilon }(\hat {\boldsymbol {x}},\hat {t}))^2)$ in (5.2):

(5.3)\begin{equation} \left({\frac{\partial }{\partial \hat{t}}}-\hat{\nabla}^{6}\right) \hat{\epsilon}(\hat{\boldsymbol{x}},\hat{t}) ={-}\hat{\nabla}^2 \hat{p}_e(\hat{\boldsymbol{x}}).\end{equation}

We can then derive analytically the film height profile in a similar manner as in the work of Tulchinsky & Gat (2016), who studied the response of elastic sheets in the context of impact mitigation. Indeed, the Green's function of (5.3) is

(5.4)\begin{equation} G_{el}(\hat{\boldsymbol{x}},\hat{t})=\frac{\mathcal{H}(\hat{t})}{(2{\rm \pi})^{{2}}} \int_{\mathbb{R}^{{2}}} \, {\rm e}^{-|\boldsymbol{k}|^{6} \hat{t}} \, {\rm e}^{{\rm i} \boldsymbol{k} \boldsymbol{\cdot} \hat{\boldsymbol{x}}} \, {\rm d}\boldsymbol{k},\end{equation}

and the solutions to a Dirac load and to a confined uniform load are given by (4.9) and (4.10), respectively, provided $G_{el}$ be replaced with $G_{cap}$. The solution to a Dirac load, $\hat {p}_e(\hat {\boldsymbol {x}})=\delta _{Dirac}(\hat {\boldsymbol {x}})$, exists and simplifies to the following Hankel transform:

(5.5a,b)\begin{equation} \hat{\epsilon}(\hat{r},\hat{t}) = \hat{t}^{1/3} \phi\left(\frac{\hat{r}}{\hat{t}^{1/6}}\right), \quad \phi(y)=\frac{1}{2{\rm \pi}}\int_0^{+\infty} J_0(uy)\frac{1- {\rm e}^{{-}u^6}}{u^3}\,{\rm d}u, \end{equation}

where $\hat {r}=\lVert {\hat {\boldsymbol {x}}}\rVert$. At the axis of symmetry $\hat {r}=0$ this gives

(5.6a)$$\begin{gather} \hat{\epsilon}(0,\hat{t})=\frac{\varGamma(2/3)}{4{\rm \pi}}\hat{t}^{1/3}, \end{gather}$$

(5.6b)$$\begin{gather}\hat{h}(0,\hat{t})=1-\left(\frac{\hat{t}}{\hat{\tau}}\right)^{1/3}, \quad \hat{\tau} = \left(\frac{4{\rm \pi}}{\varGamma(2/3)\beta_{el}}\right)^3, \end{gather}$$

with $\hat {\tau } \simeq 799/\beta _{el}^3$ the characteristic time scale of the process and $\varGamma$ Euler's gamma function ($\varGamma (1/3)\simeq 2.68, \varGamma (2/3)\simeq 1.35$). The bending pressure $\hat {p}(\hat {r},\hat {t})=\beta _{el}\hat {\nabla }^4\hat {\epsilon }(\hat {r},\hat {t})$ at the axis of symmetry is given by

(5.7)\begin{equation} \hat{p}(0,\hat{t})=\beta_{el} \frac{\varGamma(1/3)}{12{\rm \pi}} \hat{t}^{{-}1/3}.\end{equation}

The self-similar solution (5.5a,b) is only valid when the load is a Dirac distribution. However, it can be expected that it is also an intermediate asymptotic solution (Barenblatt 1996) towards which solutions to other loads converge. We verify this in figure 9(a) for a uniform and confined load: $\hat {p}_e(\hat {r}<1)=1/{\rm \pi}$ and $\hat {p}_e(\hat {r}>1)=0$. This is a property also present for the levelling scenario for elastic (Pedersen, Salez & Carlson 2021b) and capillary (Benzaquen et al. 2013, 2014) interfaces, including polymer melts (McGraw et al. 2013), when an initial local deformation is allowed to relax, as well as for gravity currents (Ball & Huppert 2019) governed by a similar equation. However, under load and contrary to the aforementioned cases, the deformation grows with time and the linear approximation always eventually breaks down so that (5.5a,b) is only valid in a hypothetical intermediate regime. In figure 9(b) we show by solving numerically the nonlinear equation (5.2) in the case of a confined uniform load that the self-similar regime is indeed reached for $\beta _{el} < 1$. The solution to the linear equation predicts that the height at the axis of symmetry becomes zero in finite time: this singularity is prevented by nonlinear effects that become significant from time $\hat {t}\simeq \hat {\tau }/10$. For larger loads, $\beta _{el} > 1$, the asymptotic self-similar regime does not have enough time to be reached before nonlinear effects appear.

Figure 9. Deformation profile of an elastic sheet supported by viscous fluid and exposed to an external load. (a) Solution of the linear evolution equation (5.3) in cylindrical coordinates for (right panel) a narrow load approaching a Dirac distribution, $p_e(\hat {r})=p_0\exp [-(\hat {r}/a)^2]/\hat {r}$ with $a=1/200$ and (left panel) a confined uniform load, $p_e(\hat {r})=p_0(1+\text{erf}[a(1-\hat {r})])$ with $a=50, \text{erf}$ the error function. In both cases $p_0$ is chosen such that $\int _{\mathbb {R}^+} 2{\rm \pi} p_e(\hat {r}) \hat {r} \,{\rm d}\hat {r}=1$. The rescaled solution converges to the self-similar solution (5.5a,b), and faster for the narrow load (right) than for the distributed one (left). (b) Evolution of the deformation at the axis of symmetry ($\hat {r}=0$) from the nonlinear evolution equation (5.2) and for a confined uniform load.

5.2. Response of an elastic sheet to a settling droplet

We return to the case of a settling droplet where the weight of the droplet acts on the sheet, set by (2.1c), so that we can expect to be able to characterize its self-similar asymptotic response without a priori knowledge of the structure of the air layer. We let $r_c=(2\delta /3)^{1/2}$ and $p_c=2{\rm \pi} /\delta$, the values obtained in the rigid case. Then $\beta _{el}=8{\rm \pi} \alpha \delta /9 h_s, t_c=32\alpha \delta ^3/9h_s^3\lambda$, and the rescaled height $\hat {h}_2(\hat {r},\hat {t})$ is governed by (5.2), with $p_2(\hat {r},\hat {t})$ in lieu of $p_e(r)$. For $\beta _{el} \lesssim 1$ the sheet reaches a self-similar evolution, until nonlinear effects appear at time $\hat {\tau }\simeq 800/\beta _{el}^3 \simeq 37 (h_s/\alpha \delta )^3$, when $h_2(0,t) \ll h_s$. For $\beta _{el} \gtrsim 1$, the sheet does not enter the self-similar regime before nonlinear effects play a role.

For $\beta _{el} \lesssim 1$ and $t< t_c\hat {\tau }$, we then expect from (5.5a,b) and (5.6) the following self-similar evolution of the elastic sheet:

(5.8a)\begin{gather} \frac{h_s-h_2(r,t)}{h_s-h_2(0,t)}=\phi\left(\frac{\hat{r}}{\hat{t}^{1/6}}\right) =\phi\left(\frac{r}{t^{1/6}} \times \left[\frac{12\alpha}{h_s^3\lambda}\right]^{1/6} \right), \end{gather}

(5.8b)\begin{gather}1-\frac{h_2(0,t)}{h_s} = \left(\frac{\hat{t}}{\hat{\tau}}\right)^{1/3} \simeq 0.2 \times \left(\frac{t}{\lambda \alpha^2}\right)^{1/3}, \end{gather}

while the pressure in the viscous film beneath the sheet at $r=0$ is found from (5.7)

(5.9)\begin{equation} p_3(0,t) = \frac{\varGamma(1/3)}{12{\rm \pi}} \frac{\beta_{el} h_s}{\alpha r_c^4} \left(\frac{t}{t_0}\right)^{1/3} \simeq 0.68 \left(\frac{\alpha}{\lambda h_s^3 t}\right)^{1/3}. \end{equation}

Figure 10 shows that, if $\beta _{el} \lesssim 1$ to ensure that the self-similar dynamics is reached, and for $t \lesssim t_0 \tau \simeq 130\lambda ^{-1}\alpha ^{-2}$ to ensure that deformations remain small, (5.8) can be verified. The results also show a third necessary condition, $\alpha \lesssim 1$, which is linked to the breakdown of the assumption of a localized load that we discuss next.

Figure 10. Deformations of an elastic sheet supported by a viscous layer for a droplet settling with $\delta =0.05$. (a) Normalized deformation of the elastic sheet, shown here for $\lambda =10^{-5}, h_s=10$, and (left panel) $\alpha =1$, (right panel) $\alpha =10$. (b) Deformation of the elastic sheet at the axis of symmetry ($r=0$) for $\lambda =10^{-5}$ and various values of $\alpha$ and $h_s$, leading to different values of $\beta _{el}$.

5.3. Numerical results

In figure 11 we show the evolution of the air layer thickness for a droplet with $\delta =0.05$ from numerical simulations of (2.1) and (5.1). For $\alpha < 1$, i.e. for droplets smaller than the elastogravity length, we observe the same behaviour of the air layer as in the solid case discussed in § 2.3. This gives a criterion for when elasticity effects can be neglected in the droplet settling process. In this scenario, we derived analytically and verified above the self-similar response of the elastic sheet.

Figure 11. Evolution of the air layer thickness for a droplet settling on an elastic sheet with $\delta =0.05$. (a) Profiles for $\alpha =1, h_s=1, \lambda =10^{-3}$ of the droplet $h_1(r,t)$ and the viscous film $h_2(r,t)$ at $t=250$. We note that the thickness of the elastic sheet is not represented to scale for clarity. (b) Profile of the air layer $H(r,t)=h_1(r,t)-h_2(r,t)$ for $\alpha =1, h_s=1, \lambda =10^{-3}$ and $t=250$, 1000, 5000 and $10^5$. (c–e) Evolution of the air layer thickness at the axis of symmetry ($r=0$), at its minimum value, and of the dimple radius for $\lambda =10^{-3}, h_s=1$ and $\alpha =0.01, 0.1, 1, 10$. (f –h) Likewise for $\alpha =1$ and various values of $\lambda$ and $h_s$ The dashed lines in (d) and (e) show the height and radius at the local minimum corresponding to the dimple neck, for $\alpha =10$ and $t\simeq 10^4$ the minimum radius is not located anymore at the neck.

For $\alpha \gtrsim 1$, we observe a significant deviation in the dynamics of the air layer as shown in figure 11(a–e). It continues to present a dimple, but the height at the axis of symmetry saturates and evolves towards a constant value. The minimum height, at the neck, continues to decrease while remaining larger than for the rigid case. Contrary to the rigid case, the radius of the neck of the dimple also monotonically increases with time. We note that this dimple is present while both the droplet and the sheet heights are increasing near the axis of symmetry, and its neck is approximately located at the inflection point of the second spatial derivative of $h_2$. Contrary to the case of a rigid solid coated by a thin viscous film, without an elastic sheet atop, this new behaviour not only depends on the property of the sheet characterized by $\alpha$, but also on the properties on the supporting viscous film height and viscosity through $h_s$ and $\lambda$ as show in figure 11(f –h). Namely, the thickness and dimple radius increase with increasing initial height and decreasing viscosity.

The continuous increase of the radius of the neck of the dimple helps to explain why the analytical results for the sheet height presented in the previous section fail when $\alpha \gtrsim 1$, even when $\beta _{el} \lesssim 1$ and $\hat {t}<\hat {\tau }$. The analysis is based on the results derived in § 5.1, which formally only applies if the load that the droplet puts on the sheet were concentrated into a single point in space. The results continue to apply asymptotically when the load is localized to a small region as compared with the lateral extent of the deformations of the sheet, which evolves as $\hat {r}_{sheet}(\hat {t})\sim 3.4 \hat {t}^{1/6}$, i.e. $r_{sheet}(t)\sim \alpha ^{-1/6} h_s^{1/2} \lambda ^{1/6} t^{1/6}$ (we define $r_{sheet}(t)$ as the radial position of the maximum of $h_2(r,t)$, see figure 10$a$). For large enough $\alpha$ the dimple radius increases as fast as the radial extent of the sheet deformations and the load exerted on the sheet cannot be assumed to be localized to a very small region anymore (figure 11$e$).

At long times $t\gtrsim \tau \simeq 37 (h_s/\alpha \delta )^3$ and regardless of the value of $\alpha$, nonlinear effects of the thin film become important and the behaviour completely changes. Specifically, the air film no longer forms a dimple but instead evolves into a rippled deformation similar to the capillary case (§ 4.4), as illustrated in figures 8(c,d) and 11(d,e), where it can be seen that the minimum thickness is not located anymore at the neck of the dimple, and that the pressure in the film below the elastic sheet is no longer negligible against the pressure in the air film above it.

6. Discussion

We have studied the settling of a droplet onto three classes of thin soft layers coating a rigid solid: a compressible elastic material, a viscous fluid, or an elastic sheet supported by a viscous film. First, when the soft layer is made of a compressible elastic material (§ 3), we found that depending on its softness the air layer separating the droplet and the elastic solid transitions from the dimple characterizing droplet settling onto a rigid wall to a layer with near-uniform thickness. This thickness is always larger than the minimum thickness of the rigid case. Second, we studied in § 4 the settling towards a thin viscous film, generalizing the analysis of Yiantsios & Davis (1990) for an infinite bath of fluid. We derived the long-term asymptotics of both the film and of the air layer separating the film and the droplet in the limit of small film deformations. The dimple shows the same structure as for the case of a rigid substrate, but its thickness is asymptotically larger by a factor $1+\xi /\delta = 1 + \sigma _1/\sigma _3$ with $\sigma _1$ and $\sigma _3$ the surface tension of the droplet and of the liquid film, respectively. When the film deformation becomes large, we observed rippled deformations which modify the drainage dynamics and in particular lead to a sudden decrease of the minimum air layer thickness. Finally, we described in § 5 a droplet settling towards an elastic sheet supported by a thin layer of viscous fluid. We found that when the droplet is smaller than the elastogravity length the dimple profile of the air layer is similar to the solid case, and derived analytically the response of the elastic sheet. For droplets of size similar to or larger than the elastogravity length, we observed numerically an increase of both the minimum thickness of the dimple and of its radius. In the three cases studied, the deformation of the surface always leads to an increase of the minimum thickness of the dimple, provided that these deformations are small. This suggests that surface deformations upon droplet settling can delay impact or coalescence of the droplet with the substrate.

We observed a shift in dynamics for large deformations of the viscous film coating a solid substrate, with or without an elastic sheet on top. The dimple evolves then towards a more complex structure which strongly affects the drainage process and, to our knowledge, had not been described before for interactions between droplets and compliant surfaces. This so-called rippled dimple occurs because of the nonlinear response of viscous films to large deformations. A better understanding of this regime could also have applications beyond settling droplets settling and impacting drops. Indeed, the response of a capillary interface under load appears for instance in the deformation of glassy materials by nanobubbles (Ren et al. 2020) and is relevant for lubricant-impregnated surfaces (Pack et al. 2017) as well as other industrial applications (Carou et al. 2009; Lunz & Howell 2018), while an elastic sheet supported by a viscous fluid can be found in various situations, from technological (Rogers, Someya & Huang 2010) and biological (Bongrand 2018) applications up to the geophysical scale (Michaut 2011).

In our study of the settling of a droplet onto a thin liquid film, we focused on very viscous films and verified that our choice of parameters justifies the assumptions required for no slip of air at the interfaces of both the droplet and film. While this assumption only alters the time scale of the dynamic on solid surfaces, where the no-slip condition continues to apply at one of the interfaces, both numerical (Yiantsios & Davis 1990; Duchemin & Josserand 2020) and experimental (Lo et al. 2017) studies have shown that slip can significantly alter the dynamics of droplet deposition on a liquid film. The analysis presented in § 4 henceforth only applies in a regime where the air layer is wide enough and the viscosity ratios large enough. It also applies when considering the presence of surfactants, which can alter the boundary conditions at the interfaces, in either the droplet or the film. It has long been known that a drop falling onto a dirty interface takes much more time to coalesce than on a clean one (Reynolds 1881; Rayleigh 1882), and indeed small amounts of surfactant have been shown to considerably delay the coalescence process (Amarouchene, Cristobal & Kellay 2001).

For a droplet settling on an elastic sheet, we rationalized the dynamics for droplets smaller than the elastogravity length and only considered bending stresses in the description of the sheet. When the deformations of the sheet become significant compared with its thickness, additional effects from stretching could appear. Including stretching in the elastic model requires a more complete description, e.g. the full Föppl–von Kármán equations (Landau & Lifshitz 1959).

Finally, in all the cases we have studied the minimum thickness of the air layer continuously decreases with time and there would eventually be direct contact between the droplet and the substrate. This usually occurs due to the influence of surface forces such as van der Walls or electrostatic interactions (Oron, Davis & Bankoff 1997; Israelachvili 2011), rarefied gas effects (Duchemin & Josserand 2012), nano-roughness (Kolinski, Mahadevan & Rubinstein 2014; Li, Vakarelski & Thoroddsen 2015) or other small scale phenomena which can otherwise be neglected for most of the drainage dynamics. We did not consider them in this study but isolated the coupling between the deformations of the droplet and of the substrate. They should be incorporated when the minimum thickness of the air layer reaches the order of a hundred nanometres or less as they lead to the rupture of the air film then and to contact (Couder et al. 2005; Chan et al. 2011). In the case of settling onto a viscous film, this supporting layer could also rupture due to surface forces (Zhang & Lister 1999; Carlson & Mahadevan 2016). Contact can occur in an axisymmetric manner (Chan et al. 2011), but Lo et al. (2017) showed experimentally that when a large drop settles on a rigid surface, symmetry eventually breaks and contacts occur earlier than for an axisymmetric situation. Axisymmetry could be preserved for deposition on thin liquid films in their experiments. It would be interesting to understand when the deposition process can become non-axisymmetric and, more generally, to investigate to which extent the observations made and conclusions drawn in this work hold for drops with larger Bond number.

This work also opens the question of possible effects of non-uniform substrates on the drop settling dynamics. Directional transport of droplets above a substrate can for instance be reached through a Leidenfrost dynamics above a textured solid (Lagubeau et al. 2011) or a Marangoni dynamics of a drop above a liquid film with a temperature gradient (Davanlou & Kumar 2015). Gradients of substrate stiffness (Style et al. 2013) or bending rigidity (Bradley et al. 2019) can also lead to the transport of a droplet. Our work gives a first minimal description of how gravitational settling dynamics of droplets can be affected by a thin compliant layer, but there are many natural extensions such as how nonlinear elastic effects, gradients in substrate properties and adhesive contact affect the flow.