Pre-impact dynamics of a droplet impinging on a deformable surface

To determine the behaviour of the liquid droplet free surface, close to the point of initial contact which is $x=y=0$, we take the following scaling

where the scales of $t$ and $f$ come from the desire to study short time behaviour and from the form of equation (5), respectively. The $O(\varepsilon^{2})$ time scale is that expected for the traversing, at a unit non-dimensional velocity, of the thin air gap which has normal width of $O(\varepsilon^{2})$. The scales (9) lead us to take asymptotic expansions of the liquid velocity components and pressure in the following form

The vertical velocity scale is order unity due to the order unity downward velocity of the droplet, then the horizontal velocity scale follows from continuity. The large pressure scale arises from seeking a balance between the liquid acceleration and the pressure gradient.

Now, for $Re\gg 1$, upon substitution of scales (9) and expansions (10) into the governing equations (3), the leading order momentum equations and continuity equation are those of unsteady potential flow

The leading order kinematic condition (6) now reduces to

while the far-field droplet behaviour reduces to

from (5).

From equations (11) it can be shown that $P_{l}$ satisfies Laplace’s equation and due to the far-field boundedness (13) and condition (12) the pressure profile on the droplet interface $P(X,T)=P_{l}(X,0,T)$ and free-surface profile $F(X,T)$ are related by

The bar denotes the Cauchy principle value integral.

In the thin gas film between the droplet and the deformable surface, we assume that the vertical length scale is an order of magnitude smaller than the horizontal length scale to capture the slenderness of the film. Also, unlike in the droplet formulation, we now need to consider the deformable surface $g$, which we suppose to have the same scale as the droplet free surface $f$ in order for the surface deformation to have a leading order influence. We then take the following scaling in terms of our previously defined small parameter $\varepsilon$,

The scales (15) lead to asymptotic expansions of the form

where the horizontal velocity scale is expected to be large compared with the vertical velocity scale in order to balance the continuity equation.

When substituting scales (15) and expansions (16) into the governing equations for the gas (4), for $Re\gg 1$ the leading order equations take the form of lubrication flow,

The inertial and acceleration terms do not appear in the left hand side of equations (17a-17b) because these terms are negligible compared to the pressure gradient term. This assumption requires $\rho_{g}/\rho_{l}\ll\varepsilon$. Also, it is assumed that there is a balance in the horizontal momentum equation between the pressure gradient and the viscous terms (Smith, Li, and Wu, 2003), which leads to the definition

This, combined with the requirement $Re\gg 1$ in the liquid, gives the condition that $\varepsilon$ must satisfy in order for our model to be valid, namely

For example, water and air give rise to the range $10^{-3}\ll\varepsilon\ll 0.27$.

The leading order kinematic conditions (6) and (7) now reduce to

and

respectively. The vertical momentum equation (17b) implies that $P_{g}(X,Y,T)=P_{g}(X,0,T)$, then integrating equation (17a) twice in $Y$ and applying conditions (20-21), with normal stress condition (8), gives

the (Reynolds) lubrication equation that helps to link all three unknown quantities $F$, $G$ and $P$.

Finally, we require another equation linking the pressure in the air film with the deformation of the surface. The model we will use here was first described by Carpenter and Garrad (1985) to model surface coatings as a elastic plate (or tensioned membrane) supported above a rigid surface by an array of springs. It is derived from a nonlinear model that includes all relevant physics which is then linearised under the assumption of longitudinal deflections being much smaller than transverse ones. An excellent derivation of this equation can be found in Alexander, Kirk, and Papageorgiou (2020) (we will ignore viscous traction in the air film in our model). This model is also used in a number of other studies on deformable surfaces (Carpenter and Garrad, 1986; Gajjar and Sibanda, 1996; Davies and Carpenter, 1997; Pruessner, 2013; Pruessner and Smith, 2015). The relevant equation takes the non-dimensional form

where the non-dimensional coefficients $e_{i}$ are $(e_{1},e_{2},e_{3},e_{4},e_{5})=(-B^{*}/R^{3}V^{2}\rho_{l},\allowbreak T_{t}^{*}/RV^{2}\rho_{l},\allowbreak-\kappa^{*}R/V^{2}\rho_{l},\allowbreak-M^{*}/R\rho_{l},\allowbreak-C^{*}/V\rho_{l})$ and $p_{s}$ is the non-dimensional relative base pressure, which is taken to be zero. We choose to use this thin-membrane type model over more simpler models as there is commonality in such membranes in nature and practical settings, such as droplet impact onto leaves Gilet and Bourouiba (2015), butterflies Li et al. (2018), umbrellas and raincoats, and relevance in many other applications of droplet impact on deformable surfaces discussed elsewhere in this paper. More simple surface deformation models can be extracted as subsets from equation (24) and one such, a Kelvin-Voigt model of viscoelasticity, is described in detail in III while the equation is considered in full generality in IV. The recent study by Pegg (2019) is likewise based on the Smith, Li, and Wu (2003) theory but with a surface equation generally different from ours; the present work is more focused on soft deformable surfaces which are found to lead to new non-intuitive results. The constants $B^{*}$, $T_{t}^{*}$, $\kappa^{*}$, $M^{*}$ and $C^{*}$ correspond to the flexural rigidity, the longitudinal tension, the stiffness, the mass density and the damping constant of the surface, respectively. In order to apply equation (23) to our scaled liquid droplet application described above, we must apply scales for $g$, $x$, $t$ and $p$ given in (15-16). This gives

where $(\tilde{e}_{1},\tilde{e}_{2},\tilde{e}_{3},\tilde{e}_{4},\tilde{e}_{5})=(\varepsilon^{-1}e_{1},\varepsilon e_{2},\varepsilon^{3}e_{3},\varepsilon^{-1}e_{4},\varepsilon e_{5})$ and $P_{s}=\varepsilon^{-1}p_{s}$. Therefore, the five non-dimensional parameters which control the system are the following

The relative size and importance of each term will be discussed in each case considered in this study. In practise, each of the dimensional structural parameters $B^{*}$, $T_{t}^{*}$, $\kappa^{*}$, $M^{*}$ and $C^{*}$ can vary dramatically depending on the situation, so therefore so can each $e_{i}$ in (25). As droplet impact is a ubiquitous phenomenon, it can occur on many different surfaces of differing elastic properties, for example leaves (Gilet and Bourouiba, 2015), thin foils (Contino et al., 2019), skin (Tagawa et al., 2013) and foodstuffs (Andrade, Skurtys, and Osorio, 2012). It is also possible to use surface engineering to modify the elastic properties of a surface to our advantage (Howland et al., 2016; Weisensee et al., 2016; Vasileiou et al., 2016; Vasileiou, Schutzius, and Poulikakos, 2017; Langley, Castrejon-Pita, and Thoroddsen, 2020). Therefore our focus will be on studying different cases of equation (24), with a very wide range of value of the parameters $e_{i}$.

The deformable surface will be considered to be initially zero and stationary. Boundary conditions will be discussed separately for each case considered. The system to be solved is the non-linear set of governing equations (14), (22) and (24), subject to the boundary condition (13), the pressure $P$ being initially zero and decaying at infinity and the relevant boundary and initial condition on $G$. The governing equations require a numerical treatment. The numerical scheme is slightly different for each case considered and so will be discussed separately for each case.

It is useful now to mention the fluid parameters used in recent experiments. Of most relevance here, we will consider the experiments performed by Langley, Castrejon-Pita, and Thoroddsen (2020) for the pre-impact behaviour of ethanol droplets impacting soft surfaces of varying stiffness. They considered a parameter regime where the Reynolds number $Re$ ranged from 1209 to 20394 and the Weber number $We$ ranged from 17 to 2825. Hence both are typically large and in particular the assumption of a quasi-inviscid liquid droplet seems valid in this regime. Also, in the present study surface tension has been ignored and the scaled surface tension forces are given by

where $\sigma$ is the surface tension coefficient and $We=\rho_{l}V^{2}R/\sigma$ is the Weber number. As the parameter $\varepsilon$ is small and $We$ is typically large, this again seems a valid assumption for the vast majority of the evolution. As the droplet approaches the surface, we do still expect the free-surface to deform to the point where high curvatures are observed, as seen in Smith, Li, and Wu (2003). This would result in a large value for $\nabla^{2}F$, and thus surface tension forces may become significant at this stage. However, for some cases in the following section the free-surface curvature will not become large at any stage and so surface tension effects remain small. The effect of surface tension is considered in Purvis and Smith (2004) for droplet impacts in our parameter regime and in Vanden-Broeck and Smith (2008) for a similar air-cushioning related problem for higher Reynolds numbers than in this work. In the present study, however, as has already been assumed, surface tension effects will be ignored. For all parameters considered in Langley, Castrejon-Pita, and Thoroddsen (2020) the Froude number $Fr=V/\sqrt{gR}$, where $g$ is the gravitational acceleration, is also large; hence gravitational effects are ignored.

The effects of compressibility in the air film are also ignored in the present study. A scaling argument performed by Mandre, Mani, and Brenner (2009), where they balanced the gas pressure gradient with the droplet deceleration, found that compressibility effects in the air can be ignored if the compressiblity parameter $\delta$ satisfies

where $p_{0}^{*}$ is the surrounding ambient pressure. In the study of Langley, Castrejon-Pita, and Thoroddsen (2020) there are impacts considered both in the compressible ($\delta<1$) and incompressible regimes ($\delta\gg 1$), so we will still attempt to draw comparisons with their work later in V. For detailed discussions of compressible gas-cushioned droplet impacts see Mani, Mandre, and Brenner (2010) and Hicks and Purvis (2013) also.

Figure 2 summarises the main model assumptions for a water droplet impact, when $p_{0}^{*}=10^{5}$ Pa.

There is a clear and rather large range of droplet velocities and radii where our model is valid and this is the gray shaded region in figure 2, where the cross shows where the droplet radius is 1 mm and the velocity is 1 m s^-1, a point well within our regime.

Figure 3 shows the evolution of the droplet free-surface and pressure profile for a flat rigid surface ($G=0$) and for a viscoelastic surface for a range of values of $\tilde{e}_{3}=\tilde{e}_{5}$. We chose to vary the parameters as $\tilde{e}_{3}=\tilde{e}_{5}$ because in practise they are dependant on each other, however different combinations are considered later on. The top two panels are the solution for a flat rigid surface which has been previously reported a number of times (Smith, Li, and Wu, 2003; Hicks and Purvis, 2017). As the droplet approaches the surface, there is an initial build-up of pressure directly beneath the droplet. The build-up of pressure acts to decelerate the falling droplet and as the gap between the droplet and the surface narrows, the pressure build-up is large enough to decelerate the droplet free-surface to rest at the center point. At this stage, the droplet begins to deform either side of the center point and eventually overtake it, resulting in an approach to touchdown in two locations. This, in turn, results in the pressure bifurcating away from the center point into two pressure peaks. The process continues until the moment of touchdown on the surface, which results in the entrapment of an air pocket, which may then subsequently form a bubble (van Dam and Le Clerc, 2004; Thoroddsen et al., 2005).

The result for the flat rigid surface in figure 3 is then compared to the result for impact onto a viscoelastic surface for decreasing magnitudes of the surface parameters $\tilde{e}_{3}=\tilde{e}_{5}$. Lowering the magnitude of the surface parameters $\tilde{e}_{3}$ and $\tilde{e}_{5}$ corresponds to lowering the surface stiffness and (viscous) damping, which results in larger surface deformations. Rows two to four in figure 3 show the solutions for the free-surface and the pressure profile evolution, along with the shape of the surface as touchdown is approached, for $\tilde{e}_{3}=\tilde{e}_{5}=-1$, -0.2 and -0.1, respectively. For a water droplet in air of radius 1 mm and impact velocity 1 m s^-1, these parameter variations correspond to a surface spring stiffness in the range $\kappa^{*}\sim 10^{9}$ - $10^{10}$ Pa m^-1 and (viscous) damping of $C^{*}\sim 10^{3}$ - $10^{4}$ Pa s. We note that the spring stiffness values here are larger than the stiffnesses of the soft solids considered in Langley, Castrejon-Pita, and Thoroddsen (2020), however the spring stiffness defined in equation (24) is a stiffness $per~{}unit~{}length$ and the length scales in this problem are typically very small. Smaller values of the dimensionless parameter $\tilde{e}_{3}$ are considered in III.2.

A number of conclusions can be drawn from the influence of surface deformability in this viscoelastic model on the pre-impact dynamics of droplet impact. First of all, there is a profound delay to touchdown. In the absence of air cushioning the droplet would impact on an undisturbed surface at $X=Y=0$ at $T=0$. The presence of air cushioning delays this touchdown into positive time. In figure 3 the solution at $T=0$ is the dash-dotted line, while the positive time solutions are the solid lines. As the droplet approaches the soft viscoelastic surface, the build-up of pressure beneath the droplet acts now to not only deform and decelerate the droplet free surface, but also push the surface away from the droplet. This results in a slower closing of the air gap between the droplet and the soft viscoelastic surface, with lower magnitudes of surface parameters $\tilde{e}_{3}$ and $\tilde{e}_{5}$ resulting in a slower closing of the air gap. In consequence there is a larger delay to touchdown as illustrated by the thick solid lines in figure 3, corresponding to the solution as touchdown is approached, where the time has increased from $T=5.84$ for a flat rigid surface to $T=22.3$ for the soft viscoelastic surface with $\tilde{e}_{3}=\tilde{e}_{5}=-0.1$. The dimensional time scale is given by

For a 1 mm water droplet in air with impact velocity 1 m s^-1, the touchdown delay incurred here is 11.4 $\mu$s.

The slower closing of the air gap also influences the build-up of pressure beneath the droplet. In figure 3 as the parameters $\tilde{e}_{3}$ and $\tilde{e}_{5}$ are decreased in magnitude, the bifurcating behaviour of the pressure profiles is still present. However, the pressure peak amplitude as touchdown is approached is lowered by the introduction of surface deformability, with larger surface deformations resulting in lower pressures.

Figure 4 highlights the link between the rate at which the air gap closes and the maximum pressure. The lower pressures seen for more deformable surfaces are found to be present at all time.

It is also interesting to note the shape of the pressure solution as touchdown is approached. Figure 5 shows the air film thickness $H=F-G$ and pressure profile solution for a range of values $\tilde{e}_{3}=\tilde{e}_{5}$. It can be seen that for surface parameters of sufficiently high magnitude, the pressure profile solution near touchdown is virtually identical to that of the flat surface solution (Smith, Li, and Wu, 2003; Hicks and Purvis, 2017), with a very sharp pressure peak located at the cusp of the air film thickness $H$. As the magnitude of the surface parameters is decreased the sharp pressure peak near touchdown still exists up until around $\tilde{e}_{3}=\tilde{e}_{5}=-0.2$ (this solution is also shown in figure 3) where a rounder pressure peak located just behind the cusp of the air film thickness $H$ is seen. As the magnitude of the surface parameters is decreased further to $\tilde{e}_{3}=\tilde{e}_{5}=-0.1$, the rounder pressure peak overtakes the sharp one near touchdown. The location of the lower, rounder pressure peak just behind the cusp of the air film thickness, as opposed to at it, could help explain the extended air gliding behaviour for droplet impacts on to softer solids seen in experiments by Langley, Castrejon-Pita, and Thoroddsen (2020). Air gliding is when a droplet skids on a thin air film as opposed to touching down; the numerical solutions of the air film thickness in figure 5(a) appear to show the onset of such behaviour.

A vitally important outcome of the dynamics described here is the area of entrapped air underneath the droplet at impact. Trapped air at impact can cause potential problems when it comes to using droplets to spray coat materials or in cooling processes. From the numerical results of pre-impact air-water-substrate interaction in figure 3 it is clear that air entrapment still occurs. What is perhaps more clear in figure 5 is that the presence of substrate deformability in the viscoelastic model leads to an increase in the area of air entrapped at impact. The framework of our model is in two dimensions, and it is realised that droplet impact is clearly a three dimensional phenomenon. Despite this, we are able to use our model to make a qualitative assessment of the size of entrapped air for variations in the substrate stiffness and (viscous) damping.

The X-wise symmetry of the solutions allows us to calculate the dimensionless area of entrapped air $B$ as

where $X_{d}$ is the positive $X$ station of minimum air film thickness and $T_{d}$ is the touchdown time. It should be stressed that, numerically, touchdown where $H=0$ is never quite realised due to the parabolic degeneracy of the lubrication equation (22). Therefore, the time of touchdown has to be pre-determined as a point when the air film thickness reduces to a very small positive value. Numerically, the smaller the grid size and time step in the numerical scheme, the smaller we can make this value. However, this has to be balanced with the increased computational cost due to the huge amount of simulations needed to be run to build a parametric picture. In the present section we therefore set this pre-determined value of air film thickness, where the droplet is considered to have reached touchdown, as $H_{min}=0.1$ (the solutions in figures 3 and 5 are plotted up until $H_{min}=0.2$, for comparison).

Figure 6 shows how the dimensionless area of entrapped air $B$ depends on the surface parameters $\tilde{e}_{3}$ and $\tilde{e}_{5}$, for the viscoelastic surface model. Clearly, for any decrease in magnitude of $\tilde{e}_{3}$ or $\tilde{e}_{5}$ there is an increase in the area of entrapped air. The dimensional entrapped air area $b^{*}$ can be related to the dimensional area $B$ by applying the vertical and horizontal length scales, which yields

where the pre-factor $B$ is, in the viscoelastic model, a function of the parameters $\tilde{e}_{3}$ and $\tilde{e}_{5}$ and its dependency is shown in figure 6 for parameters in the range [-10,-0.1]. For parameter values of magnitude lower than 0.1, the area of entrapped air continues to increase and for parameters of magnitude higher than 10 the behaviour of the droplet is identical to that of a flat rigid surface.

Equation (31) for the dimensional entrapped air area highlights the importance of the droplet velocity and radius on the size of entrapped air at impact. For a flat rigid surface, where $B$ is a given constant, droplets with a larger radius will entrap more air due to having a larger free-surface to interact with the air prior to impact and droplets with a larger velocity will entrap less air due to the droplet having less time to deform prior to impact. In the viscoelastic model, the numerical pre-factor $B$ depends on $\tilde{e}_{3}$ and $\tilde{e}_{5}$ and formulas for these parameters are given in (25). It is interesting to note that the parameter $\tilde{e}_{3}$ does not depend on the droplet radius $R$. It is therefore of interest to see how the size of entrapped air varies with droplet velocity and radius for given fluid and structural parameters (the structural parameters of concern here are the stiffness $\kappa^{*}$ and the damping $C^{*}$).

Figure 7 gives a comparison produced by results for the dimensional entrapped air area of a water droplet in air over wide ranges of the droplet radius $R$ and velocity $V$, for a flat rigid surface and soft viscoelastic surfaces. We chose to vary the droplet velocity from 0.1 m s^-1 to 2 m s^-1 and the radius from $10^{-4}$ m to $10^{-2}$ m in order to keep the parameter regime mostly contained in the valid model parameter regime given in figure 2. Figure 7(a) shows the variations of the entrapped air area $b^{*}$ as the droplet velocity and radius vary, comparing a flat rigid surface to a soft viscoelastic surface with spring stiffness $\kappa^{*}=6$ MPa m^-1 and (viscous) damping $C^{*}=20$ kPa s. For all combinations of droplet velocity and droplet radius, more air is entrapped by the soft viscoelastic surface, as is to be expected. For small droplet impact velocities, the effect of the soft solid is minimal and it behaves much like a rigid surface. As the droplet velocity increases, the effect of the soft viscoelastic surface is far more profound. From equation (31) we can see that for a flat rigid surface the area of entrapped air decreases with increased droplet velocity, whereas for the soft viscoelastic surface this decrease in entrapped air area is delayed and even halted for increased droplet velocity, for the parameters under consideration here. This can be seen more clearly in figure 7(c) where the entrapped air area is given as a function of droplet velocity for a 1 mm water droplet in air, for a variety of spring stiffness $\kappa^{*}$ and (viscous) damping $C^{*}$. This is due to higher air film pressures, induced by higher droplet velocities, causing larger surface deformation prior to impact. By considering the shape of the contour lines in figure 7(a), variations of the droplet radius have a less profound influence on the increase in entrapped air on a soft viscoelastic solid compared to a flat rigid surface. This is shown more clearly in figure 7(b), where the entrapped air area is plotted as a function of droplet radius for an impact velocity of 1 m s^-1. On a logarithmic scale, the increased amount of entrapped air due to impact onto a soft viscoelastic surface is almost independent of droplet radius.

Figure 8 shows the solution profiles for the free-surface $F$ and the pressure $P$ for three different values of the stiffness parameter $\tilde{e}_{3}$, with the other surface parameters fixed. In practice zero stiffness is not possible, hence, while all the other parameters are fixed, $\tilde{e}_{3}\rightarrow 0^{-}$ is a limiting scenario (the same applies for the other parameters in IV.2). The dashed line in each plot of figure 8(a) corresponds to the solution of the flexible surface shape $G$ as touchdown is approached. The results illustrate the outcomes of variations in the stiffness parameter $\tilde{e}_{3}$, however variations in any of the other surface parameters yield qualitatively similar results. The solution in the flat surface case is not shown here for brevity, see figure 3 for comparisons.

The aim of this section is to highlight the effect surface flexibility can have on the pre-impact phase of droplet impact, as opposed to the effect of a soft deformable surface considered in III. Figure 8 shows that much of the same conclusions can be drawn. First, increased surface flexibility leads to a delay to touchdown. This delay to touchdown acts to further decelerate the droplet free surface and results in a lower pressure build-up underneath the droplet. The pressure peak amplitude is again lower near touchdown for more flexible surfaces. Figure 9 shows the minimum air film thickness and the maximum pressure as functions of time, where it can be seen that this lower pressure buildup is present for all time. Figure 10 compares the solutions of the air film thickness and the pressure near touchdown, showing clearly the reduced pressure peak and increased air entrapment near touchdown for reductions in magnitude of the stiffness parameter $\tilde{e}_{3}$.

The main difference to highlight when comparing air cushioning behaviour of droplet impact on a flexible surface to impact on a soft viscoelastic surface is the shape of the pressure profile upon touchdown. Here, in figure 8, reductions in the magnitude of the parameter $\tilde{e}_{3}$ lead to reductions in the amplitude of the pressure peak near touchdown, but the peak remains sharp. This is unlike the finding for the viscoelastic surface, where increased deformability leads to a decrease in the pressure peak amplitude and also a rounder, softer peak. This could be due, in part, to the more abrupt approach to touchdown seen in figure 9, for a flexible surface, compared with the gentler approach seen in figure 4, for a soft viscoelastic surface. Also, the droplet free surface cusp remains relatively sharp in the flexible surface case as deformability increases, while in the viscoelastic case it becomes rounded, which is likely to play a key role in the shape of the pressure peak. However, as the droplet free-surface becomes sharp, surface tension effects will become significant and would need to be considered here in the flexible surface case. At leading order, the local touchdown behaviour here is identical to that described in Smith, Li, and Wu (2003).

Using equation (30) we are again able to make a qualitative assessment of how variations in the surface properties affect the size of entrapped air. Here, we perform a study of how individual variations in the parameters $\tilde{e}_{1}$, $\tilde{e}_{2}$, $\tilde{e}_{3}$ and $\tilde{e}_{4}$ (with $\tilde{e}_{5}=0$) alter the area of entrapped air at touchdown. The size of the air gap at which the calculation of equation (30) is performed is 0.26 for this model, which is the smallest achievable for all compared results and larger than that considered in III due to the more difficult numerical task.

Figure 11 shows how individual variations in the surface parameters $\tilde{e}_{1}$, $\tilde{e}_{2}$, $\tilde{e}_{3}$ and $\tilde{e}_{4}$ affect the entrapped air area at touchdown. Logarithmic variations of each parameter lead to a similar dependency on the entrapped air area. Each figure exhibits plateauing behaviour to the flat rigid surface value for large magnitude parameters and to the value for the solution as the parameter tends to zero, except for figure 11(d), where, in the current setting, variations of the magnitude of the mass density parameter $\tilde{e}_{4}$ much lower than 1 lead to an unbounded surface velocity $\partial G/\partial T$ at some point. From figure 11 it is clear that any increase in flexibility of the substrate leads to an increase in entrapped air.

In Langley, Castrejon-Pita, and Thoroddsen (2020) it was found, experimentally, that droplet impacts onto softer, more deformable, solids would entrap more air. They performed experiments using ultra-high-speed interferometry to capture the droplet free surface at impact, and varied the droplet velocity and the surface stiffness. We have qualitative agreement with these findings in our analytical results. We found in III and IV that reductions in surface stiffness resulted in an increase in entrapped air. A more subtle point to make from our analytical work is the non-intuitive dependency of the size of entrapped air on the droplet velocity. On a rigid flat surface, the size of entrapped air will decrease for higher velocities, whereas for impact onto a soft viscoelastic surface this reduction is delayed and even halted for increased impact velocities. This is alluded to in Langley, Castrejon-Pita, and Thoroddsen (2020) also and, although described in terms of gas compression while our work is incompressible, they show that entrapped air can increase with increased droplet velocity for impact onto soft solids.

Howland et al. (2016) focused their study on how soft deformable surfaces affect droplet splashing. They performed experiments by impacting ethanol droplets on silicone or acrylic substrates of varying stiffnesses, and found that the lower the stiffness, the less likely the droplet was to splash. This was found to be most likely due to higher stiffness substrate having sheet ejection of higher velocity, resulting in the sheet leaving the surface and breaking up, forming a corona splash. By contrast, for less stiff surfaces, the ejected sheet is of lower velocity and does initially leave the substrate, but can fall back down onto the substrate and slow down, suppressing the splash. They were unable to investigate the sheet ejection experimentally, because of the small time and length scales associated with it. Hence, they investigated it using numerical simulation, and found that lowering the surface stiffness reduced the contact pressure just before the sheet is ejected. Our results in section III and IV show clearly that reducing the surface stiffness leads to a reduction in the pressure peak just prior to impact. In particular, the results in III for droplet impact onto a viscoelastic solid show a softening of the pressure peak prior to touchdown. Although our study is solely focused on pre-impact behaviour, this reduction in pressure appears important in understanding the reduced sheet ejection speed mentioned in Howland et al. (2016), which results in potential splash suppression.

The coupled system (37) was solved numerically using a scheme identical to that outlined in III, with the new local pressure-shape relationship (37a). As is to be expected, this scheme is far faster than the one involving the Cauchy-Hilbert integral. An appropriate spatial domain here was found to be $[-20,20]$.

Figure 12 shows the time-marched solutions of ${H}$ and ${P}$. The mechanisms are essentially the same as before. The droplet is released at some suitably large negative time and begins to approach the surface. The air film thickness is initially parabolic and as the air film thickness begins to decrease, the pressure begins to rise in a single peak and results in the air film thickness deviating from the parabolic shape. The pressure then starts to have two peaks as the air film thickness approaches zero in two different locations equidistant from the center line.

The solutions of ${H}$ and ${P}$ in Figure 12 show clearly two traits of droplet impact on a very deformable surface that can be intuitively expected from the results in III and IV. Figure 12 shows solutions up until $T=60$. The longer time scales associated here are clear to see and it is expected that touchdown is not reached in finite time. We can also see that the horizontal bubble extent has significantly increased.

The computational results given in figure 12 show that the time scales involved in this system are large and suggest that perhaps touchdown in finite time is not reached on this time scale. Hence, we seek a solution at large positive time. Suppose that the scaled air thickness and pressure are rising relatively fast in spatial terms, near a region where $X=cT^{\alpha}$ say, with $\alpha$ a positive constant. The length scaling in this region then takes the form

where $\xi$ is order unity, the constants $\alpha$ and $m$ are unknown and in order for the region to be local we require that $m<\alpha$. Expansions of $H$ and $P$ are then taken in the rather general form

where $\lambda$ is an unknown constant. It is to be expected that $\lambda\leq 0$ as $H$ is not growing in time, and the expansion of $P$ is inferred from seeking a balance in equation (37a). Substitution of the expansions (40) into the governing equations (37) leads to the system

to leading order, subject to

and also $3\lambda<2\alpha+1$. The system (41) can be further reduced to an ordinary differential equation in $\hat{H}$ alone,

where $k$ is an integration constant and must be positive in order to keep $\hat{H}$ positive. Equations of the form (43) occur commonly in a number of applied mathematics problems (Kalliadasis and Chang, 1994; Braun and Fitt, 2003; Purvis and Smith, 2004) and it is expected that the solution $\hat{H}$ will tend to a non-zero constant $\alpha ck/12$ as $\xi\rightarrow-\infty$. Perturbations to this constant value occur, such that

where $\Re$ denotes the real part and $\gamma_{1}$ and $\gamma_{2}$ are complex constants. Note, if we were to take the perturbation in the form $\gamma_{1}\exp(\gamma_{2}\xi)$, where $\gamma_{1}$ and $\gamma_{2}$ are real constants, then that would lead to $\gamma_{2}<0$ and thus an exponentially growing perturbation. The complex constant $\gamma_{2}$ satisfies $\gamma_{2}^{3}=-q$, where $q=(12/\alpha ck^{3/4})^{4}$, and the complex constant $\gamma_{1}$ remains arbitrary. There are then three possible solutions for $\gamma_{2}$, and we wish to choose the ones such that $\Re[\gamma_{2}]>0$ so that we have a decaying perturbation. There are two such solutions, $\gamma_{2}=q^{1/3}(1\pm i\sqrt{3})/2$, which leads to the large negative $\xi$ asymptote taking the form

where $\gamma_{1,1}$ and $\gamma_{1,2}$ remain arbitrary, still. For large positive $\xi$ it can be readily shown that

where $\lambda_{1}=(12/\alpha c)^{1/3}$. From equation (41a) we can also show that the corresponding $\hat{P}$ asymptote is

where $\lambda_{2}=\alpha^{2}c^{2}\lambda_{1}$.

Equation (43) was then solved numerically as a boundary value problem using Newton-Raphson iterations, with the asymptotes for $\hat{H}$ (45-46) imposed at suitable large negative and positive $\xi$ values, respectively. The corresponding solution for $\hat{P}$ is then calculated from equation (41a) using central finite differences. Without loss of generality, the constants $\alpha^{2}c^{2}$, $12/(\alpha c)$ and $k$ can be normalised to unity by a division of $\hat{H}$, $\xi$, $\hat{P}$ by $\alpha ck/12$, $(\alpha c/12)^{4/3}k$, $(12^{5}\alpha c)^{1/3}/k$ in turn.

Figure 13 shows the solution of $\hat{H}$ and $\hat{P}$ for $\gamma_{1,1}=\gamma_{1,2}=1$. Variations in the value, even sign, of the parameters $\gamma_{1,1}$ and $\gamma_{1,2}$ yield identical results to that presented in figure 13.

The solutions given in figure 13 show excellent agreement locally with the solutions of the full system (37), solved numerically and presented in figure 12, at large time, especially the solution found for $\hat{P}$ which is qualitatively identical to the solution of $\bar{P}$ at large times (the large time solution to the full system is shown in the inset of figure 13, for comparison).

Now, let us move further rightwards into positive $\xi$ by considering $\xi=\mathcal{D}\bar{X}$, where $\mathcal{D}\gg 1$ and $\bar{X}$ is of order unity. We can infer expansions of $\hat{H}$ and $\hat{P}$ from their large positive $\xi$ asymptotes (46-47),

To leading order, equation (43) then yields $\mathrm{d}^{3}\hat{H}_{0}/\mathrm{d}\bar{X}^{3}=0$ and matching requires $\hat{H}_{0}\sim\lambda_{1}\bar{X}$ as $\bar{X}\rightarrow 0^{+}$. Hence

where $\mu_{1}$ is a positive constant. Now we can see that, at leading order for both $\hat{H}$ and $\hat{P}$,

which is a form resembling the far field condition ($(38\text{a},\text{b})$). In light of (39), (40a) and (49a), $H$ emerges as

just to the right of the local zone, which, as far as the $X^{2}$ requirement in ($(38\text{a},\text{b})$a) is concerned, indicates that the constants must satisfy

hence any $\lambda\leq 0$ works here. To the right of the above local region we have the balance $H_{TT}\sim 0$, which when matching to ($(38\text{a},\text{b})$a) at infinity gives the solution

which in turn matches with the far-field solution.

The most acceptable looking solution to arise from this analysis would be if $\lambda=0$, resulting in $\alpha=1$ and $m=0$, with $\xi=X-cT$ and $H$ and $P$ depending on $\xi$ alone. This nonlinear travelling wave form is described fully by (39-53), capturing the large time features of the full system (37) successfully and confirming the absence of touchdown on this time scale.

The analysis here is a similar one to that seen in Purvis and Smith (2004) and is analogous to a finite-time break up formulation Smith (1988); Smith, Li, and Wu (2003). A physical interpretation of this large time analysis of the large surface deformation system could well be the increased gliding extent observed in Langley, Castrejon-Pita, and Thoroddsen (2020) for droplet impacts onto soft solids. Gliding occurs when the droplet does not make contact with the substrate at the kink of the dimple, and instead glides on a thin air layer. In Langley, Castrejon-Pita, and Thoroddsen (2020) this is seen to occur more frequently and to a greater extent for more deformable surfaces. This could be interpreted as what is occurring in the numerical results for the reduced system in figure 12, with the increased time scales and horizontal extent of the kink being potential traits of droplet gliding in the absence of defects and localised wetting (Langley, Li, and Thoroddsen, 2017).