2 Formulation of the problem

Figure 1.

A side-on view of the problem. A spherical droplet of radius $R$ , travelling with constant speed $V$ , impacts onto a circular elastic plate of radius $L$ (thin horizontal line) housed in a solid housing of otherwise infinite extent (thick horizontal line).

The axisymmetric problem of a liquid droplet impacting onto a partly elastic substrate is formulated in non-dimensional variables. The liquid drop is taken to be spherical before impact, with radius $R$ . The impact speed $V$ is constant. At the instant of impact $(t=0)$ , the drop touches the elastic circular plate at its central point. The circular elastic plate of radius $L$ sits in an otherwise rigid substrate, and is simply supported at its edge (see figure 1). The liquid in the drop is assumed inviscid and incompressible, and the subsequent flow is irrotational and axisymmetric. The effects of the surrounding gas, gravity and surface tension are neglected; these assumptions are more formally justified below. The radius of the elastic plate, $L$ , is assumed to be much smaller than the radius of the drop, $R$ . This condition is motivated by a potential application of the present study to the design of smart surfaces. The ratio

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D716}=L/R\end{eqnarray}$$

is a small parameter in the present study. We consider the frame of reference in which the droplet is stationary and the plate moves up towards it. This is equivalent to considering the droplet impacting a stationary plate.

The early stage of the impact with a strong coupling between the deflection of the elastic plate and the liquid motion is considered here. The conditions of strong coupling are derived below. The dimensional scales are taken to be: $L$ is the length scale, $V$ is the velocity scale and $\unicode[STIX]{x1D70C}V^{2}R/L$ is the pressure scale. The typical time scale $T$ of this stage of the impact is taken as the time required for the elastic plate to be completely wetted. Geometrical considerations yield $T=L^{2}/(VR)$ . The displacements of both the plate and free surface of the liquid drop are of the order of $VT$ . All variables used below are non-dimensional.

The problem is studied by using the cylindrical coordinate system $(r,z)$ . Initially, $t=0$ , the drop is spherical and touches the flat horizontal substrate, $z=0$ , at a single point which is taken as the origin of the coordinate system (see figure 2). The initial surface of the drop is described by the equation

(2.2) $$\begin{eqnarray}z=\frac{\unicode[STIX]{x1D716}}{2}(r^{2}+z^{2}).\end{eqnarray}$$

For mathematical convenience we chose a frame of reference where the droplet is initially at rest and the substrate hits the liquid drop at unit speed from below. The position of the moving substrate is given by

(2.3) $$\begin{eqnarray}z=\unicode[STIX]{x1D716}(t-w(r,t)),\end{eqnarray}$$

where $w(r,t)$ is the deflection of the elastic part, $r<1$ , of the substrate, and $w(r,t)=0$ elsewhere. The deflection $w(r,t)$ is taken to be positive in the negative $z$ -direction.

The flow in the droplet is described by the velocity potential $\unicode[STIX]{x1D711}(r,z,t)$ . The potential satisfies Laplace’s equation in the flow region, $\unicode[STIX]{x1D6FA}(t,\unicode[STIX]{x1D716})$ , and is subject to the kinematic boundary condition,

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{z}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}_{r}\unicode[STIX]{x1D711}_{r}+\unicode[STIX]{x1D702}_{t},\end{eqnarray}$$

on the free surface of the drop, $z=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}(r,t)$ , where $\unicode[STIX]{x1D702}(r,0)=(1/2)r^{2}+O(\unicode[STIX]{x1D716})$ and $r=O(1)$ , and the body boundary condition

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{z}=1-w_{t}-\unicode[STIX]{x1D716}w_{r}\unicode[STIX]{x1D711}_{r},\end{eqnarray}$$

on the moving substrate (2.3), in the wetted region where $r<a(t,\unicode[STIX]{x1D716})$ . This wetted region is unknown in advance and should be found as part of the solution; the unknown function $a(t,\unicode[STIX]{x1D716})$ gives the radius of the contact region (neglecting the thin jet discussed later). The dynamic boundary condition on the free surface of the impacted drop is

(2.6) $$\begin{eqnarray}p(r,\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}(r,t),t)=0,\end{eqnarray}$$

where the hydrodynamic pressure $p(r,z,t)$ is given by the Bernoulli equation,

(2.7) $$\begin{eqnarray}p(r,z,t)=-\unicode[STIX]{x1D711}_{t}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D716}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}|^{2}.\end{eqnarray}$$

The plate deflection $w(r,t)$ is taken to be described by the thin plate theory:

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}w_{tt}+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FB}^{4}w=p(r,\unicode[STIX]{x1D716}(t-w(r,t)),t),\quad (r\leqslant 1),\end{eqnarray}$$

where

(2.9a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\frac{\unicode[STIX]{x1D70C}_{p}h}{\unicode[STIX]{x1D70C}L},\quad \unicode[STIX]{x1D6FD}=\frac{Eh^{3}\unicode[STIX]{x1D716}^{2}}{12(1-\unicode[STIX]{x1D708}^{2})\unicode[STIX]{x1D70C}V^{2}L^{3}},\end{eqnarray}$$

$\unicode[STIX]{x1D70C}_{p}$ is the density of the elastic plate, $h$ is the plate thickness, $E$ is the Young modulus of the plate material, and $\unicode[STIX]{x1D708}$ is the Poisson ratio. The parameter $\unicode[STIX]{x1D6FC}$ indicates the importance of the structural mass per unit area of the plate, $\unicode[STIX]{x1D70C}_{p}h$ , compared to the added mass of the liquid per unit area, which is of the order of $O(\unicode[STIX]{x1D70C}L)$ . The parameter $\unicode[STIX]{x1D6FD}$ can be considered as the dynamic rigidity of the plate. This parameter depends on the elastic characteristics of the plate, speed of impact and the size of the drop. We assume that both $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are of order $O(1)$ in the present analysis, which guarantees strong coupling between hydrodynamic loads and the plate deflection. Outside the wetted region, $r>a(t,\unicode[STIX]{x1D716})$ , the right-hand side in (2.8) is zero. We assume that the elastic plate is simply supported at $r=1$ and flat initially. Then the plate equation (2.8) is to be solved subject to the initial and boundary conditions

(2.10a,b ) $$\begin{eqnarray}\displaystyle w=0,\quad w_{rr}+\frac{\unicode[STIX]{x1D708}}{r}w_{r}=0,\quad (r=1), & & \displaystyle\end{eqnarray}$$

(2.11a,b ) $$\begin{eqnarray}\displaystyle w(r,0)=0,\quad w_{t}(r,0)=0\quad (r\leqslant 1). & & \displaystyle\end{eqnarray}$$

The formulation of the fully coupled problem is completed by the initial conditions

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}(r,z,0)=0, & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle a(0,\unicode[STIX]{x1D716})=0, & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(r,0)={\textstyle \frac{1}{2}}r^{2}+O(\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

where $r=O(1)$ and $z=O(1)$ . The pressure $p(r,z,t)$ in (2.6) and (2.7) does not account for ambient pressure, hydrostatic pressure and the pressure due to surface tension effects. The hydrostatic pressure is of the order of $\unicode[STIX]{x1D70C}gR$ , which is much smaller than the dynamic pressure scale $\unicode[STIX]{x1D70C}V^{2}R/L$ if $\unicode[STIX]{x1D716}(gR/V^{2})\ll 1$ . The surface tension can be neglected at leading order if the inverse Weber number $We^{-1}=\unicode[STIX]{x1D70E}/(\unicode[STIX]{x1D70C}V^{2}R)\ll 1$ , where $\unicode[STIX]{x1D70E}$ is the coefficient of surface tension of the liquid. For a water droplet of density $\unicode[STIX]{x1D70C}=1000~\text{kg}~\text{m}^{-3}$ , surface tension $\unicode[STIX]{x1D70E}=7.2\times 10^{-2}~\text{N}~\text{m}^{-1}$ and radius $R=5\times 10^{-3}~\text{m}$ , with $g=9.81~\text{m}~\text{s}^{-2}$ impacting an elastic circular disk of radius $L=1\times 10^{-3}~\text{m}$ at velocity $V=1~\text{m}~\text{s}^{-1}$ , we have $\unicode[STIX]{x1D70E}/(\unicode[STIX]{x1D70C}V^{2}R)<0.0144$ and $\unicode[STIX]{x1D716}(gR/V^{2})<10^{-2}$ . Hence we neglect surface tension and gravity, as both provide smaller contributions to the hydrodynamic pressure compared to that from inertia.

The problem formulated above is strongly coupled. The plate deflection depends on the hydrodynamic loads through the right-hand side of (2.8), and both the hydrodynamic pressure (2.7) and the flow in the liquid droplet depend on the elastic deflection of the plate through the body boundary condition (2.5). The flow region, $\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D716},t)$ , and the wetted area of the plate, $r<a(t,\unicode[STIX]{x1D716})$ , are unknown in advance and should be determined along with the hydrodynamic and structural characteristics of the problem.

The problem (2.2)–(2.14) can be simplified during the early stage of impact. The approximate solution can be obtained by asymptotic methods as $\unicode[STIX]{x1D716}\rightarrow 0$ . Equations (2.2)–(2.3) show that the flow region can be approximated by the upper half-space, $z>0$ , to leading order, and the boundary conditions (2.4)–(2.6) can be linearised and imposed on the plane $z=0$ . In addition, the linearised dynamic condition (2.6) can be integrated in time using the initial conditions. Below we keep the original notations of the unknown functions for their leading-order terms. At leading order, the hydrodynamic part of the problem reads

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}=0,\quad (z>0), & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}=0,\quad (z=0,r>a(t)), & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{z}=1-w_{t}(r,t),\quad (z=0,r\leqslant a(t)), & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}\rightarrow 0,\quad (r^{2}+z^{2}\rightarrow \infty ), & \displaystyle\end{eqnarray}$$

where $w(r,t)=0$ for $r>1$ . The hydrodynamic pressure is given by the linearised Bernoulli equation, $p(r,z,t)=-\unicode[STIX]{x1D711}_{t}(r,z,t)$ . The plate equation (2.8) at leading order then becomes

(2.19) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D6FC}w_{t}(r,t)+\unicode[STIX]{x1D711}(r,0,t))+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FB}^{4}w=0,\quad (r<1). & & \displaystyle\end{eqnarray}$$

Equation (2.19) is solved subject to the boundary conditions (2.10) and initial conditions (2.11). The shape of the liquid free surface, $z=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}(r,t)$ , is provided at leading order by $\unicode[STIX]{x1D702}_{t}(r,t)=\unicode[STIX]{x1D711}_{z}(r,0,t)$ , where $r>a(t)$ , which follows from the kinematic condition (2.4). Note that $\unicode[STIX]{x1D702}(r,0)\rightarrow r^{2}/2$ as $\unicode[STIX]{x1D716}\rightarrow 0$ . The equations of flow, (2.15)–(2.18), and the equation of the plate deflection (2.19) at leading order are linear but crucially still coupled.

The radius of the contact region between the liquid and substrate, $r<a(t)$ , is determined by using the condition that the vertical coordinate of the free surface, $z=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}(a(t),t)$ at the contact line $r=a(t)$ , and the vertical coordinate of the elastic substrate (2.3) are equal, namely

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D702}(a(t),t)=t-w(a(t),t).\end{eqnarray}$$

Condition (2.20) is known as the Wagner condition (Wagner Reference Wagner1932). This condition assumes that the spray jet formed at the periphery of the contact region is thin and can be neglected at leading order as $\unicode[STIX]{x1D716}\rightarrow 0$ . The Wagner condition makes the leading-order problem nonlinear. As in Korobkin (Reference Korobkin1982), we can reduce the Wagner condition (2.20) to

(2.21) $$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x03C0}/2}\sin (\unicode[STIX]{x1D703})F(a(t)\sin (\unicode[STIX]{x1D703}),t)\,\text{d}\unicode[STIX]{x1D703}=0,\end{eqnarray}$$

where $F(r,t)=t-w(r,t)-r^{2}/2$ is the distance between the initial position of the liquid free surface and the current position of the impacting elastic substrate. See Korobkin (Reference Korobkin1996) and Korobkin & Scolan (Reference Korobkin and Scolan2006) for full details of this approach for related problems. Note that the radius of the contact region $a(t)$ depends strongly on the plate deflection $w(r,t)$ .

The hydrodynamic problem (2.15)–(2.18) is solved by the method of dual integral equations (Sneddon Reference Sneddon1966) and the plate deflection is obtained by the method of normal modes (Korobkin Reference Korobkin1998; Scolan Reference Scolan2004) applied to the structural problem (2.19), (2.10), (2.11), where the radius $a(t)$ of the contact region is determined by equation (2.21).

Figure 2.

(a) A sketch of the problem away from the contact region. The light grey paraboloid shows the liquid free surface $z=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}(r,t)$ , the dark grey circle is the elastic plate with deformation $z=\unicode[STIX]{x1D716}(t-w(r,t))$ and the light grey square region indicates the infinite rigid substrate for the elastic plate. (b) A side-on view of the impact situation for time $t>0$ with plate deflection $w(r,t)$ shown. The substrate moves upwards towards the droplet, which is initially at rest, at a unit speed in the non-dimensional variables.

3 Coupled problem of hydroelastic impact

The formulated problem (2.10)–(2.11), (2.15)–(2.21) is coupled. The hydrodynamic part of the problem (2.15) subject to (2.16)–(2.18), and the structural part of the problem (2.19), subject to (2.10)–(2.11) should be solved simultaneously, together with the Wagner condition (2.21) for the unknown radius of the contact region, $a(t)$ .

3.1 Structural problem

It is convenient to introduce a new unknown function,

(3.1) $$\begin{eqnarray}q(r,t)=-\frac{1}{\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D6FC}w_{t}(r,t)+\unicode[STIX]{x1D711}(r,0,t)),\end{eqnarray}$$

where $0\leqslant r<1$ , and to rewrite equation (2.19) as

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle q_{t}=\unicode[STIX]{x1D6FB}^{4}w, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}w_{t}+\unicode[STIX]{x1D711}(r,0,t)=-\unicode[STIX]{x1D6FD}q(r,t). & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D711}(r,0,t)=0$ , where $a(t)\leqslant r<1$ , during the early stage when the elastic plate is only partly in contact with the liquid, $a(t)<1$ . Within the method of normal modes, the deflection $w(r,t)$ and auxiliary function $q(r,t)$ are sought in the forms

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle w(r,t)=\mathop{\sum }_{n=1}^{\infty }A_{n}(t)w_{n}(r), & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle q(r,t)=\mathop{\sum }_{n=1}^{\infty }k_{n}^{4}q_{n}(t)w_{n}(r), & \displaystyle\end{eqnarray}$$

where $A_{n}(t)$ and $q_{n}(t)$ are the coefficients to be determined. The functions $w_{n}(r)$ are the non-trivial bounded solutions to the homogeneous boundary value problem

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{4}w_{n}=k_{n}^{4}w_{n},\quad r<1, & & \displaystyle\end{eqnarray}$$

(3.7a,b ) $$\begin{eqnarray}\displaystyle w_{n}(1)=0,\quad w_{n}^{\prime \prime }(1)+\unicode[STIX]{x1D708}w_{n}^{\prime }(1)=0, & & \displaystyle\end{eqnarray}$$

and $k_{n}$ are the corresponding eigenvalues. The functions $w_{n}(r)$ describe the axisymmetric shapes of free vibrations of a circular simply supported plate with frequencies proportional to $k_{n}^{2}$ (Leissa Reference Leissa1969). The solutions have the form

(3.8) $$\begin{eqnarray}w_{n}(r)=\text{J}_{0}(k_{n}r)-\frac{\text{J}_{0}(k_{n})}{\text{I}_{0}(k_{n})}\text{I}_{0}(k_{n}r),\end{eqnarray}$$

where $k_{n}$ are the solutions of the equation

(3.9) $$\begin{eqnarray}\frac{\text{J}_{1}(k_{n})}{\text{J}_{0}(k_{n})}+\frac{\text{I}_{1}(k_{n})}{\text{I}_{0}(k_{n})}=\frac{2k_{n}}{1-\unicode[STIX]{x1D708}},\end{eqnarray}$$

with $n\geqslant 1$ and $k_{n+1}>k_{n}$ . Here $\text{J}_{n}(r)$ and $\text{I}_{n}(r)$ are the Bessel functions of the first kind and modified Bessel functions of the first kind of order $n$ respectively. The functions (3.8) are orthogonal but not normalised:

(3.10) $$\begin{eqnarray}\int _{0}^{1}w_{n}(r)w_{m}(r)r\,\text{d}r=U_{n}\unicode[STIX]{x1D6FF}_{nm},\end{eqnarray}$$

where $U_{n}$ are given in appendix A, and where $\unicode[STIX]{x1D6FF}_{nm}=0$ for $n\neq m$ and $\unicode[STIX]{x1D6FF}_{nn}=1$ .

3.2 Coupled problem

The equations (3.3), (3.4) and (3.5) together with initial conditions (2.11) provide the system of ordinary differential equations for the vector $\boldsymbol{q}(t)=\{q_{1}(t),q_{2}(t),\ldots \}$ ,

(3.11a,b ) $$\begin{eqnarray}\dot{\boldsymbol{q}}=\boldsymbol{A}(t),\quad \boldsymbol{q}(0)=0,\end{eqnarray}$$

where $\boldsymbol{A}(t)=\{A_{1}(t),A_{2}(t),\ldots \}$ is the vector of coefficients in the series (3.4), and overdot stands for time derivative. The boundary problem (2.15)–(2.18) for the velocity potential $\unicode[STIX]{x1D711}(r,z,t)$ and the series for the plate deflection (3.4) lead to the following decomposition of the potential:

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D711}(r,z,t)=\unicode[STIX]{x1D711}_{0}(r,z,a)-\mathop{\sum }_{n=1}^{\infty }{\dot{A}}_{n}(t)\unicode[STIX]{x1D711}_{n}(r,z,a).\end{eqnarray}$$

Here $\unicode[STIX]{x1D711}_{0}(r,z,a)$ is the potential of the flow caused by the impact of the rigid circular disk of radius $a$ with the condition on the disc, $\unicode[STIX]{x1D711}_{0,z}(r,0,a)=1$ . The potentials $\unicode[STIX]{x1D711}_{n}(r,z,a)$ satisfy (2.15) and (2.18), with the condition (2.17) becoming $\unicode[STIX]{x1D711}_{n,z}(r,0,a)=w_{n}(r)$ , for $r<a$ . Substituting (3.12) and (3.4) into the second equation of (3.3), multiplying both sides of this equation by $w_{m}(r)r$ , and integrating in $r$ from $0$ to $1$ , we find

(3.13) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FC}{\dot{A}}_{m}U_{m}+\int _{0}^{a}\unicode[STIX]{x1D711}_{0}(r,0,a)w_{m}(r)r\,\text{d}r-\mathop{\sum }_{n=0}^{\infty }{\dot{A}}_{n}(t)\int _{0}^{a}\unicode[STIX]{x1D711}_{n}(r,0,a)w_{m}(r)r\,\text{d}r\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D6FD}k_{m}^{4}q_{m}(t)U_{m},\end{eqnarray}$$

which can be written more concisely in vector form as

(3.14a,b ) $$\begin{eqnarray}\dot{\boldsymbol{A}}=-(\unicode[STIX]{x1D726}+\unicode[STIX]{x1D652}(a))^{-1}(\unicode[STIX]{x1D655}\boldsymbol{q}-\boldsymbol{g}(a)),\quad \boldsymbol{A}(0)=0.\end{eqnarray}$$

Here $\unicode[STIX]{x1D726}=\text{diag}\{\unicode[STIX]{x1D6FC}U_{1},\unicode[STIX]{x1D6FC}U_{2},\ldots \}$ and $\unicode[STIX]{x1D655}=\text{diag}\{\unicode[STIX]{x1D6FD}k_{1}^{4}U_{1},\unicode[STIX]{x1D6FD}k_{2}^{4}U_{2},\ldots \}$ are diagonal matrices. The elements of the vector $\boldsymbol{g}(a)$ are

(3.15) $$\begin{eqnarray}g_{n}(a)=-\int _{0}^{a}\unicode[STIX]{x1D711}_{0}(r,0,a)w_{n}(r)r\,\text{d}r,\end{eqnarray}$$

and the symmetric added mass matrix $\unicode[STIX]{x1D652}(a)$ has the elements

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D61E}_{nm}(a)=-\int _{0}^{a}\unicode[STIX]{x1D711}_{n}(r,0,a)w_{m}(r)r\,\text{d}r.\end{eqnarray}$$

The system of ordinary differential equations (3.11) and (3.14) is valid for $a\leqslant 1$ . For the later stages, when $a(t)>1$ and the elastic plate is completely wetted, we should change the upper limits in (3.15) and (3.16) to $1$ and set $w_{n}(r)=0$ , where $r>1$ in the boundary conditions for the potentials $\unicode[STIX]{x1D711}_{n}(r,z,a)$ . The coefficients in (3.14) depend on the radius of the wetted area, $a(t)$ , which is defined by equation (2.21). Substituting for $F(r,t)$ and the series (3.4) in (2.21), we find

(3.17) $$\begin{eqnarray}t-\frac{1}{2}a^{2}(t)\int _{0}^{\unicode[STIX]{x03C0}/2}\sin ^{3}(\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}-\mathop{\sum }_{n=1}^{\infty }A_{n}(t)\int _{0}^{\unicode[STIX]{x03C0}/2}\sin (\unicode[STIX]{x1D703})w_{n}(a\sin (\unicode[STIX]{x1D703}))\,\text{d}\unicode[STIX]{x1D703}=0,\end{eqnarray}$$

and then

(3.18) $$\begin{eqnarray}a^{2}(t)=3t-3\boldsymbol{A}(t)\boldsymbol{\cdot }\boldsymbol{Q}(a),\end{eqnarray}$$

where

(3.19) $$\begin{eqnarray}Q_{n}(a)=\int _{0}^{\unicode[STIX]{x03C0}/2}\sin (\unicode[STIX]{x1D703})w_{n}(a\sin (\unicode[STIX]{x1D703}))\,\text{d}\unicode[STIX]{x1D703}=\frac{1}{a}\int _{0}^{a}\frac{\unicode[STIX]{x1D70C}w_{n}(\unicode[STIX]{x1D70C})\,\text{d}\unicode[STIX]{x1D70C}}{\sqrt{a^{2}-\unicode[STIX]{x1D70C}^{2}}}.\end{eqnarray}$$

Note that $w_{n}(a\sin (\unicode[STIX]{x1D703}))\equiv 0$ , where $a\sin (\unicode[STIX]{x1D703})>1$ . The system (3.11), (3.14) and (3.18) can be solved numerically once the integrals (3.15), (3.16) and (3.19) are known functions of $a$ . The functions $Q_{n}(a)$ and $g_{n}(a)$ are evaluated in appendices B and C, respectively.

3.3 Hydrodynamic problem

The hydrodynamic problem (2.15)–(2.18), together with the decomposition (3.12), leads to two mixed boundary value problems. The rigid part of the decomposition (3.12) becomes

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}_{0}=0,\quad (z>0), & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{0}=0,\quad (z=0,r>a(t)), & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{0,z}=1,\quad (z=0,r\leqslant a(t)), & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{0}\rightarrow 0,\quad (r^{2}+z^{2}\rightarrow \infty ), & \displaystyle\end{eqnarray}$$

whilst the elastic part of the decomposition requires

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}_{n}=0,\quad (z>0), & \displaystyle\end{eqnarray}$$

(3.25) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{n}=0,\quad (z=0,r>a(t)), & \displaystyle\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{n,z}=w_{n}(r),\quad (z=0,r\leqslant a(t)), & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{n}\rightarrow 0,\quad (r^{2}+z^{2}\rightarrow \infty ), & \displaystyle\end{eqnarray}$$

where $n\geqslant 1$ .

These mixed boundary value problems are solved through the use of Hankel transformations. Equations (3.20) and (3.23) give

(3.28) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{0}(r,z,a)=\int _{0}^{\infty }G_{0}(\unicode[STIX]{x1D706},a)e^{-\unicode[STIX]{x1D706}z}\text{J}_{0}(\unicode[STIX]{x1D706}r)\,\text{d}\unicode[STIX]{x1D706},\end{eqnarray}$$

where $G_{0}(\unicode[STIX]{x1D706},a)$ is a new unknown function. Substituting (3.28) in the boundary conditions (3.21) and (3.22) gives us the dual integral equations

(3.29) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\infty }G_{0}(\unicode[STIX]{x1D706},a)\text{J}_{0}(\unicode[STIX]{x1D706}r)\,\text{d}\unicode[STIX]{x1D706}=0,\quad r>a & \displaystyle\end{eqnarray}$$

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{\infty }G_{0}(\unicode[STIX]{x1D706},a)\unicode[STIX]{x1D706}\text{J}_{0}(\unicode[STIX]{x1D706}r)\,\text{d}\unicode[STIX]{x1D706}=-1,\quad r\leqslant a. & \displaystyle\end{eqnarray}$$

These equations were solved by Sneddon (Reference Sneddon1966), yielding

(3.31) $$\begin{eqnarray}G_{0}(\unicode[STIX]{x1D706},a)=-\sqrt{\frac{2\unicode[STIX]{x1D706}}{\unicode[STIX]{x03C0}}}\int _{0}^{a}\sqrt{x}\text{J}_{1/2}(x\unicode[STIX]{x1D706})\,\text{d}x\int _{0}^{x}\frac{\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}}{\sqrt{x^{2}-\unicode[STIX]{x1D70C}^{2}}}.\end{eqnarray}$$

Equations (3.31) and (3.28) provide

(3.32) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{0}(r,0,a)=-\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\int _{0}^{a}\sqrt{x}\,\text{d}x\int _{0}^{x}\frac{\unicode[STIX]{x1D70C}\,\text{d}\unicode[STIX]{x1D70C}}{\sqrt{x^{2}-\unicode[STIX]{x1D70C}^{2}}}\int _{0}^{\infty }\sqrt{\unicode[STIX]{x1D706}}\text{J}_{0}(\unicode[STIX]{x1D706}r)\text{J}_{1/2}(\unicode[STIX]{x1D706}x)\,\text{d}\unicode[STIX]{x1D706}.\end{eqnarray}$$

This can be integrated directly (with the $\unicode[STIX]{x1D706}$ integral given by 6.575.1 of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007)) to give an explicit expression for the rigid velocity potential,

(3.33) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}_{0}(r,0,a)=-\frac{2}{\unicode[STIX]{x03C0}}\sqrt{a^{2}-r^{2}},\quad (r\leqslant a). & & \displaystyle\end{eqnarray}$$

The elastic part of the problem (3.24)–(3.27) is solved using the same method. The elastic terms of the velocity potential, $\unicode[STIX]{x1D711}_{n}(r,0,a)$ , are given by

(3.34) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{n}(r,0,a)=-\frac{2}{\unicode[STIX]{x03C0}}\int _{r}^{a}\frac{\text{d}x}{\sqrt{x^{2}-r^{2}}}\int _{0}^{x}\frac{\unicode[STIX]{x1D70C}w_{n}(\unicode[STIX]{x1D70C})\,\text{d}\unicode[STIX]{x1D70C}}{\sqrt{x^{2}-\unicode[STIX]{x1D70C}^{2}}}.\end{eqnarray}$$

By using (3.19) we can express (3.34) as

(3.35) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{n}(r,0,a)=-\frac{2}{\unicode[STIX]{x03C0}}\int _{r}^{a}\frac{xQ_{n}(x)\,\text{d}x}{\sqrt{x^{2}-r^{2}}}=-\unicode[STIX]{x1D6F7}_{n}(r,a),\end{eqnarray}$$

and (3.12) as

(3.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}(r,0,a)=-\frac{2}{\unicode[STIX]{x03C0}}\sqrt{a^{2}-r^{2}}+\mathop{\sum }_{n=1}^{\infty }{\dot{A}}_{n}(t)\unicode[STIX]{x1D6F7}_{n}(r,a),\quad (r\leqslant a). & & \displaystyle\end{eqnarray}$$

The added-mass elements (3.16) then become

(3.37) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D61E}_{nm}=\frac{2}{\unicode[STIX]{x03C0}}\int _{0}^{a}\left(\int _{r}^{a}\frac{xQ_{n}(x)\,\text{d}x}{\sqrt{x^{2}-r^{2}}}\right)w_{m}(r)\,\text{d}r=\frac{2}{\unicode[STIX]{x03C0}}\int _{0}^{a}x^{2}Q_{n}(x)Q_{m}(x)\,\text{d}x. & & \displaystyle\end{eqnarray}$$

The integrals (3.37) are evaluated analytically in appendix D. Integrating in (3.35) by parts, we obtain

(3.38) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{n}(r,a)=\frac{2}{\unicode[STIX]{x03C0}}Q_{n}(a)\sqrt{a^{2}-r^{2}}-\frac{2}{\unicode[STIX]{x03C0}}\int _{r}^{a}\sqrt{x^{2}-r^{2}}Q_{n}^{\prime }(x)\,\text{d}x,\end{eqnarray}$$

where $Q_{n}(a)$ is given by (3.19).

Equations (3.38) and (3.36) provide the behaviour of the velocity potential close to the contact line as $r\rightarrow a(t)$ ,

(3.39) $$\begin{eqnarray}\unicode[STIX]{x1D711}(r,0,t)=B(t)\sqrt{a(t)-r}+O((a-r)^{3/2}),\end{eqnarray}$$

where

(3.40) $$\begin{eqnarray}B(t)=\frac{2}{\unicode[STIX]{x03C0}}\sqrt{2a}\left(\mathop{\sum }_{n=1}^{\infty }{\dot{A}}_{n}(t)Q_{n}(a)-1\right).\end{eqnarray}$$

Therefore the Wagner model predicts a square-root singularity in the radial velocity of the flow in the main region, $(\unicode[STIX]{x2202}\unicode[STIX]{x1D711}/\unicode[STIX]{x2202}r)(r,0,t)$ , and in the hydrodynamic pressure, $p(r,0,t)=-\unicode[STIX]{x2202}\unicode[STIX]{x1D711}/\unicode[STIX]{x2202}t$ , at the contact line as $r\rightarrow a(t)$ . The coefficient of this singularity, $B(t)$ , depends on the plate deflection and is calculated as part of the solution in the main flow region. The Wagner model is not valid at the periphery of the contact region. In Wagner-type problems, this singularity is resolved by introducing an inner region, the so-called jet-root region, around the contact line (Howison et al. Reference Howison, Ockendon and Wilson1991). This locally resolves the singularity by the presence of a thin jet running along the substrate. This local inner solution is discussed below.

4 Jet-root and jet regions

The dynamics of the jet and its dependence on elastic oscillations of the plate play a key role in potential splashing of the droplet. We distinguish the jet-root region, which is in a small vicinity of the contact line and is characterised by large velocities and pressures, and the jet itself. The equations governing the flows in both the jet-root region and the jet region are derived by using stretched variables. For the jet-root region $r=a(t)+\unicode[STIX]{x1D716}^{2}X$ , $z=\unicode[STIX]{x1D716}(t-w(a,t))+\unicode[STIX]{x1D716}^{2}Z$ and $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D716}(\unicode[STIX]{x1D719}(X,Z,t)+{\dot{a}}X)$ are used in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ . Howison et al. (Reference Howison, Ockendon and Wilson1991) and Oliver (Reference Oliver2002) contain full details of the case for impact with a rigid substrate. We do not provide details of the asymptotic analysis here, and instead just highlight how the elastic substrate affects the solution.

In the jet-root region the singularity found in the pressure above is resolved by considering a locally two-dimensional flow problem around the contact line $r=a(t)$ subject to the far-field condition which follows from (3.39) as $X^{2}+Y^{2}\rightarrow \infty$ . Locally the jet-root flow is the same as in the problem of a rapid rigid body impact onto a water free surface. The latter problem was studied by Oliver (Reference Oliver2002). In our situation the key difference is in the $B(t)$ term of the far-field condition, which contains additional elastic terms: see equation (3.40). The important results are the initial thickness of the jet,

(4.1) $$\begin{eqnarray}H_{j}(t)=\frac{\unicode[STIX]{x03C0}}{16}\frac{B^{2}}{{\dot{a}}^{2}},\end{eqnarray}$$

and the velocity of fluid entering the jet $2{\dot{a}}(t)$ , in the static frame of reference.

Note that without the elastic deflection of the plate, $A_{j}(t)=0$ in (3.18) and (3.40), the Wagner condition requires $a(t)=\sqrt{3t}$ , giving $B=-2\sqrt{2a}/\unicode[STIX]{x03C0}$ in the non-dimensional variables. The jet thickness (4.1) is then given by $2t^{3/2}/(\sqrt{3}\unicode[STIX]{x03C0})$ . In dimensional variables, the jet thickness at the entrance to the jet is  $2R(Vt_{d}/R)^{3/2}/(\sqrt{3}\unicode[STIX]{x03C0})$ , where $t_{d}$ is the dimensional time and $Vt_{d}$ is the displacement of the droplet. The obtained results are in agreement with the results by Korobkin & Scolan (Reference Korobkin and Scolan2003) in the problem of an elliptic paraboloid impact onto the flat water surface. Note that the uniform velocity in the jet, $2{\dot{a}}(t)$ , and the constant jet thickness, $H_{j}$ , are approached exponentially quickly as $X$ increases. In particular, $\unicode[STIX]{x1D719}_{X}(X,0,t)={\dot{a}}(t)+O(\exp [-\unicode[STIX]{x03C0}X/2H_{j}])$ as $X\rightarrow +\infty$ along the plate. The pressure in the jet-root region is much greater than in the main flow region, with $p={\dot{a}}^{2}/(2\unicode[STIX]{x1D716})$ at the stagnation point. However, the hydrodynamic pressure approaches zero in the jet exponentially quickly as $X\rightarrow +\infty$ . With the elastic deformation of the plate, the jet thickness and jet velocity can be found from (4.1) once the coupled hydrodynamic-elastic problem described above is solved.

Although the curvature in the jet-root region becomes higher than in the main flow region, we still neglect the effect of surface tension in the current paper. Using the jet-root scales we find a Weber number local to the jet root given by $\unicode[STIX]{x1D70C}V^{2}R\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70E}$ , which compares with $\unicode[STIX]{x1D70C}V^{2}R/\unicode[STIX]{x1D70E}\unicode[STIX]{x1D716}$ in the main flow region considered in § 3. We note here that since the Weber number in the jet-root region is two orders of magnitude smaller, the overall effect of curvature (proportional to $1/We$ ) can become large and cause surface tension to be a leading-order effect. For instance, when considering the impact of a water droplet with a $0.5$ mm radius elastic plate at $1~\text{m}~\text{s}^{-1}$ , we find $\unicode[STIX]{x1D70C}V^{2}R\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70E}\approx 7$ . More generally, for a typical range of water droplet sizes and impact velocities the local Weber number in the jet root ranges from approximately $1$ to $450$ . Solving the jet root with strong surface tension included is a challenging and interesting problem in its own right, but is not the focus of this work. We restrict ourselves to situations where the surface tension may be considered a higher-order correction.

As the flow leaves the jet-root region it enters a long thin jet detailed here. Motivated by matching to the jet-root region, the axisymmetric flow in the jet is described in cylindrical coordinates by the shallow-water equations,

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h_{0}}{\unicode[STIX]{x2202}t}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(ruh_{0})=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}r}=0, & \displaystyle\end{eqnarray}$$

where $r>a(t)$ . Here we have introduced the leading-order jet thickness $h(r,t)$ and radial velocity $u(r,t)$ . Matching onto the jet-root region provides boundary conditions

(4.4a,b ) $$\begin{eqnarray}u(a,t)=2{\dot{a}}(t),\quad h(a,t)=H_{j}(t).\end{eqnarray}$$

We note here that the jet thickness in the dimensional variables is $L^{3}h(r,t)/R^{2}$ . Details of these equations can be found in Howison et al. (Reference Howison, Ockendon and Wilson1991) for two-dimensional problems, in Korobkin (Reference Korobkin1997) for axisymmetric problems, and in Oliver (Reference Oliver2002) for three-dimensional problems. See also Moore & Oliver (Reference Moore and Oliver2014) and Moore et al. (Reference Moore, Howison, Ockendon and Oliver2012) for more details and the case of an oblique impact. The shallow-water problem (4.2)–(4.4) is solved using the method of characteristics. We introduce the Lagrangian coordinate $\unicode[STIX]{x1D70F}$ corresponding to the time at which a given material particle enters the jet. The position of this particle at a later time $t$ is given by

(4.5) $$\begin{eqnarray}r=a(\unicode[STIX]{x1D70F})+2{\dot{a}}(\unicode[STIX]{x1D70F})(t-\unicode[STIX]{x1D70F}),\end{eqnarray}$$

and the jet thickness by

(4.6) $$\begin{eqnarray}h(r,t)=H_{j}(\unicode[STIX]{x1D70F})\left(1-\frac{2\ddot{a}(\unicode[STIX]{x1D70F})(t-\unicode[STIX]{x1D70F})}{{\dot{a}}(\unicode[STIX]{x1D70F})}\right)^{-1}\left(1+\frac{2{\dot{a}}(\unicode[STIX]{x1D70F})}{a(\unicode[STIX]{x1D70F})}(t-\unicode[STIX]{x1D70F})\right)^{-1}.\end{eqnarray}$$

This jet thickness has been derived under the assumption that the flow in the main region is described by the Wagner model of liquid impact. The Wagner model requires, in particular, that the contact region expands in time, namely that ${\dot{a}}(t)>0$ . Therefore $1+2{\dot{a}}(\unicode[STIX]{x1D70F})(t-\unicode[STIX]{x1D70F})/a(\unicode[STIX]{x1D70F})\geqslant 1$ in (4.6). The last term in (4.6) describes a decrease of the jet thickness as the particles spread further from the contact line in order to conserve mass. We note that as in (for example) Howison et al. (Reference Howison, Ockendon and Wilson1991), Oliver (Reference Oliver2002) and Howison et al. (Reference Howison, Ockendon, Oliver, Purvis and Smith2005), at leading order, equation (4.5) predicts that the splash jet extends to infinity. This infinite jet is only present for the impact of solid and liquid bodies bounded by a smooth surface before impact where ${\dot{a}}(t)\rightarrow \infty$ as $t\rightarrow 0$ . For non-smooth bodies, such as the water entry of a rigid cone at constant speed, the jet is finite in length.

It can be observed from (4.6) that the jet thickness $h(r,t)$ grows beyond all bounds if the acceleration of the contact line $\ddot{a}(\unicode[STIX]{x1D70F})$ is positive as $(2\ddot{a}(\unicode[STIX]{x1D70F})(t-\unicode[STIX]{x1D70F}))/{\dot{a}}(\unicode[STIX]{x1D70F})\rightarrow 1$ for some $t$ . This predicted unbounded growth in $h(r,t)$ violates the assumptions of our shallow-water model, as the vertical extent of the jet grows beyond the $O(\unicode[STIX]{x1D716}^{2})$ thickness we have assumed. We note that the vertical velocity of the jet at the point of unbounded growth is more singular than its height. As the critical time at which the shallow-water model breaks down is approached, the jet vertical velocity can become the same order of magnitude as the jet horizontal velocity, indicating a splash. Note that the contact line acceleration is strictly negative for a rigid substrate, where $a(t)=\sqrt{3t}$ (3.18). In this case $h(r,t)$ in (4.6) is always bounded. With the coupled elastic plate however, oscillations appear in the contact line acceleration $\ddot{a}(t)$ and it can become positive, leading to splashing being prompted in the jetting region.

5 Numerical solution

The system (3.11), (3.14) and (3.18) which describes the response of the elastic plate and the flow in the main flow region is solved numerically. Note that the right-hand sides of (3.11) and (3.14) depend on the unknown vector functions $\boldsymbol{q}$ , $\boldsymbol{A}$ and the radius of the contact region $a(t)$ , where $\text{d}a/\text{d}t>0$ during the impact stage. It is convenient numerically to take the radius $a$ as a new time-like variable and consider the time-varying functions vectors $\boldsymbol{q}$ , $\boldsymbol{A}$ as functions of $a$ . Differentiating the Wagner condition (3.18) with respect to time, we obtain

(5.1) $$\begin{eqnarray}\frac{\text{d}t}{\text{d}a}=\frac{\frac{2}{3}a+\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{Q}^{\prime }(a)}{1-\dot{\boldsymbol{A}}\boldsymbol{\cdot }\boldsymbol{Q}(a)}=G(\boldsymbol{A},\boldsymbol{q},a),\end{eqnarray}$$

where $\dot{\boldsymbol{A}}$ is given by (3.14). Note that $G(\boldsymbol{A},\boldsymbol{q},a)$ does not depend on the unknown function $t(a)$ . Next we multiply both sides of equations (3.11) and (3.14) by $\text{d}t/\text{d}a$ and use (5.1) to derive the equations for $\boldsymbol{A}(a)$ and $\boldsymbol{q}(a)$ :

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{A}}{\text{d}a}=-(\unicode[STIX]{x1D726}+\unicode[STIX]{x1D652}(a))^{-1}(\unicode[STIX]{x1D655}\boldsymbol{q}(a)-\boldsymbol{g}(a))G(\boldsymbol{A},\boldsymbol{q},a), & \displaystyle\end{eqnarray}$$

(5.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{q}}{\text{d}a}=\boldsymbol{A}(a)G(\boldsymbol{A},\boldsymbol{q},a). & \displaystyle\end{eqnarray}$$

The initial conditions for the system of ordinary differential equations (5.1)–(5.3) are

(5.4a,b ) $$\begin{eqnarray}\boldsymbol{A}=0,\quad \boldsymbol{q}=0,\quad \text{at }t=0\text{ or }a=0.\end{eqnarray}$$

The initial value problem (5.1)–(5.4) is regular at $a=0$ . In particular, $t=O(a^{2})$ . The new time-like variable allows us to avoid difficulties with starting simulations from $t=0$ , where $\text{d}a/\text{d}t\rightarrow \infty$ as $t\rightarrow 0$ (see for example Korobkin & Khabakhpasheva (Reference Korobkin and Khabakhpasheva2006) for more details). The system of ordinary differential equations (5.1)–(5.3) is truncated to $N$ modes and integrated by using the fourth-order Runge–Kutta scheme. In total $2N+1$ equations are integrated. The number of retained modes is not an issue in the present problem because all elements of the system (5.1)–(5.3) are given by analytic formulae. The step of integration, $\unicode[STIX]{x1D6FF}a$ , is related to the time step, $\unicode[STIX]{x1D6FF}t$ , by equation (5.1), $\unicode[STIX]{x1D6FF}t=G\unicode[STIX]{x1D6FF}a$ , where $G=O(a)$ for small time. To properly capture the highest $N$ th mode, the time step $\unicode[STIX]{x1D6FF}t$ is taken to be smaller than the non-dimensional period of the $N$ th mode. The natural frequency of the $N$ th mode is equal to $(\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC})^{1/2}k_{N}^{2}$ , which follows from the plate equation (2.8) and equation (3.6). The time step $\unicode[STIX]{x1D6FF}t$ is at least six times smaller than the period of the $N$ th mode if

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}t<(\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD})^{1/2}k_{N}^{-2}.\end{eqnarray}$$

Correspondingly the step $\unicode[STIX]{x1D6FF}a$ should be smaller than

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}a<(\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD})^{1/2}k_{N}^{-2}G^{-1}(\boldsymbol{A},\boldsymbol{q},a).\end{eqnarray}$$

It is seen that the upper limit of the integration step $\unicode[STIX]{x1D6FF}a$ depends on the solution through the speed of the contact radius $\text{d}a/\text{d}t=G^{-1}(\boldsymbol{A},\boldsymbol{q},a)$ . Initially the speed is very high, which makes it possible to use a relatively large step $\unicode[STIX]{x1D6FF}a$ . However, as the speed $\text{d}a/\text{d}t$ becomes small the right-hand side in (5.6) decreases, and this inequality requires a very small step of integration. In such conditions we swap between numerical stepping in contact line radius, $a$ , and stepping in time, $t$ , depending on the contact line velocity $\text{d}a/\text{d}t$ . We make this transition at approximately $\text{d}a/\text{d}t=2.5$ . If the stepping in time is used, the radius $a(t)$ is obtained from the Wagner condition, (3.18), by a root-finding method.

It was confirmed numerically that the solution converges with increasing number of modes and that the number of modes used in our calculations captures the physics of the full system with sufficient accuracy. We confirmed the convergence with number of modes by analysing in particular the behaviour of the added-mass matrix, $\unicode[STIX]{x1D652}(a)$ . It was found that as the number of modes increases the elements of the added-mass matrix rapidly decay, as will be seen in appendix D, figure 21. The convergence with number of modes was found to be fairly good – with ten modes providing very good accuracy during the impact stage with a partially wetted plate. Later, when the elastic plate is fully wetted, the method is found to have substantially better convergence – with five modes providing an adequate long-term representation of the flow and plate deflection. This stems from the elastic effects (and hence the added-mass contribution) becoming much less important compared with the rigid terms during the stage with a fully wetted plate.

Once (5.1)–(5.4) have been resolved the jet thickness (4.6) yields conditions providing the time splashing is observed. For negative acceleration $\ddot{a}(\unicode[STIX]{x1D70F})$ , we have

(5.7) $$\begin{eqnarray}1-\frac{2\ddot{a}(\unicode[STIX]{x1D70F})(t-\unicode[STIX]{x1D70F})}{{\dot{a}}(\unicode[STIX]{x1D70F})}=1+\frac{2|\ddot{a}(\unicode[STIX]{x1D70F})|(t-\unicode[STIX]{x1D70F})}{{\dot{a}}(\unicode[STIX]{x1D70F})}\geqslant 1,\end{eqnarray}$$

where $t\geqslant \unicode[STIX]{x1D70F}$ . If $\ddot{a}(\unicode[STIX]{x1D70F})>0$ for a certain $\unicode[STIX]{x1D70F}$ , then $h(r,t)\rightarrow \infty$ as $t\rightarrow t_{c}(\unicode[STIX]{x1D70F})$ , where $t_{c}(\unicode[STIX]{x1D70F})$ is determined by the equation

(5.8) $$\begin{eqnarray}\frac{2\ddot{a}(\unicode[STIX]{x1D70F})}{{\dot{a}}(\unicode[STIX]{x1D70F})}(t_{c}(\unicode[STIX]{x1D70F})-\unicode[STIX]{x1D70F})=1,\end{eqnarray}$$

and $t_{c}(\unicode[STIX]{x1D70F})>\unicode[STIX]{x1D70F}$ . For a given function $a(\unicode[STIX]{x1D70F})$ , to determine the position and time of formation of the splash we need to find the minimum value, $t_{min}$ , of the function $t_{c}(\unicode[STIX]{x1D70F})$ . To calculate $t_{min}$ and the radius $r_{s}$ at which the splashing occurs, we find the maxima of

(5.9) $$\begin{eqnarray}\frac{1}{t_{c}(\unicode[STIX]{x1D70F})}=\frac{2\ddot{a}(\unicode[STIX]{x1D70F})}{{\dot{a}}(\unicode[STIX]{x1D70F})+2\unicode[STIX]{x1D70F}\ddot{a}(\unicode[STIX]{x1D70F})},\end{eqnarray}$$

with the constraint that

(5.10) $$\begin{eqnarray}0<\frac{1}{t_{c}(\unicode[STIX]{x1D70F})}<\frac{1}{\unicode[STIX]{x1D70F}}.\end{eqnarray}$$

This approach provides a graphical way to find the time $t_{min}$ and the corresponding value of $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{min}$ , such that $t_{c}(\unicode[STIX]{x1D70F}_{min})=t_{min}$ . The radius $r_{s}$ at which splashing occurs then follows from (4.5):

(5.11) $$\begin{eqnarray}r_{s}=a(\unicode[STIX]{x1D70F}_{min})+\frac{{\dot{a}}^{2}(\unicode[STIX]{x1D70F}_{min})}{\ddot{a}(\unicode[STIX]{x1D70F}_{min})}.\end{eqnarray}$$

Although the elastic plate motion additionally affects fine details of the jet behaviour through the initial jet thickness $H_{j}$ (see (3.39) and (4.1)), its primary effect is to be responsible for triggering blow-up of the jet through the varying acceleration of the contact line.

Figure 3.

Parametric plane showing the regions of negative and unbounded contact line velocity. Pluses represent simulations which encounter contact line shrinkage, crosses represent simulations with no problems, and filled circles represent simulations with an unbounded contact line velocity.

Although the method of normal modes provides a good basis to investigate the velocity field in the droplet, plate deflection and elastic stresses, it does not allow us to accurately determine the pressure distribution over the wetted region. The pressure singularity at the contact line, coupled with poor convergence of the modal series for the plate acceleration $w_{tt}(r,t)$ , make a series representation of the pressure impractical, as discussed by Korobkin (Reference Korobkin1998). However, our main goal in this paper is to study the effect of elasticity on the spreading and splashing of a droplet. We assume that the pressure in the wetted part of the substrate does not go down to the vapour pressure due to the elastic plate vibration. If so, cavitation does not occur, and both the spreading and splashing of the droplet are not affected by the details of the pressure distribution.