3.1. Solution of the amplitude equation

Before looking in more detail at solutions, let us first observe that (3.1) can be non-dimensionalised by introducing dimensionless time $\tilde {\tau } = \sqrt {A}\tau$ and amplitude $\tilde {y} = y/|y_0|$, where $y_0$ is the amplitude at some initial time $\tau _0 = t_i - t_0$, which then leads to the one parameter equation

(3.4)\begin{equation} \frac{\textrm{d}^2\tilde{y}}{\textrm{d}\tilde{\tau}^2} + \left(1 - \frac{B}{\tilde{\tau}^2}\right) \tilde{y} = 0. \end{equation}

Clearly, this equation has two distinct limits, namely an early time limit, $\tilde {\tau } \gg \sqrt {B}$, where the time-dependent last term is negligible and the equation reduces to that of a harmonic oscillator with dimensionless angular frequency 1 (corresponding to a dimensional angular frequency $\sqrt {A}$), and a late time limit, $\tilde {\tau } \ll \sqrt {B}$, just before impact at $\tilde {\tau } = 0$ where the last term in (3.4) is dominant

(3.5)\begin{gather} \frac{\textrm{d}^2\tilde{y}}{\textrm{d}\tilde{\tau}^2} + \tilde{y} \approx 0\quad\textrm{for}\ \tilde{\tau} \gg \sqrt{B}, \end{gather}

(3.6)\begin{gather}\frac{\textrm{d}^2\tilde{y}}{\textrm{d}\tilde{\tau}^2} - \frac{B}{\tilde{\tau}^2} \tilde{y} \approx 0\quad\textrm{for}\ \tilde{\tau} \ll \sqrt{B}. \end{gather}

Allowing for complex solutions, we may write the general solution to (3.5) as

(3.7)\begin{equation} \tilde{y} \approx \exp[\textrm{i}\psi + \textrm{i}\tilde{\tau}] \quad(\textrm{for}\ \tilde{\tau} \gg \sqrt{B}), \end{equation}

where the prefactor of the exponential has been set to $1$, consistent with the definition of the dimensionless amplitude $\tilde {y}$, and $\psi$ is some constant initial phase. Also the second equation (3.6) can be solved straightforwardly by noting that one should try a power-law solution of the form $\tilde {y} \sim \tilde {\tau }^\kappa$, which when inserted into (3.6) leads to the algebraic relation $\kappa (\kappa -1) - B = 0$, which has solutions $\kappa _\pm = \tfrac {1}{2}(1\pm \sqrt {4B+1})$ with which

(3.8)\begin{equation} \tilde{y} \approx \tilde{y}_+\tilde{\tau}^{\kappa_+} + \tilde{y}_-\tilde{\tau}^{\kappa_-} \approx \tilde{y}_-\tilde{\tau}^{\kappa_-}\quad\textrm{with}\ \kappa_\pm{=}\tfrac{1}{2}(1\pm \sqrt{4B+1})\quad(\textrm{for}\ \tilde{\tau} \ll \sqrt{B}). \end{equation}

Clearly, as $\tilde {\tau }\to 0$, the solution is dominated by the second term, for which the exponent $\kappa _-$ is always negative. It is good to note that one may in fact also write down analytical solutions to (3.4) in terms of Bessel functions. For completeness, these are provided in Appendix B, but they are of little use for the purposes of this article.

In figure 3(a) we plot the real part $\tilde {y}_r$ of the dimensionless amplitude $\tilde {y}$ as a function of dimensionless time $\tilde {\tau }$ for $B=1$ (red curve). Clearly, for large $\tilde {\tau }$ the solution is oscillatory and is fitted nicely by the asymptotic solution (3.7) (dotted black curve). As expected, the solution starts to diverge when $\tilde {\tau } \approx \sqrt {B} = 1$, which may be better observed in the inset where the same data are plotted doubly logarithmically. Clearly, $\tilde {y}$ diverges as a power law, which is well fitted by the dotted black line, which corresponds to (3.8), with $\tilde {y}_- = 1.4$. In the same plots, the blue curve corresponds to the magnitude $|\tilde {y}(\tilde {\tau })|$ of the complex solution with the same boundary conditions, which in this case happens to coincide with $\tilde {y}_r$. This is not always the case, since the divergence may happen in any direction in the complex plane, depending on initial conditions, as evidenced by figure 3(b).

Figure 3.

(a) Solution of the amplitude equation for $B=1$. The main figure shows the amplitude $\tilde {y}$ as a function of time remaining to impact $\tilde {\tau }$. The red curve shows the real part of the solution $\tilde {y}_r(\tilde {\tau })$ with boundary conditions $\tilde {y}(20) = 1$, $\dot {\tilde {y}}(20) = 0$ at $\tilde {\tau }_0=20$, whereas the blue curve shows the magnitude $|\tilde {y}(\tilde {\tau })|$ of the complex solution starting from the same boundary conditions. The dotted black line corresponds to the approximate solution in the large-$\tilde {\tau }$ limit (3.5), again starting from the same point. The inset shows the same data in a doubly logarithmic plot and the dotted black line the approximate solution $\tilde {y} =\tilde {y}_-\tilde {\tau }^{\kappa _-}$ in the small-$\tilde {\tau }$ limit (3.6) with $\tilde {y}_- = 1.4$. (b) Projection onto the complex plane $\tilde {y} = \tilde {y}_r + \textrm {i}\tilde {y}_i$ of solutions of the amplitude equation starting from the same boundary conditions as in (a), but taken at different initial times $\tilde {\tau }_0 = 20 - n{\rm \pi} /3$ for $n=0,\ldots,5$. Clearly, after spending time on the unit circle (corresponding to oscillations), divergence may occur in any direction in the complex plane. The curve closest to the positive $\tilde {y}_r$-axis corresponds to the one plotted in figure (a), i.e. $\tilde {\tau }_0 = 20$.

For the purpose of our stability analysis it is sufficient to look at the magnitude $|\tilde {y}(\tilde {\tau })|$. Clearly, for $\tilde {\tau } \gg \sqrt {B}$, we find that $|\tilde {y}(\tilde {\tau })| \approx 1$ (see (3.7)), whereas close to impact ($\tilde {\tau } \ll \sqrt {B}$) we have $|\tilde {y}(\tilde {\tau })| \approx |\tilde{y}_{-}|\tilde{\tau}^{\kappa_{-}}$ (from (3.8)). By simply matching both asymptotic solutions at the cross-over point $\tilde {\tau } = \sqrt {B}$, we obtain the unknown $|\tilde{y}_{-}|$ and arrive at a continuous approximate solution $|\tilde {y}|_{appr} \approx |\tilde {y}|$

(3.9)\begin{equation} \left|\tilde{y}\right|_{appr}(\tilde{\tau}) = \begin{cases} 1 & \textrm{for}\ \tilde{\tau}\geqslant\sqrt{B},\\ \left(\dfrac{\tilde{\tau}}{\sqrt{B}}\right)^{-{1}/{2}(\sqrt{4B+1}-1)} & \textrm{for}\ \tilde{\tau}<\sqrt{B}. \end{cases} \end{equation}

To understand how the solution depends on the parameter $B$, we plot the magnitude of the dimensionless amplitude $|\tilde {y}|$ (solid lines) together with the approximate solution $|\tilde {y}|_{appr}$ (dashed lines) for several different values of $B$ in figure 4(a), all starting with $\tilde {y}(\tau _0) = 1$ and $\dot {\tilde {y}}(\tau _0) = 0$ at $\tilde {\tau }_0 = 1000$. Clearly, all curves start to diverge for $\tilde {\tau }$ below $\sqrt {B}$, but the rate at which this happens strongly depends on $B$, consistent with the above approximate solution. Note that the approximate solutions only fairly coincide with the exact ones: whereas the power-law exponent is highly accurate, one observes a shift between $|\tilde {y}|$ and $|\tilde {y}|_{appr}$ that is especially pronounced for large $B$ (inset of figure 4a). Here, it is good to note that, although a matched asymptotic expansion can be constructed, it is algebraically complex and not converging quickly such that the approximate solution of (3.9) is preferred. Also, the cross-over point $\tilde {\tau }<\sqrt {B}$ marks exactly that point in time at which the solution can start to grow, i.e. where the time-dependent squared frequency becomes negative.

Figure 4.

(a) Doubly logarithmic plot of the absolute amplitude $|\tilde {y}|$ as a function of time $\tilde {\tau }$ for 5 different values of $B = 0.01$, 0.1, 1, 10, 100 (solid lines), in a close up near to the asymptotic early-time solution ($\tilde {\tau }>\sqrt {B}$). The dashed lines represent the approximate solution (3.9), whereas the horizontal black dashed-dotted line represents a threshold value $p=1.5$ (i.e. $\log _{10} p = 0.176$). Note that the differences between the full and approximate solutions are small for this value of $p$. The inset contains the same data plotted for a larger range in the vertical coordinate to indicate that all data diverge as a power law for $\tilde {\tau }\to 0$, with the strongly $B$-dependent exponent $\kappa _- = (\sqrt {4B+1}-1)/2$. The vertical dotted lines indicate the cross-over point $\tilde {\tau }=\sqrt {B}$. (b) Doubly logarithmic plot of the absolute amplitude growth rate $\tilde {\varGamma }$ (solid lines) vs $\tilde {\tau }$. The dashed lines represent the approximate solution $\tilde {\varGamma }_{appr}$ and the dashed dotted black line corresponds to the onset value of $\tilde {\varGamma }_{appr}$ at $\tilde {\tau } = \sqrt {B}$. All quantities are dimensionless (see main text).

3.2. Absolute amplitude growth rate

A second quantity that is worth looking at is the absolute amplitude growth rate $\tilde {\varGamma }$, which we define as the rate at which the absolute amplitude grows

(3.10)\begin{equation} \tilde{\varGamma} \equiv{-} \frac{1}{|\tilde{y}|}\frac{\textrm{d}|\tilde{y}|}{\textrm{d}\tilde{\tau}}, \end{equation}

where the minus sign guarantees that $\tilde {\varGamma } \geqslant 0$ in the divergent region ($\tilde {\tau } \downarrow 0$). This definition coincides with the growth rate as it is commonly defined in stability problems for the case of exponential growth (and for which it is constant). The absolute amplitude growth rate is plotted in figure 4(b) (solid lines) for the data found in figure 4(a), where one observes that for $\tilde {\tau } \gg \sqrt {B}$ the growth rate is very small, and subsequently rapidly increases for $\tilde {\tau } \gtrsim \sqrt {B}$, until, for $\tilde {\tau }< \sqrt {B}$, the solution converges to a power law with exponent $-1$. This last feature is easily understood by computing the absolute amplitude growth rate from the approximate solution $|\tilde {y}|_{appr}$

(3.11)\begin{equation} \tilde{\varGamma}_{appr}(\tilde{\tau}) ={-}\dfrac{1}{|\tilde{y}|_{appr}}\dfrac{\textrm{d}|\tilde{y}|_{appr}}{\textrm{d}\tilde{\tau}} =\begin{cases} 0 & \textrm{for}\ \tilde{\tau}\geqslant\sqrt{B},\\ \dfrac{\sqrt{4B+1}-1}{2\tilde{\tau}} & \textrm{for}\ \tilde{\tau}<\sqrt{B}, \end{cases} \end{equation}

which indeed corroborates the observed power-law exponent. In fact, the approximate absolute amplitude growth rate $\tilde {\varGamma }_{appr}$ (dashed lines) excellently coincides with the numerical solution $\tilde {\varGamma }$ in the late regime ($\tilde {\tau }<\sqrt {B}$), which once more confirms that the power-law exponent for $\tilde {y}$ in this limit is well predicted by the approximate solution. A final feature that is worth pointing to in figure 4(b) is that initially $\tilde {\varGamma }$ appears to oscillate with a small amplitude (${\sim }10^{-4}$) for $B\geqslant 1$. Although this detail is not relevant for the stability question we are addressing, it is good to note that this is not some numerical inaccuracy, but a property of the equations, as is shown in Appendix C where we will derive a set of equations for $\tilde {\varGamma }$, which is not only elegant, but also easier to solve numerically.

3.3. Determining stability criteria

An important question we will need to address is to determine the onset time of the instability, i.e. the time at which the solution $|\tilde {y}(\tilde {\tau })|$ becomes unstable. We will now set out to define and determine this threshold value.

From the above discussion it is clear that $|\tilde {y}(\tilde {\tau })|$ is stable at early times, and starts to diverge when $\tilde {\tau }$ becomes smaller than $\sqrt {B}$. What is complicating this task, however, is that the divergence is not uniform, but strongly depends of the value of the parameter $B$. Naively, one could just set instability onset to the time below which the (approximate) dimensionless amplitude first starts to rise above $1$, i.e. at $\tilde {\tau } = \sqrt {B}$. This would, however, seriously over-predict the onset time for small values of $B$. e.g. looking at the blue curve corresponding to $B=0.1$ in figure 4(a) one observes that $\sqrt {B} \approx 0.316$, whereas the solution only surpasses the dashed-dotted black line, indicating a growth of the amplitude to just $1.5$ times its original size, at $\tilde {\tau } \approx 3.78\times 10^{-3}$, which is almost two orders of magnitude smaller than $\sqrt {B}$.

The need to distinguish these slow growing modes with small $B$ forces us to define the onset time by setting a threshold value $p$ for the absolute amplitude $|\tilde {y}|$ to determine the onset time and define it as

(3.12)\begin{equation} |\tilde{y}|(\tilde{\tau}_{ons}) = p. \end{equation}

This does more justice to the dynamics of the amplitude equation, but comes at the expense of having to introduce an additional parameter to the problem, namely, the threshold value $p$. Also, setting a threshold value to the full, numerical solution $|\tilde {y}|$ will not lead to a practical algebraic expression. However, using the approximate solution (3.9), we can directly solve $|\tilde {y}|(\tilde {\tau }_{ons}) = p$ as

(3.13)\begin{equation} \tilde{\tau}_{ons, appr} = \sqrt{B} \,p^{{-}2/(\sqrt{4B+1}-1)} = \sqrt{B}\exp\left[-\frac{2\log(p)}{\sqrt{4B+1}-1}\right]. \end{equation}

Coincidentally, at a reasonable threshold value $p=1.5$ (horizontal dashed-dotted black line in figure 4a) the discrepancy between $|\tilde {y}|$ and $|\tilde {y}|_{appr}$ is observed to be very small. Thus, for such threshold values we may expect the approximate onset times to lie close to the numerical ones.

In figure 5(a) we plot the numerically determined onset times $\tilde {\tau }_{ons}$ (solid lines) together with the approximation $\tilde {\tau }_{ons, appr}$ from the above equation (dashed lines) for different values of the threshold $p$. Clearly, the approximate solution gives an accurate description of the actual onset time, especially for intermediate threshold values, $p\approx 1.5$). For the smallest $p$ value plotted, there is a clear shift of the numerical result to slightly larger onset times.

Figure 5.

(a) Onset time $\tilde {\tau }_{ons}$ as a function of $B$ for three different threshold values $p = 1.1$, $1.5$, $2.0$, both numerically from the direct solution $\tilde {y}(\tilde {\tau })$ (solid lines) and using the approximate expression (3.13) (dashed lines). Note that the undulations that are visible, e.g. in the solution for $p=1.1$, are due to initial conditions, which are here taken to be $\tilde {y}(\tilde {\tau }_0) =1$ and $\dot {\tilde {y}}(\tilde {\tau }_0) =0$, for $\tilde {\tau }_0 = 4B$. The undulations may be suppressed taking larger values for $\tilde {\tau }_0$. The dashed dotted line represents the semi-classical result $\tilde {\tau }_{sclas} = \sqrt {B}$. (b) Multiplication factor $M$ as a function of $B$. Again, the dashed-dotted line denotes the semi-classical result $M_{sclas} = \sqrt {B}$. Also plotted are the numerically determined maximum growth factors $|\tilde {y}|_{max}$ for two representative parameter settings corresponding to wave (black dashed-dotted line) and disk impact (red dashed-dotted line). See Appendix D for details of this estimate.

Equally significant is the observation that the onset time quantitatively depends on the choice of the threshold value $p$. This can be easily understood when looking back at figure 4(a), where one observes that, especially for smallest values of $B$, the power-law exponent in the divergent region, $\kappa _- = \tfrac {1}{2}(1-\sqrt {4B+1})$, becomes vanishingly small, corresponding to large shifts in onset time for varying $p$. This is an issue that especially impacts this time-dependent analysis, but is also worthy of consideration for the classical analysis, as discussed in the next subsection.

Given that the onset time depends on the choice of the threshold value, it may be good to look for an additional quantity that may provide some insight in the speed with which the solution diverges. Naturally, one may be tempted to use the absolute amplitude growth rate $\tilde {\varGamma }$, but this quantity still depends on time in the divergent regime, and additionally is also affected by the non-dimensionalisation of the time coordinate. One can, however, define another quantity that is not affected by the above two complications, namely the multiplication factor $M = \tilde {\varGamma }\tilde {\tau }$, which is simply the product of $\tilde {\varGamma }$ and the amount of time remaining until impact, $\tilde {\tau }$. Since asymptotically, $\tilde {\varGamma }\sim 1/\tilde {\tau }$, once in the diverging regime, $M$ is not expected to depend on time and, given the excellent agreement of the approximate and actual expressions in this regime, one may simply compute $M$ as

(3.14)\begin{equation} M = \tilde{\varGamma}\,\tilde{\tau} \approx \frac{\sqrt{4B+1}-1}{2\tilde{\tau}}\,\tilde{\tau} = \frac{1}{2}(\sqrt{4B+1}-1)\quad (\textrm{for}\ \tilde{\tau}\ll \sqrt{B}). \end{equation}

This expression is plotted in figure 5(b). With the definition (3.10) of $\tilde {\varGamma }$, one may adopt $M$ as a measure for the relative increase $\varDelta |\tilde {y}|/|\tilde {y}|$ in the diverging time interval $\tilde {\tau }_{ons}$, but of course cannot be identified with it, since the true amplitude is diverging with time $\tilde {\tau }$. This implies that, even within the linear context, $M$ will be a lower bound on the actual increase in amplitude, which may even be orders of magnitude larger than $M$.

This is confirmed by numerically solving for the maximally possible amplitude that the solution can obtain, realising that the amplitude $|y(\tau )|$ can never become larger than the remaining gap width ($U_0\tau$) between the solid and the free surface. This condition leads to a maximum growth factor $|\tilde {y}|_{max} = |y|_{max}/y_0$ of the original amplitude $y_0$, which are also plotted in figure 5(b) for some typical parameter values. More details on the calculation can be found in Appendix D and in the next section, where the two curves will be discussed in some greater detail. It suffices to note that the more realistic estimate of growth $|\tilde {y}|_{max}$ is orders of magnitude larger than $M$, at the expense of explicitly depending on all relevant system parameters, including the initial amplitude $y_0$.

Finally, the divergence of the multiplication factor $M$ for large $B$ (i.e. small wavelength $\lambda$) may look puzzling at first sight, but one should realise that the amount of time in which this divergence takes place scales with $1/\sqrt {A} \sim \lambda ^{3/2}$ and consequently becomes vanishingly small for small $\lambda$.

3.4. Stability criteria and (semi-)classical stability analysis

The final subject to discuss in this section is how all of the above relates to a traditional stability analysis, where in (3.4) the factor $\tilde {\varOmega }^2 \equiv 1-B/\tilde {\tau }^2$ would be time independent. Before looking at the conditions under which (3.4) may behave classically, let us first naively neglect the time dependence of the equations and write down the results that we may classically expect. This corresponds to assuming (wrongly, as we know) that $\tilde {\varOmega }^2$ changes slowly with time, a situation which we will call this the semi-classical picture. First of all, (3.4) would become unstable as soon as $\tilde {\varOmega }^2$ becomes zero, leading to the semi-classical onset time

(3.15)\begin{equation} \tilde{\tau}_{ons,sclass} = \sqrt{B}, \end{equation}

which result we have plotted in figure 5(a) as the black dashed-dotted line. The semi-classical absolute amplitude growth rate is simply zero above $\tilde {\tau }_{ons, sclass}$ and $\sqrt {-\tilde {\varOmega }^2}$ below it, i.e.

(3.16)\begin{equation} \tilde{\varGamma}_{sclass}(\tilde{\tau}) = \begin{cases} 0 & \textrm{for}\ \tilde{\tau}\geqslant\sqrt{B}\,,\\ \dfrac{\sqrt{B-\tilde{\tau}^2}}{\tilde{\tau}} \approx \dfrac{\sqrt{B}}{\tilde{\tau}} & \textrm{for}\ \tilde{\tau}<\sqrt{B}, \end{cases} \end{equation}

where the last approximate equality holds at the very last stages ($\tilde {\tau } \ll \sqrt {B}$). Clearly, using this last expression we may write down a semi-classical multiplication factor as

(3.17)\begin{equation} M_{sclass} = \tilde{\varGamma}_{sclass}\,\tilde{\tau} = \sqrt{B-\tilde{\tau}^2} \approx \sqrt{B}. \end{equation}

It is this expression that we have plotted together with the full expression (3.14) in figure 5(b). It needs to be stressed that this semi-classical picture does not represent the result of traditional time-independent stability analysis, but what one would get if one would naively apply its results to this time-dependent case. Based on the above one could say that the full expressions converge to the semi-classical ones in the limit $B \gg 1$, and drastically diverge when $B \lesssim 1$. The latter is therefore the region where one would expect to obtain results that strongly diverge from the time-independent picture. Most specifically, for small $B$ solutions are seen to be considerably stabilised with respect to the time-independent case: the onset time is strongly suppressed (which implies that they become unstable at a much later point in time, i.e. at smaller $\tilde {\tau }$) and also the multiplication factor is much smaller (${\sim }B$ as compared with ${\sim }\sqrt {B}$ in the semi-classical case).

Finally, we ask ourselves the question under what conditions taking a time-independent approach would be justified. This is expected to be the case when $\tilde {\varOmega }^2$ changes sufficiently slowly as a function of time, such that the dynamics is only minorly influenced. This implies that the time derivative of $\tilde {\varOmega }^2$ must be sufficiently small in the region where the solution becomes unstable, i.e. close to onset time

(3.18)\begin{equation} \left.\frac{\textrm{d}(\tilde{\varOmega}^2)}{\textrm{d}\tilde{\tau}}\right|_{\tilde{\tau}_{ons}} = \left.\frac{2B}{\tilde{\tau}^3}\right|_{\tilde{\tau} = \sqrt{B}}=\frac{2}{\sqrt{B}}\ll1\quad\Rightarrow\quad B\gg4, \end{equation}

confirming that a time-independent approach may be justified for very large values of $B$.

4.1. Wave impact against a vertical wall

Starting with the case where a wave impacts a vertical wall, we first notice that since gravity acts parallel to the wall, we may set $g=0$ in the expression (3.2) for $A$, such that $A = \alpha _3k^3$ and $B=\beta _2k^2$, with which

(4.2)\begin{align} \tau_{ons,w}(k) &= \sqrt{\frac{\beta_2}{\alpha_3}}k^{{-}1/2}\exp\left[-\frac{2\log(p)}{\sqrt{4\beta_2k^2+1}-1}\right] \nonumber\\ &= \sqrt{w(\delta)(1+\delta)\frac{\rho_l}{\sigma}}\,\frac{q_0}{U_0}k^{{-}1/2}\exp \left[-\frac{2\log(p)}{\sqrt{4w(\delta)\dfrac{q_0^2}{U_0^2}k^2+1}-1}\right], \end{align}

where in the last equation we have inserted the expressions (3.2) and (3.3) for $\alpha _3$ and $\beta _2$, respectively, and where we have defined the function $w(\delta )$ as

(4.3)\begin{equation} w(\delta) = \frac{\delta}{(1+\delta)^2} \end{equation}

for notational convenience. Note that, in the limit $\delta \ll 1$, which is generally satisfied, $w(\delta ) \approx \delta$.

To plot the expressions for the onset time and multiplication factor, we take the air–water interface at $20\,^\circ \textrm {C}$ and standard atmospheric pressure as an example. This fixes the liquid and gas densities and the interfacial tension to $\rho _l = 998\ \textrm {kg}\ \textrm {m}^{-3}$, $\rho _g = 1.20\ \textrm {kg}\ \textrm {m}^{-3}$ and $\sigma = 0.073\ \textrm {N}\ \textrm {m}^{-1}$. In figure 6 we plot the onset time $\tau _{ons}$ as a function of the wavelength $\lambda = 2{\rm \pi} /k$ for different values of $q_0/U_0$, together with the semi-classical result

(4.4)\begin{equation} \tau_{ons,w,sclass}(k) = \sqrt{\frac{B}{A}} = \sqrt{\frac{\beta_2}{\alpha_3}}\frac{1}{\sqrt{k}}= \sqrt{w(\delta)(1+\delta)\frac{\rho_l}{\sigma}}\,\frac{q_0}{U_0}\frac{1}{\sqrt{k}}. \end{equation}

Clearly, whereas the semi-classical result $\tau _{ons,w,sclass}(\lambda )$ diverges for $\lambda \to \infty$, indicating that long wavelengths are always unstable, the time-dependent result $\tau _{ons,w}(\lambda )$ exhibits a maximum at some finite, marginally stable wavelength $\lambda _{marg}$, for which the liquid interface first becomes unstable as the wave approaches the wall (i.e. for $\tau \to 0$). This marginal wavelength is thus solely selected as a consequence of the time dependence of the process.

Figure 6.

Stability diagram for a wave impacting a vertical wall. Lower (main) panel, doubly logarithmic plot of the onset time $\tau _{ons}$ vs wavelength $\lambda$ for three different values of the gas flow rate parameter $q_0/U_0$ (solid lines; see legend), with $p=1.5$. The dashed lines represent the semi-classical result $\tau _{ons,sclass}$ and vertical dotted lines indicate for which $\lambda$ the value $B=1$ is reached. The dashed-dotted black curve indicates the location of the marginal onset time $\tau _{marg}$ as a function of the marginal wavelength $\lambda _{marg}$. Upper panel, multiplication factor $M$ as a function of the wavelength $\lambda$ for the same three cases as in the lower panel (solid lines), together with the semiclassical result $M_{sclass}$ (dashed lines). The dashed-dotted black line corresponds to the numerically calculated maximum growth factor $|\tilde {y}|_{max}$ (see Appendix D) corresponding to the largest gas flow rate parameter $q_0/U_0 = 0.50$ m and an impact velocity $U_0 = 6.0\ \textrm {m}\ \textrm {s}^{-1}$. Densities and interfacial tension correspond to those of water and air at $20\,^\circ \textrm {C}$ and atmospheric pressure ($\rho _l = 998\ \textrm {kg}\ \textrm {m}^{-3}$, $\rho _g = 1.20\ \textrm {kg}\ \textrm {m}^{-3}$ and $\sigma = 0.073\ \textrm {N}\ \textrm {m}^{-1}$).

To find the marginal wavelength we need to determine the location of the maximum by solving $\textrm {d}\tau _{ons,w}/\textrm {d}k = 0$, which is done in Appendix E and results in

(4.5)\begin{equation} \lambda_{marg} = 2{\rm \pi} \frac{m(\delta)}{g(p)}\,\frac{q_0}{U_0}, \end{equation}

where the functions $g(p)$ and $m(\delta )$ are defined in (E6) and (E11). Here, note that $m(\delta ) \approx \sqrt {\delta }$ in the case that $\delta \ll 1$. Reinserting the expression for $k_{marg} = 2{\rm \pi} /\lambda _{marg}$ into (4.2), provides us with an expression for the marginal onset time $\tau _{marg} = \tau _{ons,w}(k_{marg})$

(4.6)\begin{equation} \tau_{marg} = h(p)\,n(\delta)\left(\frac{\rho_l}{\sigma}\right)^{1/2}\left(\frac{q_0}{U_0}\right)^{3/2}, \end{equation}

where $h(p)$ and $n(\delta )$ are defined in (E8) and (E12). The line connecting the different $(\lambda _{marg},\tau _{marg})$-pairs for varying gas flow rate parameter $q_0/U_0$ is plotted as the dashed-dotted line in figure 6, through the marginal wavelengths and onset times of the three plotted curves, which are indicated by the asterisks.

From the above equations it follows that the marginal wavelength only depends on the gas flow rate parameter $q_0/U_0$, to which it is proportional, and on the density ratio, through $\lambda _{marg} \sim \delta ^{1/2}(q_0/U_0)$ using the small $\delta$ approximation provided in (E11). There is no surface tension dependence. In contrast, the onset time does depend on both surface tension and gas flow rate, namely as $\tau _{marg} \sim \delta ^{3/4}(\rho _l/\sigma )^{1/2}(q_0/U_0)^{3/2}$, now using the small $\delta$ approximation for $n(\delta )$ (E12). As a result, the gas velocity at the marginal onset time, $U_{g,{marg}} = U_g(\tau _{marg})$, can be estimated as $U_{g,{marg}} \sim \delta ^{-3/4}(\rho _l/\sigma )^{-1/2}(q_0/U_0)^{-1/2}$, which implies that the marginal gas speed decreases with increasing $q_0/U_0$.

The multiplication factor $M$ decreases with increasing wavelength, which can be traced back directly to the fact that $B \sim k^2$. For the threshold value $p = 1.5$ for which the stability diagram has been plotted $M$ is larger than unity for wavelengths smaller than $\lambda _{marg}$ but above this marginal value rapidly decreases to negligibly small values. This is consistent with the sharp cutoff of the onset time with increasing $\lambda$ where $\tau _{ons}$ drops by over five orders of magnitude within a decade of $\lambda$.

Turning to the numerically determined maximum growth factor $|\tilde {y}|_{max}$ corresponding to the largest gas flow rate parameter $q_0/U_0 = 0.50$ m (dashed-dotted black line in the upper panel of figure 6) we see that at the marginal onset time $\tau _{marg}$ (black star) the initial amplitude will grow by a factor $|\tilde {y}|_{max} \approx 1\times 10^5$ (${\approx }2.5$ mm for an initial amplitude of $25$ nm). This confirms that the growth rate of the marginal wavelength is sufficient to be observable. Also note that for small wavelengths (i.e. where $M$ diverges) $|\tilde {y}|_{max}$ tends to a constant large value, implying that all small wavelength disturbances may grow approximately by the same amount. It is good to realise that $\lambda _{marg}$ is the first disturbance to start growing, whereas the other modes may only start to grow later (i.e. at smaller values of $\tau$), when the surface has already been deformed with wavelength $\lambda _{marg}$.

4.2. Disk impact on a liquid surface

We subsequently turn to the impact of a disk on a liquid surface, in which case the gravitational term needs to be included in the equations. This implies that $A = \alpha _1k + \alpha _3k^3$ and $B=\beta _2k^2$, with which

(4.7)\begin{align} \tau_{ons,d}(k) &= \sqrt{\frac{\beta_2k^2}{\alpha_1k + \alpha_3k^3}}\exp\left[-\frac{2\log(p)}{\sqrt{4\beta_2k^2+1}-1}\right]\nonumber\\ &= \sqrt{\frac{\delta}{1+\delta}}\,\frac{R_0k}{2\sqrt{(1-\delta)gk + \dfrac{\sigma}{\rho_l}k^3}}\exp\left[-\frac{2\log(p)}{\sqrt{w(\delta)R_0^2k^2+1}-1}\right], \end{align}

again using (3.2) and (3.3) for $\alpha _1$, $\alpha _3$ and $\beta _2$ to arrive at the last expression, and identifying $q_0/U_0 = R_0/2$.

Again we use the example of an air–water interface at $20\,^\circ \textrm {C}$ and standard atmospheric pressure to plot the onset time as a function of wavelength in figure 7. It is compared with the semi-classical result which in this case equals

(4.8)\begin{equation} \tau_{ons,d,sclass}(k) = \sqrt{\frac{B}{A}} = \sqrt{\frac{\beta_2k^2}{\alpha_1k + \alpha_3k^3}} =\frac{1}{2} \sqrt{\frac{\delta}{1+\delta}}\,\sqrt{\frac{\rho_l}{(1-\delta)\rho_lgk + \sigma k^3}}R_0 k. \end{equation}

For disk impact, the semi-classical result is non-monotonic (as can be seen in figure 7 where it is plotted for three different values of the disk radius $R_0$) and attains a maximum for the marginal wavelength, which is calculated from the condition $\textrm {d}\tau _{ons}/\textrm {d}k = 0$, which leads directly to

(4.9)\begin{equation} \lambda_{marg,sclass} = 2{\rm \pi}\sqrt{\frac{\alpha_3}{\alpha_1}} = 2{\rm \pi}\sqrt{\frac{\sigma}{(1-\delta)\rho_l g}}. \end{equation}

This is $2{\rm \pi}$ times the capillary length, and as a consequence, independent of $q_0/U_0=R_0/2$. As may be expected from this fact, it is indeed exactly the same result as one obtains in a classical Kelvin–Helmholtz stability analysis, i.e. with a constant gas flow over the interface. The corresponding semi-classical marginal onset time is obtained by inserting $k_{marg,sclass}$ into (4.8)

(4.10)\begin{equation} \tau_{marg,sclass} = \frac{\beta_2^{1/2}}{(4\alpha_1\alpha_3)^{1/4}} = \frac{1}{4}\sqrt{\frac{2\delta}{1+\delta}}\left(\frac{\rho_l}{(1-\delta)\sigma g}\right)^{1/4}R_0. \end{equation}

In the limit of small $\delta$ we have $\tau _{marg,sclass} \sim \delta ^{1/2}(\rho _l/(\sigma g))^{1/4}R_0$. Clearly, also in this disk impact case the onset time grows with $q_0/U_0 = R_0/2$, but not as fast as in the wave impact case, where we do not have a proportionality, but a power law with exponent $3/2$ (4.6). We again estimate the gas velocity at this semi-classical marginal onset time, $U_{g,{marg,sclass}} = U_g(\tau _{marg,sclass})$, can be estimated as $U_{g,{marg,sclass}} \sim \delta ^{-1/2}(\sigma g/\rho _l)^{1/4}$, which, just as in the classical case, is found to be independent of $R_0$.

Figure 7.

Stability diagram for a disk impacting onto a liquid surface. Lower (main) panel, doubly logarithmic plot of the onset time $\tau _{ons}$ vs wavelength $\lambda$ for three different values of the disk radius $R_0 = 2q_0/U_0$ (solid lines; see legend), with $p=1.1$. The dashed lines represent the semi-classical result $\tau _{ons,sclass}$ and vertical dotted lines indicate for which $\lambda$ the value $B=1$ is reached. The vertical dashed-dotted black curve indicates the location of the semi-classical marginal wavelength $\lambda _{marg,sclass}$. Upper panel, factor $M$ as a function of the wavelength $\lambda$ for the same three cases as in the lower panel (solid lines), together with the semiclassical result $M_{sclass}$ (dashed lines). The dashed-dotted red line corresponds to the numerically calculated maximum growth factor $|\tilde {y}|_{max}$ (see Appendix D) corresponding to the middle radius $R_0 = 0.10$ m and an impact velocity $U_0 = 2.0\ \textrm {m}\ \textrm {s}^{-1}$. Densities and interfacial tension correspond to those of water and air at $20\,^\circ \textrm {C}$ and atmospheric pressure ($\rho _l = 998\ \textrm {kg}\ \textrm {m}^{-3}$, $\rho _g = 1.20\ \textrm {kg}\ \textrm {m}^{-3}$ and $\sigma = 0.073\ \textrm {N}\ \textrm {m}^{-1}$).

Above we are discussing the semi-classical limit, which is plotted as the dashed lines in figure 7. When we compare this with the full solution (solid lines), we see that for large values of $R_0$ the semi-classical results for the marginal wavelength and marginal onset time are a good approximation to the true ones, with a cutoff happening for larger wavelengths, but as soon as $R_0$ comes in the centimetre range, the cutoff starts to interfere with the maximum in the semi-classical result (represented by the dashed-dotted black line in figure 7), and the maximum shifts to smaller wavelengths, as can be clearly seen in figure 7 in the curve corresponding to $R_0 = 2.0$ cm.

It is not easy to find a general analytic expression for the marginal wavelength in the case of the disk as it involves solving a fourth order polynomial equation which can be shown ( Appendix E) to lead to a function of the form

(4.11)\begin{equation} \lambda_{marg,d} = \sqrt{w(\delta)}\,R_0\,F\left(\log(p),\frac{\alpha_1\beta_2}{\alpha_3}\right) = \sqrt{w(\delta)}\,R_0\,F\left(\log(p),\frac{w(\delta)\,(1-\delta) \rho_l g R_0^2}{4\sigma}\right), \end{equation}

where $F(\log (p),\zeta ) \sim \zeta ^{-1/2}$ for large $\zeta$ (i.e. large $R_0$) in order for the function to have the correct asymptotics leading to (4.9).

Also in this case, the multiplication factor $M$ decreases with increasing wavelength, which is again connected to $B \sim k^2$, and is even identical to the wave impact case (identifying $q_0/U_0 = R_0/2$). We observe that the multiplication factor $M$ at the location of the maximum in the $\tau _{ons}$-curve for the lowest value of $R_0$ is significantly smaller for the disk impact (figure 7) when comparing with wave impact. This can only partly be traced back to the somewhat smaller value of the threshold ($p=1.1$) used to produce the $\tau _{ons}(\lambda )$ and suggests that due to the reduced growth rate the instability may be suppressed altogether at these smaller disk sizes.

In this case, the numerically determined maximum growth factor $|\tilde {y}|_{max}$ is chosen to correspond to the middle disk size $q_0/U_0 = R_0/2 = 0.050$ m (dashed-dotted red line in the top figure 7). The initial amplitude will grow by a factor $|\tilde {y}|_{max} \approx 1\times 10^2$ at the marginal onset time, which is much less than in the wave impact case and makes observation of the marginal wavelength questionable, leading to the tentative conclusion that the influence of including time dependence for small disk sizes is not only to shift the marginal wavelength to smaller values and to retard the marginal onset time (i.e. to smaller $\tau _{marg}$), but also to suppress the growth rate of the marginally stable wavelength.

4.3. The role of the threshold value $p$

In all of the above analysis, there is a pronounced influence of the choice of the threshold value $p$ on the results presented in this section. Generally speaking, for larger values of the threshold the marginal wavelength becomes smaller. For values of $p$ between $1.1$ and $2.0$ the function $g(p)$ may be approximated accurately by $g(p) \approx (4.6(p-1))^{1/2}$, which implies a factor $3$ change of the marginal wavelength when varying $p$ in that range. Clearly, using values smaller than $p=1.1$ makes little sense, both because the approximation becomes less accurate and a small threshold value is hardly predictive for instability to occur, especially for small values of $B$. On the other limit, taking values larger than $p=2.0$ one may argue that when the amplitude has grown to two times its original size, the instability onset must have occurred in the past. From a comparison of the exact and the approximate solutions with the amplitude equation (cf. figures 4a and 5), one may infer that taking $p=1.5$ in general is a good compromise with sufficient accuracy to have predictive value.

Although the somewhat problematic task of having to choose a threshold value appears to be a feature of this particular time-dependent problem and without doubt is closely related to the fact that the equation parameters themselves are rapidly changing functions of time, it is good to realise that this problem is not completely absent in a classical stability analysis: when one obtains a marginally stable wavelength, then by definition its growth rate is exactly zero for the corresponding setting of the system parameters. So, it is a priori unclear by how much one needs to detune the system parameters to actually observe the instability in an experiment, which may be related to the amount of time one has to observe the instability (especially true for convective instabilities), or the extent to which one is able to control fluctuations in the undisturbed set-up. A classical (although very nonlinear) example that comes to mind is the onset of turbulence in pipe flow.

4.4. Advective terms

Finally, we take two steps back in our analysis. The first step back is to where we needed to introduce a moving reference frame in order to have a valid base state for our analysis, which was moving with a velocity of magnitude $U_l(\tau )$ in the negative $x$-direction (2.5). This implies that we now need to transform the result of our analysis back into the laboratory frame, namely by propagating it with a velocity $U_l(\tau )$ in the positive $x$-direction.

The second step back goes to the point where in order to obtain the simple form of the amplitude equation (2.16), a transformation of the original time-dependent amplitude $\varepsilon (\tau )$ into a new amplitude $\bar {\varepsilon }(\tau )$ was introduced by means of a phase factor. To obtain the original amplitude we need to transform our result $\bar {\varepsilon }(\tau ) ({=}y(\tau ))$ back

(4.12)\begin{equation} \varepsilon(\tau) = \bar{\varepsilon}(\tau)\exp[-\textrm{i}kH(\tau)], \end{equation}

with $H(\tau )$ as defined in (2.15). Using $U_g(\tau ) = q_0/(U_0\tau )$ we can compute this function as

(4.13)\begin{equation} H(\tau) = \frac{2\delta}{1-\delta^2}\int_{-\infty}^t U_g(\tau')\,\textrm{d}t' ={-} \frac{2\delta}{1-\delta^2}\frac{q_0}{U_0}\left[\log(\tau) - \log(\infty)\right]. \end{equation}

Although strictly speaking the term $\log (\infty )$ is divergent, we may fix it at some large value $\log (\tau _s)$ to verify that it only contributes some arbitrary constant phase factor that is immaterial for the result we are after, and further neglect it. Inserting the above expression in (4.12), we find that it represents an oscillatory factor, with a frequency and wave speed that diverge as $\tau$ approaches zero, or quantitatively

(4.14)\begin{equation} c_H(\tau) = \frac{\omega(\tau)}{k} = \frac{H(\tau)}{\tau} = \frac{2\delta}{1-\delta^2}\frac{q_0}{U_0}\frac{\log(\tau)}{\tau}. \end{equation}

Now, we have identified two typical velocities that transform any unstable wave pattern, it will be advected with a velocity $U_l(\tau )$ and modulated with a wave speed $c_H(\tau )$. To estimate their magnitude we will now compare these speeds with the other reference speeds that are at our disposal, namely the gas speed $U_g(\tau )$ and the approach speed of the impacting wave or disk $U_0$, thereby focussing on the marginal onset time $\tau _{marg}$.

Starting with $U_l(\tau )$, we find that its ratio with the gas velocity is constant

(4.15)\begin{equation} \frac{U_l(\tau)}{U_g(\tau)} = \frac{\delta}{1-\delta} \approx \delta \quad(\textrm{for}\ \delta\ll1), \end{equation}

and small, since the gas to liquid density ratio $\delta$ is small. For $c_H(\tau )$ we find similarly

(4.16)\begin{equation} \frac{c_H(\tau)}{U_g(\tau)} = \frac{2\delta}{1-\delta^2}\log(\tau) \approx 2\delta \log(\tau_{marg})\quad(\textrm{for}\ \delta\ll1\ \textrm{and}\ \tau=\tau_{marg}). \end{equation}

Clearly, also in this case, since the dependence on $\tau _{marg}$ is logarithmic, the density ratio dominates the result, and we may conclude that generally the intrinsic velocities that disturb the wave pattern are one to three orders of magnitude smaller than the velocity that creates it, the gas velocity $U_g(\tau )$.

It is, however, good to note that, whereas $U_l(\tau )$ advects the pattern in the same direction as the gas velocity, the wave with speed $c_H(\tau )$ propagates in the opposite direction, since $\log (\tau _{marg}) < 0$. Nevertheless, all are diverging for $\tau \to 0$. Therefore it is appropriate to also compare them with the velocity of approach, $U_0$, at the time of marginal onset. First, we turn to the wave impact case, and compare the advective velocity $U_l(\tau )$ with $U_0$, for which

(4.17)\begin{equation} \frac{U_l(\tau_{marg})}{U_0} = \frac{\delta}{(1-\delta)n(\delta)h(p)}\sqrt{\frac{\sigma}{\rho_lU_0q_0}} \approx \delta^{1/4}\sqrt{\frac{\sigma}{\rho_lU_0q_0}} \quad(\textrm{for}\ \delta\ll1), \end{equation}

where we have used expression (4.6) and its small $\delta$ approximation, and for simplicity took $h(p) \approx 1$. Note that $\delta ^{1/4}$ is not very small (0.2–0.3), but that the argument of the square root is to be interpreted as an inverse Weber number (${We} = \rho _lU_0q_0/\sigma$), where for all but the very smallest impacts, ${We}$ is of the order $100$ or more commonly even larger. That is, usually, the advection velocity $U_l(\tau _{marg})$ is at least one to two orders of magnitude smaller than the impact velocity $U_0$. A similar conclusion may be drawn for the wave speed $c_H(\tau )$ for which

(4.18)\begin{equation} \frac{c_H(\tau_{marg})}{U_0} = \frac{2\delta\log(\tau_{marg})}{(1-\delta^2)n(\delta)h(p)}\sqrt{\frac{\sigma}{\rho_lU_0q_0}} \approx 2\delta^{1/4} \log(\tau_{marg})\sqrt{\frac{\sigma}{\rho_lU_0q_0}} \quad(\textrm{for}\ \delta\ll1), \end{equation}

which is, however, a bit larger than the previous ratio due to the factor $2 \log (\tau _{marg})$, but still typically at least an order of magnitude smaller than one in realistic cases.

Turning to disk impact, we compare the advective velocity $U_l(\tau )$ with $U_0$ at the semi-classical marginal onset time $\tau =\tau _{marg,sclass}$ (4.10)

(4.19)\begin{equation} \frac{U_l(\tau_{marg,sclass})}{U_0} = \frac{\sqrt{2}\,\delta}{(1-\delta)^{3/4}(1+\delta)^{1/2}w(\delta)^{1/2}} \frac{\left(\sigma g /\rho_l\right)^{1/4}}{U_0} \approx \sqrt{2\,\delta} \frac{U_{cap}}{U_0} \quad(\textrm{for}\ \delta\ll1), \end{equation}

where $U_{cap} \equiv (\sigma g /\rho _l)^{1/4}$ is the capillary velocity scale well known from the minimum propagation speed of waves on a liquid interface, to which it is proportional. For the water–air interface we find that $U_{cap} \approx 0.164\ \textrm {m}\ \textrm {s}^{-1}$, and therefore we may conclude that also in this case under typical impact conditions where one would expect instabilities to be of importance ($U_0 > 1.0\ \textrm {m}\ \textrm {s}^{-1}$), the advective velocity is at least an order of magnitude smaller than $U_0$. A similar conclusion can be drawn for the wave speed.

In conclusion, the advection and wave speeds within the emerging patterns are expected to be small: they are a factor $\delta$ smaller than the diverging gas speed and at marginal instability onset time they are at least an order of magnitude smaller than the approach speed $U_0$.

5.1. The influence of gas viscosity

To characterise what happens at the (undisturbed and flat) gas–liquid interface due to the presence of gas viscosity, we imagine solving for a laminar flow that, far away from the interface, is identical to the gas flow we studied in the earlier sections of this work, $U_g(\tau ) = q_0/(U_0\tau )$. We somewhat simplify the problem by assuming that the liquid remains stationary during the development of this boundary layer, which implies that the expected small motions within the liquid are neglected. This turns the interface into a flat and rigid boundary, and we aim at solving the laminar boundary layer flow that forms above this interface. Due to translational invariance we may assume that this flow is described by a flow field of the form $\boldsymbol {u} = u(z,\tau )\boldsymbol {e}_x$, with boundary conditions $u(0,\tau )=0$ and $u(\infty,\tau ) = U_g(\tau )$, as depicted in figure 8(a). Here, $\nu$ represents the kinematic viscosity of the gas and $\boldsymbol {e}_x$ denotes the unit vector in the $x$-direction.

Figure 8.

(a) Schematic of the situation calculated in Appendix F, where a boundary layer of thickness $\delta _{BL}$ is forming above a flat (liquid or solid) interface due to a gas flow of velocity $U_g(\tau )$ that diverges for $\tau \to 0$ in the horizontal ($x$) direction. The velocity in the boundary layer $u(z,\tau )$ is assumed to be parallel to the interface and only to depend on $z$ and $\tau$. (b) Solution of the dimensionless boundary layer profile $Y(\eta ) = u(z,\tau )/U_g(\tau )$ as a function of $\eta = z/\sqrt {\nu \tau }$. The horizontal red dashed line at $\eta _{BL}$ indicates the size of the boundary layer. (c) Twice the boundary layer thickness $2\delta _{BL} = 2\sqrt {\nu \tau }$ together with the gap width $d = U_0\tau$, both as a function of $\tau$. The intersection of the two curves (black asterisk) indicates the time $\tau$ beyond which the system is dominated by gas viscosity. For these curves, $\nu = 15.1\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$ corresponding to air at $20\,^\circ \textrm {C}$ and atmospheric pressure, and $U_0 = 1.0\ \textrm {m}\ \textrm {s}^{-1}$.

If the gas velocity would be constant, the problem would correspond to the famous textbook problem of the impulsively started plate, in a moving reference frame in which the plate is at rest, however, where the analytic solution is an error function. In fact, for the above form of $U_g(\tau )$, we can find an analytic, self-similar solution for the time-dependent problem as well, when we start at some very far initial time $\tau _0$ that does not influence the system at the time $\tau$ of interest any longer

(5.1)\begin{equation} \frac{u(z,\tau)}{U_g(\tau)} = \frac{1}{2}\sqrt{\rm \pi}\,\left(\frac{z}{\sqrt{\nu\tau}}\right)\exp \left[\frac{z^2}{4\nu\tau}\right]\textrm{erfc}\left[\frac{z}{2\sqrt{\nu\tau}}\right] , \end{equation}

which expression is derived in Appendix F and plotted in figure 8(b). Most important for the purpose of this section, however, is that the boundary layer thickness $\delta _{BL}(\tau )$ is given by

(5.2)\begin{equation} \delta_{BL}(\tau) = 13.93 \sqrt{\nu\tau}, \end{equation}

where the prefactor $13.93$ is determined from the point at which $u(z,\tau )=0.99U_g(\tau )$, or, equivalently, $Y(\eta )=0.99$, indicated by the horizontal dashed red line in figure 8(b). Note that approaching the impact moment (i.e. for $\tau \to 0$), the boundary layer decreases in size, which sounds counterintuitive at first sight, for those familiar with the behaviour of viscous boundary layer development for steady asymptotic conditions where the boundary layer thickness is always an increasing function of time. It can, however, be understood intuitively by realising that boundary layer formation needs to keep up with the divergence of the outer flow $U_g(\tau )$, where for equal subsequent steps in time, the steps in velocity diverge as well, which makes the boundary layer thickness shrink in time.

However, even if the boundary layer thickness decreases in time, it does so as the square root of $\tau$, and as a consequence is at some point in time overtaken by the thickness of the gap between the liquid and the solid, which decreases linearly in time, as $d(\tau ) = U_0\tau$. In figure 8(c), we compare twice the boundary layer thickness $\delta _{BL}(\tau )$ with the gap thickness $d(\tau )$ and see that there is the expected intersection point at

(5.3)\begin{equation} \tau_{visc} = 776 \frac{\nu}{U_0^2}, \end{equation}

where for all $\tau <\tau _{visc}$ the viscous boundary layers at the two sides of the gap will touch and a viscous Poiseuille flow will start to form with a rapidly increasing pressure head and decreasing velocity that will prevent the type of instabilities discussed in the previous section from occurring. On the other hand, if an instability develops at a time $\tau$ before $\tau _{visc}$, the fact of its growth is not likely to be hindered substantially by the action of gas viscosity.

A short comment is in place since in the above calculation a no-slip boundary condition has also been assumed on the gas–liquid interface, whereas a no-shear boundary condition would be more appropriate. This is, however, justified as long as the dynamic viscosity in the liquid phase is much larger than that in the gas phase, which is satisfied for many gas–liquid combinations, e.g. air and water at atmospheric conditions.

5.2. The influence of gas compressibility

As the gas flow speed $U_g(\tau )$ is diverging, at some point in time it will get close to sonic values, i.e. it will become equal to the speed of sound in the gas, $c_g$. Taking the Mach number ${Ma}$ of order unity as the threshold, ${Ma} = U_g/c_g = 1$ we find

(5.4)\begin{equation} \tau_{comp} = \frac{q_0}{U_0\,c_g}. \end{equation}

For all $\tau <\tau _{comp}$, the gas velocity in the gap becomes supersonic, and the flow will become choked which will cause the pressure in the air pocket (or centre of the disk) to rise, and simultaneously prevent an instability as from the previous section to occur. As in the above discussed case, if an instability arises for $\tau >\tau _{comp}$ the subsequent choking of the flow may influence but not prevent the development of the instability.

Now let us estimate, for a moderate impact speed $U_0 = 1.0\ \textrm {m}\ \textrm {s}^{-1}$, what the influence of gas viscosity and gas compressibility is on the stability diagram for the wave impact case (subsection 4.1), where we have used the properties of water and air at atmospheric pressure and a temperature of $20\,^\circ \textrm {C}$. In figure 9(a), we again plot the onset time as a function of the wavelength, but we now add the time $\tau _{visc}$ beyond which viscous effects prevent the development of a Kelvin–Helmholtz instability (dashed horizontal black line), together with the times $\tau _{comp}$ beyond which this happens due to choked flow (dotted horizontal lines). Since $\tau _{comp}$ depends on $q_0/U_0$ there are three of those. Clearly, $\tau _{comp}$ is smaller than $\tau _{visc}$ and the maximum $\tau _extrm{marg}$ of the $\tau _{ons}(\lambda )$-curves for all values of $q_0/U_0$, such that we can conclude that gas compressibility is not an issue for the parameter settings of figure 9(a). This is, however, not the case for gas viscosity: for the lowest value of $q_0/U_0$ (${=}0.02$ m) the entire instability onset curve lies within the red shaded area, and the occurrence of an instability is therefore out of the question. For the intermediate value ($q_0/U_0 = 0.10$ m), only the maximum itself lies outside the red shaded area, and it seems doubtful that an instability will occur. Only for the highest value, $q_0/U_0 = 0.50$ m, a substantial part of the instability onset curve lies above the $\tau _{visc}$-threshold and an instability is expected to be observed.

Figure 9.

(a) Doubly logarithmic plot of the onset time $\tau _{ons}$ vs wavelength $\lambda$ for the case of a wave impacting a vertical wall, identical to the one in figure 6, i.e. for the same three different values of the gas flow rate parameter $q_0/U_0$ (solid lines) and $p=1.5$. The horizontal dashed black line represents the threshold time $\tau _{visc}$ beyond which the flow in the gap becomes dominated by viscosity (red shaded area), and the horizontal dotted lines indicate the threshold times $\tau _{comp}$ beyond which the flow is choked, for the three values of $q_0/U_0$, respectively. Here, we also fixed $U_0 = 1.0\ \textrm {m}\ \textrm {s}^{-1}$. (b) Phase diagram, stating in what regions of the control parameter space of the wave impact problem, spanned by the length scale $q_0/U_0$ and velocity scale $U_0$, instabilities may occur, or viscosity or compressibility effects prevent the occurrence of an instability. The dotted black square indicates the region where experimental situations are likely to occur. See text for further explanation. Densities, interfacial tension, gas viscosity and gas speed of sound correspond to those of water and air at $20\,^\circ \textrm {C}$ and atmospheric pressure ($\rho _l = 998\ \textrm {kg}\ \textrm {m}^{-3}$, $\rho _g = 1.20\ \textrm {kg}\ \textrm {m}^{-3}$, $\sigma = 0.073\ \textrm {N}\ \textrm {m}^{-1}$, $\nu = 1.51\times 10^{-5}\ \textrm {m}^2\ \textrm {s}^{-1}$ and $c_g = 343\ \textrm {m}\ \textrm {s}^{-1}$).

Finally, we explore in what regions of the direct control parameter space of the wave impact problem, spanned by the length scale $q_0/U_0$ and velocity scale $U_0$, Kelvin–Helmholtz instabilities may occur and where viscosity or compressibility effects prevent the occurrence of such an instability. To that end, for each point in this parameter space, we look at the ordering of the three time scales $\tau _{marg}$, $\tau _{visc}$ and $\tau _{comp}$, provided by (4.6), (5.3) and (5.4), respectively. e.g. if $\tau _{marg} > \tau _{visc},\tau _{comp}$, then instability is expected to occur. If, however, $\tau _{visc} > \tau _{marg},\tau _{comp}$, viscous effects prevent the development of an instability, etc. The resulting phase diagram is plotted in figure 9(b). Each of the lines is computed by equating two of the three representative time scales, and all three intersect in a triple point, the coordinates of which can be computed as

(5.5)\begin{equation} \left(U_0,\frac{q_0}{U_0}\right) = \left(h(p)\,n(\delta)\,\sqrt{776\,\frac{\rho_l\nu c_g^3}{\sigma}} ,\frac{\sigma}{\rho_l \,(h(p)\,n(\delta)\,c_g)^2}\right). \end{equation}

The lines separate the regions where instabilities are expected, where they are impeded by viscosity, and where they are impeded by gas compressibility. The most interesting area is indicated by the dotted black line in the form of a rectangle, which provides the bounds of situations that may practically occur. That is, for water and air, the competition of gas viscosity and Kelvin–Helmholtz instability dominate the scene, whereas gas compressibility is expected to play a very minor role, simply because it is restricted to the simultaneous incidence of very small $q_0/U_0$ and very large $U_0$, highly unlikely to be found in an actual experimental situation. However, examining (5.5), one observes that for a (e.g. cryogenic) gas with a lower speed of sound the triple point would move towards or even into the area of interest.

The general shape of the phase diagram 9(b) is understood intuitively when one realises that $q_0/U_0$ is the main parameter regulating the instability onset time. That is, for large $q_0/U_0$ instabilities occur before viscosity and compressibility effects become relevant. For small $q_0/U_0$ viscosity and compressibility do become important before instability onset, and, clearly, for small values of the impact speed $U_0$ viscosity dominates, whereas for large impact speeds compressibility is the determining factor. From a more quantitative point of view, it is quite remarkable that gas viscosity is of very substantial importance, and even capable of stabilising interfaces at the considerable length scale of ${\sim }0.10$ m, whereas gas compressibility turns out to play a minor role.

A final word of caution is appropriate because continuum theory and with that all of the above reasoning breaks down at the Knudsen limit, i.e. when the distance between the impacting solid and liquid becomes of the order of the mean free path $\ell$. For an ideal gas, $\ell = k_BT/(\sqrt {2}p\varSigma )$ with $k_B$ Boltzmann's constant and $\varSigma$ the collisional cross-sectional area. For air at $20\,^\circ \textrm {C}$ and atmospheric pressure, $\ell \approx 70$ nm, leading to Knudsen times $\tau _K = \ell /U_0$ which are significantly smaller than the other time scales and are situated in the lower left corner of figure 9(b) and are therefore irrelevant. However, the Knudsen limit may become an issue at lower ambient pressures.