## Abstract

The evolution of the combined solid–fluid motion when a solid body undergoes a skimming impact with (and rebounds from) a shallow liquid layer is investigated. A model is derived coupling the motion of the body to the fluid dynamics of the liquid layer. This predicts that the lift on the body induced by the pressure in the liquid layer is sufficient to entirely retard its incident downward motion before causing the body to rise out of the liquid. Water exit is predicted at a finite scaled time. Analysis for the small-time behaviour immediately after touchdown, and also as water exit is approached, shows close agreement with numerical prediction.

## 1. Introduction

The childhood pastime of stone skipping or *ducks and drakes* (as it is known in the UK) involves throwing a smooth stone almost horizontally in an attempt to produce a series of oblique skimming impacts with water, each separated by a phase in which the stone rebounds into the air. Experimental studies of stone skipping phenomena have been conducted by Clanet *et al.* (2004) and Rosellini *et al.* (2005) in particular. At each impact with the water, the pressure on the stone produces a positive lift overall on the stone, carrying it out of the water. If the effect of air resistance is neglected, then the stone is carried along the classical parabolic trajectory in the air under the influence of gravity to another skimming impact with the liquid.

In the present work the original motivation arose in aircraft-icing applications, in which ice crystals (often of a long thin rod-like shape) impact the leading edge of an aerofoil that has become coated in a thin liquid layer as the result of flying through clouds; the interactions between a solid body and a shallow liquid layer are examined here during one of these skimming impacts. In the current model, the vertical height and the rotation of the body are free to respond to liquid pressure-induced lift and moment of inertia, respectively. The horizontal retardation of the body owing to the liquid drag is shown to be a secondary effect. Using *inter alia* a flow structure based on that of Smith & Ellis (2010) and the momentum jump conditions of Tuck & Dixon (1989) a series of nonlinear ordinary differential equations are derived governing the vertical height of the body, its rotation and the position of the leading edge of the interactions between the liquid and solid body. The inherent nonlinearity in the model allows the subsequent rebound phase of the impact to be captured.

The steady-flow study by Tuck & Dixon (1989) produced a relationship for surf-skimming interactions based upon momentum conservation between the leading-edge height, the trailing-edge height and the ambient liquid layer depth. This was confirmed by studying the inner flow close to the leading edge of the wetted region between the body and the liquid. A momentum jump at the leading edge of the body is also present in the multi-body interaction studies of Bowles & Smith (2000) and Smith & Ellis (2010). An important study by Howison *et al.* (2004) incorporates a tangential velocity into the slamming theory for the deep water impacts of Wagner (1932) and Howison *et al.* (1991, 2002), and the shallow water impacts of Korobkin (1995, 1999). Mention should also be made here of related studies by Oliver (2002), Hewitt *et al.* (submitted) and Ellis *et al.* (submitted). Throughout the study of Howison *et al.* (2004) the vertical velocity of the solid body is assumed to remain constant even after the body enters the liquid and so the rebound phase of the skimmer motion is excluded. The asymptotic theory for the slamming problem has been recently extended to second-order accuracy (Oliver 2007), while boundary integral solutions for the pressure, free-surface profiles and load on the plate for this problem have been presented by Iafrati & Korobkin (2004, 2008). High-speed droplet impacts with water have been considered by Purvis & Smith (2005) and Howison *et al.* (2005), with the former study also incorporating impacts at incidence. At slower speeds the flow under a fixed surf board or an inclined sluice gate is considered by Binder & Vanden-Broeck (2005, 2007), in a regime in which surface tension has a significant effect on the free-surface evolution. Interactions between the liquid layer and body have been incorporated in a number of previous impact studies, including in the form of an elastic response in the body (Korobkin & Khabakhpasheva 1999) and a trapped gas cavity in the region of impact (Korobkin 1996). The latter study is notable in that a system of ordinary differential equations governing the boundaries of the contact patch between fluid and solid and the response of the body to the impact is again derived. Also of note in this area are the recent studies by Pak & Lauga (2010) and Mandre & Mahadevan (2010).

A model for the skimming impact of a solid body and a shallow liquid layer is derived in §2; the general theory here excludes deep water. Section 3 focuses this general theory, initially to the case of a parabolic body geometry and latterly to the case of a flat plate. Section 4 describes numerical results and a solution technique for the flat plate case. Also considered are the behaviour in the wake and at small times immediately after impact. The behaviour of the body as water exit is approached is considered in §5. Finally, further discussion and future extensions of the work are presented in §6.

## 2. Model development

The oblique impact of an ice crystal hitting a liquid coating on an aircraft surface motivates the investigation of the motion of a solid body as it undergoes a skimming impact with an otherwise still liquid layer of finite depth and interacts with the liquid flow (see figure 1). A solid body with mass kg is considered as it initially impacts upon a horizontal liquid layer. Throughout the motion the drag and lift on the solid body from the air are neglected; so before the body enters the liquid, gravity is the only force acting upon the body. The fluid–solid interactions in this problem are complicated and, therefore, our initial focus considers only the problem in two dimensions. However, there is no *a priori* impediment in principle to extending the subsequent analysis to three dimensions. The motion is studied in the Cartesian plane, where is parallel to the initial liquid layer surface and points vertically upwards; the bars, where used, denote dimensional quantities. In this coordinate frame, at the point of initial touchdown the body is assumed to be moving with velocity m s^{−1}, where for skimming impact we expect and .

As the body is carried into the liquid by its initial downward momentum, extra hydrodynamic forces start to act upon the body in addition to gravity. These are the drag slowing the body in the horizontal direction and the lift acting in the vertical direction, which both arise as integrals of the fluid-dynamical pressure over the wetted body surface.

### (a) Model development and non-dimensionalization

The motion is considered in a frame of reference centred on the horizontal coordinate of the centre of mass of the body. The vertical position of the centre of mass is denoted . We choose a length scale based on the distance from the centre of the body to its trailing edge, so that . For a body with uniform density, this means the total dimensional distance from the leading edge to the trailing edge is roughly . The velocity is non-dimensionalized using the horizontal speed , while the convective timescale and pressure scale are also used. Here is the uniform liquid density, while is the liquid viscosity. When applied to the two-dimensional Navier–Stokes equations, the non-dimensional velocity ** u**=(

*u*,

*v*) and pressure

*p*are related through 2.1where

*t*is the non-dimensional time. Here the Reynolds number is , and the Froude number is , where

*g*is the acceleration owing to gravity. Their inverses measure the significance of viscous and gravitational effects, respectively. If the body is smooth, for example with length m and incident speed m s

^{−1}, then

*Re*=190 000 and

*Fr*=51, indicating that effects owing to viscosity and gravity do not affect the liquid flow.

The interaction between the liquid layer and the solid body occurs over the region where the two are in contact. For instance, the conservation of momentum in the horizontal and vertical directions implies2.2aand2.2brespectively, where the integration contour follows the wetted portion of the body. These equations balance the component of the hydrodynamic forces on the body (the integral of the pressure over the wetted region of the body) against the product of the mass and the acceleration of the body. The quantity is the non-dimensional mass. We note that in two dimensions the liquid density has dimensions of kg m^{−2}. Similarly, conservation of angular momentum implies2.3where again the integration is carried out around the wetted portion of the body. The non-dimensional moment of inertia of the body, *i*, is related to the dimensional moment of inertia .

Next, we assume that the solid body lands on a shallow liquid layer, of depth , where and the depth of the liquid layer is assumed to be small compared with the horizontal extent of the solid body (see figure 1). In this case the horizontal distance *x* and time *t* are of *O*(1) and the vertical height scales as *y*=*ϵY* . In the wetted region *x*_{1}≤*x*≤*x*_{0} say, the liquid layer thickness *h*(*x*,*t*) depends on the position *Y* =*Y*_{m}(*t*) of the centre of mass of the body, the angle of rotation *θ*(*t*) of the body and the thickness of the body, which is 2*T*(*x*) say. Thus,2.4Here *Y*_{m} and *h* are measured from the bed where *Y* =0. In front of the body (*x*<*x*_{1}) and in the wake (*x*>*x*_{0}), the function *Y* =*h*(*x*,*t*) denotes the free surface of the liquid layer. The quantities *x*_{0}, *x*_{1}, *Y*_{m} and *θ* are unknown functions of time *t*. A third non-dimensional group can be constructed: the Weber number , whose inverse measures the importance of surface tension on the free surface. Here *σ* is the surface tension coefficient and for the above parameters *We*=170 000, indicating that surface tension nominally has negligible effect on the evolution of the liquid free surface.

The large difference between the horizontal and vertical speeds is exploited by seeking solutions with the form2.5where *x*, *Y* , *u* and *V* are *O*(1). If *ϵ* is small, but also the Froude number *Fr* and the reduced Reynolds number *ϵ*^{2}*Re* are large, then the effects of gravity and viscosity are negligible in the problem. In this case, the vertical momentum equation within equation (2.1) immediately implies *p*=*p*(*x*,*t*). The flow is also irrotational, as a consequence of which the scaling implies ∂*u*/∂*Y* =0, and the horizontal momentum conservation equation in (2.1) simplifies to give2.6Additionally, this means *u*=*u*(*x*,*t*). Therefore, if the mass conservation equation is integrated from the base at *Y* =0 (where a no-penetration boundary condition is applied) to the body (where a kinematic boundary condition is applied), then2.7

A key feature of the model is the pressure jump at *x*=*x*_{1}+, the leading edge of the wetted region (see figure 1). Here, following Tuck & Dixon (1989), the jump condition2.8with2.9applies where *h*(*x*_{1},*t*) is the position of the body at *x*=*x*_{1}. Similar pressure jump conditions are employed in other problems at the leading edges of solid bodies interacting with steady or unsteady fluid flows (Bowles & Smith 2000; Smith & Ellis 2010).

Concerning the trailing edge, if we subtract off the atmospheric pressure from the pressure in the liquid and assume the pressure varies continuously at the trailing edge *x*=*x*_{0}, then *p*(*x*_{0},*t*)=0. Similarly, the free surface *h*(*x*,*t*) is continuous across *x*=*x*_{0}. So in the wake region (*x*>*x*_{0}),2.10

The body is in contact with the liquid layer only in the region *x*_{1}≤*x*≤*x*_{0} (see figure 1). Outside this region the pressure on the body is zero as in equation (2.10). Therefore, when the variables are scaled to account for the shallow liquid layer the vertical and angular momentum equations (2.2b) and (2.3) imply2.11aand2.11brespectively, where is the scaled mass of the body and is the scaled moment of inertia. The horizontal momentum conservation equation (2.2a) is unaltered. In the case in which the body responds to interactions with the liquid layer both *M* and *I* are generally *O*(1). We subsequently neglect the effect of gravity in the vertical momentum conservation equation, which means we must consider the case in which *ϵFr* is large.

The above system governs the behaviour of a solid body during a skimming impact with a shallow liquid layer. The main unknowns in this system are *h*, *u*, *p*, *Y*_{m}, *θ*, *x*_{1} and *x*_{0}, which are expected to be of order unity as the result of the above scaling. Our aim is to investigate the behaviour of a solid with horizontal speed 1, which is initially falling (d*Y*_{m}/d*t*(0)<0) and touches the liquid layer surface at time *t*=0, with *Y*_{m}(0) such that the body just touches the surface at one point where *x*_{1}=*x*_{0}. We expect to be given the entry conditions for the body involving its angle of rotation *θ*(0) (which can be either positive or negative) and its angular velocity d*θ*/d*t*(0) (which again can be either positive or negative).

### (b) Dynamics for a thin body

If the vertical extent of the body is much less than that of the liquid layer, such that the body is relatively thin, then2.12where the terms with tildes are small. The trailing edge of the body is taken to be sharp rather than smooth and the liquid free surface detaches from the body there (see figure 1). Thus the distance between the trailing edge and the centre of mass *x*_{m} is fixed. Also the general leading-edge position *x*_{1} is of *O*(1). The body shape equation (2.4) implies2.13

Substitution of equation (2.12) into equations (2.7) and (2.6) in the region *x*_{1}≤*x*≤*x*_{0} gives2.14aand2.14bAt the known trailing edge *x*_{0}, the equi-pressure condition implies , while in the wake region *x*>*x*_{0} we find2.15At the unknown leading edge *x*_{1} the jump condition (2.8) yields2.16with2.17Also the vertical momentum equation (2.11a) simplifies to2.18

The integral term in equation (2.2a) is small for small and therefore the body moves with a constant horizontal velocity to leading order. Finally, equation (2.11b) simplifies to become2.19

## 3. Dynamics for a body with parabolic thickness

The liquid in the wetted region [*x*_{1},*x*_{0}] is of primary interest. If the thickness there is at most quadratic in *x*, say with constants *A*, *B* and *C*, then from equation (2.13)3.1is also quadratic. In this region (2.14) indicates that is at most linear in *x*. Therefore, we can write3.2where *γ*_{n} for *n*=0, 1, 2 and 3 are unknown functions of time, which can be determined by substituting the definition of pressure into the remaining governing equations as follows.

The jump condition (2.16) with equation (2.17) implies3.3

Next, integration of equation (2.14b) with equation (3.1) gives3.4where *D*(*t*) collects together terms which are just functions of time or are constant, and represents the fluid velocity at *x*=0. If this velocity profile is evaluated at *x*=*x*_{1}, and substituted into the jump condition (2.16), then a second equation for leading-edge pressure is obtained, with the form3.5

Next the pressure (3.2) and velocity profiles (3.4) are substituted into the liquid momentum conservation equation (2.14a). Matching the coefficients of the terms which are constant, and linear and quadratic in *x*, gives three expressions for *γ*_{1}, *γ*_{2} and *γ*_{3}, with the forms
3.6a
3.6b
and
3.6crespectively. At the trailing edge the equi-pressure (Kutta) condition implies3.7which allows us to find *γ*_{0} using (3.6*a*–*c*). In the wake region *x*>*x*_{0}, the free-surface height and liquid velocity are governed by3.8

Finally, momentum conservation on the body requires3.9aand3.9b

### (a) The straight flat body

For the remainder of the paper, we concentrate on a flat plate occupying [−*x*_{0},*x*_{0}]. In this case, the parameters representing the shape of the body *A*=*B*=*C*=0, and the horizontal position of the centre of mass of a uniform body *x*_{m}=0 also. The equations governing the body motion simplify significantly. The equation for the trailing-edge pressure (3.7) and equation (3.6c) for *γ*_{3} are unchanged, but are included here for completeness. At this stage, we also drop the tilde indicating small aspect ratio body.

The trailing-edge pressure condition (3.7) implies3.10a

The two expressions resulting from the jump condition at the leading edge imply3.10bfrom equation (3.3), and3.10cfrom equation (3.5).

The equations (3.6) for *γ*_{1}, *γ*_{2} and *γ*_{3} imply
3.10d
3.10e
and
3.10f

The linear and angular momentum conservation equations (3.9*a*–*b*) become3.10gand3.10hrespectively.

Equations (3.10*a*–*h*) constitute a system of eight equations for the quantities *Y* , *θ*, *x*_{1}, *D* and *γ*_{n} for *n*=0, 1, 2 and 3, which are all unknown functions of time. The explicit dependence on *x* is eliminated from this system. The remaining unknown functions of time *t* evolve as the solid body skims along the liquid surface.

Before analysing the behaviour of this system of equations, it is worth comparing the above modelling approach with the interesting ad hoc models of Clanet *et al.* (2004) and Rosellini *et al.* (2005), which describe stone skipping in deep water. These earlier models incorporate a lift force on the body which is proportional to a dynamic pressure response , multiplied by the wetted area of the body. This is essentially equivalent to equation (2.2), but under the assumption that the pressure beneath the body is uniform. This approach allows the body to react and rebound from the liquid. However, the assumed lack of spatial variations in the liquid pressure prevents changes in the angle of attack over the course of an impact, which result from changes in angular momentum. This effect is included in our current model rationally (see, for instance, equation (3.10h)). Similar views hold for the interesting study of Hewitt *et al.* (submitted), kindly pointed out to us by Dr James Oliver after submission of the present work and having some overlap with ours but differences also, for example with regard to a full account of water exit as in §5 below.

## 4. Skimming impact results

Solutions were obtained by two methods and their results agreed. One method is now described. The system of equations (3.10) contains second-order derivatives of both vertical position *Y* and the angle of the body *θ*. In many problems it is convenient to eliminate these second derivatives by introducing new variables and writing the system purely in terms of first derivatives. In this case, we introduce the vertical velocity *V* and angular velocity *ω* satisfying4.1aand4.1bAnd so the momentum equations take the form4.1cand4.1da pair of first-order ordinary differential equations. Similarly, substituting for *V* in equation (3.10b) gives4.1eNext we use equation (3.10e) to find a first-order ordinary differential equation for *x*_{1}, with the form4.1fThis gives a set of six ordinary differential equations for the quantities *Y* , *V* , *θ*, *ω*, *D* and *x*_{1}. However, the values of the *γ* variables remain to be found.

We now seek an algebraic system which we can use to find the values of the *γ* variables. Firstly, the wake condition (3.10c) gives one relationship between the *γ* variables. Secondly, starting with equation (3.10h), we can use equation (3.10a) to eliminate d^{2}*θ*/d*t*^{2}, leaving4.2aThirdly, we can substitute from equations (3.10g) and (4.1b) into equation (3.10d) to generate another relationship between the *γ* variables, with the form4.2bFinally, we can eliminate d*x*_{1}/d*t* from equations (3.10e) and (3.10f), and, after some manipulation, we find4.2cThe equations (3.10c) and (4.2*a*–*c*) form a linear system of equations to find the values of the *γ* variables providing the values of *Y* , *V* , *θ*, *ω*, *D* and *x*_{1} are known. Together this system of first-order ordinary differential equations and linear algebraic equations for the *γ* variables leads to a numerical scheme to find the evolution of the skimmer. Further manipulation leads to the complete elimination of the *γ* variables from the system, but this generates much more complicated differential equations and does not provide any additional insight into the behaviour of the system.

The system is solved using a seventh and eighth-order Runge–Kutta adaptive step-size integration starting from the small-time solution determined in §4*a* below. Numerical solutions are presented in figure 2*a* for a case in which *M*=2, *I*=1, *V*_{0}=−1 and *θ*_{0}=−2. The vertical height *Y* of the centre of mass of the plate is shown by the thick solid line. The height of the centre of mass initially falls, as the downward momentum of the plate carries it into the liquid. Subsequently, the plate begins to feel a positive lift from the increased liquid pressure in the wetted region, which acts to slow the descent of the plate and from *t*=2 onwards the plate starts to rise out of the water. The length of the wetted region tends to zero at *t*=5.385, indicating that the plate is once again clear of the water or, strictly, approaching such a state. Also shown in figure 2*a* are the values of *θ*, *x*_{1} and *D*. The position of the leading edge *x*_{1} starts close to the trailing edge *x*_{0} as the plate initially touches down at just one point. The position of the leading edge then moves in the negative direction as more and more of the plate becomes submerged, before it reaches a maximum. After this point the position of the leading edge tends back towards the fixed trailing-edge position again as the plate exits the liquid. While the leading-edge position *x*_{1} is to the right of the centre of mass, there is an anticlockwise moment on the plate owing to the pressure being positive only over the wetted region. This acts to reduce the angle *θ* of the plate. The corresponding values of the *γ* variables are shown in figure 2*b*.

The evolution of the leading-edge position *x*_{1} throughout an impact is in figure 3. The basic configuration in which *M*=2 and *I*=1 with an initial impact velocity *V* =−1 and plate rotation *θ*=−2 is shown by the thick solid line in each part of the figure. Also shown in figure 3*a* is the leading-edge evolution for an impact velocity *V* =−2 and −3, while the other parameters are unchanged. These results show that if a plate touches down with a greater initial impact velocity, then it penetrates deeper into the water, which makes good physical sense. However, the case *V* =−3 illustrates that there is an upper limit on the impact velocity with regard to rebounding, as in this case the leading-edge position *x*_{1} enters the range *x*<−1 for a portion of the evolution. This corresponds to a leading-edge position off the plate and is therefore unphysical. Instead the front of the plate becomes submerged and the theory breaks down.

Figure 3*b* shows the leading-edge evolution for *θ*=−1 and −2 and figure 3*c* shows the evolution for *M*=1 and 3, in addition to the basic configuration. As the initial plate rotation is decreased the penetration depth required to move the leading-edge position towards the front of the plate falls and again the front of the plate becomes submerged in the *θ*=−1 case when *x*_{1}=−1. With increases in plate mass there is a reduction in the deceleration on the plate upon entering the water, allowing deeper penetration. Again, this eventually causes the front of the plate to become submerged and the theory to break down.

The pressures generated in the interaction region are shown in figure 4. Figure 4*a* shows pressure profiles along the wetted region of the body at time increments of 0.5, starting at the point of initial touchdown. In all profiles the pressure tends to zero at the trailing edge of the wetted region, satisfying the Kutta condition. For the standard impact configuration, the pressures along the body are largest at the leading edge of the interactions between the body and the liquid layer. The pressures on the body grow as the leading-edge position moves away from the trailing-edge position. However, the maximum pressure on the plate occurs before the leading-edge position reaches its maximum separation from the trailing edge, and subsequently the pressures along the body start to fall. Later in the impact, the pressures along the body are much smaller than just after entry, owing to the change of direction of the body. Indeed the pressures along the body appear to tend to zero as water exit is approached. The leading-edge pressure , shown in figure 4*b* as a function of time over the course of the impact, increases rapidly from zero as the body enters the liquid layer, with the maximum pressure on the body being found at the leading edge. As water exit is approached, the leading-edge pressure tends to zero. The inset in figure 4*b* illustrates the very final moments before water exit and shows that there is a sudden reduction in the leading-edge pressure at this stage. This is considered further in §5.

The position of the contact region between the body and the liquid layer is given by equation (3.1) and the velocity is given by equation (3.4), in the region *x*_{1}≤*x*≤*x*_{0}. The height of the plate is shown in the left-hand pane of figure 5*a* at time increments of 0.5, starting at the point of initial impact. This shows the plate initially descending into the water layer before the overall positive pressure subsequently causes the plate to rise out of the water again. Velocity profiles at the same time increments are shown in the left-hand pane of figure 5*b*.

### (a) Small-time solutions

The configuration of the plate upon impact is given by the centre of mass *Y*_{0}, the vertical velocity *V*_{0}, the angle of the plate *θ*_{0}, the horizontal position of the trailing edge *x*_{0} and *ω*_{0} the angular velocity of the plate. Given this initial configuration, analytical solutions for the plate evolution at small times after entry are sought, using asymptotic expansions of the form
4.3a
4.3b
4.3c
and
4.3d
For small times after touchdown, the leading-edge contact point will move in the negative *x* direction as more of the plate becomes submerged, and therefore we expect . The presence of the initial plate vertical velocity *V*_{0} and angular velocity *ω*_{0} as coefficients of the linear terms in *Y* and *θ*, respectively, is an immediate consequence of the Taylor expansion of these quantities at small time. In conjunction with these expansions we need additional expansions for each of the *γ* variables and these take the form
4.4a
4.4b
4.4c
and
4.4d
The goal of this section is to determine additional information about the unknown coefficients in these expansions which are not known from the initial configuration, such as and *D*_{1}.

If we substitute these expansions into the wake equation (3.10a), then we find4.5aand4.5bgiving a relationship between the coefficients of the *γ* expansions at each order.

The leading-edge pressure term . When the conditions (4.5) from the wake equation are used together with the above expansions (4.3) and (4.4), we find4.6This indicates that the pressure at the leading edge is small when the plate enters the liquid for the first time. Therefore, when we substitute into equation (3.10b), we findThe leading-edge pressure on the left-hand side is *O*(*t*), and this must be the same size as the term on the right-hand side. This means *Y*_{0}+*x*_{0}*θ*_{0}=0, which is a consequence of the geometrical configuration of a plate entering an initially undisturbed liquid layer. Therefore, the remaining terms are related through4.7Similarly, when we substitute into equation (3.10c)This time the requirement for the right-hand side to be at most *O*(*t*) requires4.8Therefore, when we match terms at *O*(*t*), we find4.9

Equation (3.10d) impliesMatching terms at leading order implies4.10

The remaining equations guarantee the absence of terms in *Y* and *θ* of terms with size *O*(*t*^{2}) and *O*(*t*^{3}).

Equation (3.10e) implies *γ*_{2(0)}=*ω*_{0}, and therefore *γ*_{2(0)} can be eliminated from the above expressions. Equation (4.8) determines *D*_{0}, while equations (4.7) and (4.9)–(4.10) can be thought of as a system of three equations to determine the value of *D*_{1}, *γ*_{1(0)}, and , when we are given the initial plate height *Y*_{0}, vertical velocity *V*_{0}, angle *θ*_{0} and angular velocity *ω*_{0}.

For instance is given by4.11where *x*_{0} is the position of the trailing edge of the plate. Having found we can then find *D*_{1} and *γ*_{1(0)}, giving additional terms in the expansions.

The small-time asymptotes for the leading-edge position are shown in figure 3 in addition to the numerical solutions. In figure 3*a*,*b*, the linear asymptotes vary with changes in impact velocity and initial angle, respectively, through equation (4.11). In each case, the linear asymptotic theory is in excellent agreement with the numerical profiles. For the variationsQ1 in plate mass shown in figure 3*c*, both the linear asymptotes and the numerical results are in excellent agreement for small times and both initially appear independent of the plate mass. This can be seen with the lack of dependence upon the mass in equation (4.11).

### (b) Wake behaviour

The height of the free surface and the fluid velocity perturbation in the wake are governed by equation (3.8). Using the method of characteristics we find4.12a function of *x*−*t* determined by *U*(*x*_{0}−*t*)=*u*(*x*_{0},*t*), the velocity profile (3.4) evaluated at the trailing edge. Similarly4.13with *H* determined by . Here, *h*=*h*(*x*_{0},*t*) is the free-surface height evaluated just downstream of the trailing edge.

The free-surface height and velocity perturbation profiles are shown in the right-hand panes of figure 5*a*,*b*, respectively. As required, the velocity profiles are continuous across the trailing edge of the plate. At the rear of the wake the liquid velocity profile tends to the undisturbed velocity. Immediately prior to the end of the wake the velocity perturbation is positive, causing a build-up of liquid and a rise in the free-surface height locally at this point.

The stone-skipping experiments of Clanet *et al.* (2004) and Rosellini *et al.* (2005) are conducted in deep water, which makes direct quantitative comparison with the current study of shallow water skimming difficult. However, beyond the obvious similarity that both the experiments and the current theory show the body rebounding from the liquid, the experimental photographs through the wake have many features in common with the free-surface profiles shown in figure 5*a*. In both cases, after an initial impact phase, the global free-surface minimum detaches from the trailing edge of the body and moves downstream into the wake. Further downstream in the wake from this minimum, the theory predicts and experiments show that the free-surface profile recovers rapidly to a height above the undisturbed free surface. The jump in free-surface height at the rear of the wake in the profiles in figure 5*a* can then be interpreted as the splash present at the rear of the wake in the experiments.

## 5. Finite-time exit from the liquid

The numerical results suggest the plate rebounds after its initial touchdown and finally exits or, rather, approaches exiting the liquid completely after a finite time. To analyse the behaviour close to the exit time *t*_{e}, we investigate expansions of the key quantities. In contrast to the small-time behaviour, the numerical results suggest that both *γ*_{0} and *γ*_{1} become unbounded as (see figure 2*b*). Therefore, we propose expansions of these quantities with the form5.1aand5.1bwhere the coefficients of the expansions are to be determined and the exact functional form of is unknown at this stage, beyond the fact that we expect singular behaviour as . The numerical results suggest that the values of *γ*_{2} and *γ*_{3} remain finite as and therefore we consider expansions5.1cand5.1din this limit. The numerics suggest that the leading edge of the plate moves towards the trailing edge of the plate with a finite velocity as , which motivates an initial expansion for the leading-edge position in which the largest non-constant term is linear in (*t*_{e}−*t*), i.e.5.2where *x*_{0} is the known fixed position of the trailing edge.

The choice of the same singular function in both the expansions of *γ*_{0} and *γ*_{1} is motivated by considering the wake condition (3.10a). As the trailing-edge pressure is zero throughout the motion, the leading-order terms imply the same singular function is present in both the expansions of *γ*_{0} and *γ*_{1}, and that5.3Similarly, matching the coefficients of the constant terms implies5.4

Unlike the small-time behaviour, as the plate approaches water exit the configuration of the plate is unknown and the motion of the plate as it leaves the water is of interest. In order to match the asymptotic behaviour close to water exit, we assume the vertical position and velocity, and the angular rotation and velocity at the exact moment of water exit are given by the final value computed by the numerics. As these quantities all remain bounded at large time, regular expansions with the form5.5aand5.5bare proposed for the vertical height and plate angle, with the final plate height (velocity) on water exit being denoted and the final plate angular rotation (velocity) on water exit being denoted . The coefficients of the linear terms in each expansion are immediate from the Taylor expansion of each quantity about the exit time.

Next we consider the evolution of *D* based on equation (3.10d). As the quantity *γ*_{1} is unbounded as , we expect a second term in this equation also to become unbounded. The numerics suggest that d*D*/d*t* is unbounded for large time, which is consistent with the expansions (5.5) for *Y* and *θ*. Matching terms at leading order suggests an expansion for *D* with the form5.6The quantity *D* remains bounded for large time, so we expect to be some function which is between constant and linear in size. This implies that, while is singular, it is smaller than *O*((*t*_{e}−*t*)^{−1}).

It is next expedient to consider the leading-edge pressure , as water exit is approached. The numerical results suggest that the leading-edge pressure is small as water exit is approached. Applying the wake conditions (5.3) and (5.4) to the leading-edge pressure condition shows why this is the case as the largest term remaining in the expansion has the form5.7where the first term is some function between constant and linear in size.

Having determined the largest term in the pressure expansion, we next consider which terms in equations (3.10b) and (3.10c) have equal size. At small time in equation (3.10b), we found geometrical considerations and the undisturbed liquid free surface implied that *Y* +*x*_{1}*θ* was small. However, at large time this quantity is *O*(1) and we are forced to conclude that 1−d*x*_{1}/d*t* is therefore small. This requires that in equation (5.2) and we need to consider higher order terms in the expansion of *x*_{1}. If we consider an expansion with the form5.8where *g* is some function between linear and quadratic. A non-regular expansion is required to balance the leading-order pressure (5.7). Taken together, this implies the leading-order terms in equation (3.10b) are given by5.9

Next equating the right-hand sides of equations (3.10b) and (3.10c) implies5.10which implies the term in square brackets is small at large time. This gives a relationship for *D*_{0} with the form5.11To leading order this implies5.12

From equations (5.9) and (5.12), *f* satisfies the ordinary differential equation5.13Therefore,5.14for some constant of integration *α*. The result (5.14) is as found in a different context by Smith & Ellis (2010). Further integration can be conducted to find *f* itself, which can be written in terms of elliptical integrals. However, this is ultimately not instructive. The numerical predictions for the leading-edge position *x*_{1} are compared with the linear asymptotic theory for finite-time exit from the liquid in figure 3. The asymptotes provide good agreement with the numerical solutions, particularly for lower impact velocities, initial plate angles and plate masses. The agreement is not quite as close as that predicted in the small-time case, but for the finite-time exit problem there is a clearly delicate additional term (see equation (5.14)) in the expansion between the linear and quadratic terms which is not present at small time.

Having determined the form of , the leading-order pressure behaves as5.15as . This behaviour close to water exit can be seen in the leading-edge pressure evolution shown in the inset of figure 4*b*, where the pressure profile tends to zero quite suddenly as the plate leaves the liquid layer. The value of the constant of integration *α* can be estimated graphically by considering the function5.16as . For the standard configuration of plate impact described above, a value for *α*=6.0×10^{−6} is estimated, indicating very rapid changes in the evolution very close to water exit.

## 6. Conclusions

A consistent model has been presented to explain the interactions between a rigid plate and a shallow liquid layer, as the plate undergoes an oblique skimming impact with the water layer and then rebounds from it. The plate is subject to changes in vertical motion and rotation in response to the pressure in the liquid layer, in addition to the major horizontal motion. The induced motion causes the plate to rebound from the liquid into the air within finite time, providing the leading edge of the plate/liquid interaction region is confined to the length of the plate. Physically the plate rebounds from the water because the pressure induced in the liquid layer is broadly positive, which causes deceleration of the downward motion of the body and then upward acceleration, ultimately leading to water exit. After rebounding back into the air, the body may be carried to a subsequent skimming impact with the liquid layer under the influence of gravity. Asymptotic behaviours for small time and close to water exit have been determined and show close agreement with numerical results.

The model in §3 allows impacts to be considered in cases in which the body has some prescribed thickness, rather than just being a plate. Additionally, if rather than having a sharp end point the body has a smoothly varying body shape, then the trailing-edge position *x*_{0} can move. Subsequent work could allow for the free movement of the trailing edge, and in this case the model would be extended to incorporate an additional momentum jump condition at the trailing edge. This allows solutions with splash jets both in front and behind the body, as in Howison *et al.* (2004). Other factors have been neglected or shown to be small in creating this simplified model. For heavier bodies, gravity must be included if the mass of the skimmer . Drag on the body from the water is shown to be negligible in the problem described. However, over a series of impacts a skimming stone will eventually lose horizontal momentum and fail to rebound from the liquid as higher order effects dissipate energy. Interactions with the air may play a role in some circumstances with drag in the air slowing the horizontal motion. It is also known that in a narrow gap between a solid and a liquid free surface significant additional pressure variations may be present in the gas as impact is approached owing to air-cushioning. In normal impacts this can delay the moment of impact (Smith *et al.* 2003), while in oblique impacts the gas pressure can accelerate the liquid free surface towards the body (Smith *et al.* 2006; Hicks & Purvis 2010). Finally, the current study is primarily focused on skimming impacts with a shallow liquid layer, which arises for instance in the industrial context of ice crystal impingement onto a liquid-coated aircraft surface. However, the extension of this model to capture solid–fluid interactions in oblique impacts with deep water remains an open problem.

## Acknowledgements

The authors would like to thank Roger Gent (AeroTex UK) and Andrew Ellis for their helpful early discussions. Thanks are also due to Alexander Korobkin and James Oliver for their insights, to colleagues in the icing and impact group, and to EPSRC (RA) and UCL (RAIS) for their support to P.D.H.

- Received June 14, 2010.
- Accepted July 22, 2010.

- © 2010 The Royal Society