2. The skimming model

The model used here is essentially that of Hicks & Smith (Reference Hicks and Smith2010). The body of mass $m_D$ and length $2L_D$ is assumed to be thin and to impact upon the water obliquely at a small angle (say $\phi$, to be defined explicitly later) such that its horizontal velocity is $U_D$, with the subscript $D$ denoting a dimensional quantity. The trailing edge impacts first. The assumption that the trailing edge touches down first may, for some specific underbody shapes $F$, rule out some ranges of $\theta$ (usually small values). The configuration is sketched in figure 1 in non-dimensional terms with the vertical scale stretched.

Figure 1. In non-dimensional terms, a schematic of the thin body with a sharp trailing edge skimming on a liquid layer. The solid arrows indicate the flow direction in a frame of reference in which the body does not appear to move horizontally. Its centre of mass is at height $d_1 + d_2 Y$ and varies with time. The leading contact position, $x_1$, also varies with time whilst the trailing edge of the body and of the wetted (pressured) region, $x_0$, remains fixed. The subscript L refers to the liquid layer. (The figure is not to scale and the small angles of incidence and inclination are accentuated.)

We work in a coordinate frame ($x$ horizontal, $y$ vertically upward) moving horizontally with the horizontal velocity component of the centre of mass of the body, so that the water layer flows with unit speed horizontally towards the station $x = -1$ of the leading edge of the body, where the Cartesian coordinates $\boldsymbol {x}$, the time $t$, the fluid velocity vector $\boldsymbol {U}$ and the pressure variation $P$ are non-dimensionalised with respect to $L_D, L_D /U_D, U_D, \rho _D U_D^2$, respectively. Here, $\rho _D$ is the constant fluid (water) density $\rho _w$ referred to in the Introduction, the atmosphere (air) is taken to be dynamically negligible during contact of the body with the water surface and the atmospheric pressure is taken as zero without loss of generality. The constant unit approach speed over the time scales of current interest is due to the relative smallness of the drag forces acting on the body surface; the horizontal component of velocity thereby remains almost constant (and much greater than the vertical component), at least until much longer time scales may come into play.

There are several additional features of the liquid-layer–body skimming interaction displayed in figure 1. Many of these are explained further throughout the paper; however, we provide a brief overview of the salient points here. Within the frame of reference centred on the horizontal motion of the body, the oncoming water is on the left in the figure and the wake on the right. (We emphasise that figure 1 is a sketch and is not to scale. In particular, the wake shape drawn is schematic; precise shapes are given numerically in figure 5 of Hicks & Smith Reference Hicks and Smith2010.) The origin of coordinates used here is on the wall ($y = 0$) directly below the centre of mass. Only the underbody shape is directly influential within the skimming interaction and it is defined by a curve $F(x)$ within an $O(d_2)$ variation in shape. In scaled quantities, with the base at $y = 0$, the centre of mass has $d_1 y = d_1 + d_2 Y(t)$ and the water layer height has $d_1 y = d_1 + d_2 h(x,t)$. Here, $d_1$ is the undisturbed water layer thickness while $d_2$ is the typical disturbance to the thickness, of order $d_1$ or smaller; Appendix A includes discussion of the case where $d_2$ is comparable to $d_1$ whereas the majority of the paper is for the case $d_2 \ll d_1$. By comparison, the overbody dynamics is negligible such that the overbody can take the form of an arbitrarily smooth curve (within the limits of scalings below), although the shape of the overbody does influence the moment of inertia and hence the rotation of the body. Ahead of the contact point, $x_1(t)$, the flow can reverse leading to a thin splash jet that may extend ahead of the body. Furthermore, surrounding the moving position $x_1(t)$, there is a ‘square’ Euler region (with its $x$-wise extremities marked by long dashed lines in figure 1) that remains close to the contact point and describes turnover, upstream of which the oncoming water layer is undisturbed to leading order. This effective contact point is sometimes called the turnover point or jet root; herein we will use ‘contact point’ to mean ‘effective contact point’. The main region of interest (where the liquid-layer–body interaction occurs) is the so-called wetted region or pressured region that is given by ($[x_1(t), x_0]$) where $-1 < x_1(t) < x \leq x_0 = 1$. Of note, the pressure $P$ is zero (atmospheric) on all the free surfaces that lie outside the wetted region. Since the splash jet on the body is one of those free surfaces, the ‘wetted region’ defined by $x_1(t) < x < 1$ might be described better as the ‘pressured region’ but we shall keep to the conventional term of wetted region below. Finally, the schematic diagram in figure 1 is given in scaled terms as explained next.

In the regime of current interest, the governing equations (as in Hicks & Smith Reference Hicks and Smith2010) can be linearised using the following reasoning and related scales. The water layer is taken to be thin in the sense that its non-dimensional thickness $d_1$ is small compared with the body length of unity. The underbody shape variation is with respect to $x$ only and is of order $d_2$ which is also small such that $d_2$ is small relative to $d_1$, while the angle of attack $d_2 \theta$ and the angle of inclination $\phi$ are likewise taken to be of order $d_2$. So the scaled angle $\theta$ is of order unity. Thus $d_2 \ll d_1 \ll 1$. In consequence, we can expect the majority of the water flow under the body to be only slightly disturbed from its oncoming state of unit velocity. We observe that, strictly, the derivation of the governing equations and associated conditions here is based on the nonlinear thin-layer formulation described by Hicks & Smith (Reference Hicks and Smith2010), using much of the flow structure in Tuck & Dixon (Reference Tuck and Dixon1989), but followed by linearisation as in the former paper: the derivation is summarised in Appendix A.

For the fluid flow, then, the governing equations are the linearised shallow-water equations. Gravity and capillary effects are supposed negligible by virtue of large Froude and Weber numbers while viscous effects can also be neglected for large Reynolds numbers provided the flow remains attached, particularly along the underbody and along the bottom of the water layer. The corresponding asymptotic expansion of the flow solution has the form

(2.1a,b)\begin{equation} \boldsymbol{U} = (1 + d_2 u / d_1, d_2 v) + \cdots,\quad P = d_2 p /d_1 + \cdots, \end{equation}

with $\boldsymbol {x} = (x, d_1 y)$ for $t$ of order unity. Substitution into the Navier–Stokes or Euler equations (here, to re-emphasise, viscous effects are supposed negligible) then yields the horizontal momentum equation

(2.2a)\begin{equation} u_t + u_x ={-} p_x, \end{equation}

while the vertical momentum shows that $p = p(x, t)$ is independent of $y$. The vorticity here is zero. The continuity equation coupled with the kinematic boundary condition at the surface of the water $y = 1 + d_2 h(x,t) /d_1$ and the lack of vorticity at leading order yields the balance

(2.2b)\begin{equation} h_t + h_x + u_x = 0. \end{equation}

Here, $h$ is zero at the incident wall-layer height. Equations (2.2a) and (2.2b) hold for $x_1 < x < x_0$, the unknown moving location $x = x_1(t)$ being the upstream contact point and the trailing-edge location $x_0$ being fixed at unity. The contact conditions involve unknown jumps of $u$, $p$ from their oncoming values of zero such that

(2.2c)$$\begin{gather} p + (1 - x_1^\prime)u = 0 \quad \text{at } x = x_1^+, \end{gather}$$

(2.2d)$$\begin{gather}u + (1 - x_1^\prime)h = 0 \quad \text{at } x = x_1^+, \end{gather}$$

from Tuck & Dixon (Reference Tuck and Dixon1989) and Hicks & Smith (Reference Hicks and Smith2010) (the prime notation indicates differentiation with respect to time here and throughout). We observe the present model for the fluid flow is in addition the same as that of Case 3 in Howison et al. (Reference Howison, Ockendon and Oliver2004). A thin jet of water is also provoked along the underbody ahead of the contact point. The requirement at the sharp trailing edge is

(2.2e)\begin{equation} p = 0 \quad \text{at } x = 1, \end{equation}

which follows from the Kutta condition.

As regards the body motion, the underbody surface moves according to the relation (holding for the present small displacements)

(2.3a)\begin{equation} h(x, t) = Y(t) + x \theta (t) + F(x), \end{equation}

where $Y(t)$ is the perturbation height of the centre of mass, $\theta (t)$ is the finite scaled rotation angle (the real unscaled angle, to repeat, is $d_2 \theta$, as noted earlier, which is small) and $F(x)$ denotes the fixed shape of the underbody. For ${\rm d}^2F(x)/{{\rm d}x}^2 < 0$ the underbody is concave and conversely the underbody is convex for ${\rm d}^2F(x)/{{\rm d}x}^2 > 0$. The centre of mass is taken to be at $x = 0$ midway between the leading and trailing edges of the body, as a central example. We should remark that the $x$-wise variation of the underbody shape is supposed to be $O(d_2)$ typically in a successfully completed skip, whereas the $x$-wise variation in body thickness plays virtually no direct role in the current interaction. Only the wetted (pressured) part of the underbody is subjected to the fluid flow pressure of (2.1a,b) whereas the unwetted part suffers merely the air pressure of zero. Hence the body-motion equations are

(2.3b)$$\begin{gather} M Y^{\prime \prime} = \int_{x_1}^1 p(x,t)\,{{\rm d}x}, \end{gather}$$

(2.3c)$$\begin{gather}I \theta^{\prime \prime} = \int_{x_1}^1 x p(x, t)\,{{\rm d}x}, \end{gather}$$

describing the vertical and rotational motions in turn. The working here is of course in the non-inertial frame of reference with the origin remaining in the same vertical position but moving horizontally at the horizontal speed of the body. The scaled mass $M$ and moment of inertia $I$ are related to their dimensional counterparts $m_D$, $i_D$ by means of $M = m_D d_1 /( \rho _D L_D^2), I = i_D d_1 /( \rho _D L_D^4)$, respectively. The acceleration/force balance for the horizontal motion in contrast confirms that the body velocity in the $x$-wise direction remains constant over the current time scales.

There is an interesting question of the relation between the present approach and previous models, especially concerning certain papers, e.g. Wagner (Reference Wagner1931), Tuck & Dixon (Reference Tuck and Dixon1989), Howison et al. (Reference Howison, Ockendon and Oliver2004), Korobkin (Reference Korobkin2004) and Hicks & Smith (Reference Hicks and Smith2010). Our approach is based on Tuck & Dixon (Reference Tuck and Dixon1989) as extended to unsteady interaction by Hicks & Smith (Reference Hicks and Smith2010) and we work totally within that framework of shallow water as described in Appendix A. That framework includes an exact form of the conditions holding at contact points and those conditions enable us to proceed readily to the present conditions (2.2c) and (2.2d). The current investigation is then over the body length scale on which full interaction takes place. It is recognised that once linearisation is applied Wagner properties hold on a smaller length scale buried inside the Euler region but on the other hand those properties have negligible effect on the all-important jump conditions (2.2c) and (2.2d) at leading order. We should repeat also our main concern is with a freely moving body and corresponding novel issues.

In terms of a problem statement, the skimming problem involves solving the system (2.2) and (2.3) for the unknowns $u(x, t), p(x, t), h(x, t), x_1(t), Y(t), \theta (t)$ subject to prescribed values of the initial height $Y(0)(= Y_0\ {\rm say})$, vertical velocity $V = V_0 = Y^\prime (0)$, initial angle $\theta (0)$ ($= \theta _0$ say) and rotation velocity $\omega = \omega _0 = \theta ^\prime (0)$ at time zero. The initial touchdown being at the trailing edge implies that $Y_0 = - \theta_0,\ {\rm with}\ \theta_0$ being negative, and at the trailing edge the perturbation height $h$ is initially zero, while $x_1$ is initially unity. The overall problem is nonlinear because of the conditions in (2.2c) and (2.2d). This, to repeat, is in essence the model of Hicks & Smith (Reference Hicks and Smith2010) who presented solutions of the system for a variety of parameter values. The special case of a flat underbody can be treated more analytically as shown by Hicks & Smith (Reference Hicks and Smith2010), leading to a nonlinear system of ordinary differential and algebraic equations. In the current study, however, a somewhat broader argument based on pressure and velocity as in the present section will be preferred since it holds for more general body shapes. Our concern in the remainder of this paper is with the skimming response as the scaled mass and moment of inertia are increased.

3. Numerical solutions for increasing mass

Numerical solutions of (2.2) and (2.3) were obtained using a seventh- and eighth-order Runge–Kutta adaptive step-size integration (Fehlberg Reference Fehlberg1968; Palmer & Smith Reference Palmer and Smith2020). Importantly, there are three geometric conditions that are obtained at leading order from the small-time analysis (Hicks & Smith Reference Hicks and Smith2010; Palmer & Smith Reference Palmer and Smith2020) that must be met to initialise the numerical solution. Firstly, $Y_0 + x_0 \theta _0 = 0$ must be fulfilled. This is a constraint on the position of the body as it enters the initially undisturbed water layer. Secondly, $D_0 = x_0^2 \omega _0/2 + x_0 (V_0 + \theta _0 + a_0)$ must also be met; $D_0$ is the initial fluid velocity at $x = 0$ and is here related to the initial trailing-edge position of the body, angular momentum of the body, velocity of the body, inclination of the body and the constant term of the underbody shape function. Thirdly, there is a condition for the position of the contact point at small time given by

(3.1)\begin{equation} x_1^\prime(0) = \frac{-3(x_0 \omega_0 + V_0) + \sqrt{9(x_0 \omega_0 + V_0)^2 + 8 \theta_0 (x_0 \omega_0 + V_0)}}{4 \theta_0}, \end{equation}

as in Hicks & Smith (Reference Hicks and Smith2010) and Palmer & Smith (Reference Palmer and Smith2020). The limit case of a flat underbody,

(3.2)\begin{equation} F(x) = 0, \end{equation}

is used in all the computational results in this section.

We first focus on the following representative initial conditions for each case: $Y_0 = 4, V_0 = -1$, scaled angle $\theta _0 = -4, \omega _0 = 0$ with $D_0 = x_0^2 \omega _0/2 + x_0 (V_0 + \theta _0)$ and $M = 3I$, with varied $I = 1,2,3,4$. It is noteworthy that physical sense requires the limitation $I < M$ because of the definitions of $M$ and $I$ in terms of integrals.

The initial $Y_0$ and $\theta _0$ values here are each slightly larger than the cases shown in Hicks & Smith (Reference Hicks and Smith2010) to allow for larger values of $M$ and $I$ before flooding (where fluid is able to run over the top of the body, Smith & Liu Reference Smith and Liu2017) and possibly sinking. Results for $M = 3I$ are presented since they capture the typical dynamics of the skimming motion. Physically, the moment of inertia for a thin flat body with uniform density can be at most $M$ (Liu Reference Liu2017). However, in many physical applications, the body spins around its vertical axis which stabilises its contact angle with water via gyroscopic effect. This effect can be modelled in two dimensions by increasing the body's moment of inertia to be greater than what is seemingly possible for a flat body (Liu Reference Liu2017). When $M = 2I$ (for comparison), the body can complete successful skims for larger values of $I$, whilst larger ratios cause the body to sink for smaller values of $I$. In each of these cases, however, the same dynamics is typically seen (as indicated by the analysis in § 4).

The evolution of the leading contact position of the wetted region $x_1$ throughout an impact is shown in figure 2 to highlight solution trends as $M$ and $I$ increase. The same qualitative behaviour is seen in each case with the leading contact position being initially close to the trailing edge at touchdown and immediately moving in the negative $x$-direction as the body descends into the water layer. Eventually a maximum wetted area is reached, after which the body tends to rebound, if $M$ and $I$ are not too large for the given initial conditions, with the leading contact position moving downstream and returning to the trailing-edge position as water exit is approached. The physical limitation of the rebound is that $x_1 \in [-1,1 ]$ since $x = - 1$ marks the leading edge of the body itself.

Figure 2. Evolution of $x_1$ as a function of time for varying $M = 3I$, $I = 1,2,3,4$. The linear behaviours at small times and for exit (dashed lines) are shown. See table 1 for the other parameter values.

Table 1. Table of control parameters for the numerical investigations.

The time to body exit increases with increasing $M$ and $I$, as does the extent of body wetting as measured from $x_1$ to $x_0 (=1)$. When $I = 4, M = 12$ (indeed, when $I > 3, M > 9$) the body is predicted to flood and quite likely sink (Smith & Liu Reference Smith and Liu2017) as shown by the entire underbody extent becoming wetted ($x_1(t) < -1$): the modelling of the skimming motion is then no longer considered to be physical and thus the numerical solution is terminated. Overall, as a check, the small time response in (3.1) holds well throughout the body entry, with the solution for the body's trajectory remaining closer to the asymptote for longer as the mass increases. The exit solution in the computations, however, as $x_1$ returns to unity, at time $t_e$ say, only matches the analytical exit form

(3.3)\begin{equation} {x}_1-1 ={-} (t_e - t) + O((t_e - t)/ \ln(t_e - t)), \end{equation}

over shorter time frames with increasing $M$ and $I$ (a phenomenon that involves a delicate effect of the logarithmic dependence as analysed in § 5). This dynamics is noticeably different from that examined in Palmer & Smith (Reference Palmer and Smith2020) (in brief, for larger mass: the skimming motion is longer and deeper as measured by $x_1$ vs time, the body undergoes a super-elastic response (see figure 4), and the small-time entry solution holds for longer). In addition, we have found no examples where $dx_1/dt$ ever exceeds unity. This tends to suggest that the contact point never recedes in the laboratory frame of reference.

Now we consider the contact-point pressure at $x_1$ presented in figure 3. In general, the pressure initially increases to a maximum and after that it falls. As exit is approached this decline becomes slower until the body is on the verge of exit, at which point there is rapid change in the pressure. This occurs for all values of $M$ and $I$ and is in keeping with the trends found in Hicks & Smith (Reference Hicks and Smith2010) and Palmer & Smith (Reference Palmer and Smith2020). As $M$ and $I$ increase the qualitative shape of the pressure profile changes, with pressure increasing for a longer period of time as the body descends further into the liquid layer. For the largest $M$ and $I$ for which rebound occurs $(M = 9,\ I = 3)$, a peak is reached, after which the pressure falls as before; however, there is a quantitative change in the profile and a larger positive pressure is sustained, declining more slowly. The consequences of these differences in pressure are seen in figures 4 and 5, regarding the vertical and rotational progression of the body, respectively. When $I = 4,\ M = 12$ the body floods and the pressure solution terminates at the time when $x_1(t) < -1$. Up until this juncture the trend in the curves for increasing mass continues.

Figure 3. Evolution of the pressure on the body at the moving contact point $x_1$ as a function of time for varying $M = 3I$, $I = 1,2,3,4$. See table 1 for the other parameter values. The solution curves here terminate at the times when either exiting or flooding occurs.

Figure 4. (a) Evolution of $Y$ as a function of time. (b) Evolution of $V$ as a function of time. Each plot is for varying $M = 3I$, $I = 1,2,3,4$. See table 1 for the other parameter values. The solution curves here terminate at the times when either exiting or flooding occurs.

Figure 5. (a) Evolution of $\theta$ as a function of time. (b) Evolution of $\omega$ as a function of time. Each plot is for varying $M = 3I$, $I = 1,2,3,4$. See table 1 for the other parameter values. The solution curves here terminate at the times when either exiting or flooding occurs.

In figure 4, the vertical height of the body's centre of mass, $Y$, initially falls linearly as the downward momentum carries the body into the liquid layer. The body then begins to rise out of the water as the pressure on the wetted region increases, see figure 3. Meanwhile, the (initially negative) vertical velocity increases throughout the majority of the skimming motion even as the pressure falls away until some point is reached where the velocity becomes near constant. When $M$ and $I$ are comparatively small, the body leaves the water at a lower height than it entered. However, as $M$ and $I$ increase there is a change in the dynamics and a super-elastic response, as previously shown in Liu (Reference Liu2017) and Palmer & Smith (Reference Palmer and Smith2020), occurs with the body leaving the water layer at a far greater height and magnitude of velocity than it entered. This super-elastic response is a result of the nonlinear deformation of the liquid layer where a relatively small portion of the large horizontal kinetic energy of the liquid layer is converted into (relatively small) vertical kinetic energy of the body (through the body's angle of inclination and the developing splash jet); see also Hewitt et al. (Reference Hewitt, Balmforth and McElwaine2011). In the laboratory frame, while the body can rise with a vertical speed/angle greater than that with which it impacts, its horizontal momentum must be converted into vertical momentum by means of the fluid flow and so the body loses some horizontal velocity, albeit a small amount of higher order.

These differences in the entry and exit height of the body can be understood when considered in parallel with the rotational behaviour of the body, as shown in figure 5. In each case here there is initially an anti-clockwise moment on the body due to a positive pressure over the wetted region as it progresses through the water layer. This causes the body's angle of inclination to increase throughout its skimming motion and the angular velocity to simultaneously increase rapidly. Again, for the lower $M$ and $I$ cases, the trend is consistent with the change in angle being monotonic and the angular velocity tending towards a constant as exit is approached. By contrast, for increased $M$ and $I$ the body is more resistive to the rotational torque and as such does not rotate as significantly. Indeed, the body's rotation reverses, becoming clockwise at some point, and thus a larger negative angle of trajectory is maintained.

Physically, in the smaller $M$ and $I$ cases, the body is descending into the water layer, undergoing a significant rotation which acts to lift the rear of the body upwards. This causes the leading contact position of the wetted region to retreat towards the trailing edge, producing a rapid decrease in the wetted region towards body exit. A swift change in pressure then arises during the short exit period with the vertical and angular accelerations diminishing toward zero. This helps to show why the exit asymptote only holds well within a short time frame.

In the larger $M$ and $I$ cases on the other hand the close geometrical relation between $Y$ and $\theta$ also contributes to the super elastic response and significantly large $Y,V$ values upon exit. Note that there is a larger disturbance of the water layer, producing an isolated wave like response in the water under the body. This movement of the water layer and the lack of substantial rotation causes contact to be maintained with the body. As such, a larger positive pressure persists across the body (figure 3) as the water layer deforms and continues to push against the body. The delayed exit and slower decline in pressure then leads to the super-elastic response, pushing the body out of the water with greater velocity and height. This is in accordance with observations from empirical studies (Hewitt et al. Reference Hewitt, Balmforth and McElwaine2011).

Finally, regarding the case when $I = 4,\ M = 12$ in figures 4 and 5, the curves again terminate when $x_1(t) < -1$, having shown the same overall trend as mass increases. For $I > 4, M > 12$ it is expected that the body trajectory will become more ballistic. As such, each of the solution curves in figures 2 and 3 will become increasingly linear (with corresponding simplifications in figures 4, 5), terminating at earlier times for $x_1(t) < -1$.

4. Analysis for large mass

Given the general trends of delayed exiting (or flooding) as $M$ and $I$ increase in the numerical results above, along with higher pressures and slower rotation, the detailed account below returns to consider the original interaction equations of § 2 when the body becomes heavier in the sense of increased $M$ and $I$. There are two temporal stages then. These are initially studied for a general underbody with fixed shape

(4.1)\begin{equation} F(x) = O(1), \end{equation}

with their implications for the flat underbody of § 3 to be discussed later.

4.1. Stage 1 – the majority of the skim

This stage covers most of the body's descent and ascent (if there is no sinking) in the water after entry when a distinguished time scale of $O(M^{1/3})$ brings in the lift and moment. As such, the orders of magnitude indicate that $x_1(t) - x_0 = O(1)$, with the time scaling $t = M^{1/3} t^*$. The particular scaling here follows from the expectation that only a relatively narrow band of initial body velocities can lead to skimming instead of flooding. As such, the height $Y$ and the angle $\theta$ remain equal and opposite constant values to a first approximation (following the small-time condition of $Y_0 + x_0 \theta _0 = 0$ with $x_0$ set equal to unity) but subject to small perturbations of unknown order $E$ say over the long time scale of the complete skim. It follows from (2.3a) that the depth $h$ is simply an $O(1)$ function of $x$ to leading order with explicit time dependence arising first in the order-$E$ perturbation, while the conditions at the trailing edge and moving contact imply that the flow velocity $u$ and the pressure $p$ have the same format as $h$. Hence both of the body-motion balances in (2.3b) and (2.3c) indicate that $M E(\sigma )^{-2}$ must be $O(1)$ to avoid a ballistic motion leading to flooding. Here, $\sigma$ is the unknown time scale. We anticipate further that the $E$-perturbations should also bring in the time derivatives in $h, u, p$, which suggests the balance $E \sim (\sigma )^{-1}$. Combining the two balances here yields $M(\sigma )^{-3}$ being $O(1)$ and therefore the time scaling of $M^{1/3}$. Moreover it is inferred that the wetting speed scales as $x_1^\prime = M^{-1/3}\,{\rm d} x_1/{\rm d} t^*$ in this stage.

With the given underbody function $F$ taken as $O(1)$ in general, we expect the variations in $p, h, u, x_1, Y, \theta$ to be as follows:

(4.2a)$$\begin{gather} (p, h, u, x_1) = ({p}_0, {h}_0, {u}_0, \bar{x}_1) + M^{{-}1/3}({p}_1, {h}_1, {u}_1, \overline{\bar{x}}_1) + \cdots, \end{gather}$$

(4.2b)$$\begin{gather}(Y, \theta) = ({-}s,s) + M^{{-}1/3}({Y}_0(t^*), {\theta}_0(t^*)) + \cdots. \end{gather}$$

The $O(1)$ constant $s$ is expected to be negative, corresponding to the centre of mass being above the water, while the small perturbations in $Y$, $\theta$ occurring over the present long time scale indicate the narrow band of impact velocities (of order $M^{-2/3}$) involved here.

We substitute the expansion into (2.2) and (2.3). At leading order, the unknown terms having subscript zero together with $\bar {x}_1$ are seen to satisfy the reduced equations

(4.3a)$$\begin{gather} {{h}_{0x}} + {{u}_{0x}} = 0, \end{gather}$$

(4.3b)$$\begin{gather}{{u}_{0x}} ={-}{{p}_{0x}}, \end{gather}$$

(4.4)\begin{gather} {h}_0 ={-}s + sx -{F}(x). \end{gather}

Thus, $h_0$ only depends on the $x$ value. The boundary conditions at leading order become

(4.5a)$$\begin{gather} \text{at the trailing edge}: \quad {p}_0(1, t^*) = 0, \end{gather}$$

(4.5b)$$\begin{gather}\text{at the contact point}: \quad {p}_0(\bar{x}_1,t^*) ={-}{u}_0(\bar{x}_1,t^*) = {h}_0(\bar{x}_1,t^*). \end{gather}$$

Hence, we have the simple solution

(4.6a)\begin{equation} p_0 = h_0 ={-} u_0 ={-} s + sx - F(x), \end{equation}

along with the body-motion balances, which now require

(4.6b)$$\begin{gather} {Y}_0^{\prime \prime} = \int_{\bar{x}_1}^1 {h}_0(x, t^*)\,{{\rm d}x}, \end{gather}$$

(4.6c)$$\begin{gather}I^* {\theta}_0^{\prime \prime} = \int_{\bar{x}_1}^1 x{h}_0(x, t^*)\,{{\rm d}x}, \end{gather}$$

where the $O(1)$ positive constant $I^* = I/M$ and $0< I^*<1$. This restriction is in keeping with the limitation on $I$ given after (3.2). The integral contributions in (4.6b) and (4.6c) and the perturbation form of (4.2a) and (4.2b) are the aspects that distinguish the present time scale. However, we are one equation short since, although (4.6a) determines $h_0$, (4.6b) and (4.6c) yield only two equations for three unknown functions of $t^*$, namely ${Y}_0, {\theta }_0$ and $\bar {x}_1$. We continue to the next order to derive the third equation.

At the next order explicit time dependence appears in the form

(4.7a)$$\begin{gather} {{h}_{0t^*}} + {{h}_{1x}} + {{u}_{1x}} = 0, \end{gather}$$

(4.7b)$$\begin{gather}{{u}_{0t^*}} + {{u}_{1x}} ={-}{{p}_{1x}}, \end{gather}$$

(4.8)\begin{gather} {h}_1 = {Y}_0(t^*) + {\theta}_0(t^*), \end{gather}

and although $h_0$ and $u_0$ are quasi-steady effects $Y_0$ and $\theta _0$ turn out not to be. The relevant boundary conditions are now as follows:

(4.9a)$$\begin{gather} \text{at the trailing edge}: \quad {p}_1 = 0; \end{gather}$$

(4.9b)$$\begin{gather}\text{at the contact point}: \quad {p}_1 + \overline{\bar{x}}_1 {p}_0^\prime(\bar{x}_1) ={-}{u}_1 - \overline{\bar{x}}_1 {u}_0^\prime(\bar{x}_1) + \bar{x}_1^\prime {u}_0(\bar{x}_1) , \end{gather}$$

(4.9c)$$\begin{gather}\text{at the contact point}: \quad {u}_1 + \overline{\bar{x}}_1 {u}_0^\prime(\bar{x}_1) ={-}{h}_1 - \overline{\bar{x}}_1 {h}_0^\prime(\bar{x}_1) + \bar{x}_1^\prime {h}_0(\bar{x}_1) . \end{gather}$$

Here, we have left out $t^*$ explicitly from $h_0, u_0$ in the view of (4.6a). The reduced equations, (4.7a) and (4.7b), imply that

(4.10a)$$\begin{gather} {u}_1(x, t^*) ={-}{h}_1(x, t^*) + U_1 (t^*), \end{gather}$$

(4.10b)$$\begin{gather}{p}_1(x, t^*) = {h}_1(x, t^*) + P_1(t^*), \end{gather}$$

where the flow velocity contribution $U_1(t^*)$ is to be found but the Kutta condition determines the pressure contribution here as $P_1(t^*) = -{h}_1(1, t^*)$. Then, when we impose (4.9b), (4.9c) and also use (4.8) we find much cancellation takes place and leads to the two equations

(4.11)$$\begin{gather} U_1(t^*) = (s - s \bar{x}_1 + F(\bar{x}_1))\bar{x}_1^\prime + Y_0 + \theta_0, \end{gather}$$

(4.12)$$\begin{gather}U_1(t^*) ={-}(s - s \bar{x}_1 + F(\bar{x}_1))\bar{x}_1^\prime, \end{gather}$$

for $U_1(t^*)$ and $\bar {x}_1^\prime (t^*)$ effectively. Hence we obtain the following:

(4.13)\begin{equation} {Y}_0(t^*) + {\theta}_0(t^*) = 2 \bar{x}_1^\prime(t^*) \lbrace -s + s\bar{x}_1(t^*) - F(\bar{x}_1)\rbrace, \end{equation}

as the third equation for $Y_0, \theta _0, \bar {x}_1$. Thus, the nonlinear problem in this stage reduces to solving (4.6b), (4.6c) and (4.13).

Now focusing on the basic case of a flat underbody surface where $F$ is zero, the three equations reduce to

(4.14a)$$\begin{gather} \frac{Y_0^{\prime \prime}(t^*)}{s} ={-} \frac{1}{2} (\bar{x}_1(t^*)-1)^2, \end{gather}$$

(4.14b)$$\begin{gather}\frac{I^* \theta_0^{\prime \prime}(t^*)}{s} ={-} \frac{1}{6} (\bar{x}_1(t^*)-1)^2 (2\bar{x}_1(t^*)+1), \end{gather}$$

(4.14c)$$\begin{gather}\frac{{Y}_0(t^*) + {\theta}_0(t^*)}{s} = 2\bar{x}_1^\prime(t^*) (\bar{x}_1(t^*)-1). \end{gather}$$

Combining these, a nonlinear ordinary differential equation (ODE) is obtained for the function $\bar {x}_1(t^*)$, namely

(4.15)\begin{equation} \left( \bar{x}_1^\prime (\bar{x}_1 - 1) \right)^{\prime \prime} = \left(\bar{x}_1 - 1\right) ( a_0 + a_1 \bar{x}_1 + a_2 \bar{x}_1^2 ), \end{equation}

where $a_0 = (3I^* + 1)/(12 I^*), a_1 = (1 - 3I^*)/(12I^*)$ and $a_2 = -1/(6I^*)$ are known $O(1)$ constants. For small times $t^*$ the solution $\bar {x}_1$ must behave as $1 - Bt^* + \cdots$ in order to match with the initiation form found in (3.1). In the present context $B$ is an arbitrary positive constant.

To simplify we put $\bar {x}_1 - 1 = - (b_1/b_2) Q^{1/2}, t^* = (-2b_1)^{-1/3} T$ with $b_1 = a_1 + 2a_2 < 0$ and $b_2 = a_2 < 0$. This leaves us with a nonlinear ODE for the dependence of $Q$ on $T$, namely

(4.16a)\begin{equation} d^3 Q/dT^3 ={-}Q + Q^{3/2}. \end{equation}

The initial condition to match the trend emerging from the small-time behaviour is

(4.16b)\begin{equation} Q \sim C T^2,\quad \text{for small } T, \end{equation}

with the initiation constant $C$ being a positive constant related to $B$ (and therefore to $\bar {x}_1^\prime (0))$. Thus we have a one-parameter system (4.16a) and (4.16b) to solve. The solution is plotted in figure 6(a) for varying $C$. Higher values of the initiation constant $C$ lead to a blowup of $Q$ within a finite scaled time, implying subsequent flooding and possible sinking of the body, whereas lower values point to $Q$ returning to zero within a finite time. The latter indicates an approach to a full rebound, i.e. an exit phenomenon. Notably the separatrix is at $C_{sep} = 0.1017179$. The mathematical and physical meaning of $Q$ (and its integral shown in figure 6b) as a function of scaled time will be explored further in § 4.2 as regards the exit and in § 5 where comparisons with the full model are discussed.

Figure 6. Large-mass analysis for stage 1. (a) Plots of $Q$ vs $T$ with varying values of constant $C$. (b) Plots of the integral of $Q$ with respect to $T$ vs $T$ with varying values of constant $C$. There is a separatrix for $C_{sep} = 0.1017179$. (c) Solutions to (4.14a)–(4.14c) for stage 1 perturbation values (4.2b) with $Y_0(0) = \theta _0(0) = \omega _0(0) = 0$ and $V_0(0) = -0.4973$ to produce a near maximal skim where $\bar {x}_1 \sim -1$ at the point of rebound.

Properties for small and large initiation constant $C$ are described in Appendix B. For $C$ less than some value $C_{max} < C_{sep}$ (notably $C_{max}$ is not the maximum value of $C$ overall, rather the largest value of $C$ for which a successful skim may occur) the most important behaviour to study is the response near liftoff as the contact-related function $Q \rightarrow 0^+$ at a finite time $T = T_{0^-}$ say. (Here, $C_{max}$ is the maximum value of $C$ that allows liftoff to ensue.) The response is linear in the sense that $Q \sim \alpha _1(F_0 - F)$ where the positive constant $\alpha _1$ can be derived numerically from the above solutions. However, it follows that, near a finite scaled time $t^* = t_{0^-}^*$, the behaviour is irregular in terms of $\bar {x}_1$, specifically

(4.17a)\begin{equation} \bar{x}_1 \sim 1-\alpha_2^{1/2}(t_0^*-t^*)^{1/2}, \end{equation}

where $\alpha _2 = (b_1/b_2)^2(-2b_1)^{1/3}\alpha _1$, and so the scaled contact-point velocity is

(4.17b)\begin{equation} \frac{{\rm d}\bar{x}_1}{{\rm d} t} = M^{{-}1/3} \frac{{\rm d}\bar{x}_1}{{\rm d} t^*} \sim \frac{1}{2} M^{{-}1/3} \alpha_2^{1/2} (t_0^* - t^*)^{{-}1/2} \quad \text{as } t^* \rightarrow t_{0^-}^*. \end{equation}

Therefore the relatively slow evolution of $\bar {x}_1$ during the above stage quickens significantly now and indeed it reaches $O(1)$ when $(t_0^* - t^*)$ reduces to the order $M^{-2/3}$. This defines the final exit stage.

4.2. Stage 2 – exiting from the water

The final stage has the contact point approaching the trailing edge closely and has the spatial and temporal scales $x = 1 + M^{-1/3} \hat {x}$, $t = M^{1/3}t_0^* + M^{-1/3} \hat {t}$, respectively. The solution now expands as

(4.18a)$$\begin{gather} (p, h, u, F) = M^{{-}1/3} (\hat{p}, \hat{h}, \hat{u}, \hat{F}) + \cdots, \end{gather}$$

(4.18b)$$\begin{gather}(Y, \theta) = ({-}s,s) + M^{{-}1/3} (\hat{Y}(\hat{t}),\hat{\theta}(\hat{t}))+\cdots, \end{gather}$$

(4.18c)$$\begin{gather}x_1 = 1 + M^{{-}1/3} \hat{x}_1 +\cdots, \end{gather}$$

this being the form inferred from the ending of the previous stage. The contact-point velocity here is of order unity.

Substituting (4.18a)–(4.18c) into the full problem we find $\hat {Y}(\hat {t}),\hat {\theta }(\hat {t})$ are constants, given by the matching at large negative $\hat {t}$, and $\hat {h} = \beta + s \hat {x}$. Here, the constant $\beta = \hat {Y} + \hat {\theta }$ is known. Working through, we obtain the following equation for $Z = \hat {x}_1(\hat {t})$:

(4.19)\begin{equation} Z Z^{\prime \prime} (\beta + sZ) +2s Z{Z^{\prime}}^2 - 2sZZ^\prime + \beta (Z^\prime - 1)^2 = 0. \end{equation}

Putting $\beta = sB$ (noting that we expect $\beta > 0, s < 0$ and so the constant $B$ is negative), we write $Z = -BR$, $\hat {t} = -B\bar{T}$ (so that $\bar{T}$ marches forward in time) to give us a nonlinear equation for $R(\bar{T})$, namely

(4.20a)\begin{equation} R R^{\prime \prime} (1 - R) - 2 R{R^{\prime}}^2 + 2RR^\prime + (R^\prime - 1)^2 = 0, \end{equation}

which is free of parameters. The boundary condition here is

(4.20b)\begin{equation} R \sim{-}|\bar{T}|^{1/2} \quad \text{as } \bar{T} \rightarrow -\infty, \end{equation}

from matching with the end of the previous stage in (4.17a). The solution is shown in figure 7. The match (4.20b) is observed at large negative $\bar{T}$ values in the figure while the approach to final exit is highlighted in the inset.

Figure 7. Large-mass analysis for stage 2. Graphs of $R$ (blue) and its derivative (red) vs scaled time $\bar{T}$, along with asymptotes (dashed) at large negative $\bar{T}$ values, from (4.20b) and its derivative, and asymptotes at near-exit times, from (4.23) and its derivative. The inset shows a close-up view very near exit.

At final exit, as $\bar{T}$ tends to some finite value $\bar{T}_E$, the contact-point position $R$ tends to zero. Therefore, to examine the behaviour of the system as exit is approached we integrate (4.20a) to

(4.21)\begin{equation} R^2(1 - R^\prime) + R(R^\prime - 2) = (\bar{T}_E - \bar{T}) ,\end{equation}

on the assumption that $R^\prime R$ tends to zero; the subsequent working shows the assumption is valid. As the contact point recedes, $R$ rapidly tends to zero and we can neglect the term with a factor $R^2$ on the left-hand side at leading order. Thus, transforming back to $Z$ and $\hat {t}$ the near-exit solution is governed by:

(4.22)\begin{equation} Z(Z^\prime - 2) = (\hat{t}_E - \hat{t}). \end{equation}

The solution to (4.22) can be found implicitly, showing the final exit behaviour to be

(4.23)\begin{equation} Z = \hat{x}_1 \sim{-}(\hat{t}_E - \hat{t}) \lbrace 1 + (\ln(\hat{t}_E - \hat{t}))^{{-}1}\rbrace + \cdots.\end{equation}

The implicit solutions of (4.22) and (4.23) are shown in figure 8. In particular, the comparison in figure 8 and the result (4.23) agree with and confirm the exit response shown in (3.3). The logarithmic form only appears in the present large-mass scenario when stage 2 is encountered, very near to the final exit, a feature that explains the last-gasp behaviour seen in the numerical results of the full system (2.2)–(2.3) in figures throughout § 3. See also figure 7 and its inset.

Figure 8. Exact solution of the local exit equation (4.22) (solid). The final asymptote (4.23) is also shown (dotted).

5. Comparisons

The analysis above demonstrates several interesting scales and trends present in the skimming process of a large-mass body that would benefit from further exploration and discussion. In § 4.1, it was shown that the skimming process as described by (2.2) and (2.3) reduces to a one-parameter ODE for a body of large mass (4.16a) and (4.16b). As shown in figures 6(a) and 6(b) the behaviour of the one-parameter system qualitatively changes with the constant $C$, from the initial condition (4.16b), with higher values ($C > C_{sep}$) resulting in a finite time blowup. In practice, prior to such blowup, flooding occurs (with $x_1 < -1$ indicating that the liquid layer completely covers and submerges the body). Flooding also applies for $C_{max} < C < C_{sep}$. Of particular interest here are the cases where $C < C_{max}$ and how the dynamics of the skimming motion is affected by the stage 1 time and body wetting scales (in § 4.1 the time scaling is shown to be $M^{1/3}$ whilst the extent of body wetting $x_1 - x_0$ remains $O(1)$). Indeed, from (4.16a) and (4.16b), it is possible to find initial conditions that cause bodies of varying mass to undergo skimming motions of near maximal wetting (i.e. the position of the contact point $x_1$ getting as close as possible to $-1$ at the body leading edge yet still completing a full skim).

Using the definitions of $Q, a_0, a_1, a_2, b_1$ and $b_2$ (presented between (4.15) and (4.16a)), this ‘maximal’ skim occurs when $x = x_{LE} = -1$ with the equivalent value for $Q$ given by $Q_{LE} = (16/9) (1 + I^*)^{-2}$. Following from (4.16b) it is possible to find the value $C_{max} < C_{sep}$ such that the maximum value of $Q$ equals $Q_{LE}$, representing a maximal skim in the full system. For example, in figure 6(a) the scaled moment of inertia is $I^* = 1/3$, and therefore we seek the value of $C_{max}$ that leads to the maximum value of $Q$ being $Q_{LE} = 1$ (we recall that in figure 6(a) the $Q$ solution was seen to have been able to overshoot $Q_{LE}$ for $C_{max} < C < C_{sep}$, but this results in the body flooding).

Once the value of $C_{max}$ has been obtained for a given $I^*$, the small-time behaviour in (4.16b) gives that

(5.1)\begin{equation} x_1 - 1 \sim \frac{3+3I^*}{2} \left(\frac{1+I^*}{2I^*}\right)^{1/3} C_{max}^{1/2} t^*. \end{equation}

Inserting this value into (4.15) provides a relationship between the initial conditions and $C_{max}$ as follows:

(5.2)\begin{equation} \frac{Y_0(t^*) + \theta_0(t^*)}{s} \sim 2 \left(\frac{3+3I^*}{2}\right)^2 \left(\frac{1+I^*}{2I^*}\right)^{2/3} C_{max} t^*. \end{equation}

Therefore, given two of the three initial conditions $s, Y_0(t^*)$ or $\theta _0(t^*)$ as known or chosen, the third may be calculated to produce a near maximal skimming motion. Additionally, $V_0(0)$ and $\omega _0(0)$ may be freely chosen to fulfil this relationship since, from (4.2b), $V_0(t^*) = M^{-2/3} Y_0(t^*)$ and $\omega _0(t^*) = M^{-2/3} \theta _0(t^*)$ may be substituted into the above. It should be noted that this method can be used to produce the asymptotic large-mass solution for any two given initial conditions and maximum $x_1$ value as desired.

In practice, numerically solving the system (4.16a) and (4.16b) is quick and thus iterative methods can be used to find $C_{max}$ for a given value of $I^*$. Using a seventh-/eighth-order time step adaptive Runge–Kutta method with small error tolerances we found that for $I^* = 1/3$, $C_{max} = 0.097899$. For $\omega _0 = 0$ and $s = -4$ we thus expect for maximal skim the value

(5.3)\begin{equation} V_0 \sim 2^{(11/3)} C_{max} s ={-}0.4973 M^{{-}2/3}. \end{equation}

The solutions to (4.14a)–(4.14c) are presented in figure 6(c) for the above parameters and are plotted against the scaled time $t^*$.

The near maximal underbody wetting is clearly seen in the $x_1$ profile. The profiles corresponding to $Y(t^*)$ and $V(t^*)$ show that the super-elastic response as previously seen in § 3 is captured by the reduced form in § 4. The super-elastic behaviour is enhanced with increased mass since this leads to more body wetting and a longer, sustained contact with the water layer. In addition the increased mass results show the interesting rotational dynamics from § 3 with the body's angle of attack increasing in the clockwise direction after rebound. Together with the trend in body velocity, the increasing clockwise body angle helps the body maintain contact with the liquid layer as the trailing edge rotates into the liquid layer. These results are indicative of an isolated wave and water jet building under the body as it rebounds, the size of which increases with the body's mass, and yields the sustained interaction with the water sheet.

Next, we make comparisons with the numerical findings for the full system (2.2)–(2.3). The value on the right-hand side in (5.3) is found to be marginally too negative for maximal skimming according to the numerical results over the range of $M/3$ between 10 and 50. Instead the value $V_0 = -0.4968 M^{-2/3}$ produces the desired maximal results over that range. This very small discrepancy concerning the difference between the asymptotic and the full solutions is likely due to the use of the small-time solution (from which the full system deviates before rebound) in determining (5.1), the limitation of numerical accuracies in both the iterative method to find $C_{max}$ and the numerical solution, or possibly the finite $M$-range mentioned. An $O(1)$ time shift also clearly plays a role. However, inserting the value of $C_{max}$ for $V_0= -0.4968$ into (4.16a) and (4.16b) we obtain a maximum value of $Q = 0.9975$, which is very close to maximal wetting indeed. Thus, presented in figures 9 to 11 are the results of (2.2)–(2.3) with the above initial conditions and $V_0 = -0.4968 M^{-2/3}$.

Figure 9. Evolution of $x_1$ vs $t^* = M^{-1/3}t$ for $M = 3I$ and varying $I = 10,20,30,40,50$. The initial conditions have been chosen to achieve the near maximal wetting of the underbody. The asymptotic solution is shown in each case (black curves).

Figure 10. (a) Value of $(Y-Y_0)M^{1/3}$ vs $t^*$. (b) Value of $VM^{2/3}$ vs $t^*$. Each plot is for $M = 3I$ and varying $I = 10,20,30,40,50$. The initial conditions have been chosen to achieve the near maximal wetting of the underbody. The asymptotic solution is shown (black curves). Also shown are $Y_0(t^*), V_0(t^*)$ from (4.2b), for comparison.

Figure 11. (a) Value of $(\theta -\theta _0)M^{1/3}$ vs $t^*$. (b) Value of $\omega M^{2/3}$ vs $t^*$. Each plot is for $M = 3I$ and varying $I = 10,20,30,40,50$. The initial conditions have been chosen to achieve the near maximal wetting of the underbody. The asymptotic solution is shown in each case (black curves). Also shown are $\theta _0(t^*), \omega _0(t^*)$ from (4.2b), for comparison.

The specific comparisons in figures 9 to 11 concern $x_1, Y, V, \theta$ and $\omega$ vs the scaled time $t^*$ for the maximal cases. Plotted in figure 10(a) in the form $M^{1/3} (Y - Y(0))$ are solutions from the full system along with the asymptote $Y_0(t^*)$ from § 4, for comparison purposes, and similarly in figures 10(b) and 11(a,b). The trends holding and the scales showing appear to be in keeping with the analytical predictions of the previous section. This is subject to the qualifications described in the previous paragraph of course and we note that the approach of $\theta$ and $\omega$ to the asymptote is relatively slow for increasing mass. Here, the analysis implies that the acceleration $Y^{\prime \prime }$ in height must be positive since it can be shown to be proportional to $Q$, while the angular acceleration $\theta ^{\prime \prime }$ is initially positive and later negative, and both analytical aspects are reflected in the full results in figures 10 and 11; likewise the analytical integral shown in figure 6(b) is proportional to the velocity $Y^\prime$ and this asymptotic trend agrees with the full trend of figure 10.