Skimming impact of a thin heavy body on a shallow liquid layer

Abstract

This study addresses the question of whether a thin, relatively heavy solid body with a smooth under-surface can skim on a shallow layer of liquid (for example water), i.e. impact on the layer and rebound from it. The body impacts obliquely onto the liquid layer with the trailing edge of the underbody making the initial contact. The wetted region then spreads along the underbody and eventually either retracts, generating a rebound, or continues to the leading edge of the body and possibly leads to the body sinking. The present inviscid study involves numerical investigations for increased mass ($M$, in scaled terms) and moment of inertia ($I$, proportional to the mass) together with an asymptotic analysis of the influential parameters and dynamics at different stages of the skimming motion. Comparisons between the asymptotic analysis and numerical results show close agreement as the body mass becomes large. A major finding is that, for a given impact angle of the underbody relative to the liquid surface, only a narrow band of initial conditions is found to allow the heavy-body skim to take place. This band includes reduced impact velocities of the body vertically and rotationally, both decreasing like $M^{-2/3}$, while the associated total time of the skim from entry to exit is found to increase like $M^{1/3}$ typically. Increased mass thereby enhances the super-elastic behaviour of the skim.

1. Introduction

The issue of whether a body of comparatively large mass and smooth under-surface (underbody) can skim across a layer of water or other liquid is considered here. There are many interesting areas of application. An industrial one concerns aircraft icing where the body or particle of concern is an ice shard or ice crystal and the water layer present is situated on the surface of a wing or fuselage or inside an engine intake (Gent, Dart & Cansdale 2000; Mason, Strapp & Chow 2006; Purvis & Smith 2016; Smith et al. 2019). The skimming or sinking of the ice particle here influences heat transfer and hence the creation of further ice with a consequent risk to aircraft safety. Another area concerns landings of airborne vehicles and other bodies on a sea or a lake where again the matter of safety is a prime issue (Von Kármán 1929; Wagner 1931). Meteor impacts can also lead to skimming. The skimming of an object across a water surface is also probably familiar to many readers from the recreational game of stone skipping or skimming where a thin stone is thrown with a significant horizontal velocity in an attempt to cause the stone to bounce off the water surface. Ideally, the stone descends into the water with a small downward velocity, impacts and then is subjected to pressures on the underbody from the water layer which produce a positive lift and push the stone back out of the water, followed by further such bounces thereafter. Here, the onset of the skimming motion is modelled and, as such, the air effect is assumed to be dynamically negligible when compared with the liquid layer's interaction with and effect on the body.

The behaviour of the combined liquid-layer–body motion is particularly sensitive when the body has relatively large mass since a ballistic body motion would be expected for a wide range of impact conditions. The ballistic effect is due to the body's large mass, multiplied by acceleration, dominating the overall interaction under many conditions of incident velocity and body attitude: we have in mind here a thin body at small incidence rather than a bluff body with high inertia or a spear entering lengthways. The ballistic effect reduces the influence of the liquid layer, causing the body's trajectory and motion to be largely unaffected by the water layer flow. As a result, the body would simply sink as a classical projectile to the bottom of the layer. Thus, this raises the question of whether a heavy body may skim and, if so, how? That is the central question here.

A balance in terms of orders of magnitude between the vertical mass acceleration of the body and the lift due to water-flow pressure suggests a critical parameter of approximately

(1.1)\begin{equation} \frac{\rho_b}{\rho_w} \left( \frac{H}{L_b}\right)^2, \end{equation}

where $\rho _b$ is the body density, $\rho _w$ is the water density, $H$ is the water depth and $L_b$ is the body length. The parameter (cf. body mass per unit horizontal distance relative to water mass) is based on the typical time scale being that of the water transit under the body and on the body mass being approximately the body density multiplied by the water volume beneath the underbody. The latter factor leads to the squared contribution displayed above. For shallow water, the depth ratio ($H/L_b$) will be somewhat less than unity and, for a ‘heavy’ body, the density ratio ($\rho _b/\rho _w$) is anticipated to be somewhere between 2 and 10 say (this range includes various metals for example). As a result, the above critical parameter is nominally large. Then again, experience indicates that the details of the initial impact conditions (angle, orientation) and the body shape can play a very significant role in determining the trajectory of the body on the water surface. This further suggests quantification of the skim evolution would be desirable.

In general the topic of skimming, particularly in the present scenario of a thin body with a sharp trailing edge on water with finite depth, has been subject to numerous investigations: for example Green (1935, 1936), Tuck & Dixon (1989), Hicks & Smith (2010), Liu & Smith (2014) and Palmer & Smith (2020). In such cases the water is assumed to detach smoothly from the body's sharp trailing edge, whereas ahead of the body a splash jet can also be emitted (see also Hewitt, Balmforth & McElwaine 2011). One modelling challenge here is to determine the contact region between the body and the water, thereby obtaining the resulting lift force on the body, the lift being mainly due to the water-flow pressure. The early mathematical treatments of a solid object slamming on an undisturbed free surface of water were pioneered in Wagner (1931), Wagner (1932) and Von Kármán (1929). There are also more recent very interesting works by Edge (1968), Watanabe (1986), Cointe & Armand (1987), Greenhow (1987, 1988), Scolan & Korobkin (2001), Korobkin (2004), Howison, Ockendon & Oliver (2004) and Khabakhpasheva & Korobkin (2013). The range of phenomena present may be complex and there are many parameters potentially involved including the body shape and the initial conditions at impact. It is notable, though, that most studies impose a prescribed body motion whereas the present work allows free movement of the body – this free movement is felt to play a crucial part as far as the question of heavy-body skimming is concerned.

The present investigation takes the Hicks & Smith (2010) interaction theory as its basis. This assumes two-dimensional motion of the thin body and the water, the latter being treated as inviscid and incompressible. The water layer is also thin. At impact the body's inclination to the water surface is taken to be small and the body velocity is taken as almost but not quite parallel to that surface. The effects of the major interaction parameters and of different body shapes of concern are to be described by means of nonlinear analysis and computation. Solutions for the resulting fluid–body interaction including the thin-layer flow dynamics are described by Hicks & Smith (2010) and more recently by Liu & Smith (2014) for reduced inclinations, by Smith & Liu (2017) for sinking bodies, by Palmer & Smith (2020) for impacts of increased-thickness bodies and by Liu & Smith (2021) for smooth-body impacts. We consider herein whether the theory predicts skimming to take place with a body of comparatively large mass and if so under what conditions.

Section 2 describes the framework for the present study and the interaction structure including splash jets at the unknown moving front of the wetted (pressured) surface on the underbody. The work considers shallow water, with the liquid layer being thin relative to the body length, and variations in the underbody shape significantly smaller than the liquid-layer thickness. (This places the fluid flow part of the current study in an intersectional regime between the Wagner (1931) theory and the Tuck & Dixon (1989) and Korobkin (2004) theory.) Section 3 then addresses computational solutions of the water–body interaction, showing results for the body motion, the water flow, the pressures and the moving front for increasingly large values of the appropriate scaled mass and moment of inertia. Here super-elastic behaviour is encountered: the body rebounds with a velocity greater than that with which it originally hit the liquid surface. Next, § 4 presents an analysis to identify the trends in the interactive properties for enhanced mass of the body. Section 5 describes comparisons between the computational and analytical results, while § 6 gives the conclusion and further comments including practical application and the overall predictions for heavy-body impacts and rebounds.

2. The skimming model

The model used here is essentially that of Hicks & Smith (2010). 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 2010.) 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 2010) 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 (2010), using much of the flow structure in Tuck & Dixon (1989), 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 (1989) and Hicks & Smith (2010) (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. (2004). 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 (1931), Tuck & Dixon (1989), Howison et al. (2004), Korobkin (2004) and Hicks & Smith (2010). Our approach is based on Tuck & Dixon (1989) as extended to unsteady interaction by Hicks & Smith (2010) 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 (2010) 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 (2010), 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 1968; Palmer & Smith 2020). Importantly, there are three geometric conditions that are obtained at leading order from the small-time analysis (Hicks & Smith 2010; Palmer & Smith 2020) 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 (2010) and Palmer & Smith (2020). 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 (2010) 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 2017) 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 2017). 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 2017). 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 2017) 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 (2020) (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 (2010) and Palmer & Smith (2020). 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 (2017) and Palmer & Smith (2020), 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. (2011). 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. 2011).

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.

6. Conclusion

The theory of oblique skimming impact on a shallow liquid layer has been used here to help understand the behaviour for a body of enlarged mass. The trailing edge of the body initially descends into the water layer with the wetted (pressured) region extending from this point, a feature that was consistent across all cases. For cases where the body did not simply remain submerged, the wetted region later receded back towards this point. If the mass is not too large then, depending on the initial conditions, the body may complete a successful skimming motion where the body is able to rebound and exit the liquid layer with the maximal leading contact point of the wetted region remaining strictly within the confines of the horizontal extent of the body. If, however, the mass is too large, the body will become flooded and continue to descend into the liquid layer. In either case, the body undergoes changes in its vertical position, vertical velocity, inclination and angular velocity throughout the motion in response to the changes in the underbody pressures. When the body completes a successful skim and exits the liquid layer into the air, the body may continue under the influence of gravity to a subsequent skimming impact with the liquid layer. The initial conditions for such a successful skim lie in a narrow band of impact velocities decreasing as the $-$2/3 power of body mass and producing a skim time scale that is enhanced as the 1/3 power of mass.

In this work two distinct phases of the skimming motion were found, i.e. the majority of the motion and the exit stage, for which the analysis showed how the dynamics of the body trajectory changed dramatically with mass. Qualitative agreement exists with the flat plate analysis of Hicks & Smith (2010) for smaller mass. For larger mass a super-elastic response due to the vertical motion and rotational dynamics of the body is seen. Indeed, as noted above, the analysis for large mass confirms the existence of a small band of velocities and initial inclinations that may lead to successful skimming motions; conditions for near maximal skimming motions are readily calculable from the reduced-order equations. (Additionally, the logarithmic effect at exit as seen in the previous literature (Hicks & Smith 2010; Liu & Smith 2014; Palmer & Smith 2020) is shown to also be relevant in the large-mass regime and is of similar form to that in Palmer & Smith (2020), see also our (3.3). An enhanced formula capturing this effect for the final stage prior to exit is provided in implicit form in figures 7 and 8.)

The asymptotic results provide physical understanding of the skimming interaction and of the solution trends with regard to the underbody pressure, rotational dynamics and vertical motion. These results help to add overall clarity to the process that is governed by the complex coupled nonlinear system (2.2) and (2.3). In general the asymptotics and reduced-order analysis hold well in each case for the majority and for the exit parts of the skimming motion. In practical terms the assumption of large values for the scaled mass $M$ leads to the ratio of the body density to the liquid density being supposed large (exactly how large depends on the overbody shape as well as the underbody shape). Stone and many metals in particular give density ratios in the range 2–10 if the liquid involved is water.

To highlight the wider scope of the results presented in this paper, figure 12 presents a comparison for increasing mass between the full solution and asymptotic solution for a flat body from (4.14a)–(4.14c) and for a convex underbody of $O(1)$ thickness where

(6.1)\begin{equation} F(x) ={-}(1 - x^2). \end{equation}

(Equations (4.6a)–(4.6c) and (4.13) are used for the curved body case (6.1).) The $O(1)$ thickness here is relative to the original scales of § 2 of course. In the curved body case the initial angle, say $\hat {\theta }_0$, between the underbody at the trailing edge and the water layer, is equal to the initial angle for the flat plate case to account for the thickness of the body near the trailing edge. Overall, the same general trends (shallower and shorter skimming motions) are seen for the curved body with increasing mass and the asymptotic solution holds reasonably well. These results also indicate that bodies with some curvature can complete successful skimming motions for a wider range of velocities, including for conditions that would lead to the sinking of equivalent flat bodies.

Figure 12. Evolution of $x_1$ vs $t^* = M^{-1/3}t$ for $M = 3I$ and varying $I = 50, 100$ (blue and red curves, respectively). Dashed curves relates to the flat underbody case ($F = 0$) and the dotted curves relate to a curved underbody case ($F = (1 - x^2)$). The asymptotic solutions are shown in each case (black curves). The initial conditions for both cases are the same and have been chosen to achieve the near maximal wetting of the underbody for the flat body case.

For subsequent work, it would be interesting to investigate fully the interplay between larger mass and increased body curvature. This would include seeking to understand how these two properties affect the presence of different stages in the skimming motion and affect the overall skimming dynamics. Moreover, introducing further flexibility into the modelling set-up would be of interest, in particular by allowing the trailing edge to move freely, or by changing the underlying modelling assumptions by introducing a new scaling that voids (2.1a,b). In particular, a moving trailing edge, as noted in Hicks & Smith (2010), would allow for solutions with splash jets both in front and behind the body, as in Howison et al. (2004) and Liu & Smith (2021). Given the significant splash jets present ahead of the body in the work presented here, it would be very interesting to see how these develop at the trailing edge of the body also for increasing mass.

Further extensions to this work include the incorporation of gravity effects for larger bodies or the consideration of a series of impacts for variously shaped bodies, extending the work presented in Liu & Smith (2014). Exploring fully the relation between the present shallow-water regime and the deeper-water models of Wagner (1931), Korobkin (2004) and Howison et al. (2004), is also of interest, especially the properties of fluid–body interaction for the regime where the body length and water depth are comparable. Perhaps the greater challenge, however, is that associated with having the streamwise variation in underbody height become comparable to the incident liquid-layer thickness (see figure 1 and the description in Appendix A), where the fluid flow part of the interaction becomes nonlinear.