1. Introduction
Free-surface problems in fluid dynamics are usually challenging. The case of incompressible, inviscid flow, in situations where surface tension may also be neglected, is of particular importance for analysing ocean waves. If a body of inviscid fluid is disturbed, at least some fluid elements will be set on trajectories which will take them outside the initial boundary. Once a region of fluid starts to escape from the bulk, it may form a slender jet. As the jet extends, pressure gradients reduce, and the fluid elements which enter the jet move ballistically.
In many problems, such as ocean waves and splashes of objects dropped into the sea, gravity will play an important role. However, it is of interest to understand the solution of this problem in contexts where there is no gravity, for two reasons. Firstly, a clear understanding of the case without gravity is a foundation for understanding cases where gravity is important. Secondly, in astrophysical processes, bodies of fluid could be subjected to strong tidal forces during a transient encounter, before freezing into an exotically extended shape. A possible example is the Oumuamua object (Meech et al. Reference Meech, Weryk and Micheli2017), an asteroid size object of interstellar origin, which appeared to have a highly elongated structure.
In this paper we consider the simplest possible non-trivial case. There is an infinitely deep body of incompressible, inviscid fluid, which is initially stationary, with a flat interface, and the complementary half-space is assumed occupied by a gas of negligible density. The interface has zero surface tension. The system is subject to an impulsive disturbance, such that at
$t=0$
, the velocity is instantaneously
${\boldsymbol{u}}={\boldsymbol{\nabla }}\phi$
, where
$\phi$
is a velocity potential. Only the two-dimensional version of this problem will be addressed here. The natural questions are: From which points on the surface (if any) do the jets emerge? What is the velocity of the tip of the jet? How does the shape of the jet evolve? How much material is ejected in the jet? Or, as a possible alternative, does all of the fluid eventually find its way into the jet?
There is a substantial literature which considers jets arising from various types of shock or singularity which might be imposed on the velocity field. Longuet-Higgins (Reference Longuet-Higgins1983) argued that, when a standing wave reaches its maximum amplitude, the velocity and acceleration become infinite at the finite-angle cusp which was predicted by Penney & Price (Reference Penney and Price1952), causing a narrow jet to shoot upwards. In axisymmetric cases these jets have a conical envelope, with a cone angle
$\gamma$
which decreases as a function of time in a manner predicted in Longuet-Higgins (Reference Longuet-Higgins1976), and the flow field inside the cone is hyperbolic. Longuet-Higgins (Reference Longuet-Higgins1983) also presents evidence that the same jet structure describes jets projected upwards from bubbles as they break through a liquid surface. Recent works developing these themes include Gordillo & Blanco-Rodriguez (Reference Gordillo and Blanco–Rodríguez2023), Basak, Farsoiya & Dasgupta (Reference Basak, Farsoiya and Dasgupta2021), Kayal, Basak & Dasgupta (Reference Kayal, Basak and Dasgupta2022) and Kayal & Dasgupta (Reference Kayal and Dasgupta2023). Jets can also be created by focussing circularly symmetric surface waves to a point; see Zeff et al. (Reference Zeff, Kleber, Fineberg and Lathrop2000), Gordillo, Onuki & Tagawa (Reference Gordillo, Onuki and Tagawa2020) and McAllister et al. (Reference McAllister, Draycott, Davey, Yang, Adcock, Liao and van den Bremer2022). McAllister et al. (Reference McAllister, Draycott, Davey, Yang, Adcock, Liao and van den Bremer2022) discusses the evolution of the jet observed in their experiment in terms of the Longuet-Higgins theory. Other violent events which can trigger the formation of a jet include splashback from dropped objects impacting at high Reynolds number (see, e.g. O’Brien & Hodgins Reference O’Brien and Hodgins1995; Gekle et al. Reference Gekle, Gordillo, van der Meer and Lohse2009; Gekle & Gordillo Reference Gekle and Gordillo2010), or waves impacting a wall (Cooker & Peregrine Reference Cooker and Peregrine1992; Cooker Reference Cooker2002). These latter problems have multiple complications, involving gravity, vorticity generated by shear at the surface of the projectile, cavitation and turbulent dissipation.
In this work we adopt a distinct and complementary approach, examining the long-time evolution of jets which arise from a generic and non-singular initial condition. The existing literature does not appear to present a satisfying account of this foundational question. An exactly solvable case, that of a circular body subjected to a simple extensional flow, was discovered by Dirichlet (Reference Dirichlet1860), who showed that the circle deforms into an ellipse with increasing aspect ratio. Further developments of this approach are discussed in Longuet-Higgins (Reference Longuet-Higgins1972). The more general case which we consider was previously discussed by Longuet-Higgins (Reference Longuet-Higgins2001a
), who made the intriguing suggestion that the jet tip may have a cusped profile, described by a ‘canonical form’, which is somewhat reminiscent of the normal form for a cusp in ‘catastrophe theory’ (reviewed in Poston & Stewart Reference Poston and Stewart1996). We find that the form of the jet profile is different from that proposed by Longuet-Higgins, the width profile being a hyperbola, rather than a 3/2 power-law cusp. Some closely related ideas are presented in further works by Longuet-Higgins and collaborators, see Longuet-Higgins & Dommermuth (Reference Longuet-Higgins and Dommermuth2001a
,
Reference Longuet-Higginsb
) and Longuet-Higgins (Reference Longuet-Higgins2001b
). A paper by Longuet-Higgins & Dommermuth (Reference Longuet-Higgins and Dommermuth2001b
) considers the same initial conditions as our study, but with gravity included. Their simulations are initially similar to ours, but the ejected jet falls back into the bulk of the fluid before the long-time regime is established (see also Kayal & Dasgupta Reference Kayal and Dasgupta2023). Despite the fact that both cases may be described as ‘hyperbolic’ our jet has a very different structure from that considered in Longuet-Higgins (Reference Longuet-Higgins1983). Our jet, illustrated schematically in figure 1, has a cross-section which is approximated by a hyperbola. There is, therefore, no well-defined cone angle
$\gamma$
, which is used to characterise the jet in Longuet-Higgins (Reference Longuet-Higgins1983).
Figure 1. Schematic illustration of our conclusions about the form of the jet, as anticipated in (1.4)–(1.7). At large times, the tip of the jet extends with a constant speed. The profile of the base of the jet approaches a curve which is independent of time, apart from a downward displacement proportional to
$\ln t$
. The body of the jet has a hyperbolic form:
$y\Delta x=2C$
, where
$C$
is a constant.
We remark that there are yet more phenomena in fluid flow which are described as ‘jets’. A well-known review (Eggers & Villermaux Reference Eggers and Villermaux2008) defines a ‘jet’ as a ‘collimated stream of matter having a more-or-less columnar shape’. We shall use the term jet despite the flow being tapering rather than columnar. We are also aware of two other approaches to the problem of the motion of a two-dimensional perfect fluid with a free surface. The first of these, as exemplified by Kuznetsov et al. (Reference Kuznetsov, Spector and Zakharov1994), Zakharov et al. (Reference Zakharov, Kuznetsov and Dyachenko1996) and Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Zakharov1996), uses a Hamiltonian formalism with the surface displacement and surface potential as conjugate fields, combined with a conformal mapping to the half-plane. However, this approach is perturbative in small displacements and angles. Examples of the second approach are Karabut et al. (Reference Karabut and Zhuravleva2014) and Karabut et al. (Reference Karabut and Zhuravleva2024). These papers give exact solutions contrived to obey an auxiliary equation which keeps the free surface moving with a constant velocity. Neither of these approaches were useful for our problem, because we consider a non-perturbative flow of a rather general form.
For potential flow,
${\boldsymbol{u}}={\boldsymbol{\nabla }}\phi$
. If the fluid is incompressible, we have
${\nabla} ^2\phi =0$
. If the flow is inviscid, the potential and pressure
$p$
are related by the time-dependent Bernoulli equation
\begin{equation} \frac {\partial \phi }{\partial t}+\frac {1}{2}({\boldsymbol{\nabla }}\phi )^2+\frac {1}{\rho }p=\varPhi (t) . \end{equation}
In the following, we take
$\rho =1$
, where
$\rho$
is fluid density. We assume that the initial boundary is
$y=0$
, and that the pressure is
$p=0$
at the surface. The fluid is disturbed by an impulse at
$t=0$
, which will be assumed to be periodic in
$x$
. To avoid unnecessary complications, we discuss the following simple form for the initial velocity potential:
\begin{equation} \phi (x,y,0)=\cos (x)\exp (y)+\epsilon \cos (2x)\exp (2y), \end{equation}
and we assume
$|\epsilon | \lt 1/4$
. The initial velocity field is
\begin{equation} {\boldsymbol{u}}=\left (-\sin (x)\exp (y)-2\epsilon \sin (2x)\exp (2y),\cos (x)\exp (y) +2\epsilon \cos (2x)\exp (2y)\right )\! . \end{equation}
Note that there is a converging fixed point of the
$x$
-coordinate of fluid elements on the surface at
$x=0$
, although the surface itself moves in the
$y$
-direction, initially rising with speed
$u_y=1+2\epsilon$
at that point. Similarly there is a diverging fixed point of the surface flow at
$x=\pm \pi$
where the initial vertical speed
$u_y=-1+2\epsilon$
is downward. (The condition
$\epsilon \lt 1/4$
was imposed to guarantee that the surface flow is diverging at
$x=\pm \pi$
, see (2.9)).
We expect that the upward, converging fixed point of the surface flow becomes the tip of a jet moving vertically upwards. As material in the vicinity of this point escapes from the surface, the velocity gradients decrease, and (1.1) indicates that pressure gradients approach zero, so that the escaping fluid elements move ballistically (i.e. with constant velocity) at large time. The occurrence of ballistic motion in jets has been remarked upon by other authors: see for example Gekle and Gordillo (Reference Gekle and Gordillo2010). In particular, the speed of the tip is expected to approach a constant value
$v^\ast$
as
$t\to \infty$
. The diverging fixed point of the surface flow is expected to become the lowest point of a trough in the surface. It is a non-trivial task to determine the form of the jet, which in this paper will be described by estimating the surface profile
$y_{\textit{s}}(x,t)$
.
We were not able to find an exact solution to determining the surface profile. Figure 1 is a schematic illustration of the form of the jet at different times. After a brief initial transient, the position of the tip is approximated by
\begin{equation} y_{\textit{s}}(0,t)\sim v^\ast t+y_0, \end{equation}
where both the tip speed
$v^\ast$
and the offset
$y_0$
must be determined numerically (by fitting the empirically determined height of the jet at large
$t$
). The principal results of our analysis are a theory for the width of the jet, and for the depth of the trough. We argue that, as a function of height
$y$
above the original surface, the full width of the jet is
\begin{equation} \Delta x\sim \frac {2C}{y}, \end{equation}
for some constant
$C$
. This relation applies at sufficiently large times, with
$y$
satisfying
$y\ll v^\ast t$
and
$y\gg 1$
. At the base of the jet, we argue that the lowest point of the trough is, asymptotically as
$t\to \infty$
\begin{equation} y_{\textit{s}}(\pm \pi )\sim B-D\ln (t), \end{equation}
for some constants
$B$
, and
$D \gt 0$
. We also argue that, as
$t\to \infty$
, the base of the jet retains its form, so that
\begin{equation} y_{\textit{s}}(x,t)\sim F(x)+B-D\ln (t), \end{equation}
for some function
$F(x)$
, satisfying
$F(\pm \pi )=0$
and
$F(x)\to \infty$
as
$x\to 0$
. This approximation is accurate in the limit as
$t\to \infty$
, except in a small (and shrinking) interval surrounding the tip of the jet at
$x=0$
.
We obtain formulae for the profile of the jet at large times which depend upon the initial conditions via just two parameters, namely the total kinetic energy
$E$
, and the tip speed of the jet,
$v^\ast$
. The kinetic energy
$E$
is determined exactly from the initial impulse. The tip speed
$v^\ast$
depends upon the initial motion of the surface, and had to be determined by performing a numerical integration of the flow until the point at which the tip speed is no longer varying (figure 3 shows that this time interval is quite short).
Section 2 will discuss some general principles which form the basis for our analysis. Section 3 will present the argument which supports (1.5). In § 4 we derive an expression for the surface profile
$y_s(x,t)$
valid in the vicinity of the jet tip, and use conservation of energy to relate the tip velocity
$v^\ast$
to the coefficient
$C$
. Section 5 discusses the form of the base of the jet, justifying (1.6) and (1.7). Our investigations were informed by numerical simulations of the flow, based upon the method of Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976). We had to make significant modifications to their method, which are described in Appendix A. Section 6 uses an interpolation to combine our two approximate solutions (which describe, respectively, the jet and the base region), showing excellent agreement with numerical simulations. This final section also presents a brief conclusion.
2. Basic principles
2.1. Equation of motion for boundary
Analysing the motion of free surfaces in inertia-dominated flows is a challenging problem. A very significant advance was introduced in Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976), who gave equations of motion for fluid elements on the boundary, and for values of the velocity potential on the boundary. These equations are to be supplemented with a procedure for finding a harmonic function which represents the velocity potential in the interior of the region. In Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976), a boundary integral method is used to determine this harmonic function. In this paper, we use analytical approximations instead (primarily, using the ‘slender-body’ concept).
Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976) considered a two-dimensional flow, in a region with Cartesian coordinates
$(x,y)$
. They show that the equations of motion for the boundary are
\begin{equation} \frac {\textrm{D}x}{\textrm{D}t}=\frac {\partial \phi }{\partial x} \ ,\ \ \ \frac {\textrm{D}y}{\textrm{D}t}=\frac {\partial \phi }{\partial y} \ ,\ \ \ \frac {\textrm{D}\phi }{\textrm{D}t}=\frac {1}{2}({\boldsymbol{\nabla }}\phi )^2 . \end{equation}
The potential
$\phi (x,y,t)$
is a harmonic function, which is determined by the evolution of
$\phi$
on the boundary. (It is necessary to consider the potential in the interior region in order to determine the component of
${\boldsymbol{\nabla }}\phi$
normal to the boundary, which is required in (2.1)). Their paper also introduced a conformal transformation from the
$z=x+\textrm {i}y$
variable, to
$w=u+\textrm {i}v$
\begin{equation} w=u+\textrm {i}v=r\exp (\textrm {i}\theta )=\exp [-\textrm {i}(x+\textrm {i}y)]=\exp (-\textrm {i}z) . \end{equation}
The potential
$\phi$
remains a harmonic function when expressed in terms of the
$(u,v)$
variables, but the equations of motion for the boundary are transformed. We shall require them in both Cartesian and polar coordinates of the
$w$
-plane. The polar coordinate form was given in Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976)
\begin{equation} \frac {\textrm{D}r}{\textrm{D}t}=r^2\frac {\partial \phi }{\partial r} \ ,\ \ \ \frac {\textrm{D}\theta }{\textrm{D}t}=\frac {\partial \phi }{\partial \theta }, \end{equation}
and for evolution of potential on the boundary
\begin{equation} \frac {\textrm{D}\phi }{\textrm{D}t}=\frac {1}{2}\left [\left (r\frac {\partial \phi }{\partial r}\right )^2 +\left (\frac {\partial \phi }{\partial \theta }\right )^2\right ]\! . \end{equation}
In our calculations, we found it convenient to use Cartesian coordinates of the
$w$
-plane. In terms of the
$(u,v)$
variables, (2.3), (2.4) transform to
\begin{align} \frac {\textrm{D}u}{\textrm{D}t} &=(u^2+v^2)\frac {\partial \phi }{\partial u} ,\nonumber \\ \frac {\textrm{D}v}{\textrm{D}t} &=(u^2+v^2)\frac {\partial \phi }{\partial v} , \end{align}
\begin{align} \frac {\textrm{D}\phi }{\textrm{D}t} &=\frac {(u^2+v^2)}{2} \left [\left (\frac {\partial \phi }{\partial u}\right )^2+\left (\frac {\partial \phi }{\partial v}\right )^2\right ]\! . \end{align}
In the
$w$
-plane, the initial condition is a unit circle centred on the origin, and the potential at
$t=0$
is
$\phi (u,v)=\textit{Re}(w+\epsilon w^2)$
.
2.2. Global conservation laws
The initial kinetic energy
\begin{equation} E=\frac {1}{2}\int _{-\pi }^\pi \textrm {d}x\int _{-\infty }^0\textrm {d}y\ [{\boldsymbol{\nabla }}\phi ]^2=\frac {\pi }{2} [1+2\epsilon ^2], \end{equation}
is conserved, as is the volume, implying
\begin{equation} {\mathcal V}=\int _{-\pi }^\pi \textrm {d}x\ y_{\textit{s}}(x,t)=0, \end{equation}
where
$y_{\textit{s}}(x,t)$
is the surface height, initially
$y_{\textit{s}}(x,0)=0$
. We do not use conservation of momentum and angular momentum explicitly.
2.3. A local conservation principle
Consider the motion of material points on the surface. Because the surface has constant pressure, there is no acceleration of fluid elements along lines tangent to the surface. Pressure gradients in the direction normal to the boundary can accelerate fluid elements normal to the boundary. Changes of the relative speed of two elements on the surface can, therefore, only occur if the surface is curved, so that the normals point in different directions.
Consider the relevance of these observations to the motion of points on the surface which start in the vicinity of the unstable fixed point of the surface flow, at
$x=\pm \pi$
. A point which is displaced from this fixed point by a distance
$\delta x$
has horizontal speed (to leading order in
$\delta x$
)
\begin{equation} u_x=(1-4\epsilon )\delta x . \end{equation}
The point continues to move with the same horizontal speed until the slope of the boundary is significantly different from zero.
2.4. Velocity gradient
While not essential to our arguments, it is informative to consider the velocity gradient matrix
\begin{equation} \boldsymbol {A}=\left (\begin{array}{cc} \partial ^2_{xx}\phi & \partial ^2_{xy}\phi \cr \partial ^2_{xy}\phi & \partial ^2_{yy}\phi \end{array}\right )\! . \end{equation}
In two-dimensional potential flow the velocity gradient matrix is a traceless symmetric matrix. It follows that the eigenvectors are orthogonal, and the eigenvalues are both real and that their sum is zero.
The equation of motion for
$\boldsymbol {A}$
is (in the case of an inviscid fluid)
\begin{equation} \frac {\textrm{D}\boldsymbol {A}}{\textrm{D}t}+\boldsymbol {A}^2=-\frac {1}{\rho }\boldsymbol {H}, \end{equation}
where
$\boldsymbol {H}$
is the Hessian matrix of the pressure (with elements
$H_{\textit{ij}}=\partial ^2_{\textit{ij}}p$
) (see Cantwell Reference Cantwell1992).
At long times, we expect that most pairs of material points become increasingly widely separated, implying that pressure gradients decrease. We should, therefore, expect that the pressure Hessian approaches zero at large time. If
$\boldsymbol {H}$
is negligible, (2.11) indicates that
$A_{\textit{ij}}\sim 1/t$
as
$t\to \infty$
. This, in turn, would indicate that we might write
\begin{equation} \phi (x,y,t)\sim \frac {1}{t}\psi (x,y), \end{equation}
in the limit as
$t\to \infty$
, where
$\psi (x,y)$
is a harmonic function. Approximations of this form will figure in our analysis.
2.5. Physical principles of the solution
Immediately after the impulse is applied at
$t=0$
, the pressure in the fluid can be determined. (For completeness, the calculation is discussed in Appendix B). We find that the pressure at depth increases instantaneously, implying an upward acceleration of fluid elements close to the surface. This must be a transient effect otherwise it would accelerate them to the extent that energy would not be conserved. The pressure field immediately following the initial impulse can, therefore, be thought of as a short-lived pulse which imparts a secondary impulse to the fluid elements. We were not able to make a precise statement about this secondary impulse, so that the initial response of the system must be determined by a numerical calculation.
Once the jet starts to extend out of the body of the fluid, all fluid elements in the jet end up close to the surface. As a consequence, pressure gradients must rapidly become negligible, implying that all fluid elements in the jet move ballistically. This means that, in the long-time limit, an element which reaches a height
$y$
above the fluid surface at time
$t$
has a vertical speed given approximately by
\begin{equation} u_y\sim \frac {y-y_0}{t}, \end{equation}
for some choice of offset,
$y_0$
. Figure 2(a) shows the numerically computed vertical speed of different fluid elements on the surface as a function of time, for
$\epsilon =0$
. The elements close to the tip (top curve) rapidly approach a constant vertical speed. Those elements which are close to the trough have a speed which decays slowly, as
$1/t$
, in accord with (1.7). Figure 2(b) shows the form of the surface at three values of the time
$t$
, showing that the surface has undergone a significant deformation by the time the tip reaches its maximum speed. This suggests that it would be difficult to make an analytical estimate for
$v^\ast$
and therefore we use the numerically determined value of
$v^\ast =1.84$
for
$\epsilon =0$
. Figure 3 shows the time dependence of the speed of the tip (at
$x=0$
) for different values of
$\epsilon$
, and how its steady long term speed,
$v^\ast$
, depends upon
$\epsilon$
.
Figure 2. (a) Numerically computed vertical speed of 11 elements initially evenly spaced along the surface (
$y=0$
) from
$x=-\pi$
(initial speed
$-1$
) to
$x=0$
(initial speed
$1$
). These data show that the transient pressure pulse discussed in § 2.5 is short lived. (b) Fluid surface at
$t=0.5$
,
$t=1$
and
$t=1.5$
. Throughout,
$\epsilon =0$
.
Figure 3. (a) Tip velocity as a function of time for different values of
$\epsilon$
as discussed in § 2.5. (b) Constant tip velocity (post initial acceleration) as a function of
$\epsilon$
(
$\epsilon \geqslant 0$
).
Section 3 will describe how (2.5) and (2.6) can be used, together with (2.13), to deduce the hyperbolic form of the body of the jet, described by (1.5). Here, we use the equations of motion for material points on the boundary, expressed in the
$(u,v)$
coordinate system. In this coordinate system, the curve representing the surface extends very rapidly along the positive
$u$
axis as
$t\to \infty$
: if (as expected)
$y\sim v^\ast t$
, then
$r\sim \exp (v^\ast t)$
. This suggests that the closed curve representing the surface has a much greater extension along the
$u$
-axis than along the
$v$
-axis. This would have the benefit of allowing us to use a ‘slender-body’ approach to determining the potential.
Our results indicate that the base of the jet does not get narrower as time increases, rather it approaches a time-independent profile, which simply sinks lower as time increases, in accord with (1.7). Figure 4 shows evidence of this from our numerical model.
Figure 4. Fluid boundary at
$t=0,1,\ldots ,8$
from numeric calculations (with
$\epsilon =0$
), vertically shifted so that
$y=0$
at
$x=\pm \pi$
at all times.
3. Form of body of jet
In this section we discuss the long-time behaviour of the jet. It is assumed that the jet has already started to escape from the surface, so that the form of the boundary in the
$w$
-plane is extended along the
$u$
-axis, as discussed in § 2.5. The form of (1.1) suggests considering solutions of the form (2.12), because the time derivative and the gradient-squared term both lead to factors of
$t^{-2}$
. Furthermore, because we assume that the boundary is slender, being extended in the
$u$
-direction, we can approximate the solution of the Laplace equation by a function which is restricted to having a quadratic dependence upon
$v$
\begin{equation} \phi (u,v,t)\sim \frac {1}{t}\left [g(u)-\frac {1}{2}g^{\prime \prime}(u)v^2\right ]\!, \end{equation}
for some general smooth function
$g(u)$
.
Now we introduce the pivotal assumption, namely that at large time, the material points on the boundary move tangentially along this curve, which is based on observations from our numerical work. This curve is defined by a function
$f$
\begin{equation} v=f(u) . \end{equation}
Now using (2.5) to determine a differential equation for
$f(u)$
\begin{equation} f^{\prime}= \displaystyle\frac{\displaystyle\frac{\textrm{D}v}{\textrm{D}t}}{\displaystyle\frac{\textrm{D}u}{\textrm{D}t}} = \displaystyle\frac{\displaystyle\frac{\partial \phi }{\partial v}}{\displaystyle\frac{\partial \phi }{\partial u}}= - \displaystyle\frac {g^{\prime \prime}(u)}{g^{\prime}(u)}f, \end{equation}
and integrating gives
\begin{equation} f=\frac {C}{g^{\prime}} . \end{equation}
In order to determine the form of the jet, we need to identify the function
$g(u)$
. Our ‘slender-body’ assumption means
$v$
is negligible relative to
$u$
, so we consider
$v=0$
. We can infer that, close to
$x=0$
and when
$t$
is large,
$\phi \sim y^2/2t$
, because this gives a vertical velocity
$u_y=\partial \phi /\partial y=y/t$
in accord with (2.13). A more refined estimate of
$\phi$
can be obtained from substituting
$\phi (u,0,t)= ({1}/{t})g(u)$
into (2.6), leading to
\begin{equation} g(u)=\frac {1}{2}(ug^\prime (u))^2 . \end{equation}
Hence
\begin{equation} g(u)=\frac {(\ln (\beta u))^2}{2}, \end{equation}
where
$\beta$
is a constant of integration. Hence
\begin{equation} f(u)=C\frac {u}{\ln (\beta u)} . \end{equation}
Using the fact that the profile must be symmetric about
$x=0$
, and mapping back to the
$(x,y)$
-plane using
$x = -\arctan (f(u)/u)$
and
$y = \ln (r) = \ln (\sqrt {u^2 + f(u)^2})$
gives, for large
$u$
\begin{align} x &=\pm \arctan \left (\frac {C}{\ln (\beta u)}\right ) \approx \pm \frac {C}{\ln (\beta u)}, \end{align}
\begin{align} y &= \ln (u) + \ln \left ( \sqrt {1 + \frac {C^2}{\ln (\beta u)^2}} \right ) \approx \ln (u) . \end{align}
We conclude that, when
$u\gg 1$
(implying that
$y$
is large)
\begin{equation} x \sim \pm \frac {C}{y+\ln (\beta )}; \end{equation}
a hyperbola, close to the
$x$
-axis. Note that, because we assumed that the trajectories are tangent to the boundary curve at any instant in time,
$\beta$
may be time-dependent. Equation (3.10) is then consistent with (1.7). The fact that we were able to find a solution validates our assumption that the velocity is tangent to the surface. Comparison with (1.6) indicates that we should choose
$\beta$
so that
$|x|$
becomes large when
$y\to B-D\ln t$
, that is as
$t\to \infty$
the form of the boundary at large
$y$
is
\begin{equation} x \sim \pm \frac {C}{y-B+D\ln t}, \end{equation}
where the constants,
$B$
,
$C$
,
$D$
, are still to be determined.
Figure 5 illustrates the numerically computed surface in the
$(u,v)$
plane, showing evolution from a unit circle at
$t=0$
to a curve which is extended along the positive
$u$
-axis. The large aspect ratio of this curve justifies the use of a ‘slender-body’ approach. This figure was produced by calculating the boundary in
$(x,y)$
coordinates using the methods described in Appendix A, and using (2.2) to transform to the
$(u,v)$
plane. Solving (2.5), (2.6) numerically presents a challenge, because of the explosive growth of the boundary curve when expressed in these coordinates. We were able to do this to just beyond
$t=1$
using a different approach to the calculations than our principal method. Up to the point of failure the two methods show excellent agreement. Details are provided in Appendix A.
Figure 5. The fluid surface is represented by a closed curve in
$(u,v)$
space, which starts from a unit circle at
$t=0$
, and which rapidly extends along the positive
$u$
-axis panel (a) is at
$t=0$
, panel (b) at
$t=1.6$
and panel (c) at
$t=6$
. The large aspect ratio as
$t\to \infty$
validates the use of a slender-body approximation.
4. Relation between tip speed and width parameter
In § 3 it was argued that the width of the jet is, to leading order as
$t\to \infty$
,
$\Delta x\sim 2C/y$
for some coefficient
$C$
, which is independent of time. The tip advances with a speed
$v^*$
. The values of
$C$
and
$v^*$
are related by conservation of energy.
The kinetic energy is given by (2.7). Most of the kinetic energy is contained in the region close to the tip, where the vertical speed is nearly independent of
$x$
within the jet. If the width of the jet is
$\Delta x(y)$
, then the energy is approximated by
\begin{equation} E=\frac {1}{2}\int ^{v^\ast t+y_0}_{y_{{b}}} \textrm {d}y\ \Delta x(y) [u_y(y)]^2, \end{equation}
where
$u_y(y)$
is the vertical speed at height
$y$
, and
$y_{{b}}$
a value close to the base of the jet.
We know that
$u_y(y)=(y-y_0)/t$
(see (2.13)). Close to the tip of the jet,
$y_{\textit{s}}$
may be assumed to be a quadratic function of
$x$
, and as we move away from the tip region,
$y_{\textit{s}}$
may be approximated by a hyperbola. We can, therefore, replace the profile described by (3.11)) by a function of the form
\begin{equation} y(x,t)\sim Y+\frac {C}{\sqrt {x^2+x_0^2}}, \end{equation}
where
$Y=B-D\ln t$
and where
$x_0(t)$
will be chosen so that
$y=y_0+v^\ast t$
when
$x=0$
. To leading order as
$t\to \infty$
, we have
$x_0\sim C/(v^\ast t)$
.
Note that the energy integral (4.1) is dominated by contributions from the tip, rather than the base, because the velocity approaches zero at the base. So, setting
$y^{\prime}=y-y_0$
, (4.2) implies that, as
$t\to \infty$
,
$x\sim C\sqrt {v^{\ast 2}t^2-{y'}^2}/y^{\prime}v^\ast t$
. Hence, noting that
$\Delta x=2x(y)$
, as
$t\to \infty$
, (4.1) is approximated by
\begin{equation} E=C\int _{y_{{b}}-y_0}^{v^\ast t}\textrm {d}y^{\prime}\ \frac {y^{\prime}}{v^\ast t^3}\sqrt {(v^{\ast }t)^2-y'^2}\sim \frac {C v^{\ast 2}}{3} . \end{equation}
Knowing
$E$
, we can determine
$C$
in terms of
$v^\ast$
\begin{equation} C=\frac {3\pi (1+2\epsilon ^2)}{2v^{\ast 2}} . \end{equation}
Using a precise expression for
$x_0$
, for which (4.2) gives
$y=y_0+v^\ast t$
when
$x=0$
, we obtain our best approximation for the form of the jet
\begin{equation} y(x,t) \sim B-D\ln t+\frac {C}{\sqrt {x^2+\frac {C^2}{(y_0+v^\ast t-B+D\ln t)^2}}} . \end{equation}
Figure 6 compares the numerically computed tip profile with the prediction of (4.5) using (4.4) and (5.3), for the case
$\epsilon =0$
. The fitted parameters were
$v^\ast =1.84$
,
$y_0=-0.34$
,
$B=\ln (b)=-0.798$
(see § 5.2).
5. Modelling the base of the jet
5.1. Lowest level of trough
First let us consider the height of the lowest point,
$x=\pm \pi$
, at time
$t$
. This will be estimated by assuming that at base of the jet has a fixed shape, and simply sinks in the limit as
$t\to \infty$
. The sinking of the lowest point is then determined using conservation of volume.
Using (4.2) to describe the shape of the tip of the jet, the volume of the tip above
$y=y_{{b}}$
(with
$y_{{b}}\ll v^\ast t$
), is
\begin{align} {\mathcal V}_+ &= \int _{y_{{b}}}^{v^\ast t+y_0}\textrm {d}y\ \Delta x(y) \sim 2C \int _{y_{{b}}-y_0}^{v^\ast t}\textrm {d}y^{\prime}\ \frac {\sqrt {v^{\ast 2}t^2-y'^2}}{y'v^\ast t} \nonumber \\[2pt] &\quad \sim 2 C \left [-\ln (\varGamma )+\ln \left (1+\sqrt {1-\varGamma ^2}\right )-\sqrt {1-\varGamma ^2}\right ]\!, \\[9pt] \nonumber \end{align}
where
\begin{equation} \varGamma =\frac {y_{{b}}-y_0}{v^\ast t} . \end{equation}
Note that, in the limit as
$t\to \infty$
, the second and third terms of (5.1) tend to
$\ln (2) - 1$
, so
$V_+ \sim -2C\ln (\varGamma )$
, i.e. the dependence upon
$t$
is logarithmic, with coefficient
$2C$
. The increase in
${\mathcal V}_+$
is compensated by lowering the base region, for which the profile is described by (1.7), except in a narrowing interval in the vicinity of the tip, at
$x=0$
. Accounting for the fact that the width of the region is
$2\pi$
, the contribution to the area change from the base has a logarithmic term,
${\mathcal V}_-=-2\pi D\ln t$
, where
$D$
is the coefficient in (1.7). Hence conservation of volume indicates that
\begin{equation} D=\frac {C}{\pi }=\frac {3(1+2\epsilon ^2)}{2v^{*2}} . \end{equation}
Figure 7 shows the numerically determined level of the base of the trough, compared with (1.6), using the fitted parameters.
5.2. Shape of base of trough
Consider the form of the jet when expressed in terms of the
$(u,v)$
variables (defined by (2.2)). The bottom of the trough at
$x=\pm \pi$
corresponds to a point in the
$(u,v)$
plane, with coordinate which approaches the origin as
$t\to \infty$
: taking the exponential of (1.7) suggests that the position
$(u^\ast (t),0)$
corresponding to the bottom of the trough has
\begin{equation} u^\ast (t)=-b t^{-D}; \end{equation}
(with
$B=\ln b$
). The requirement that volume is preserved, at all times, determines how the height of the base region, as expressed in (1.6), varies as a function of time. This imposes constraints on both
$D$
and
$b$
.
$D$
is determined from (5.3) while value of
$b$
needs to be determined numerically, and we take this as 0.45, which it takes at the largest value of
$t$
(approximately
$t=8$
) for which our numerical solution was stable.
Let us consider the consequences of the surface corresponding to a curve in the
$(u,v)$
plane scaling in the same manner as (5.4). Define scaled coordinates
$(U,V)$
, a transformed time variable
$\tau$
, and a transformed potential
$\varPsi (U,V)$
by writing
\begin{equation} (U,V)=t^{D}(u,v) \ ,\ \ \ \tau =\ln (t) \ ,\ \ \ \phi (u,v,t)=\frac {1}{t}\varPsi (U,V,\tau ) . \end{equation}
In these variables, the Longuet-Higgins and Cokelet equations of motion, (2.5), (2.6) become
\begin{align} \frac {\textrm{D}U}{\textrm{D}\tau} &= (U^2+V^2)\frac {\partial \varPsi }{\partial U}+DU ,\nonumber \\ \frac {\textrm{D}V}{\textrm{D}\tau } &= (U^2+V^2)\frac {\partial \varPsi }{\partial V}+DV , \end{align}
and
\begin{equation} \frac {\textrm{D}\varPsi }{\textrm{D}\tau }=\frac {(U^2+V^2)}{2} \left [\left (\frac {\partial \varPsi }{\partial U}\right )^2+\left (\frac {\partial \varPsi }{\partial V}\right )^2\right ] +\varPsi . \end{equation}
We hypothesise that the boundary approaches a fixed curve when expressed in terms of the
$(U,V)$
variables. Note that, if
$\varPsi (U,V,\tau )$
is independent of
$\tau$
, then (5.6) describes motion with a fixed velocity field in the
$(U,V)$
plane. The boundary curve might be represented by a curve representing the unstable manifold attached to an unstable fixed point at
$U^\ast =-b$
, which corresponds to the position of the unstable fixed point of the surface flow at
$x=\pm \pi$
. Note that (5.6) and (5.7) are contrived to conform to the scaling of the solution as
$t\to \infty$
: results obtained from them will not be valid for small
$t$
.
The curve in
$(U,V)$
space represents a family of curves in
$(x,y)$
space, corresponding to different values of
$t$
. Consider how the value of
$y$
changes for fixed
$x$
, as
$t$
is varied. According to (2.2), fixing
$x$
corresponds to fixing the polar angle
$\theta$
in the
$(u,v)$
plane, and hence also in the
$(U,V)$
plane. The ray at polar angle
$\theta$
intersects the curve at a point
$(U(\theta ),V(\theta ))$
, with radius
$R(\theta )=\sqrt {U^2+V^2}$
. The corresponding value of
$y$
is
$y=\ln (t^{-D}R(\theta ))$
. It follows that the values of
$y$
are shifted by the same amount, independent of
$x$
, as
$t$
increases. That is, the hypothesis that there is a fixed curve in
$(U,V)$
space is equivalent to the hypothesis that the boundary surface sinks in the
$(x,y)$
plane, without changing its form.
The unstable fixed point at
$U^\ast$
corresponds to the unstable fixed point of the surface flow at
$x=\pm \pi$
. We can expand
$\varPsi (U,V)$
about this fixed point, and determine coefficients by comparing with the unstable fixed point of the surface flow. Noting that
$\varPsi$
should be a harmonic function, to leading order
\begin{equation} \varPsi (U,V)=\frac {\mu }{2}(U^2-V^2)+\nu U, \end{equation}
where
$\mu$
and
$\nu$
are constants to be determined. To determine the fixed point
$U^\ast$
= −
$b$
, set
$\textrm {D}U/\textrm {D}\tau =0$
to obtain
\begin{equation} \mu b^2-\nu b+D=0 . \end{equation}
The equation of motion for
$V$
in the vicinity of the fixed point satisfies
\begin{equation} \frac {\textrm{D}V}{\textrm{D}\tau }\sim \big [D-\mu b^2\big ]V, \end{equation}
so that
\begin{equation} V(\tau )\sim V_0 \exp [\lambda \tau ]=V_0 t^\lambda \ ,\ \ \ \lambda =D-\mu b^2 . \end{equation}
The principle discussed in § 2.3 above implies that material points on the surface move with an approximately fixed velocity while they are still in the vicinity of the unstable fixed point of the surface flow. If their initial displacement from this fixed point is
$\delta x_0$
, the displacement after time
$t$
is
\begin{equation} \delta x(t)=\delta x_0+ (1-4\epsilon )\delta x_0 t . \end{equation}
The value of
$\lambda$
in (5.11) can be determined by comparing (5.11) and (5.12). Because
$\delta x=\pi -\theta \sim -V/U$
close to the negative
$U$
axis, (5.11) and (5.12) are only consistent as
$t\to \infty$
if
\begin{equation} \lambda =D-\mu b^2=1 . \end{equation}
We have assumed that
$b$
and
$D$
are known. Equations (5.9) and (5.13) then enable us to determine
$\mu$
and
$\nu$
. Close to the fixed point,
$U^\ast =-b$
, the equation of motion for
$U$
is approximated by
\begin{equation} \frac {\textrm{D}U}{\textrm{D}\tau }\sim \frac {D}{b}V^2 . \end{equation}
Using (5.11) with
$\lambda =1$
and as
$V=0$
at
$U=-b$
, then the curve in the
$(U,V)$
plane is approximated by
\begin{equation} V_{{b}}(U)\sim \pm \sqrt {\frac {2b}{D}(U+b)} \ ,\ \ \ U+b\gt 0 . \end{equation}
Figure 8 shows the numerically determined surface contour in
$(U,V)$
space, at
$t=2,4,8$
, using the value of
$D$
predicted by (5.3). The curves are nearly coincident, consistent with the expectation that they should approach a limiting curve as
$t\to \infty$
. The figure also shows the analytic approximation to this limiting curve, (5.15). Close to the critical point,
$(U,V)=(U^\ast ,0)$
, the numerical curves approach this parabolic approximation as
$t$
increases. Away from the critical point, there is a an increasing discrepancy, which is already significant when
$U=0$
. It would be possible in principle to continue the expansion about the fixed point to higher orders in
$V$
, but we believe that this leading-order solution is sufficient to understand the physical principles.
Figure 8. Numerically determined surface contour in
$(U,V)$
space at
$t=2$
(blue),
$t=4$
(green) and
$t=8$
(orange), using the value of
$D$
predicted by (5.3), and
$b=0.45$
. The figure also shows the analytic approximation (purple) to the limiting curve in (5.15). These data are for
$\epsilon =0$
.
Equation (5.15) expresses the form of the boundary close to the trough in the
$(U,V)$
plane. We can use (2.2) and (5.5) to transform this to a curve in the original Cartesian coordinates. The result, illustrated in figure 9, is in good agreement with numerical simulations of the boundary. The Cartesian form of the equation is not particularly informative of itself, however, so is not reproduced here.
6. Discussion
6.1. Synthesis
We have obtained two approximations for the form of the surface, one valid in the jet ((4.5), (4.4)), and the other describing the base of the jet (5.15). Figures 6 and 9 show that these approximations work well, close to
$x=0$
and
$x=\pm \pi$
, respectively. Figure 10 compares the numerically computed profile with an interpolation between these two expressions
\begin{equation} y(x,t)=w(x)\,y_{\textit{j}}(x,t)+[1-w(x)]\,y_{{b}}(x,t), \end{equation}
where
$y_{\textrm{j}}(x,t)$
and
$y_{{b}}(x,t)$
are the profiles for the jet and base regions, (4.5) and (5.15) respectively, and where the interpolating weight function was
\begin{equation} w(x)=\left \{ \begin{array}{cc} -2x^2/\pi ^2+1 & 0 \leqslant |x| \lt \pi /2, \cr 2x^2/\pi ^2-4|x|/\pi +2 & \pi /2 \leqslant |x| \leqslant \pi. \cr \end{array} \right . \ \end{equation}
(The interpolation function is arbitrary: this example changes smoothly from
$w(0)=1$
to
$w(\pi )=0$
, with
$w(\pi /2)=1/2$
.) There is excellent agreement between the interpolated theoretical function and the numerically determined surface contour, for
$t=4$
and
$t=8$
.
Figure 10. Panels (a) and (b) compare the numerically determined surface contours in
$(x,y)$
space, at
$t=4,8$
respectively, with theory ((4.5), (4.4) describing the jet and (5.15), describing the base, interpolated according to (6.1) and (6.2)). Zoomed in images of the jet tips are shown at the upper right of each panel. These data are for
$\epsilon =0$
.
6.2. Conclusions
Understanding the ejection of jets from the surface of a simple fluid is a fundamental and surprisingly challenging problem. Earlier works, such as Longuet-Higgins (Reference Longuet-Higgins1983), Zeff et al. (Reference Zeff, Kleber, Fineberg and Lathrop2000), Gordillo et al. (Reference Gordillo, Onuki and Tagawa2020), McAllister et al. (Reference McAllister, Draycott, Davey, Yang, Adcock, Liao and van den Bremer2022) and Kayal & Dasgupta (Reference Kayal and Dasgupta2023) have emphasised situations where the jet is initiated at a singularity. However, the production of ballistic jets from inertial flows must be a ubiquitous phenomenon. In this paper we considered what is arguably the simplest extension of that model: an initially flat surface of a deep body of liquid in two dimensions is subject to a spatially periodic, impulsive disturbance. Even for this case, it is difficult to produce numerical results which are reliable at large time, and it appears to be impossible to find a precise asymptotic solution which describes all regions of the motion. However, we were able to show, in § 3, that, away from the base and tip, the jet profile is hyperbolic. This is at variance with an earlier suggestion Longuet-Higgins (Reference Longuet-Higgins2001a ) that the profile is a cusp.
We also found that the base region moves downwards without changing its form, and that the depth is logarithmic in time. There is no separatrix, dividing material elements which will eventually enter the jet from those which remain in the bulk. Eventually, every fluid element, at any initial depth, becomes part of the jet.
The speed
$v^\ast$
of the tip of the jet must be determined numerically, but once this has been determined, we can predict the evolution of the jet at large time. In particular, the parameters
$C$
determining the width of the jet, and
$D$
describing the depth of the trough, are given by explicit formulae, (4.4) and (5.3). The parameter
$b$
introduced in (5.4) is not determined by an explicit formula, but it can be determined by the requirement that volume is conserved.
The argument in Longuet-Higgins (Reference Longuet-Higgins2001a
) uses a somewhat different approach: seeking a coincidence of two contours (the surface is a zero contour of both the pressure
$p$
and its material derivative,
$\textrm {D}p/\textrm {D}t$
). This is implemented by developing the velocity potential as a Taylor series. His approach suffers from the difficulty that the method develops a velocity potential which has unbounded growth at infinity, and which cannot be matched to any physical boundary condition.
There remain many open questions about the ejection of jets from bodies of inviscid fluid. In our case the point of origin of the jet and its direction (
$x=0$
, upwards) were obvious, but in cases with a more general initial condition, the number, origin, direction, speed and relative mass of the jets are hard to determine with certainty. Extensions to three dimensions are also non-trivial.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical technique
A.1. Basic method
The velocity potential at time
$t$
will be expressed in terms of Fourier coefficients,
$\alpha _n(t)$
\begin{equation} \phi (x,y,t)=\sum _{n=0}^{N_c} \alpha _n(t)\cos (nx) \exp (ny), \end{equation}
so the velocity components are
\begin{align} u_x(x,y,t)=\frac {\partial \phi }{\partial x} &= -\sum _n n\,\alpha _n(t)\sin (nx)\exp (ny) ,\nonumber \\ u_y(x,y,t) =\frac {\partial \phi }{\partial y} &= \sum _n n\,\alpha _n(t)\cos (nx)\exp (ny) . \end{align}
There are
$M$
material points on the boundary, with coordinates and potentials
$X_i(t)$
,
$Y_i(t)$
,
$\varPhi _i(t)$
. Initially, these points are taken to be uniformly spaced on the
$x$
-axis between
$-\pi$
and
$\pi$
. The initial impulse given to these material points is per (1.2), so that when
$\epsilon =0$
,
$\alpha _n(0)=\delta _{n,1}$
in (A1). Values used for results presented were
$N_c=18$
(the number of cosine Fourier components) and
$M=401$
(the number of material points).
Because (A1) is automatically a solution, there is no need to solve the Laplace equation. The numerical method used to evolve the motion of the boundary (both the positions and potentials of the material points on the boundary) is based on the approach of Longuet-Higgins & Cokelet (Reference Longuet-Higgins and Cokelet1976), in which the boundary evolves according to (2.1). Given that the coefficients
$\alpha _n$
are known at some initial time
$t$
, the new state of the boundary (coordinates and potentials) is given by
\begin{align} X_i^\prime &= X_i+\left (\frac {\partial \phi }{\partial x}\right )_i\delta t \ , \ \ \ \ Y_i^\prime =Y_i+\left (\frac {\partial \phi }{\partial y}\right )_i\delta t, \nonumber \\ \varPhi _i^\prime &= \varPhi _i+\frac {1}{2}\left (\left (\frac {\partial \phi }{\partial x}\right )_i^2+\left (\frac {\partial \phi }{\partial y}\right )_i^2\right )\delta t, \end{align}
for
$i=1,\ldots ,M$
where (A2) have been used to evaluate the velocity components on the boundary points
$X_i,Y_i$
. The evolved data (A3) are then used to determine the new Fourier representation at time
$t+\delta t$
\begin{equation} \varphi _i \big( \big\{\alpha _n^\prime \big\} \big)=\sum _{n=0}^{N_c} \alpha _n^\prime \cos \left(nX_i^\prime \right) \exp \left(nY_i^\prime \right)\!, \end{equation}
where the new Fourier coefficients
$\alpha _n^\prime$
are determined by a least-squares fit of
$\varphi _i$
to the evolved data
$\varPhi _i^\prime$
, that is by numerical minimisation of the functional
\begin{equation} f \big( \big\{\alpha _n^\prime \big\} \big)\equiv \sum _{i=1}^M\big(\varphi _i\big(\big\{\alpha _n^\prime \big\}\big)-\varPhi _i^\prime \big)^2 \ ,\ \ \end{equation}
where the earlier values
$\alpha _n$
were used by the minimisation routine as the initial guess. Once the new coefficients
$\alpha _n^\prime$
are determined, the substitutions
$\alpha _n^\prime \to \alpha _n$
and
$(X_i^\prime ,Y_i^\prime )\to (X_i,Y_i)$
in (A1) and (A2) give the new boundary point potential and velocities. Iteration of this process gives a numerical approximation for the boundary height function
$y_s(x,t)$
, although this method was found to become unstable at around
$t\approx 2$
when the jet has extended significantly beyond the body of fluid, and an extension to our method described in the following sub-section was employed to circumvent this problem.
Note that, although we are primarily concerned with the evolution of the free surface of the fluid, the numerical method employed gives a velocity potential at each time step that applies to the whole body of the fluid and this also may be used to determine positions, velocities and other characteristics of fluid elements at any depth below the surface.
A.2. Two-region parametrisation of boundary
The method for propagating the boundary requires information about the solution of the Laplace equation in the interior region. After the jet extends a significant distance above the surface, it is easier to give an accurate representation of
$\phi (x,y,t)$
if the boundary and the interior are divided into two regions, which will be referred to as the base and the jet. One benefit of this scheme is that it extends the time of numerical stability, considerably beyond
$t=2$
.
The interior potential is determined by a different method for the base (
$M_1$
points) and the jet (
$M_2$
points). The two sets are defined by overlapping threshold values for
$Y$
, namely
$Y_1$
,
$Y_2$
, with
$Y_2\gt Y_1$
. The jet region is defined by
$Y\gt Y_1$
and the base by
$Y\lt Y_2$
. We found the optimal values of
$Y_1$
and
$Y_2$
to be
$0.01$
and
$3.5$
, respectively, and these are used in the results presented. The predicted increments of
$X_i$
,
$Y_i$
and
$\varPhi _i$
are different for the two regions, and the actual increment that is applied is an interpolation
\begin{equation} \delta \varPhi _i=G\left (\frac {Y_i-Y_1}{Y_2-Y_1}\right )\delta \varPhi ^{(2)}_i +\left [ 1-G\left (\frac {Y_i-Y_1}{Y_2-Y_1}\right ) \right ]\delta \varPhi ^{(1)}_i, \end{equation}
where
$\delta \varPhi ^{(1)}_i$
and
$\delta \varPhi ^{(2)}_i$
are the increments predicted using the base and jet potentials, respectively. (Analogous interpolations are used in the process of determining the
$X_i$
and
$Y_i$
variables.) For the interpolation function
$G(x)$
, we used
\begin{equation} G(x)=\left \{ \begin{array}{cc} 0 & x\lt 0 ,\cr 2x^2 & 0\lt x\lt 1/2, \cr 4x-2x^2-1 & 1/2\lt x\lt 1,\cr 1 & x\gt 1 .\end{array} \right . \ \end{equation}
In the initial stages of the motion, the maximum height of the boundary,
$Y_{{max}}$
, is less than the first threshold,
$Y_1$
, so all material points on the boundary move according to the base potential expressed in the Fourier series, as described in the previous section. When
$Y_{{max}}$
starts to exceed
$Y_1$
those points with
$Y_i\gt Y_1$
contribute to the jet potential (described below), and are mixed in according to (A6).
In the jet region, at any given time, we assume that the field on the central axis,
$x=0$
, is given by a function
$F_0(y)$
, which can be expressed as
\begin{equation} F_0(y)=\sum _{k=0}^{N} \phi _k \left (\frac {y-Y_1}{Y_{ {max}}-Y_1}\right )^k . \end{equation}
By symmetry, away from the axis
$\phi (x,y)$
is an even function, which we expand as a series
\begin{equation} \phi (x,y)=\sum _{n=0}^{(N-1)/2}F_n(y)x^{2n} . \end{equation}
Imposing that
$\boldsymbol{\nabla} ^2\phi =0$
, we find
$F_1(y)=-F_0^{\prime \prime}(y)/2$
,
$F_2(y)=-F_1^{\prime \prime}(y)/12=F_0^{\prime \prime \prime \prime}(y)/ 24$
, etc. so that
\begin{equation} \phi (x,y)\approx F_0(y)-\frac {1}{2}F_0^{\prime \prime}(y)x^2+\frac {1}{24}F_0^{\prime \prime \prime \prime}(y)x^4 + \cdots . \end{equation}
We do a least-squares fit of (A10) to the given values of
$X_i$
,
$Y_i$
and
$\varPhi _i$
to determine the ‘jet’ coefficients
$\{\phi _k\}$
. This fit includes only those boundary points with
$Y_i\gt Y_1$
.
Here, the coefficients
$\phi _k$
may become independent of time in the long-time limit. Because the velocity profile becomes linear and we expect
$F_0(y)$
is well approximated by a quadratic, so that the
$\phi _k$
become small for
$k\geqslant 2$
. In fact we did not find this to be the case, apparently due to the least-squares fit being ill conditioned. In practice, we obtained the most consistent numerical results when we included terms up to
$\phi _{9}$
i.e.
$N=9$
. After approximately
$t=9$
the method became unstable. The numerical work was coded with Maple, and the high-precision arithmetic used by that package may have been beneficial.
Finally, when
$Y_{{max}}$
exceeds the second threshold,
$Y_2$
, the least-squares fit for the base (A5) includes only those points with
$Y_i\lt Y_2$
.
A.3. Alternative Delauney/functional method to determine the boundary in the
$(u,v)$
plane
We used an alternative method in order to validate our numerical methods. In the
$(u,v)$
plane our fluid is initially bounded by a unit circle. Here, we also used a completely different approach to solving Laplace’s equation at each time step by finding a field
$\phi (x,y)$
which minimises the functional
\begin{equation} S=\frac {1}{2}\int \textrm {d}x \int \textrm {d}y\ ({\boldsymbol{\nabla }}\phi )^2 . \end{equation}
To find a tractable form of
$S$
we use a Delaunay triangulation of the region containing the fluid (initially a unit circle), defining
$\phi (x,y)$
by its value on the vertices of each triangle,
$\phi _i$
(internal vertex, free to vary) and
$\psi _i$
(boundary vertex, fixed). The functional is then approximated by the quadratic form
\begin{equation} S=\frac {1}{2}\sum _{\textit{ij}, \textit{int}} A_{\textit{ij}} \phi _i \phi _j +\sum _{i, {\textit{int}}, j, {\textit{bdy}}} A_{\textit{ij}}\phi _i\psi _j +\frac {1}{2}\sum _{\textit{ij}, {\textit{bdy}}}A_{\textit{ij}}\psi _i\psi _j, \end{equation}
where the coefficients
$A_{\textit{ij}}$
are sums of contributions from each triangle. Minimising
$S$
is equivalent to solving a system of linear equations,
$\sum _{j, {int}} A_{\textit{ij}} \phi _j+\sum _{j, {bdy}}A_{\textit{ij}}\psi _j=0$
, to determine the internal potentials
$\phi _j$
. These then allow us to calculate gradients for the boundary points, which can then be advanced by one time step. The process is then repeated, beginning with a new Delaunay triangulation of the new fluid-filled region, for each time step.
As mentioned earlier, the
$u$
values in the jet tip region grow extremely rapidly in the
$(u,v)$
plane, and the calculation exceeds the capabilities of a standard computer after a relatively short time. We managed to project the boundary to just over
$t=1$
, much less than the time managed by our principal method described earlier in this appendix. Nevertheless, the boundary predicted by these two very different methods agrees, up to the point where the Delauney/functional method fails. Figure (11) shows the boundaries up to
$t=1$
as predicted by both numeric methods with the Delaunay/functional method results transformed back to the
$(x,y)$
plane.
Note that this second numerical method also provides details of position and velocity of fluid elements relating to discretised points below the free surface, although again, our primary interest is in the free-surface profile evolution.
Figure 11. Boundary at intervals up to
$t=1$
as determined numerically by the ‘Fourier series’ method (red diamonds) and the ‘Delaunay/functional’ method (black lines).
Appendix B. Initial pressure pulse
Here, we determine the initial pressure
$p(x,y,0)$
. Throughout the fluid the pressure satisfies
\begin{equation} \frac {\partial \phi }{\partial t}+\frac {1}{2}(\boldsymbol{\nabla }\phi )^2+p(x,y,t)=\varPhi (t), \end{equation}
and the potential is
\begin{equation} \phi (x,y,t)=\sum _{n=0}^\infty \alpha _n(t)\cos (nx)\exp (ny), \end{equation}
with the only non-zero coefficients at
$t=0$
being
$\alpha _1(0)=1$
,
$\alpha _2(0)=\epsilon$
, and possibly
$\alpha _0$
. At
$t=0$
\begin{equation} (\boldsymbol{\nabla }\phi )^2=\exp (2y)+4\epsilon \cos (x)\exp (3y)+4\epsilon ^2\exp (4y) . \end{equation}
Because
$p=0$
at the surface (initially
$y=0$
), (B2) and (B3) can only be consistent with (B1) if
$\dot \alpha _1(0)$
is non-zero (while
$\dot \alpha _n(0)=0$
for
$n\geqslant 2$
). Hence we find
\begin{equation} p(x,y,0)=\frac {1}{2}\big [1+4\epsilon ^2-\exp (2y)-4\epsilon ^2\exp (4y)\big ] +2\epsilon \exp (y)\cos (x)\left [1-\exp (2y)\right ]\! . \end{equation}
After the initial impulse which sets the fluid in motion, the pulse of pressure initially accelerates all of the surface elements upwards, as is evident from figures 2 and 3. However, we were not able to find any relation between the tip velocities in figure 3 and gradients of the initial pressure field.
Open access
