2.1. Flow in the lamella

The horizontal flow far from the substrate is denoted by $u_{\textit{outer}}$ (the outer flow). Similar to the 3-D axisymmetric case, this is assumed to be a hyperbolic flow and is given (in 2-D Cartesian coordinates) by the expression

(2.1) \begin{align} u_{\textit{outer}}(x,t)=\frac {x}{t+t_0}, \end{align}

where $t_0$ is a model parameter which describes the onset time of the flow, prior to the formation of the RL structure. The assumption of this hyperbolic flow profile in the lamella is justified by earlier simulations of 2-D droplet impact (Ó Náraigh & Mairal Reference Ó Náraigh and Mairal2023).

A boundary layer forms close to the substrate, in which the fluid velocity transitions from $u_{\textit{outer}}$ to zero, across a distance $h_{bl}$ , the boundary-layer thickness. From boundary-layer theory (White et al. Reference Frank2011), the boundary-layer thickness scales as $h_{bl}\propto \sqrt {\nu x/u_{\textit{outer}}}$ , where $\nu =\mu /\rho$ is the kinematic viscosity of the liquid. Given the functional form of $u_{\textit{outer}}$ , the $x$ -dependence cancels out, to give $h_{bl}\propto \sqrt {\nu (t+t_0)}$ . We allow for an additional degree of freedom in the problem by taking $h_{bl}\propto \sqrt {\nu (t+t_1)}$ , where $t_1$ describes the time of formation of the boundary layer (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010). Following boundary-layer theory (White et al. Reference Frank2011), we describe the vertical variation of the flow using a transition function, such that

(2.2) \begin{align} u(x,z,t)=u_{\textit{outer}}(x,t)F(z), \end{align}

where $F(z=0)=0$ and $F(z)=1$ for $z\geq h_{bl}$ . For simplicity, and following Bustamante & Ó Náraigh (Reference Bustamante and Ó Náraigh2025) and Eggers et al. (Reference Eggers, Fontelos, Josserand and Zaleski2010) for the 3-D axisymmetric case, we use a step function for the transition function, such that

(2.3) \begin{align} F(z)=\begin{cases} 0,& z\lt h_{bl},\\ 1,& z\gt h_{bl}.\end{cases} \end{align}

We have used the notation $u(x,z,t)$ for the flow in the horizontal direction; the corresponding flow in the vertical direction is denoted by $w(x,z,t)$ . We use the incompressibility condition $\partial _x u+\partial _z w=0$ to obtain an expression for $w(x,z,t)$

(2.4) \begin{align} w(x,z,t)=\begin{cases} 0,& z\lt h_{bl},\\ -\dfrac {z-h_{bl}}{t+t_0},&z\gt h_{bl}.\end{cases} \end{align}

2.2. Lamella height

To model the lamella height $h(x,t)$ , we assume the kinematic condition, such that the interface $h(x,t)$ moves with the flow. Hence, the height of the lamella $h(x,t)$ satisfies

(2.5) \begin{align} \frac {\partial h}{\partial t}+u |_{z=h}\frac {\partial h}{\partial x}=w |_{z=h}. \end{align}

In the remainder of this section, we assume we are in phase 1 of the motion, such that $h(R,t)\gt h_{bl}$ , and such that $(u,w)_{z=h}\neq 0$ . We discuss phase 2 of the motion, wherein $h(R,t)$ reaches a limiting value set by the boundary layer, and $(u,w)_{z=h}=0$ in § 4. Thus, (2.5) becomes

(2.6) \begin{align} \frac {\partial h}{\partial t}+\frac {x}{t+t_0}\frac {\partial h}{\partial x}=w |_{z=h}=-\frac {h-h_{bl}}{t+t_0}. \end{align}

Equation (2.6) has an exact solution (valid in phase 1) by the method of characteristics

(2.7) \begin{align} h(x,t)=\frac {\tau +t_0}{t+t_0}f \!\left (x\frac {\tau +t_0}{t+t_0}\right )+ \underbrace {\frac {2}{3}\left [h_{bl}(t)\frac {t+t_1}{t+t_0}-h_{bl}(\tau )\frac {\tau +t_1}{t+t_0}\right ]}_{=h_{\textit{PI}}(t)}\qquad t\geq \tau , \end{align}

where $f$ is set by the initial conditions. At late times, $f$ is taken to be a constant function, which is the so-called remote asymptotic solution (Roisman Reference Roisman2009). Here, we use the parameter $\tau$ to denote the initial time $t=\tau$ at which point the RL structure seen in figure 1 forms. As in the 3-D axisymmetric case (Marengo et al. Reference Marengo, Antonini, Roisman and Tropea2011), the contribution $h_{\textit{PI}}(t)$ is the viscous correction, it is positive definite and, as such, represents a thickening of the lamella compared with the inviscid scenario.

2.3. Mass conservation

We derive an Ordinary Differential Equation (ODE) for the volume of the rim in what follows. We omit the redundant dimension out of the plane of the page, which carries with it a factor of $\lambda$ , this being the length of the idealised cylindrical droplet. Hence, we work with the rim area $V$ , the lamella area $2\int _0^R \mathrm{d} x\int _0^{h(x,t)}\mathrm{d} z$ and the total area $V_{\textit{tot}}$ , such that

(2.8) \begin{align} V_{\textit{tot}}=V+2\int _0^R \mathrm{d} x\int _0^{h(x,t)}\mathrm{d} z. \end{align}

We use $\mathrm{d} V_{\textit{tot}}/\mathrm{d} t=0$ and (2.6) to get

(2.9) \begin{align} \frac {\mathrm{d} V}{\mathrm{d} t}=2 \!\left [\frac {R}{t+t_0}\left (1-\frac {h_{bl}}{h}\right )-\dot R\right ]h, \end{align}

where $R$ denotes the greatest extent of the lamella at time $t$ , and where we henceforth use the simplified notation

(2.10) \begin{align} h\equiv h(t)\equiv h(R,t). \end{align}

By mass conservation, (2.9) has exact solution

(2.11) \begin{align} V=V_{\textit{tot}}-2g \! \left (\!R\frac {\tau +t_0}{t+t_0}\right )-2 \! \left (\!R\frac {\tau +t_0}{t+t_0}\right ) \! h_{\textit{PI}}(t), \end{align}

where $g(s)=\int _0^s f\!(s)\,\mathrm{d} s$ .

2.4. Momentum

In a similar way, we derive an ODE for the momentum of the rim. We work with one half of the rim momentum $P$ , and one half of the lamella momentum $\rho \int _0^R \mathrm{d} x\int _0^{h(x,t)}\mathrm{d} z$ , such that the total momentum in (say) the right-hand half of the RL structure in figure 1 is

(2.12) \begin{align} P_{\textit{tot}}=P+\rho \lambda \int _0^R \mathrm{d} x\int _0^h u \,\mathrm{d} z, \end{align}

where again, $\lambda$ is the length of the idealised cylindrical droplet in the plane of the page. When $h\gt h_{bl}$ (phase 1) this becomes

(2.13) \begin{align} P_{\textit{tot}}=P+\underbrace {\rho \lambda \int _0^R \frac {x}{t+t_0}\left [h(x,t)-h_{bl}(t)\right ]\mathrm{d} x}_{=P_{\textit{lam}}}. \end{align}

By direct differentiation, we have

(2.14) \begin{align} \frac {\mathrm{d} P_{\textit{lam}}}{\mathrm{d} t} =\rho \lambda \left [\frac {R}{t+t_0}\left [h(R,t)-h_{bl}(t)\right ]\left (\dot R-\frac {R}{t+t_0}\right )-\dfrac {1}{4}\frac {h_{bl} R^2}{(t+t_0)(t+t_1)}\right ]\!. \end{align}

By Newton’s second law, we have

(2.15) \begin{align} \frac {\mathrm{d} P_{\textit{tot}}}{\mathrm{d} t}=F_{\textit{ext}}=-\lambda \gamma (1-\cos \vartheta _a). \end{align}

Here, $\vartheta _a$ is the advancing contact angle, assumed for the present purposes to be constant. Hence, combining (2.14) and (2.15) and using $h\equiv h(R,t)$ we have

(2.16) \begin{align} -\lambda \gamma (1-\cos \vartheta _a)-\frac {R}{t+t_0}\rho \lambda \left [h-h_{bl}(t)\right ]\left (\dot R-\frac {R}{t+t_0}\right )+\frac {1}{4}\frac {\rho \lambda h_{bl} R^2}{(t+t_0)(t+t_1)}=\frac {\mathrm{d} P}{\mathrm{d} t}. \end{align}

We use $\mathrm{d} P/\mathrm{d} t=\rho \lambda U (\mathrm{d} V_{1/2}/\mathrm{d} t)+\rho \lambda V_{1/2}(\mathrm{d} U/\mathrm{d} t)$ ( $V_{1/2}$ is half the rim area) and $U\equiv \dot R$ to re-write the momentum equation in a final form

(2.17) \begin{align} \rho V\frac {\mathrm{d} U}{\mathrm{d} t}=2\rho \big \{\left [h-h_{bl}(t)\right ]\left (u_0-U\right )^2+U^2h_{bl}\big \}+\frac {1}{2}\frac {\rho h_{bl} R^2}{(t+t_0)(t+t_1)}-2\gamma (1-\cos \vartheta _a). \end{align}

Here, we have divided out by the redundant length scale $\lambda$ , and we have introduced the notation $u_0\equiv R/(t+t_0)$ .

2.5. Summary of the model equations

We gather up the model equations here in one place. We emphasise that the equations are for phase 1 of the motion, wherein $h(t)\gt h_{bl}(t)$

(2.18a) \begin{align} \frac {\mathrm{d} V}{\mathrm{d} t}&=2\big [ u_0\left (h-h_{bl}\right )- \textit{Uh}\big ], \\[-9pt] \nonumber \end{align}

(2.18b) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}&=U, \\[-9pt] \nonumber \end{align}

(2.18c) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}&=2\big [(u_0-U)^2\left (h-h_{bl}\right )+U^2h_{bl}\big ]\nonumber \\&\quad +\frac {1}{2}\frac {h_{bl} R^2}{(t+t_0)(t+t_1)}-2\frac {\gamma }{\rho }(1-\cos \vartheta _a), \\[9pt] \nonumber \end{align}

where $h$ is given by (2.10), and $u_0=R/(t+t_0)$ . We analyse (2.18) in the following sections.

3.1. Exact solution

As in the 3-D axisymmetric model, we introduce the velocity defect $\varDelta =u_0-U$ . (3.1) can then be recast in terms of an ODE for $\varDelta$

(3.2) \begin{align} \frac {\mathrm{d} \Delta }{\mathrm{d} t}+\frac {\varDelta }{t+t_0}=-\frac {2h}{V}\big (\varDelta ^2-c^2\big ), \end{align}

where

(3.3) \begin{align} c= \sqrt { \frac {\gamma }{\rho h}(1-\cos \vartheta _a)} \end{align}

is the Taylor–Culick speed ( $c$ is time-dependent, through $h$ ). We further recast (3.2) as an ODE for $V \kern-1pt \Delta$

(3.4) \begin{align} \frac {\mathrm{d} }{\mathrm{d} t}(V\kern-1.5pt\Delta )+\frac {V \kern-1pt \Delta }{t+t_0}=2hc^2=2(\gamma /\rho )(1-\cos \vartheta _a). \end{align}

Thus, there is an exact solution for $V \kern-1pt \Delta$

(3.5) \begin{align} V \kern-1pt \Delta =V_{\textit{init}}\varDelta _{\textit{init}}\left (\frac {\tau +t_0}{t+t_0}\right )+\frac {\gamma }{\rho }\left (1-\cos \vartheta _a\right )\left [t+t_0-\frac {(\tau +t_0)^2}{t+t_0}\right ]\!, \end{align}

where $V_{\textit{init}}$ denotes the initial condition, i.e. $V$ evaluated at $t=\tau$ , and similarly for $\varDelta _{\textit{init}}$ . We recall the definition of $\varDelta$ , to obtain

(3.6) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}-\frac {R}{t+t_0}=-\frac {V_{\textit{init}}}{V}\varDelta _{\textit{init}}\left (\frac {\tau +t_0}{t+t_0}\right )-\frac {\gamma }{\rho }\frac {\left (1-\cos \vartheta _a\right )}{V}\left [t+t_0-\frac {(\tau +t_0)^2}{t+t_0}\right ]\!. \end{align}

This can be further re-written as

(3.7) \begin{align} \frac {\mathrm{d} }{\mathrm{d} t}\left (\frac {R}{t+t_0}\right )=-\frac {V_{\textit{init}}}{V}\frac {\varDelta _{\textit{init}}}{t+t_0}\left (\frac {\tau +t_0}{t+t_0}\right )-\frac {\gamma }{\rho }\frac {\left (1-\cos \vartheta _a\right )}{V}\left [1-\left (\frac {\tau +t_0}{t+t_0}\right )^2\right ]\!. \end{align}

To make further progress, we take $f$ to be a constant in (2.7) for $h(R,t)\equiv h$ . This approximation is valid at late times and is pursued here for mathematical convenience. The foregoing analysis is still valid without this approximation, but would become intractable. Hence, we obtain

(3.8) \begin{align} h(R,t)=\left (\frac {\tau +t_0}{t+t_0}\right )h_{\textit{init}}. \end{align}

We further assume that the lamella has flattened into a rectangular shape, hence $V=V_{\textit{tot}}-2R[(\tau +t_0)/(t+t_0)]h_{\textit{init}}$ , an assumption discussed in more detail below (§ 3.3). We also group together the terms $R[(\tau +t_0)/(t+t_0)]$ as $\eta$ . Thus, we have

(3.9) \begin{align} \frac {\mathrm{d}\eta }{\mathrm{d} t}=-\varDelta _{\textit{init}}\frac {V_{\textit{init}}}{V_{\textit{tot}}-2\eta h_{\textit{init}}}\left (\frac {\tau +t_0}{t+t_0}\right )^2-\frac {\gamma (\tau +t_0)}{\rho }\frac {\left (1-\cos \vartheta _a\right )}{V_{\textit{tot}}-2\eta h_{\textit{init}}}\left [1-\left (\frac {\tau +t_0}{t+t_0}\right )^2\right ]\!. \end{align}

This is a separable ODE. Following the steps outlined in Appendix B, (3.9) can be solved to yield

(3.10) \begin{align} Y=X- \sqrt {\epsilon ^2 X^2+A (X^3+X-2X^2)+B \epsilon (X^2-X)},\qquad X\geq 1, \end{align}

where

(3.11a) \begin{align} X&= \frac {t+t_0}{\tau +t_0}, \\[-9pt] \nonumber \end{align}

(3.11b) \begin{align} Y&=\frac {R}{V_{\textit{tot}}/(2h_{\textit{init}})}, \\[-9pt] \nonumber \end{align}

(3.11c) \begin{align} A&=\frac {c_*^2(\tau +t_0)^2}{V_{\textit{tot}}}, \\[-9pt] \nonumber \end{align}

(3.11d) \begin{align} c_*^2&=4(h_{\textit{init}}/V_{\textit{tot}})(\gamma /\rho )(1-\cos \vartheta _a), \\[-9pt] \nonumber \end{align}

(3.11e) \begin{align} B&=\frac {4\varDelta _{\textit{init}}(\tau +t_0)h_{\textit{init}}}{V_{\textit{tot}}}, \\[-9pt] \nonumber \end{align}

(3.11f) \begin{align} \epsilon &=\frac {V_{\textit{init}}}{V_{\textit{tot}}}. \\[9pt] \nonumber \end{align}

Equation (3.10) is the required exact solution of the RL model. A sample solution at ${\textit{We}}=100$ is shown in figure 2.

Figure 2.

Representative plot of the analytical solution (3.10), at ${\textit{We}}=100$ . The parameters (3.11) have been selected using the methodology in § 3.3, below.

3.2. Maximum spreading radius

The maximum spreading occurs at $\mathrm{d} R/\mathrm{d} t=0$ , hence $\mathrm{d} Y/\mathrm{d} X=0$ . Based on (3.10), this occurs when

(3.12) \begin{align}& 4 [\epsilon ^2 X^2+A (X^3+X-2X^2 )+ \epsilon \kern-1pt B (X^2-X)]\nonumber\\&\quad = [3\textit{AX}^2+2X (\epsilon ^2+\epsilon \kern-1pt B-2A )+(A-\epsilon \kern-1pt B)]^2. \end{align}

Assuming the maximum spreading occurs at large value of $t$ (hence, $X$ ), this gives $X\approx 4/(9A)$ , a result which will be verified a posteriori. Correspondingly,

(3.13) \begin{align} Y_{\textit{max}}\approx (4/9)A^{-1}-\sqrt {A (4/9)^3 A^{-3}}=\frac {4}{27}A^{-1}. \end{align}

By using the explicit definition of $A$ from (3.11) we have (details in Appendix B)

(3.14) \begin{align} \frac {R_{\textit{max}}}{R_0}\approx \frac {1}{27} \frac {V_{\textit{tot}}^2}{ h_{\textit{init}} U_0^2 (\tau +t_0)^2 R_0 }\left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}R_0}\right )\frac {\textit{We}}{1-\cos \vartheta _a}. \end{align}

In order for the large- $X$ approximation to be valid, we require $A$ to be small, and hence ${\textit{We}}$ to be large. Thus, (3.14) is valid in the large Weber number limit. This scaling behaviour ( $\mathcal{R}_{\textit{max}}/R_0\sim {\textit{We}}$ , for ${\textit{We}}\gg 1$ ) is exactly what is seen in previous published work on 2-D droplet impact (Ó Náraigh & Mairal Reference Ó Náraigh and Mairal2023), and is consistent with the energy-budget analysis (1.5).

3.3. Parameter estimation

In order to make a systematic comparison between (3.14) and both the correlation (1.5) and numerical simulations, it is necessary to estimate the pre-factors in (3.14).

Following Amirfazli et al. (Reference Amirfazli, Bustamante, Hu and Ó Náraigh2024) for the RL model in 3-D axisymmetry (inviscid case), we estimate $h_{\textit{init}}/(\tau +t_0)$ as $U_0$ , hence

(3.15) \begin{align} \frac {R_{\textit{max}}}{R_0}\sim \frac {1}{27} \frac {V_{\textit{tot}}^2}{ h_{\textit{init}}^3 R_0 }\left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}R_0}\right )\frac {\textit{We}}{1-\cos \vartheta _a},\qquad {\textit{We}}\gg 1. \end{align}

We estimate $h_{\textit{init}}$ using dimensional analysis. Hence, we assume that at the onset of the RL phase, the rim is small compared with the lamella, hence $V_{\textit{init}}=0$ . We crudely approximate the lamella as a rectangle of sides of length $2R_{\textit{init}}$ and $h_{\textit{init}}$ . On grounds of dimensional analysis, we can say that an estimation based on a more complex lamella shape will yield similar results, although the prefactors in the foregoing analysis will change.

An energy balance at this instant gives

(3.16) \begin{align} \text{Kinetic energy of sheet}+\text{Surface energy of sheet}=\frac {1}{2}\rho V_{\textit{tot}}U_0^2+\gamma (2\pi R_0). \end{align}

In more detail, this expression is:

(3.17) \begin{align}& \rho \int _0^{R_{\textit{init}}}\int _0^{h_{\textit{init}}}\left [\left (\frac {x}{\tau +t_0}\right )^2+\left (\frac {z}{\tau +t_0}\right )^2 \right ]\mathrm{d} x\,\mathrm{d} z+[2R_{\textit{init}}(1-\cos \vartheta _a)+2h_{\textit{init}}]\gamma \nonumber \\&\quad = \gamma (2\pi R_0)+\frac {1}{2}\rho V_{\textit{tot}}U_0^2. \end{align}

We evaluate the integral, and use the fact that $V_{\textit{tot}}=2R_{\textit{init}}h_{\textit{init}}=\pi R_0^2$ , hence $R_{\textit{init}}=V_{\textit{tot}}/(2h_{\textit{init}})$ . We further divide both sides by $\rho U_0^2 R_0^2$ . This gives

(3.18) \begin{align}& \frac {1}{3} \frac {1}{R_0^2 U_0^2(\tau +t_0)^2}\left [ \left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}}\right )^{\! 3} h_{\textit{init}} +\left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}}\right )\! h_{\textit{init}}^3\right ] \nonumber\\&\quad +\frac {1}{R_0{\textit{We}}}\left [ 2\left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}}\right )\left (1-\cos \vartheta _a\right )+2h_{\textit{init}}\right ]=\frac {2\pi }{\textit{We}}+\frac {1}{2}\pi .\end{align}

We multiply both sides by $h_{\textit{init}}^2$ and go over to the large- ${\textit{We}}$ limit. This gives

(3.19) \begin{align} \frac {1}{3} \underbrace {\frac {h_{\textit{init}}^2}{U_0^2(\tau +t_0)^2}}_{=1}\frac {1}{R_0^2}\left [ \left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}}\right )\!^3 h_{\textit{init}} +\left (\frac {V_{\textit{tot}}}{2h_{\textit{init}}}\right )\! h_{\textit{init}}^3\right ]\sim \frac {1}{2} \pi h_{\textit{init}}^2,\qquad {\textit{We}}\gg 1. \end{align}

After effecting cancellations, this simplifies to

(3.20) \begin{align} \frac {V_{\textit{tot}}^2}{8}\sim h_{\textit{init}}^4,\qquad {\textit{We}}\gg 1. \end{align}

We go back up to (3.15), with $h_{\textit{init}}^4\sim V_{\textit{tot}}^2/8$ for large ${\textit{We}}$ , hence

(3.21) \begin{align} \frac {R_{\textit{max}}}{R_0}\sim \frac {4\pi }{27} \frac {\textit{We}}{1-\cos \vartheta _a},\qquad {\textit{We}} \gg 1. \end{align}

To two significant figures, this gives $R_{\textit{max}}/R_0\sim 0.47\,{\textit{We}}/(1-\cos \vartheta _a)$ , valid for ${\textit{We}} \gg 1$ .

Figure 3.

Comparison between the direct numerical simulation data of Wu et al. (Reference Wu, Gui, Yang, Tu and Jiang2021) and Ó Náraigh & Mairal (Reference Ó Náraigh and Mairal2023) and the theoretical models of droplet spreading. Solid line: energy-budget analysis ((1.2)), with $b=0.5$ . Dashed line: the RL model (3.10)–(3.11). The parameters in (3.11) have been selected using the methodology in § 3.3.

3.4. Comparisons

We compare the result (3.21) with other results from the literature. First, we look at the energy budget (1.2) in the limit of large ${\textit{We}}$ and ${\textit{Re}}=\infty$ , which gives

(3.22) \begin{align} \frac {R_{\textit{max}}}{R_0}\sim \frac {\pi (1-b)}{4} \frac {\textit{We}}{1-\cos \vartheta _a}. \end{align}

Following the energy-budget results of Wildeman et al. (Reference Wildeman, Visser, Sun and Lohse2016) (albeit for 3-D axisymmetric droplets), we take $b=1/2$ as per the ‘head loss’ argument therein. To two significant figures, this gives $R_{\textit{max}}/R_0\sim 0.39\, {\textit{We}}/(1-\cos \vartheta _a)$ , to be compared with $R_{\textit{max}}/R_0\sim 0.47\,{\textit{We}}/(1-\cos \vartheta _a)$ for the RL model. Since the aim of this section is to derive a scaling law with an $O(1)$ prefactor, the discrepancy between the prefactors here can be considered small. Furthermore, we have compared both models with Direct Numerical Simulation (DNS) (figure 3), and there is good agreement between the DNS and the models, with little to be said in favour of accepting one prefactor over another.

4.1. Recasting of equations for the RL system

We recall (2.18c ) for the momentum

(4.1) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}=2\big [(u_0-U)^2\left (h-h_{bl}\right )+U^2h_{bl}\big ]+\frac {1}{2}\frac {h_{bl} R^2}{(t+t_0)(t+t_1)}-2\frac {\gamma }{\rho }(1-\cos \vartheta _a). \end{align}

We re-write

(4.2) \begin{align} (u_0-U)^2(h-h_{bl})+U^2h_{bl} = \left [ u_0\left (\!1-\frac {h_{bl}}{h}\right )-U\right ]^2h + u_0^2 h_{bl} \left (\!1-\frac {h_{bl}}{h}\right )\!. \end{align}

We introduce $\overline {u}=u_0[1-(h_{bl}/h)]$ . Thus, (4.1) becomes

(4.3) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}=2\big (\overline {u}-U\big )^2h+\underbrace {2u_0^2 h_{bl}\left (\!1-\frac {h_{bl}}{h}\right )+\frac {1}{2}\frac {h_{bl} R^2}{(t+t_0)(t+t_1)}}-2\frac {\gamma }{\rho }(1-\cos \vartheta _a). \end{align}

In a previous work on the analogous 3-D axisymmetric RL model (Bustamante & Ó Náraigh Reference Bustamante and Ó Náraigh2025), a term analogous to the one with the underbrace in (4.3) was omitted, on the basis that it was negligibly small (as justified by Gordillo et al. Reference Gordillo, Riboux and Quintero2019) and also, led to a mathematically tractable set of equations. We follow the same approach here, and we therefore work with a modified momentum equation

(4.4) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}=2 [ (\overline {u}-U)^2 -c^2]h, \end{align}

where $c^2=[\gamma /(\rho h)](1-\cos \vartheta _a)$ , and where $c$ is the Taylor–Culick speed. Therefore, for the purposes of this section, we examine the following RL model:

(4.5a) \begin{align} \frac {\mathrm{d} V}{\mathrm{d} t}&=2\left (\overline {u}-U\right )h, \\[-9pt] \nonumber \end{align}

(4.5b) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}&=U, \\[-9pt] \nonumber \end{align}

(4.5c) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}&=2 [ (\overline {u}-U)^2 -c^2]h. \\[9pt] \nonumber \end{align}

4.2. Initial conditions

Initial conditions apply at the onset of phase 1. The model has built-in an implied set of initial conditions that give rise to rim generation. At time $t=R_0/U_0=\tau$ , we take the rim volume to be zero, as in the analogous work on 3-D axisymmetric RL models (Amirfazli et al. Reference Amirfazli, Bustamante, Hu and Ó Náraigh2024; Bustamante & Ó Náraigh Reference Bustamante and Ó Náraigh2025). This gives a natural way to parametrise the rim generation phenomenon, since, by taking $V(\tau )=0$ , we obtain (from (4.5c ))

(4.6) \begin{align} (\overline {u}-U)^2h - \frac {\gamma (1-\cos \vartheta _a)}{\rho } =0 \qquad \text{at }t=\tau . \end{align}

Hence

(4.7) \begin{align} U(t=\tau )=\overline {u}-\sqrt { \frac {\gamma (1-\cos \vartheta _a)}{\rho h}}, \end{align}

or

(4.8) \begin{align} U(t=\tau )=\frac {R_{\textit{init}}}{\tau +t_0}\left (1-\frac {h_{bl}}{h_{\textit{init}}}\right )-c(\tau ). \end{align}

The values $R_{\textit{init}}$ and $h_{\textit{init}}$ describe the initial lamella, i.e. just prior to the formation of the rim. For the inviscid case (§ 3), it is important to describe carefully the dependence of these parameters on ${\textit{We}}$ , as this has implications for the value of the prefactor in $R_{\textit{max}}/R_0\sim {\textit{We}}$ for ${\textit{We}} \gg 1$ . For the viscous case, this seems less important. For instance, Roisman et al. (Reference Roisman, Rioboo and Tropea2002) determine $R_{\textit{init}}$ and $h_{\textit{init}}$ using an energy-budget analysis, whereas Eggers et al. (Reference Eggers, Fontelos, Josserand and Zaleski2010) use a geometric argument. In these studies, the final dependence of the spreading radius on ${\textit{Re}}$ and ${\textit{We}}$ is insensitive to the approach used to set the initial conditions. We use a similar geometric argument to Eggers et al. (Reference Eggers, Fontelos, Josserand and Zaleski2010) here: we assume that the droplet assumes a ‘pancake’ shape at time $\tau$ , with initial height $h_{\textit{init}}=R_0/2$ . This gives $R_{\textit{init}}=V_{\textit{tot}}/(2h_{\textit{init}})$ , where $V_{\textit{tot}}=\pi R_0^2$ is the total (conserved) area.

4.3. Phase 2 of the motion

In phase 2 of the motion, the boundary layer reaches the height of the lamella. In the present simplified mathematical model, all flow inside the lamella then stops and (2.5) for the interface height reduces to $\partial h/\partial t=0$ , hence $h=\text{Const.}$ . We use $t_*$ to denote the onset time of phase 2; $t_*$ is obtained from (2.7) by solving

(4.9) \begin{align} h(t_*)=h_{bl}(t_*). \end{align}

Here, we work with the remote asymptotic solution, such that $f(\boldsymbol{\cdot })=\text{Const.}$ in (2.7). In this case, (4.5) simplify greatly

(4.10a) \begin{align} \frac {\mathrm{d} V}{\mathrm{d} t}&= -2U h_*, \\[-9pt] \nonumber \end{align}

(4.10b) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}&= U, \\[-9pt] \nonumber \end{align}

(4.10c) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}&=2\left [U^2h_* -\frac {\gamma }{\rho }(1-\cos \vartheta _a)\right ]\!. \\[9pt] \nonumber \end{align}

Equation (4.10) have the simple fixed-point solution

(4.11) \begin{align} U=-c_*,\qquad V=2Uh_*, \end{align}

which remains valid until $R=R(t_*)-c_* t$ touches down to zero (here, $R(t_*)$ is the value of $R(t)$ at the onset of phase 2). Also,

(4.12) \begin{align} c_*=\sqrt { \frac {\gamma }{\rho h_*}(1-\cos \vartheta _a)} \end{align}

is the Taylor–Culick retraction speed. We note in passing the use of $c_*$ for a different speed in § 3; we rely on the reader’s judgment to distinguish these, based on the context.

The sharp jump between phases 1 and 2 is an artefact of the model, due to the fact that $\mathrm{d} h/\mathrm{d} t$ changes to zero abruptly in this phase. Thereafter, $R(t)$ decreases, meaning that rim retraction begins immediately in phase 2, if it is has not done so already in phase 1. Nevertheless, phase 2 of the model does capture accurately the competing effects of intertia and surface tension in rim retraction. Indeed, (4.10) is identical in form to the model put forward by Taylor (Reference Taylor1959) to describe the retraction of a liquid sheet (once the necessary changes in the problem set-up have been accomplished). The more realistic radially symmetric version of the model is also presented in the work by Taylor (Reference Taylor1959), and by Culick (Reference Culick1960). An exact solution of the radially symmetric model is presented by Bustamante & Ó Náraigh (Reference Bustamante and Ó Náraigh2025).

4.4. Analysis of the equations of motion using inequalities

Unlike in the inviscid case, the viscous RL model (4.5) does not admit an exact solution. However, we can characterise the solution using differential inequalities, similar to what was done in the 3-D axisymmetric case by Bustamante & Ó Náraigh (Reference Bustamante and Ó Náraigh2025). We summarise the approach here, and refer the interested reader to Appendix C for the lengthy calculations which underpin this summary. We start with the velocity defect

(4.13) \begin{align} \varDelta =u_0\left (1-\frac {h_{bl}}{h}\right )-U. \end{align}

Hence, (4.1) become

(4.14a) \begin{align} \frac {\mathrm{d} V}{\mathrm{d} t}&=2\Delta h, \\[-9pt] \nonumber \end{align}

(4.14b) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}&=U, \\[-9pt] \nonumber \end{align}

(4.14c) \begin{align} V\frac {\mathrm{d} U}{\mathrm{d} t}&=2 (\varDelta ^2-c^2 )h. \\[9pt] \nonumber \end{align}

By direct computation, we get

(4.15) \begin{align} \frac {\mathrm{d}\Delta }{\mathrm{d} t}+\frac {\varDelta }{t+t_0}\left (1-\frac {h_{bl}}{h}\right )= -2 (\varDelta ^2-c^2 )(h/V)-\frac {u_0}{t+t_0}\frac {h_{bl}}{h}\underbrace {\left [{1-\frac {h_{bl}}{h}}+\frac {1}{2}\frac {t+t_0}{t+t_1}\right ]}_{=\varPhi (t)\geq 0}\!. \end{align}

Since each term on the right-hand side is now of definite sign, we can conclude

(4.16) \begin{align} \frac {\mathrm{d}\Delta }{\mathrm{d} t}+\frac {\varDelta }{t+t_0}\left (1-\frac {h_{bl}}{h}\right )\leq 2c^2(h/V). \end{align}

The appearance here of a differential inequality linear in $\varDelta$ enables us to apply Gronwall’s inequality. We apply this inequality to (4.16) with $V(t=\tau )=0$ to obtain

(4.17) \begin{align} \varDelta \leq (\gamma /\rho )(1-\cos \vartheta _a)(h/V)I_h(t), \end{align}

where

(4.18) \begin{align} I_h(t)=\int _\tau ^t \frac {\mathrm{d} t}{h}. \end{align}

The full details of this calculation are provided in Appendix C. We use $\varDelta =[R/(t+t_0)] [1-(h_{bl}/h)]-(\mathrm{d} R/\mathrm{d} t)$ to re-write inequality (4.17) as

(4.19) \begin{align} \frac {\mathrm{d} R}{\mathrm{d} t}-\frac {R}{t+t_0}\left (1-\frac {h_{bl}}{h}\right )\geq -2(\gamma /\rho )(1-\cos \vartheta _a)(h/V)I_h(t). \end{align}

This can be re-written as

(4.20) \begin{align} \frac {\mathrm{d} }{\mathrm{d} t}(\textit{Rh})\geq -2(\gamma /\rho )(1-\cos \vartheta _a)(h^2/V)I_h(t). \end{align}

We note that $V=V_{\textit{tot}}-2\textit{Rh}$ and identify $\eta = \textit{Rh}$ . Hence, we re-write inequality (4.20) as

(4.21) \begin{align} \frac {\mathrm{d} \eta }{\mathrm{d} t}\geq -2(\gamma /\rho )(1-\cos \vartheta _a)\frac {h^2 I_h(t)}{V_{\textit{tot}}-2\eta }. \end{align}

This can be solved to give

(4.22) \begin{align} V_{\textit{tot}}\eta -\eta ^2\geq V_{\textit{tot}}\eta _{\textit{init}}-\eta _{\textit{init}}^2 -(\gamma /\rho )(1-\cos \vartheta _a)\underbrace {\left [2\int _\tau ^t h^2 I_h(t)\mathrm{d} t\right ]}_{=\Delta G(t)} \! . \end{align}

Inequality (4.22) is a quadratic inequality. Critical points occur at $\eta _{\textit{cr},\pm }(t)$ , where

(4.23) \begin{align} \eta _{\textit{cr},\pm }(t)=\frac {1}{2}V_{\textit{tot}}\pm [(\gamma /\rho )(1-\cos \vartheta _a)]^{1/2}[\Delta G(t)]^{1/2} \end{align}

(we write $\eta _{\textit{cr},\pm }(t)$ to indicate that the critical points are time-dependent). Since $\eta (t)=R(t)h(t)$ , we have $R(t)h(t)\leq \eta _{\textit{cr},+}(t)$ . Similarly, with the negative sign chosen, we get $R(t)h(t)\geq \eta _{\textit{cr},-}(t)$ . Since this is true for all $t$ , we must have $R(t_*)h(t_*)\geq \eta _{\textit{cr},-}(t_*)$ , hence

(4.24) \begin{align} R(t_*)\geq \frac {\eta _{\textit{cr}-}(t_*)}{h(t_*)}\mathop {=}\limits ^{\text{(4.23)}} \frac {1}{2}\frac {V_{\textit{tot}}}{h(t_*)}-\underbrace { [(\gamma /\rho )(1-\cos \vartheta _a)]^{1/2}\frac {[\Delta G(t_*)]^{1/2}}{h(t_*)}}. \end{align}

We apply the following observations to inequality (4.24):

1. (i) since the term with the underbrace is proportional to $\gamma ^{1/2}$ , and $\gamma \propto {\textit{We}}^{-1}$ , the right-hand side of the inequality is positive, for sufficiently large ${\textit{We}}$ ;

2. (ii) by definition, we have $\max _t R(t)\geq R(t_*)$ . Here, the maximum $\max _t ({\cdots} )$ is taken over all times for which the RL model applies ( $t\geq \tau$ ).

Putting these observations together, we have

(4.25) \begin{align} \max _t R(t)\geq R(t_*)\geq \frac {1}{2}\frac {V_{\textit{tot}}}{h_*}-\big [(\gamma /\rho )(1-\cos \vartheta _a)\big ]^{1/2}\frac {[\Delta G(t_*)]^{1/2}}{h_*}\gt 0. \end{align}

Here, the last inequality is true for ${\textit{We}}$ sufficiently large; without this, the string of inequalities would be vacuous.

For an upper bound, we reuse the argument due to Bustamante & Ó Náraigh (Reference Bustamante and Ó Náraigh2025). We have

(4.26) \begin{align} 2\underbrace {R_{\textit{max}}}_{=R(t_{\textit{max}})}h(t_{\textit{max}})=V_{\textit{tot}}-V(t_{\textit{max}})\leq V_{\textit{tot}}. \end{align}

Hence,

(4.27) \begin{align} R_{\textit{max}}\leq \frac {V_{\textit{tot}}}{2 h(t_{\textit{max}})}. \end{align}

If the maximum occurs in phase 1, then $h(t)$ is monotone decreasing, hence $h(t_{\textit{max}})\geq h_*$ , hence $1/h(t_{\textit{max}})\leq 1/h_*$ and

(4.28) \begin{align} R_{\textit{max}}\leq \frac {V_{\textit{tot}}}{2 h_*}. \end{align}

If the maximum occurs in phase 2, then $h(t)$ is a constant function and equal to $h_*$ , so it is still the case that $R_{\textit{max}}\leq V_{\textit{tot}}/(2h_*)$ . Hence, in both cases, inequality (4.28) is true. Hence, by a sandwich result, we have

(4.29) \begin{align} \frac {\frac {1}{2}V_{\textit{tot}} - \left [(\gamma /\rho )(1-\cos \vartheta _a)\right ]^{1/2}[\Delta G(t_*)]^{1/2}}{h_*} \leq \max _t R(t)\leq \frac {V_{\textit{tot}}}{2 h_*}. \end{align}

4.5. Evaluation of the bounds

To make further progress, we study $h_*$ in more detail, to extract its ${\textit{Re}}$ -dependence. For this purpose, we solve (2.6) numerically. We assume the remote asymptotic approximation such that $h(x,t)$ no longer depends on space, and such that (2.6) becomes an ODE, specifically

(4.30) \begin{align} \frac {\mathrm{d} h}{\mathrm{d} t}=-\frac {h-h_{bl}}{t+t_0},\qquad \tau \ll t\lt t_*,\qquad h_{bl}=\xi \sqrt {\nu (t+t_1)} \end{align}

(as we are using the remote asymptotic approximation, the inequality $\tau \ll t$ is applied here). As per the previous discussion in this section, we use the initial condition $h(\tau )=R_0/2$ . Equation (4.30) can be solved explicitly; however, it is just as straightforward to solve it numerically in Matlab. We use the stiff solver ode15s, which readily handles the transition which occurs when $\mathrm{d} h/\mathrm{d} t$ jumps abruptly to zero at $t=t_*$ . We furthermore solve (4.30) using dimensionless variables, such that all length scales are normalised by $R_0$ , and all velocity scales by $U_0$ . The fluid density $\rho$ sets the scale for mass. Equation (4.30) can then be non-dimensionalised by formally replacing $\nu$ by ${\textit{Re}}^{-1}$ , and by implicitly taking $h\rightarrow h/R_0$ , $t\rightarrow t U_0/R_0$ and similarly for $\tau$ , $t_0$ and $t_1$ . The dimensionless parameter $\xi$ is unchanged in this process.

Figure 4.

Numerical analysis of (4.30) showing the asymptotic dependence of $h_*$ , $t_*$ and $\Delta G(t_*)$ on ${\textit{Re}}$ , for large values of ${\textit{Re}}$ .

We assume that $\tau =U_0/R_0$ marks the onset of the RL phase (or, $\tau =1$ in dimensionless variables). Hence, (4.30) contains three free parameters: $\xi$ , $t_0$ and $t_1$ . We choose these parameters via nonlinear optimisation, with the objective function to be described in what follows. For the time being, we report these values here

(4.31) \begin{align} \xi =0.6908,\qquad t_0=2.0,\qquad t_1 =2.8907. \end{align}

Results are shown in figure 4. From the figure, we find

(4.32) \begin{align} h_*=k_h {\textit{Re}}^{-1/3},\qquad t_*=k_t {\textit{Re}}^{1/3},\qquad {\textit{Re}} \gg 1, \end{align}

and

(4.33) \begin{align} \Delta G =k_G {\textit{Re}}^{1/3},\qquad {\textit{Re}} \gg 1, \end{align}

where $k_h$ , $k_t$ and $k_G$ are given in table 1. We apply the results (4.32)–(4.33) to the bound (4.29). This gives the sandwich result

(4.34) \begin{align} k_1{\textit{Re}}^{1/3}-k_2 (1-\cos \vartheta _a)^{1/2} ({\textit{Re}}/{\textit{We}})^{1/2} \leq \frac {\max R(t)}{R_0}\leq k_1{\textit{Re}}^{1/3} ,\qquad {\textit{Re}} \gg 1, \end{align}

where $k_1=\pi /(2 k_h)$ , and $k_2=k_G^{1/2}/k_h$ .

Table 1.

Numerical values for the parameters $k_h$ , $k_t$ and $k_G$ , selected according to the methodology given below in § 5.3.

4.6. Correction due to the rim

The results so far have involved expressions for $R_{\textit{max}}$ , the maximum extent of the lamella. However, what is of interest is $\mathcal{R}_{\textit{max}}$ , being the maximum spreading radius of the RL structure. The maximum spreading radius can be written as $\mathcal{R}_{\textit{max}} = R_{\textit{max}} + 2a$ , where $2a$ is the footprint of the rim, and where $a$ can be estimated from figure 1 using a circular-segment construction as

(4.35) \begin{align} a=\sin \vartheta _a\sqrt {\frac {2V}{2\vartheta _a-\sin (\vartheta _a)}} ,\end{align}

where $V$ is the area of the rim (for an in-depth presentation of this construction, see Amirfazli et al. (Reference Amirfazli, Bustamante, Hu and Ó Náraigh2024)). Since $V\leq V_{\textit{tot}}$ , the correction to $\mathcal{R}_{\textit{max}}$ does not exceed a term trivially proportional to ${\textit{Re}}^0$ and hence, can be ignored in the remaining calculations. For the avoidance of doubt, in the remainder of this article we use $R_{\textit{max}}$ as the quantity of interest.

5.1. Numerical implementation

We solve (5.1) numerically using a rectangular domain of size $L_x\times L_z$ as shown in figure 5. We take $L_x=0.046\,\mathrm{m}$ , $L_z=0.012\,\mathrm{m}$ and the initial droplet size $R_0=0.003\mathrm{m}$ . To simulate droplet impact, the droplet is seeded with an initial velocity $\boldsymbol{u}=(0,0,-1)\,\mathrm{m}\,\mathrm{s}^{-1}$ . We discretise the domain using a simple block mesh in OpenFOAM with $2250$ gridpoints in the $x$ -direction and $750$ gridpoints in the $z$ -direction. A mesh-refinement study is described below. No grading of the mesh is applied in the $x$ -direction. However, simple grading is applied in the $z$ -direction such that the grid spacing $\Delta z$ of the cell closest to the bottom wall ( $z=0$ ) is half that of the cell at the top wall. This is to better capture the flow inside the lamella at late times. To capture the droplet at maximum spreading, some simulations require a larger value of $L_x$ . In these cases, $L_x$ is increased while keeping the resolution $\Delta x=L_x/n_x$ the same.

We use the PIMPLE algorithm (a combination of Pressure Implicit with Splitting of Operator (PISO) and Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithms) to solve for the pressure and implement the incompressibility condition (5.1c). This is a combination of the PISO and SIMPLE algorithms, as described in Greenshields (Reference Greenshields2025) (Chapter 5 therein). In this context, the Navier–Stokes equations are solved explicitly, with 1 outer corrector and 3 inner correctors per time step. Full details of the OpenFOAM case have been posted in the accompanying GitHub repository (Ó Náraigh Reference Ó Náraigh2025). We have also checked the OpenFOAM simulation results against a standard diffuse-interface method (Ó Náraigh & Mairal Reference Ó Náraigh and Mairal2023) and both approaches yield the same results. A simulation campaign using the diffuse-interface method is reported below, showing the same behaviour as the other campaigns, which use the OpenFOAM VOF approach.

Boundary and initial conditions are described in figure 5. We use the ‘constantContactAngle’ BC on the scalar field $\alpha$ to enforce an equilibrium contact angle at the triple point. We also use a Navier slip BC on the velocity field at the bottom wall. A standard no-slip condition is applied at the top wall. We describe the slip BC in more detail below. The third dimension into the plane of the page is chosen to be very small so as to make the simulation quasi-two-dimensional. BCs of type ‘empty’ are therefore applied to the corresponding faces. A truly 2-D simulation could yield slightly different results.

We solve (5.1) with standard transport properties for an air–water system, as summarised in table 2. Based on these properties, we have ${\textit{We}}=41.67$ and ${\textit{Re}}=3000$ . In order to study the effect of the Weber number and the Reynolds number on the maximum spreading radius, we will vary $\gamma$ and $\nu$ (water) systematically in what follows.

Table 2.

Transport properties used in the model.

5.2. Slip BC

In a 2-D setting, and using the coordinate system in figure 5, the Navier slip boundary at the bottom wall is expressed by

(5.2) \begin{align} u_x+\lambda _{\textit{s}}\frac {\partial u_x}{\partial z}=0,\qquad z=0, \end{align}

where $u_x$ denotes the velocity component in the $x$ -direction and $\lambda _{\textit{s}}$ is the slip length. If we denote $u_0$ as the tangential velocity ( $=u_x$ ) on the lower boundary and $u_1$ as the tangential velocity at the adjacent cell centre, then the discretised version of (5.2) reads

(5.3) \begin{align} u_0+\lambda _{\textit{s}}\frac {\left (u_0-u_1\right )}{d/2}=0, \end{align}

where $d$ is the size of the cell adjacent to the lower wall. Re-arranging gives

(5.4) \begin{align} u_0= (1-f)u_1,\qquad f=\frac {d}{d+2\lambda }. \end{align}

The parameter $f$ is inputted into OpenFOAM and describes the slip ( $f$ corresponds to ‘valueFraction’ in the BC of type ‘partialSlip’). A value $f=1$ corresponds to no slip and a value $f=0$ corresponds to perfect slip.

As noted by Legendre & Maglio (Reference Legendre and Maglio2015), the VOF method applied to moving contact lines exhibits poor convergence unless a slip condition with constant slip length is used. For this reason, we use a slip length $\lambda _{\textit{s}}=8.32\times 10^{-6}\,\mathrm{m}$ (hence, $(\lambda _{\textit{s}}/R_0)\times 100\lt 0.28\,\%$ , i.e. much less than the initial droplet radius). We show a mesh-refinement study at constant $\lambda _{\textit{s}}$ (hence, varying $f$ , see table 3) in figure 6. This figure shows that the mesh of size $2250\times 750$ cells is sufficient to produce a converged simulation, and hence accurately to describe the spreading behaviour. This amounts to a large number of cells for a 2-D simulation (over 1.6 million), this large number of cells is needed to accurately resolve the thin lamella structure which forms at late times. This can be seen e.g. in figure 5, which overlays annotation on a snapshot of the simulation, to produce a schematic diagram. Since the simulations are two-dimensional, the present approach of achieving converged simulations simply by adding more cells throughout the computational domain is valid. For equivalent 3-D simulations, a more computationally efficient approach would be required, e.g. local grid refinement near the bottom wall, or adaptive mesh refinement.

Table 3.

Mesh refinement study at constant slip length, $\lambda _{\textit{s}}=8.32\times 10^{-6}\,\mathrm{m}$ .

Figure 6.

Mesh refinement study at constant slip length, $\lambda _{\textit{s}}=8.32\times 10^{-6}\,\mathrm{m}$ , showing the triple point $R(t)$ as a function of time. Slip length: $\lambda _{\textit{s}}=8.32\times 10^{-6}\,\mathrm{m}$ . Other parameters as in tables 2 and 3.

We remark that a constant contact-angle model is used throughout the simulations, which rather unphysically forces the contact angle to its equilibrium value even when the triple-point velocity is non-zero. However, this simplified approach is sufficient to validate the bounds (4.34), and hence, to achieve one of the main aims of the paper. To maintain continuity with the other parts of the paper, we use the symbol $\vartheta _a$ to denote this constant contact angle. The implications of having a dynamic contact-angle model are addressed in § 6.

5.3. Results

In a first campaign of simulations (CAMPN1), we vary $\gamma$ (hence ${\textit{We}}$ ) while keeping ${\textit{Re}}$ and $\vartheta _a=90^\circ$ fixed, and we show the results in figure 7. The bounds are quite wide in general. However, the distance between the bounds and the value of $\beta _{\textit{max}}$ computed from simulations decreases with increasing ${\textit{We}}$ . This is consistent with the fact that via a sandwich result, both bounds imply that $\beta _{\textit{max}}\sim k_1{\textit{Re}}^{1/5}$ as ${\textit{We}}\rightarrow \infty$ , and for large but finite values of ${\textit{Re}}$ . In short, the bounds sharpen as ${\textit{We}}$ increases.

The bounds are more useful in practice for determining a correlation for $\beta _{\textit{max}}$ as a function of ${\textit{Re}}$ and ${\textit{We}}$ , which we proceed to do now. For this purpose, we perform further campaigns of simulations, starting with CAMPN2, where we fix ${\textit{We}}$ and ${\textit{Re}}$ as per table 2 and vary $\vartheta _a$ through a range $\vartheta _a\in [40^\circ ,130^\circ ]$ . The results from this campaign, together with CAMPN1 at varying Weber number all collapse onto a single curve (figure 8). For large values of $X=[{\textit{We}}/(1-\cos \vartheta _a)]^{1/2}$ , the variable $Y=\beta _{\textit{max}} X/{\textit{Re}}^{1/3}$ depends linearly on $X$ , as evidenced by the trend line in the figure. Given the putative correlation $\beta _{\textit{max}}=k_1 {\textit{Re}}^{1/3}$ valid at sufficiently high Weber number, we identify the linear relationship between $X$ and $Y$ as $Y=k_1 X + \text{Const.}$ (the constant may be ${\textit{Re}}$ -dependent), and we fit the value $k_1=1.24$ .

Figure 7.

Blue curve: dependence of $\beta _{\textit{max}}=\mathcal{R}_{\textit{max}}/R_0$ on Weber number, at fixed Reynolds number and contact angle $\vartheta _a=90^\circ$ (CAMPN1). Apart from the varying Weber number (through the parameter $\gamma$ ), the parameters are as given in table 2. Black solid lines: the bounds (4.34). Black dotted lines: the correlation (5.5).

Figure 8.

Collapse of results on to a single curve. Circles: varying ${\textit{We}}$ at fixed $\vartheta _a$ (CAMPN1). Squares: varying $\vartheta _a$ at fixed ${\textit{We}}$ (CAMPN2). All other parameters as in table 2. Black solid line: the lower bound given by (4.34). Broken solid line: trend line, giving the slope $Y=k_1X$ and $k_1=1.24$ .

With this value in mind, we revisit the model (4.30) and perform nonlinear least-squares fitting over the parameters $(\xi ,t_0,t_1)$ such that $\pi /(2k_h)=k_1$ is forced to take the value $1.24$ . This then provides us with the given parameter values for $\xi$ , $t_0$ and $t_1$ in (4.31). With these parameter values, we confirm that the DNS results fall comfortably inside the bounds (4.34) implied by the RL model.

In another campaign of simulations (CAMPN3), we vary $\nu$ (hence ${\textit{Re}}$ ) while keeping ${\textit{We}}$ and $\vartheta _a=90^\circ$ fixed. We are able to fit the resulting values of $\beta _{\textit{max}}$ to a curve of the form

(5.5) \begin{align} \beta _{\textit{max}}=k_1 {\textit{Re}}^{1/3}-k_0 (1-\cos \vartheta _a)^{1/2}({\textit{Re}}/{\textit{We}})^{1/2}. \end{align}

The value of $k_1$ is fixed from CAMPN1–CAMPN2 as $k_1=1.24$ ; we therefore view $k_0$ as a single fitting parameter, which is estimated based solely on the results from CAMPN3 ( $k_0=0.80$ ). Crucially, since $k_0\lt k_2$ , the fitted model (5.5) lies within the bounds (4.34) implied by the RL model.

Figure 9.

Observed values of $\mathcal{R}_{\textit{max}}$ (horizontal axis) compared with predicted values of $\mathcal{R}_{\textit{max}}$ using the correlation (5.5). The prior DNS refers to figure 5 in (Shin & Juric Reference Shin and Juric2009), figure 4 in (Gupta & Kumar Reference Gupta and Kumar2011) and figure 3 in Wu & Cao (Reference Wu and Cao2017).

To visualise the results, we plot the values of $\beta _{\textit{max}}$ from CAMPN3 on the horizontal axis, and the predicted values of $\beta _{\textit{max}}$ (from (5.5)) on the vertical axis. Crucially, we superimpose on this plot the actual versus predicted values from CAMPN1–CAMPN2, and the result is a near collapse of the data onto a straight line of slope $45^\circ$ . Results from further simulation campaigns (CAMPN4–CAMPN6) are also plotted on the same graph (figure 9), together with values for the maximum spreading radius computed in prior works by Shin & Juric (Reference Shin and Juric2009), Gupta & Kumar (Reference Gupta and Kumar2011) and Wu & Cao (Reference Wu and Cao2017). A full synopsis of the data emanating from all campaigns is given in tables 4 and 5. Notably, Campaign 5 uses the diffuse-interface method (DIM) for the droplet-impact simulations (see Ó Náraigh & Mairal (Reference Ó Náraigh and Mairal2023) for the method); the results are consistent with the other campaigns which use the VOF method in OpenFOAM.

Table 4.

Summary of the results of the simulation campaigns. The asterisk indicates a simulation where the predicted value of $\mathcal{R}_{\textit{max}}$ is negative, corresponding to a low Weber number where the correlation (5.5) does not apply.

Table 5.

Summary of the results of the simulation campaigns (continued). The asterisk indicates a simulation where the predicted value of $\mathcal{R}_{\textit{max}}$ is negative, corresponding to a low Weber number where the correlation (5.5) does not apply. CAMPN5 uses the DIM for the numerical simulations.

In summary, only the results from CAMPN1–CAMPN3 are used to fit the model to the data; thereafter, the model is used predictively, to explain the results of CAMPN4–CAMPN6. The absence of a theoretical foundation for obtaining the values of $k_1$ and $k_0$ is a limitation of the present approach. Nevertheless, the overall result is the validation of the universal scaling behaviour (1.2), in the RL regime of droplet spreading.