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 (Reference Yiantsios and Davis1990) 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 Reference Batchelor1967), $\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 Reference Yiantsios and Davis1990):

(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. Reference Alnæs, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015) 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 (Reference Yiantsios and Davis1990): $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 Reference Derjaguin and Kussakov1939; Frankel & Mysels Reference Frankel and Mysels1962). 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 Reference Duchemin and Josserand2020). 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 (Reference Yiantsios and Davis1990), 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.

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 (Reference Yiantsios and Davis1990) 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 Reference Yiantsios and Davis1991).

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. Reference Pedersen, Ren, Wang, Carlson and Salez2021a). 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 Reference Barenblatt1996) towards which solutions to any confined load would converge, similarly to the levelling scenario (Benzaquen, Salez & Raphaël Reference Benzaquen, Salez and Raphaël2013; Benzaquen et al. Reference Benzaquen, Fowler, Jubin, Salez, Dalnoki-Veress and Raphaël2014) 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 (Reference Duchemin and Josserand2020) 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. Reference Clasohm, Connor, Vinogradova and Horn2005; Ajaev, Tsekov & Vinogradova Reference Ajaev, Tsekov and Vinogradova2007, Reference Ajaev, Tsekov and Vinogradova2008), and has been referred to as rippled deformation, or as a wimple (Chan et al. Reference Chan, Klaseboer and Manica2011). 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 (Reference Duchemin and Josserand2020) 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.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 (Reference Tulchinsky and Gat2016), 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 Reference Barenblatt1996) 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 Reference Pedersen, Salez and Carlson2021b) and capillary (Benzaquen et al. Reference Benzaquen, Salez and Raphaël2013, Reference Benzaquen, Fowler, Jubin, Salez, Dalnoki-Veress and Raphaël2014) interfaces, including polymer melts (McGraw et al. Reference McGraw, Salez, Bäumchen, Raphaël and Dalnoki-Veress2013), when an initial local deformation is allowed to relax, as well as for gravity currents (Ball & Huppert Reference Ball and Huppert2019) 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.