2.1. Governing equations

Figure 3 shows a cylindrical container of radius $\hat {R}_0$ containing fluid whose surface is perturbed in the form of a Bessel mode. The displacement of the perturbed surface is measured from the flat undisturbed level and is given initially by $\hat {\eta }(\hat {r},0) = a_0{\textrm{J}}_0(k\hat {r}) = a_0{\textrm{J}}_0(({l_q}/{\hat {R}_0})\hat {r})$ where ${\textrm{J}}_0$ is the zeroth-order Bessel function of first kind and $l_q$ is explained below (2.6a,b). Surface oscillations arise due to the initial perturbation, and for a sufficiently large value of $\epsilon \equiv a_0({l_q}/{\hat {R}_0})$ (the precise value depends on the non-dimensional number $\sigma$), a jet is formed at the axis of symmetry towards the end of a surface oscillation. We study the IVP corresponding to this surface perturbation under the inviscid, irrotational approximation. Here we assume the liquid to be water and the gas above it to be air. Due to the large density ratio between the two, we may safely ignore the dynamic effect of air assuming that it exerts negligible pressure on the water surface. Thus the interface is treated as a free surface with a pressure jump across it due to surface tension.

Figure 3.

Three-dimensional visualisation of surface perturbation (grey). Surface perturbation imposed at $\hat {t}=0$ leads to a jet at the end of one oscillation for sufficiently large $\epsilon$. The coordinate system is centred at the axis of symmetry and $\hat {z}$ is measured from the unperturbed free surface.

The Laplace equation along with the nonlinear boundary conditions (kinematic boundary condition (2.2) and the Bernoulli equation with surface tension (2.3)) at the free surface and the mass conservation condition (2.4) govern the perturbation velocity potential and the perturbed interface. Equations for these are

(2.1)\begin{gather} \hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}\hat{\phi} = 0, \end{gather}

(2.2)\begin{gather}\frac{\textrm{D}\hat{\eta}}{\textrm{D}\hat{t}} - \left(\hat{\boldsymbol{\nabla}}\hat{\phi}\boldsymbol{\cdot}\boldsymbol{e}_{\hat{z}}\right)_{\hat{z}=\hat{\eta}}=0, \end{gather}

(2.3)\begin{gather}\left(\frac{\partial\hat{\phi}}{\partial \hat{t}} + \frac{1}{2}|\hat{\boldsymbol{\nabla}}\hat{\phi}|^2 + g\hat{z} + \frac{T}{\rho}\hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{n}\right)_{\hat{z}=\hat{\eta}} = 0, \end{gather}

(2.4)\begin{gather}\int_{0}^{\hat{R}_0}\hat{r}\hat{\eta}(\hat{r},\hat{t}) \,\textrm{d}\hat{r} = 0 , \end{gather}

where $T$ is the coefficient of surface tension, $\boldsymbol {e}_{\hat {z}}$ and $\boldsymbol {n}$ are unit normals to the unperturbed and perturbed free surface, respectively, and $\hat {\boldsymbol {\nabla }}$ is the gradient operator to be expressed in axisymmetric, cylindrical coordinates. We propose to solve the nonlinear IVP for (2.1), (2.2) and (2.3) up to $O(\epsilon ^2)$ in cylindrical coordinates with an imposed surface perturbation (see figure 3 and (2.6a,b)). Appropriate boundary conditions are imposed on the velocity potential $\hat {\phi }$ at the axis of symmetry ($\hat {r}=0$) and the wall of the container at $\hat {r}=\hat {R}_0$ (no-penetration). In addition, we impose symmetry conditions on the free-surface at $\hat {r}=0$ as well as the free-end edge boundary condition at $\hat {r}=\hat {R}_0$ through (2.5a,b);

(2.5a,b)\begin{gather} \left(\frac{\partial \hat{\phi}}{\partial \hat{r}}\right)_{\hat{r} = 0} = \left(\frac{\partial \hat{\phi}}{\partial \hat{r}}\right)_{\hat{r} = \hat{R}_0} = 0,\quad \left(\frac{\partial\hat{\eta}}{\partial\hat{r}}\right)_{\hat{r}=0} = \left(\frac{\partial\hat{\eta}}{\partial\hat{r}}\right)_{\hat{r}=\hat{R}_0}=0, \end{gather}

(2.6a,b)\begin{gather}\hat{\eta}(\hat{r},0) = a_0{\textrm{J}}_0\left(\frac{l_q}{\hat{R}_0}\hat{r}\right),\quad \frac{\partial\hat{\eta}}{\partial \hat{t}}(\hat{r},0) = 0, \quad \hat{\phi}(\hat{r},\hat{z},0) = 0 \end{gather}

In (2.6a,b), $l_q$ belongs to the set of (non-trivial) zeroes of ${\textrm{J}}_1(\cdot )$ with $q=1,2,3\ldots$, i.e. ${\textrm{J}}_1(l_q)=0$. This choice ensures that the initial condition satisfies (2.5a,b). For ease of algebra, we operate on (2.3) by ${\textrm {D}}/{\textrm {D}\hat {t}}$ replacing $\hat {\eta }$ in the gravity term with derivatives of $\hat {\phi }$ using (2.2) (Penney et al. Reference Penney, Price, Martin, Moyce, Penney, Price and Thornhill1952). In further calculations, we replace the kinematic boundary condition (2.2) with this equation. The resulting equations are then written in cylindrical axisymmetric coordinates ((1.1)–(1.6) in accompanying supplementary material available at https://doi.org/10.1017/jfm.2020.851) and non-dimensionalised using

(2.7a–e)\begin{equation} r = \frac{l_q\hat{r}}{\hat{R}_0}, \quad z = \frac{l_q\hat{z}}{\hat{R}_0}, \quad t = {\left( \frac{gl_q}{\hat{R}_0} \right)}^{1/2}\hat{t}, \quad \eta = \frac{l_q\hat{\eta}}{\hat{R}_0} \quad \text{and}\quad \phi = {\left( \frac{l_q^3}{g\hat{R}_0^3} \right)}^{1/2} \hat{\phi}. \end{equation}

This leads to the following set:

(2.8)\begin{align} &\frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\frac{\partial\phi}{\partial r} + \frac{\partial^2\phi}{\partial z^2} = 0, \end{align}

(2.9)\begin{align} & \left(\frac{\partial^{2} \phi}{\partial t^{2}} + \frac{\partial \phi}{\partial z} + \left\{ \frac{\partial}{\partial t} + \frac{1}{2}\left( \frac{\partial \phi}{\partial r} \right) \frac{\partial}{\partial r} + \frac{1}{2}\left( \frac{\partial \phi}{\partial z} \right) \frac{\partial}{\partial z} \right\} {|\boldsymbol{\nabla}\phi|}^{2} \right )_{z=\eta} \nonumber\\ &\quad -\sigma\left\lbrace\left[\frac{\partial}{\partial t} + \left(\frac{\partial\phi}{\partial r}\right)_{z=\eta}\frac{\partial}{\partial r}\right]\left(\frac{\dfrac{\partial^2\eta}{\partial r^2}}{\left\{1 + \left(\dfrac{\partial \eta}{\partial r}\right)^2\right\}^{3/2}} + \dfrac{1}{r}\dfrac{\dfrac{\partial\eta}{\partial r}}{\left\{1 + \left(\dfrac{\partial \eta}{\partial r}\right)^2\right\}^{1/2}}\right)\right\rbrace = 0, \end{align}

(2.10)\begin{gather} \left(\frac{\partial \phi}{\partial t} + \frac{1}{2}{ {| \boldsymbol{\nabla} \phi |}^2} + z\right)_{z=\eta} - \sigma\left(\frac{\dfrac{\partial^2\eta}{\partial r^2}}{\left\{1 + \left(\dfrac{\partial \eta}{\partial r}\right)^2 \right\}^{3/2}} + \frac{1}{r}\frac{\dfrac{\partial\eta}{\partial r}}{\left\{1 + \left(\dfrac{\partial \eta}{\partial r}\right)^2\right\}^{1/2}}\right) = 0, \end{gather}

(2.11)\begin{gather}\int_{0}^{l_q} r \eta(r,t) \,\textrm{d}r = 0, \end{gather}

(2.12a,b)\begin{gather} \left(\frac{\partial \phi}{\partial r}\right)_{r = 0} = \left(\frac{\partial \phi}{\partial r}\right)_{r = l_q} = 0,\quad \left(\frac{\partial\eta}{\partial r}\right)_{r=0} = \left(\frac{\partial\eta}{\partial r}\right)_{r=l_q}=0 , \end{gather}

(2.13a–c)\begin{gather} \text{and}\quad \eta(r,0) = \epsilon{\textrm{J}}_0(r),\quad \frac{\partial\eta}{\partial t}(r,0) = 0, \quad\phi(r,z,0) = 0, \end{gather}

Equations (2.8)–(2.13a–c) contain two non-dimensional groups viz. the nonlinearity parameter $\epsilon \equiv a_0({l_q}/{\hat {R}_0})$ and $\sigma \equiv {T({l_q}/{\hat {R}_0})^2}/{\rho g}$ which measures the relative strength of surface tension to gravity and is the inverse of Bond number. Note that $\epsilon$ appears in the initial condition while $\sigma$ appears in the boundary condition. It is important to remember in further analysis that $q$ is a fixed integer serving as an index for the primary mode. In this study $q$ has been chosen to be $5$, to remain consistent with the deep water approximation treated here. We treat $\epsilon$ as a small parameter with $\sigma = O(1)$ expanding $\phi$, $\eta$ and $t$ in (2.8)–(2.12a,b) as

(2.14)\begin{gather} \phi(r,z,t) = 0 + \epsilon\phi_1(r,z,t) + \epsilon^2\phi_2(r,z,t) + \epsilon^3\phi_3(r,z,t) + O(\epsilon^4), \end{gather}

(2.15)\begin{gather}\eta(r,t) = 0 + \epsilon\eta_1(r,t) + \epsilon^2\eta_2(r,t) + \epsilon^3\eta_3(r,t) + O(\epsilon^4), \end{gather}

(2.16)\begin{gather}\tau = t\left[1 + \epsilon^2 \varOmega_2 + O(\epsilon^3)\right]. \end{gather}

The expansion of $t$ in (2.16) is a standard Lindstedt–Poincaré technique and is necessary in order to take into account the nonlinear dependence of the oscillation frequency on the amplitude $\epsilon$ of the Bessel modes. Substituting (2.14)–(2.16) into (2.8)–(2.13a–c) and using Taylor series expansions about $z=0$, we may generically write at any $O(\epsilon ^i)$ ($i=1,2,3,\ldots$)

(2.17)\begin{gather} \frac{\partial^2\phi_i}{\partial r^2} + \frac{1}{r}\frac{\partial\phi_i}{\partial r} + \frac{\partial^2\phi_i}{\partial z^2} = 0, \end{gather}

(2.18)\begin{gather}\left(\frac{\partial^2\phi_i}{\partial \tau^2} + \frac{\partial \phi_i}{\partial z}\right)_{z=0} - \sigma\left[ \left(\frac{\partial^3 \eta_i}{\partial \tau \partial r^2}\right) + \frac{1}{r} \left(\frac{\partial ^2 \eta_i}{\partial \tau \partial r} \right) \right] = G_i(r,\tau), \end{gather}

(2.19)\begin{gather}\left(\frac{\partial \phi_i}{\partial \tau}\right)_{z=0} + \eta_i - \sigma\left[\left(\frac{\partial^2 \eta_i}{\partial r^2}\right) + \frac{1}{r} \left(\frac{\partial \eta_i}{\partial r}\right)\right] = F_i(r,\tau), \end{gather}

(2.20)\begin{gather}\int_{0}^{l_q} r \eta_i(r,t)\,\textrm{d}r = 0, \end{gather}

(2.21a,b)\begin{gather} \left(\frac{\partial \phi_i}{\partial r}\right)_{r = 0} = \left(\frac{\partial \phi_i}{\partial r}\right)_{r = l_q} = 0,\quad \left(\frac{\partial\eta_i}{\partial r}\right)_{r=0} = \left(\frac{\partial\eta_i}{\partial r}\right)_{r=l_q}=0, \end{gather}

(2.22a,b)\begin{gather} \eta_i(r,0) = {\textrm{J}}_0(r)\delta_{1i}, \quad \frac{\partial\eta_i}{\partial \tau}(r,0) = \phi_i(r,z,0)=0, \end{gather}

where $\delta _{ij}$ is the Kronecker delta. As is usual in perturbative methods, the only difference at different orders in (2.18) and (2.19) are the expressions $F_i$ and $G_i$ on the right-hand side. Here $F_1(r,\tau ) = G_1(r,\tau ) = 0$ while $F_2(r,\tau ),G_2(r,\tau ), F_3(r,\tau ), G_3(r,\tau )$ have lengthy expressions provided in appendix A. Our task is to now solve (2.17)–(2.22a,b) up to $i=2$ (second order). Note that for determining the solution up to $i=2$ completely, it is necessary to also obtain the equations at $i=3$.

Observe that $\phi (r,z,\tau )=p(\tau ){\textrm{J}}_0(\alpha r)\exp [\alpha z]$ satisfies the axisymmetric Laplace equation (2.17) for any (real) $\alpha$ and $p(\tau )$. In addition, if we choose $\eta (r,\tau ) = a(\tau ){\textrm{J}}_0(l_j({r}/{l_q}))$ and $\phi (r,z,\tau ) = p(\tau ){\textrm{J}}_0(l_j({r}/{l_q}))\exp (l_j({z}/{l_q}))$, the identity $\int _{0}^{l_q}{\textrm{J}}_0(l_j ({r}/{l_q})) r \,\textrm {d}r = 0$ and the equation ${\textrm{J}}_1(l_j)=0$ $\forall j \in \mathbb {Z}^{+}$ ensure that (2.20) and (2.21a,b) are also satisfied. Thus the general solution for $\phi _i$ and $\eta _i$ satisfying (2.17) and conditions (2.20) and (2.21a,b) may be posed as an eigenfunction expansion (Dini series, Mack Reference Mack1962)

(2.23)\begin{gather} \phi_i(r,z,\tau) = \sum_{j = 0}^{\infty} p_{i}^{(\,j)}(\tau) {\textrm{J}}_0\left(\alpha_j r\right) \exp\left(\alpha_j z\right), \end{gather}

(2.24)\begin{gather}\eta_i(r,\tau) = \sum_{j = 0}^{\infty} a_{i}^{(\,j)}(\tau) {\textrm{J}}_0\left(\alpha_j r\right)\quad \text{with} \ \alpha_j \equiv \frac{l_j}{l_q}. \end{gather}

Substituting (2.23) and (2.24) into (2.18) and (2.19) and using orthogonality relations among Bessel functions (see supplementary material), we obtain the following equations governing the coefficients $p_i^{m}(\tau )$ and $a_i^{m}(\tau )$ at any order $O(\epsilon ^i)$:

(2.25)\begin{gather} \frac{\textrm{d}^2p_{i}^{(m)}}{\textrm{d}\tau^2} + \omega_{m}^2p_{i}^{(\,j)}(\tau) = \left(1 + \sigma\alpha_m^2\right)\mathbb{G}_i^{(m)}(\tau) - \sigma\alpha_m^2\frac{\textrm{d} \mathbb{F}_i^{(m)}}{\textrm{d}\tau}, \end{gather}

(2.26)\begin{gather}a_i^{(m)}(\tau) = \frac{1}{1 + \sigma\alpha_m^2}\left(\mathbb{F}_i^{(m)}(\tau) - \frac{\textrm{d}p_i^{(m)}}{\textrm{d}\tau}\right),\quad m=0,1,2\ldots\,. \end{gather}

Here $\omega _m^2 \equiv \alpha _m(1 + \sigma \alpha _m^2)$ is the deep-water linear dispersion relation for capillary-gravity waves written in scaled variables. $\mathbb {F}_i^{(m)}(\tau )$ and $\mathbb {G}_i^{(m)}(\tau )$ in (2.25) and (2.26) are defined as

(2.27)\begin{equation} \left.\begin{gathered} \mathbb{G}_i^{(m)}(\tau) \equiv\frac{2}{l_q^2{\textrm{J}}_0^2\left(l_m\right)}\int_{0}^{l_q} \textrm{d}r r {\textrm{J}}_0\left(\alpha_m r\right) G_i(r,\tau), \\ \mathbb{F}_i^{(m)}(\tau) \equiv \frac{2}{l_q^2{\textrm{J}}_0^2\left(l_m\right)}\int_{0}^{l_q} \textrm{d}r r {\textrm{J}}_0\left(\alpha_m r\right) F_i(r,\tau). \end{gathered}\right\}\end{equation}

Equation (2.25) may now be solved for $p_i^{(m)}(\tau )$ (with suitable initial conditions) and the resultant solution determines $a_i^{(m)}(\tau )$ using (2.26). With $m=0,1,2,3\ldots$ and $i=1,2,3\ldots \,$, the initial conditions from (2.22a,b) are

(2.28a,b)\begin{gather} p_i^{(m)}(0) = 0,\quad \frac{\textrm{d}p_i}{\textrm{d}\tau}^{(m)}(0) = \mathbb{F}_i^{(m)}(0) - \left( 1 + \sigma\alpha_m^2 \right) \delta_{mq}\delta_{i1}, \end{gather}

(2.29a,b)\begin{gather}a_i^{(m)}(0) = \delta_{mq}\delta_{i1},\quad\frac{\textrm{d}a_i^{(m)}}{\textrm{d}\tau}(0) = 0. \end{gather}

We determine $p_i^{m}(\tau )$ and $a_i^{m}(\tau )$ using (2.25) and (2.26), respectively, at various orders using initial conditions (2.28a,b) and (2.29a,b).

2.2. Linear solution: $O(\epsilon )$

At $O(\epsilon )$, the solution is elementary. It is shown in the accompanying supplementary material that at linear order, the only non-zero term in the eigenfunction expansion in (2.23) and (2.24) is for $j=q$, the primary mode. With $\omega _q^2 = 1 + \sigma$, this leads to,

(2.30a,b)\begin{equation} \phi_1(r,z,\tau) ={-}\omega_q\sin(\omega_q\tau){\textrm{J}}_0(r)\exp\left(z\right) ,\quad \eta_1(r,\tau) = \cos(\omega_q\tau){\textrm{J}}_0(r). \end{equation}

2.3. Nonlinear corrections: $O(\epsilon ^2)$

At $O(\epsilon ^2)$, the calculation is algebraically tedious and quite lengthy, but otherwise standard. We outline the key results below, referring the interested reader to the accompanying supplementary material for the detailed derivation. Using (2.23) and (2.24), we obtain at $O(\epsilon ^2)$,

(2.31)\begin{align} \eta_{2}(r,\tau) &= \sum_{k = 1}^{\infty}\left[\zeta_1^{(k)}\cos(\omega_{k}\tau) + \zeta_2^{(k)}\cos(2\omega_{q}\tau) + \zeta_3^{(k)} \right]{\textrm{J}}_0(\alpha_k r) , \end{align}

(2.32)\begin{align} \phi_2(r,z,\tau) &= \frac{{\textrm{J}}_0^2(l_q)}{2}\sin\left(2\omega_q\tau\right) \nonumber\\ &\quad + \sqrt{1+\sigma}\sum_{k=1}^{\infty}\left[\xi^{(k)}_1\sin\left(\omega_k\tau\right) + \xi^{(k)}_2\sin\left(2\omega_q\tau\right)\right]{\textrm{J}}_0 \left(\alpha_k r\right)\exp\left(\alpha_k z\right), \end{align}

where expressions for $\zeta _1^{(k)},\zeta _2^{(k)},\zeta _3^{(k)}$ and $\xi _1^{(k)}, \xi _2^{(k)}$ are provided in appendix B. Note the presence of terms in (2.31) and (2.32), which are pure functions of space and time, respectively. This implies that there is a nonlinear correction to the (temporal) mean interface level, note the third term in (2.31) whose time average is non-zero. Similarly there is a nonlinear correction in $\phi$ whose spatial average is non-zero. An important thing to note in the nonlinear solution (2.31) and (2.32) is the presence of the entire spectrum (and not just the $q$th primary mode) at this order, which implies a nonlinear transfer of energy to the entire spectrum. This transfer of energy will be seen to be crucial to jetting.

Despite having expressions for $\phi (r,z,\tau )$ and $\eta (r,\tau )$ up to $O(\epsilon ^2)$, our calculation is incomplete, as a nonlinear correction to the dispersion relation $\omega_m^2 = \alpha_m(1+ \alpha_m^2\sigma)$$[1 + O(\epsilon^2)]$ is expected at this order, through $\varOmega _2$ in (2.16). For determining this, it is necessary to proceed to the third order ($O(\epsilon ^3)$) where elimination of resonant forcing of the primary mode (subscript $q$) will determine this correction. We emphasise that we do not determine the $O(\epsilon ^3)$ corrections to $\phi (r,z,\tau )$ and $\eta (r,\tau )$. A minimal $O(\epsilon ^3)$ calculation is carried out for determining $\varOmega _2$. This is done by obtaining the equations governing $p_3^{(m)}(\tau )$ and setting the coefficient of the resonant forcing term (for the primary mode) to zero. In this way, we obtain after some lengthy algebraic manipulations (see supplementary material)

(2.33)\begin{align} \varOmega_2 &= \frac{1}{2{\textrm{J}}_0^2(l_q)}\left[ \sum_{l = 1}^{\infty} \left[ \left\{\left[\alpha_l\left(2 + \sigma\right) - \alpha_l^2\right]\xi_2^{(l)} + 2\sigma\zeta_2^{(l)}\right\}{{I}}_{0 - l, 0- q, 0 - q}\right.\right.\nonumber\\ &\quad \left. - \left\{ 2\left(1 + \sigma\right)\alpha_l\xi_2^{(l)} + {\alpha_l^3\sigma}\zeta_2^{(l)} - 2\alpha_l^3\sigma\zeta_3^{(l)} \right\}{{I}}_{0 - q, 1 - q , 1 - l}\right] \nonumber\\ &\quad - \frac{\left(1 + \sigma\right)}{2}{{I}}_{0 - q, 0- q,0 - q, 0- q} + \frac{\left(1 - 4\sigma + \sigma^2\right)}{4\left(1 + \sigma\right)}{{I}}_{0 - q,0 - q, 1 - q, 1- q} \nonumber\\ &\quad \left. + \frac{3\left(1 + 3\sigma + \sigma^2 \right)}{4\left(1 + \sigma\right)}{{I}}_{0 - q, 1 - q, 1- q, 2 - q} \right]. \end{align}

Note that as shorthand notation in (2.33) we have used

(2.34)\begin{equation} {{I}}_{\nu_1- m_1, \nu_2- m_2,\ldots} \equiv \int_{0}^{1} \textrm{d}\tilde{r} \,\tilde{r}\, {\textrm{J}}_{\nu_1}\left(l_{m_1}\tilde{r}\right){\textrm{J}}_{\nu_2}\left(l_{m_2}\tilde{r}\right)\ldots \end{equation}

where $\tilde {r} \equiv r/l_q$ and $\nu _i,m_i = 0,1,2,\ldots$, e.g.

(2.35)\begin{equation} {{I}}_{0- m,1- q,1- q} \equiv \int_{0}^{1} \,\textrm{d}\tilde{r}\,\tilde{r}\,{\textrm{J}}_{0}\left(l_m\tilde{r}\right){\textrm{J}}_{1}^2\left(l_q\tilde{r}\right). \end{equation}

With $\sigma =1.086$ and parameters corresponding to air–water (see table 3), the amplitude correction to the dispersion relation for the primary mode $q=5$ evaluates to $\varOmega _2 \approx -0.0332$, using expression (2.33).

2.4. Summary of weakly nonlinear theory

We may now combine results up to $O(\epsilon ^2)$ to predict

(2.36)\begin{align} \eta(r,\tau) &= \epsilon{\textrm{J}}_0(r)\cos\left(\omega_q \tau\right) + \epsilon^2\sum_{k=1}^{\infty}\left[\zeta_1^{(k)}\cos\left(\omega_k\tau\right) + \zeta_2^{(k)}\cos\left(2\omega_q\tau\right) + \zeta_3^{(k)}\right]{\textrm{J}}_0\left(\alpha_k r\right) \nonumber\\ &\quad + O(\epsilon^3), \end{align}

(2.37)\begin{align} \phi(r,z,\tau) &={-}\epsilon \omega_q {\textrm{J}}_0(r)\exp\left(z\right)\sin\left(\omega_q\tau\right) + \epsilon^2 \frac{{\textrm{J}}_0^2(l_q)}{2}\sin\left(2\omega_q\tau\right) \nonumber\\ &\quad + \epsilon^2\sqrt{1+\sigma}\sum_{k=1}^{\infty}\left[ \xi^{(k)}_1\sin\left(\omega_k\tau\right) + \xi^{(k)}_2\sin\left(2\omega_q\tau\right) \right] {\textrm{J}}_0\left(\alpha_k r\right)\exp\left(\alpha_k z\right) + O(\epsilon^3) \end{align}

(2.38)\begin{align} \text{and}\quad \omega_k^2 &= \alpha_k\left(1+\sigma\alpha_k^2\right)\left[1+\epsilon^2\varOmega_2 + O(\epsilon^3)\right], \end{align}

where $\varOmega _2$ is given by (2.33). It may be checked that (2.36) and (2.37) satisfy the initial conditions (2.13a–c). In particular $\zeta _1^{(k)} + \zeta _2^{(k)} + \zeta _3^{(k)} =0$ for all $k=1,2\ldots$ which ensures that $\eta (r,0) = \epsilon {\textrm{J}}_0(r)$. Other initial conditions in (2.13a–c) involving the time derivative of $\eta (r,t)$ and on $\phi (r,z,0)$ are trivially verified. The singularity of the expressions for $\xi _1^{(k)}, \xi _2^{(k)}$, etc. at a possible $\omega _k = 2\omega _q$ is discussed further in § 6. As (2.36) and (2.37) are an infinite series, we truncate them for numerical evaluation in later sections. For any $r$, $z$ and $\tau$, it may be expected that the infinite sums in (2.36) and (2.37) have progressively smaller contributions at higher values of $k$. This is seen in figures 4(a) and 4(b) where we have evaluated the coefficients $\zeta _1^{(k)},\zeta _2^{(k)},\zeta _3^{(k)}$ and $\xi _1^{(k)},\xi _2^{(k)}$ for the primary mode $q=5$ as a function of $k$ in (2.36) and (2.37). It is seen that the peak values are at around $k=10$ which is expected, as this is the second multiple of $q=5$ and our calculation is accurate up to $O(\epsilon ^2)$. Thus for $q=5$, one needs to retain at least 10 terms while numerically evaluating the infinite series (2.36) and (2.37). For high accuracy, we have retained the first $85$ terms in our expansions. All the integrals that appear in $\zeta _1^{(k)}, \zeta _2^{(k)}\ldots$ have been evaluated numerically using MATLAB (MAT 2018).

Figure 4.

Here, $\zeta _1^{(k)},\zeta _2^{(k)},\zeta _3^{(k)}$ and $\xi _1^{(k)},\xi _2^{(k)}$ for $q=5$ as a function of $k$. After $k=12$, the values are quite small permitting truncation of the infinite series in (2.36) and (2.37).

5.1. Nonlinear redistribution of potential energy

We have seen that while the weakly nonlinear prediction is able to describe dimple inception and jet overshoot, it does not resolve the formation of the pronounced dimple during the second half of the oscillation or the thinning of the neck of the jet which causes droplet pinchoff. Here, we use modal analysis to understand the reasons why this is so. The interface $\hat {\eta }$ may be expressed at any time $\hat {t}_0$ using the Dini series in (5.2) with coefficients $\hat {H}(l_j;\hat {t}_0)$ (related to the Hankel transform). Here $\hat {H}(l_j;\hat {t}_0)$ may be evaluated from (5.3) (Bagrov et al. Reference Bagrov, Belov, Zadorozhyni and Trifonov2012) using $\hat {\eta }(\hat {r},\hat {t}_0)$ obtained either from theory or simulations,

(5.2)\begin{gather} \hat{\eta}(\hat{r},\hat{t}_0) = \sum_{j=1}^{\infty}\hat{H}(l_j;\hat{t}_0){\textrm{J}}_0\left(l_{j} \frac{\hat{r}}{\hat{R}_0}\right), \end{gather}

(5.3)\begin{gather}\hat{H}(l_j;\hat{t}_0)\equiv \frac{1}{\tfrac{1}{2}\hat{R}_0^2 {\textrm{J}}_0^2(l_j)}\int_{0}^{\hat{R}_0}\hat{r}\hat{\eta}(\hat{r},\hat{t}_0){\textrm{J}}_0\left(l_j\frac{\hat{r}}{\hat{R}_0}\right)\textrm{d}\hat{r}, \quad j=1,2,3\ldots\,. \end{gather}

In (5.2) and (5.3), $|\hat {H}(l_j;\hat {t}_0)|$ may be thought of as a measure of the potential energy in mode $j$ at time $\hat {t}_0$. In this study, we inject the entire initial energy as potential energy into the primary mode ($q=5$). Nonlinearity subsequently causes a part of this energy to be redistributed to modes which did not have any energy to start with. It will be seen here that this redistribution plays a crucial role in dimple formation events. The integral in (5.3) is evaluated numerically for each $j$. Figures 15(a) and 15(b) show $|\hat {H}(l_j;\hat {t}_0)|$ at the inception of dimple formation and later when a pronounced dimple forms, respectively. It is seen that at dimple inception, a large number of modes up to $j=10$ (second harmonic) share a sizeable fraction of the initial energy. At this time instant, the weakly nonlinear theory estimates the energy contained in the higher modes accurately leading to the good agreement for the dimple shape seen in the inset. At later time $\hat {t}_{0}=0.051$ s (figure 15b), when the dimple becomes more pronounced (narrower and sharper), many higher modes including noticeable peaks at $j=15$ and $20$ are seen. The potential energy for modes up to $j=10$ is well predicted by the theory, but for larger $j$ the theory suffers a significant error; this in turn leads to an unresolved dimple as seen in the inset of figure 15(b).

Figure 15.

(a) Here $\hat {H}(l_j;\hat {t}_0)$ versus $j$ ($\epsilon =0.5$) at dimple inception, $\hat {t}_0=0.02$ s, see inset for interface. At this instant, the dimple is captured by theory as significant energy transfer has occurred up to $j=10$ only. (b) Time $\hat {t}_0=0.051$ s. Dimple is sharper and narrower with significant surface energy in modes $>j=10$ also, which is not captured by the $O(\epsilon ^2)$ theory (c) Time $\hat {t}_0=0.051$ s. Interface shape obtained from (5.2) with $\hat {H}(l_j;0.051)$ obtained from simulations up to $j=10$ (pink symbols) and $j=20$ (blue symbols) while the remaining are from theory (2.36).

In order to analyse the unresolved dimple in figure 15(b) further, we perform a simple numerical test. In figure 15(c), we compare predictions from (5.2) where the first $N$ values of $\hat {H}(l_j;\hat {t}_0)$ (i.e. $j=1$ to $N$) have been obtained from simulations while the remainder ($\,j > N$) are from theory (2.36). It is seen clearly that in order to obtain the pronounced dimple seen in the inset of figure 15(b), one has to estimate $\hat {H}(l_j;\hat {t}_0)$ accurately up to $j=20$. Correcting $\hat {H}(l_j;\hat {t}_0)$ up to $j=10$ only in figure 15(c) is seen not to capture the shape of the dimple. This is consistent with the figure 15(b) where noticeable peaks are observed up to $j=20$.

Similar conclusions may be drawn from figure 16 where $\epsilon =1.2$. It is seen in figure 16(a) that at $\hat {t}_0=0.02$ s, dimple initiation is well described by the weakly nonlinear model as significant energy transfer has occurred up to $j=10$ only. However, pronounced dimple formation is not well resolved later at $\hat {t}_0=0.05$ s in figure 16(b). The spread of the energy becomes more delocalised now with small but noticeable potential energy in modes up to $j=50$ in figure 16(b), see lower inset.

Figure 16.

(a) Here $\hat {H}(l_j;\hat {t}_0)$ versus $j$ at first dimple formation ($\epsilon =1.2$), $\hat {t}_0=0.02$ s. (a) Significant energy transfer has occurred only up to modes $\leq j=10$. (b) A pronounced dimple with significant potential energy in modes $>j=10$, $\hat {t}_0=0.05$ s.

We have seen that the dimple becomes narrower with energy transfer to an increasingly broader set of modes, as $\epsilon$ is increased. The widest pronounced dimple seen in the present analysis is ${\approx }4.2$ mm wide (for $\epsilon =0.5$) and thus well below the air–water capillary-gravity wavelength of ${\approx }17$ mm. Consequently, these dimples can be treated as pure capillary waves, despite the primary perturbation being a much longer capillary-gravity wave. Our weakly nonlinear model is unable to capture the time evolution of the dimple due to its inability to resolve energy transfer beyond $j=10$. Despite this limitation, at the end of one oscillation when the jet is close to its maximum height, the local shape of the jet around $\hat {r}=0$ is well approximated by the nonlinear theory (see figure 11a–e). In order to understand this apparent paradoxical behaviour, $\hat {H}(l_j;\hat {t}_0)$ versus $j$ is plotted close to the instant of maximum jet height, in figures 17(a) and 17(b). It is seen that the maximum height of the jet is well described by the theory but the ‘shoulders’ of the jet are not. It may thus be inferred that in order to obtain the maximum height of the jet accurately, lower mode contributions ($\,j=1\text {--}10$) are important which our $O(\epsilon ^2)$ analysis achieves quite accurately. In contrast, in order to get the entire jet shape accurately, the ‘shoulder’ also needs to be resolved. This requires accurate estimation of higher-mode coefficients which require corrections presumably at $O(\epsilon ^3)$ and $O(\epsilon ^4)$. To sum up what we learn from this analysis so far, a short tabular summary of the relative strengths and weaknesses of the nonlinear model is provided in table 4.

Figure 17.

Here $\hat {H}(l_j;\hat {t}_0)$ versus $j$ when the jet is at its maximum height: (a) $\epsilon =1.2$; (b) $\epsilon =1.3$. As seen from the jet shape in the insets, the nonlinear model estimates the maximum height of the jet well but not the ‘shoulder’ of the jet, instead developing a broader base. It is concluded that in order to predict maximum height, one needs to estimate potential energy in modes up to $j=10$ only which the theory does well. However, for resolving the jet ‘shoulder’, one has to estimate the potential energy in modes $>j=10$ accurately which the $O(\epsilon ^2)$ theory cannot capture. Here $\hat {t}_0=0.067$ s.

Table 4.

Strengths and weaknesses of the $O(\epsilon ^2)$ solution to the IVP.