Elastic sphere skipping phenomena

Prior research has reported an upper bound on the impact angle of below which rigid spheres (density ρ s ) will skip on water (density ρ w )7,8,10,21. We have found that elastic spheres skip at much larger values of β o , raising the question of how the elastic response enables this enhanced skipping behaviour. To investigate the mechanisms underlying elastic sphere skipping, we film the water impact of custom-made elastomeric spheres with a high-speed camera viewing from the side. Rigid and elastic spheres having nearly the same radius R and density ρ s are shown experimentally impacting the water in Fig. 1a,b. Each sphere strikes with approximately the same speed U o and impact angle β o , but the elastic sphere has a shear modulus G that is four orders of magnitude smaller than that of the rigid sphere material. Within a few milliseconds after impact, the elastic sphere deforms dramatically and rides along the front of a cavity on the air–water interface before lifting off the surface. By the time the elastic sphere is two diameters above the surface (t≈25 ms), the rigid sphere has plunged nearly the same distance below it. The elastic sphere evidently experiences a larger upward vertical force from the water.

Figure 1: Elasticity alters sphere skipping dynamics. (a) High-speed images of an oblique impact show a rigid sphere carving an air cavity into the water as it dives below the surface (U o =24.3 m s−1; β o =29.6°; G=5.66 × 105 kPa; R=25.8 mm; and ρ*=ρ s /ρ w =0.959). Scale bar, 40 mm. (b) A highly compliant elastic sphere significantly deforms upon impact and skips off the surface (U o =22.0 m s−1; β o =32.0°; G=12.3 kPa; R=26.2 mm; and ρ*=0.937). (c) The deformed sphere resembles a disk-shaped stone oriented at a dynamic attack angle α. An inertia-dominated hydrodynamic force F acts on the flattened face, which moves with velocity described by U B and β B . Scale bar, 40 mm. (d) A surprising consequence of the interaction of sphere vibrations with the liquid interface is the formation of nested air cavities (that is, a matryoshka cavity). (e) A rigid sphere can skip if the impact angle does not exceed , leaving a smooth, asymmetric cavity on the surface (β o =17.3° for sphere shown). Full size image

The extreme sphere deformation is more evident in Fig. 1c, which shows that the water-contacting surface assumes the shape of a disk with a larger radius than that of the undeformed sphere. The disk is oriented at an attack angle α that, unlike for skipping stones4, changes in time throughout the impact. Large amplitude oscillations excited by the impact can persist in the sphere after lifting off the surface (Fig. 1b and Supplementary Movie 1) or even while the sphere is still in contact with the water (Fig. 1d and Supplementary Movie 2). In the latter case, the sphere vibrations form a group of nested cavities, or a so-called matryoshka cavity29, named after Russian nesting dolls (Fig. 1d). This phenomenon is in contrast to rigid sphere skipping, for which the cavity is asymmetric, but smooth (Fig. 1e).

Sphere deformation modes

To examine the sphere deformations more thoroughly, we implement a fully coupled numerical finite-element model in Abaqus (ref. 30) (see Methods section). Figure 2a shows the results of a numerical simulation carried out with the same sphere material properties (R, G and ρ s ) and impact conditions (U o and β o ) as the experimental test shown in Fig. 2b. The numerics reveal an elastic wave propagating around the circumference of the sphere in a counter-clockwise direction. We classify this type of wave propagation, depicted in the line drawing of Fig. 2c, as vibration mode 1−. In some cases this elastic wave impacts the air–water interface at time t w , thus initiating a matryoshka cavity (as seen numerically and experimentally in Fig. 2a,b). While tempting to attribute these kinematics to rigid body rotation, we find that the elastic wave propagation time t w is typically much smaller than the measured period of rigid rotation (see Methods and Supplementary Movie 2).

Figure 2: Fluid-structure coupling through sphere deformation modes. (a) An Abaqus numerical model captures the skipping behaviour and reveals the local relative strain in the sphere (colour contours). (b) Under the same conditions as simulated in (a) an experiment shows the formation of a nested cavity (or so-called matryoshka cavity), as also shown by the model. (a,c) The model reveals an elastic wave that propagates in a counter-clockwise direction around the sphere (classified as deformation mode 1−). At time t w , the elastic wave strikes the air–water interface. Two other deformation modes are observed; (d) mode 2: the sphere assumes an ellipsoidal shape with oscillating major and minor axes (Supplementary Movie 1); (e) mode 1+: an elastic wave propagates in the clockwise direction. (f) The attack angle α evolves in time in response to the elastic wave propagation (deformation mode 1− is pictured). The purple lines are experimental measurements of α taken at 1.33 ms intervals. Full size image

The Abaqus numerical model predicts two additional vibration modes, generally occurring with increasing impact speed and/or decreasing shear modulus. In mode 2 (Fig. 2d, Supplementary Movie 1), the sphere assumes an ellipsoidal shape with oscillating major and minor axes. In mode 1+ (Fig. 2e) an elastic wave again propagates around the circumference of the sphere, but in the clockwise direction. Finally, we observe that the attack angle α of the deformed water-contacting face evolves as a result of the elastic wave propagation (Fig. 2a,f and Supplementary Movie 2).

Skip-enhancing mechanisms

Based on our experimental observations and numerical simulations (Figs 1 and 2), we hypothesize that the elastic response of the sphere enhances skipping through two mechanisms: (1) by taking on the shape of a flat disk with an increased wetted area; and (2) by acquiring a favourable attack angle.

To connect the suggested skip-enhancing mechanisms to the vertical force acting on the sphere, we propose an analytical model of the coupled fluid-structure interaction. The sphere is idealized as an incompressible, neo-Hookean hyperelastic solid31 with shear modulus G, radius R and density ρ s . During water impact, the sphere deforms into a disk-shaped ellipsoid inclined at attack angle α(t) to the water surface (Fig. 3a). A set of equations can then be written for the motion of the centre of mass (COM) and sphere deformation in terms of general forces and tractions acting on the body (see Methods section).

Figure 3: Description of sphere deformation during impact. (a) Our analytical model assumes the deformed sphere takes an ellipsoidal shape with principal stretches λ 1 , λ 2 and λ 3 =1/λ 1 λ 2 aligned with the body-fixed {m 1 , m 2 , m 3 } coordinates, respectively. The {m 1 , m 2 , m 3 } coordinate system is inclined at attack angle α relative to the free-surface and has its origin at the sphere’s centre of mass (d 1 ,d 2 ). To model the hydrodynamic force F, we represent the sphere as a circular disk with radius λ eq R and assume that F acts at the point B on the water-contacting face, which moves with velocity described by U B and β B . (b) The attack angle α is measured from simulations using our Abaqus numerical model for spheres undergoing mode 1− deformations, with properties R=25.4 mm, ρ*=1.05 and varying shear modulus; sample α measurements are shown in Fig. 4b–d. For each β o , the attack angle is collapsed by the dimensionless time , which is proportional to the speed of mode 1− elastic waves propagating through the characteristic distance R. These attack angle data are fed into our analytical model (see Methods section). (c) The simulations also show the dependence of λ max on , where λ max is the maximum value of λ eq achieved during impact (individual marker shapes same as for (b)). The numerical results agree with predictions from our analytical model (the zoomed region shows the nine different line styles corresponding to analytical results for nine different spheres: dashed lines R=20.1 mm; solid lines R=26.2 mm; dash-dot lines R=48.8 mm; line width indicates G: thin G=12.3 kPa; middle G=28.5 kPa; and thick G=97.2 kPa). In the limit of small G, we find (valid for shallow β o ). For steeper β o , we find for decreasing . Full size image

To couple the sphere response to the fluid loading, we extend an existing hydrodynamic force model for circular disk-shaped skipping stones4 by approximating the water-contacting face as a circular disk with radius λ eq R, where λ eq approximates the sphere deformation (see Fig. 3a and Methods section). The force is modelled as

where U B and β B describe the velocity of the water-contacting face and the wetted area S w is proportional to (λ eq R)2. The equation for the vertical motion of the COM is then

where d 2 is the vertical coordinate of the COM and g is gravity. We can now predict how the hypothesized mechanisms relate to skipping. First, a larger area S w increases the magnitude of the hydrodynamic force term in equation 2. For a neo-Hookean material, S w increases with decreasing shear modulus G for a given applied compressive stress31. Second, a smaller attack angle increases cos α thereby increasing the vertical force component that lifts the sphere off the surface. We note that the governing equations for the sphere deformation predict a steady-state solution in which the attack angle evolves as in response to a circumferential wave propagating in mode 1− (see Supplementary Note 1). While the sphere deformation during impact is not steady-state, we nonetheless expect α to be governed by the speed of elastic waves propagating in the sphere through the characteristic distance R such that . Therefore, we predict that a more compliant sphere (smaller G) will assume a smaller rate of change of α, thereby increasing the upward vertical force that enables skipping.

Numerical simulations verify the expected dependence of these two mechanisms on G. Measurements from the Abaqus results show that the rate of change of the attack angle scales as for mode 1− deformations (Fig. 3b). Additionally, we find that the maximum value of λ eq achieved during impact, which we call λ max , increases with a decrease in the dimensionless term (Fig. 3c), which is the ratio of material stiffness to hydrodynamic pressure. Therefore, a smaller G yields a larger stretch and larger wetted area, as well as a smaller rate of change of α, as predicted.

To confirm that these mechanisms indeed enhance skipping, we perform experiments and simulations over a range of impact conditions and sphere properties and measure the minimum impact speed required to skip U min , as a function of G (Fig. 4). Above a certain value of G (≈103−104 kPa, depending on β o ), we recover the rigid skipping regime, in which U min is independent of shear modulus but is very sensitive to β o . For rigid spheres impacting above , where ρ*=ρ s /ρ w , prior research suggests spheres may broach (that is, become completely submerged before exiting), but not skip8 (Fig. 4f). For stiffness values below the rigid regime, the elasticity of the sphere becomes important and U min decreases monotonically with decreasing shear modulus. Our analytical model accurately predicts the experimental and numerical results in this regime. The minimum speed is also much less sensitive to β o in the elastic skipping regime and as a result we observe skipping at impact angles nearly three times larger than predicted for rigid spheres (Fig. 4 and Supplementary Fig. 1).

Figure 4: The effect of material stiffness on skipping behaviour. (a) Below a threshold stiffness (G≈103−104 kPa, depending on β o ), the elastic response of the sphere affects skipping. In this elastic regime, the minimum speed required to skip U min decreases monotonically with decreasing shear modulus G, as shown by both experimental and numerical results (triangle and star markers, respectively). In the limit of small G we expect U min ∝G1/5, which is confirmed by the numerics. As stiffness increases above G≈10 kPa, U min deviates from this relation as larger G augments the rate of change of α, thereby reducing the upward vertical force component. Our analytical model also captures this change in behaviour (solid coloured lines). For shallow impact angles, U min becomes insensitive to shear modulus for G>≈103 kPa, indicating that rigid sphere skipping behaviour is recovered. The transition between the elastic and rigid skipping regime occurs at larger G as β o increases. In the rigid skipping regime, U min is very sensitive to β o (also evident in Supplementary Fig. 1). The lower bound of the rigid skipping regime is inferred from the dark triangle symbols, which occur for experiments with . The coloured triangle markers at G=5.66 × 105 kPa result from experiments where the colour gradient on the marker indicates the uncertainty in . The upper bound on the rigid regime corresponds to ; prior literature suggests that for spheres may broach (that is, completely immerse before exiting), but not skip 8 . (b–d) Numerical simulations show that increasing G results in a larger rate of change of α (purple lines) in the elastic skipping regime. Each image shown is 6 ms after impact and the interval between α measurements is 2 ms. (e) As shear modulus increases into the rigid regime, the sphere deformation is negligible (black outline is undeformed sphere contour). (f) We have observed broaching for rigid spheres impacting with . Sphere properties for data in (a): rigid sphere experiments, R=25.8 mm, ρ*=0.959, G=5.66 × 105 kPa; all other data markers, R=26.2 mm, ρ*=0.937 for G≤12.3 kPa and ρ*=1.03 for G>12.3 kPa. The purple error bars are characteristic for experimental data. The numerical error bars represent ±1/2 of the difference in U o between the skipping and non-skipping cases used to compute U min ; error bars are offset for clarity. Full size image

While our results show that reducing shear modulus has the predicted effect on wetted area and attack angle (Fig. 3b,c), it is unclear whether one or both of these mechanisms are responsible for the observed improved skipping performance. To isolate the mechanisms we consider the limiting case of small G, for which and thus cos α≈1 over typical impact timescales. The expected dependence of U min on G in this limit can be rationalized by scaling analysis of equation 2 (see Methods and Supplementary Note 2), which gives

where t c is the collision time (that is, time in which the sphere is in contact with the water). For threshold skipping cases, we expect the characteristic acceleration to be small compared with gravity and thus, to first order, . Furthermore, in the small G limit our numerical modelling shows that (Fig. 3c). Applying this dependence and solving for U min , we find . Figure 4a shows that U min approaches the G1/5 relation in the limit of small G, indicating that the only mechanism by which reducing shear modulus enhances skipping in this limit is through the increased wetted area. However, for larger G (>≈10 kPa), U min deviates from the G1/5 relation as the coupling between shear modulus and α becomes important (Fig. 4a–d). As stiffness continues to increase, ultimately the amplitude of the deformations (affecting both S w and α) become negligible and the sphere is effectively rigid (Fig. 4a,e). Consequently, we conclude that decreasing shear modulus below this rigid boundary causes an increase in the upward vertical force that promotes skipping through both of the hypothesized mechanisms, save for the limit of small G where decreasing shear modulus only affects lift by increasing S w .

Skipping regimes

As the impact events devolve from clear skipping to water entry, we observe a transitional regime characterized by a matryoshka cavity, in which the sphere still skips (Fig. 1d; Supplementary Movie 2). The matryoshka phenomenon occurs when the total contact time of the sphere with water t c is longer than the wave time t w associated with mode 1− elastic wave propagation, such that t c /t w >1. We define t w as the time from impact until the circumferential elastic wave strikes the air–water interface (Fig. 2a,b). Experiments over a range of sphere properties and impact conditions reveal that t c /t w is governed by the impact angle β o and the ratio (Fig. 5). The dependence on these terms can be rationalized by scaling analysis with U o replacing U min in equation 3. When t c ≈t w , the characteristic sphere acceleration is much greater than g such that and we find (see Methods section). Furthermore, based on the propagation speed of mode 1− elastic waves, we expect , which is confirmed experimentally (Supplementary Fig. 2). Combining the scaling for each time gives . We can now examine the evolution of the timescale ratio in the vicinity of the transitional regime in the limit of shallow (β o →0) and steep impact angles, making use of the relationship between λ max and shown in Fig. 3c. Here, we consider steep impact angles to be , where is the maximum impact angle at which we have observed elastic sphere skipping . In the shallow β o limit, t c /t w ≈1 occurs at values of , for which λ max →1 (Fig. 3c); thus, we anticipate for shallow angles. For steep β o , the transitional regime occurs for , for which (Fig. 3c), and we expect . These limiting relations capture the evolution of t c /t w observed experimentally (Fig. 5) and provide insight into the differences observed at different impact angles. We see that for steeper β o , the sphere deformation has a larger effect on the collision time, which gets manifested as a higher sensitivity of t c /t w to .

Figure 5: Matryoshka cavities depend on two timescales. When conditions produce impacts in which the collision time t c is longer than the elastic wave propagation time t w (that is, t c /t w >1) matryoshka cavities form. Experiments show that the timescale ratio depends on and β o , with seemingly minimal dependence on R (marker size indicates R; marker shapes indicate same G as for Fig. 3b). The four coloured patches result from calculations using our analytical model with the same material properties as for the experiments. Varying R in the model over the range experimentally tested results in the spread in the patches, which is small for t c /t w <1. For shallow impact angles, scaling analysis predicts (dashed line), while for steep impact angles we expect (dash-dot line). These limiting trends capture the general evolution shown by the experimental data. Characteristic error bars are shown. Full size image

Based on our findings regarding the transitional skipping regime, we hypothesize that the same dimensionless parameters (that is, β o and ) govern all deformation modes and associated skipping behaviour. An empirical regime diagram indeed shows that these parameters classify all observed skip types (Fig. 6). Mode 1+ is promoted by shallow β o , large U o and/or small G. As stiffness becomes large relative to hydrodynamic pressure, the vibration type traverses the mode 2 and mode 1− regimes. Our analytical model correctly predicts the boundary between the mode 1− skipping (t c /t w <1) and mode 1− transitional (t c /t w >1) regimes (marked by red line on Fig. 6). A small, radius-dependent overlap region exists between the transitional regime and water entry. Finally, we have observed that the vibration mode associated with the transitional and water entry regimes is exclusively mode 1−.

Figure 6: Classification of impact phenomena for single and multiple skip events. Experiments show that all vibration modes and accompanying skipping behaviour are classified by the ratio of material stiffness to hydrodynamic pressure and β o (triangles: mode 1+; crosses: mode 2; black circles: mode 1− skipping; grey circles: mode 1− transitional; x-marks: water entry). The experimental data shown are for all three elastic sphere radii. The coloured regions provide a visual boundary of the skipping (blue), transitional (yellow) and entry (red) regimes. Our analytical model accurately predicts the boundary between mode 1− skipping and transitional events; the red line is found from the intersection of the analytical results with t c /t w =1 in Fig. 5. Given initial conditions of the first skip, our analytical model also predicts successive impact types in a multiple skip trajectory (blue squares). A multiple skip experiment validates these predictions within the limits of our test facility (purple markers; R=26.2 mm, G=12.3 kPa, ρ*=0.937). Images from the multiple skip experiment are shown in Fig. 8. Characteristic error bars are shown in grey. Full size image

Skipping sphere rebound

We quantify the rebound characteristics of skipping spheres by computing the normal and tangential coefficients of restitution, and , where U e is the exit velocity and n and t refer to components normal and tangential to the free-surface, respectively. Here, we find some interesting similarities between skipping elastic spheres and liquid droplets bouncing on inclined liquid layers28. Following the work on bouncing liquid droplets, we plot and as a function of (Fig. 7), which is equivalent to the normal Weber number with shear modulus G replacing the Laplace pressure σ/R d , where R d is the droplet radius and σ is surface tension of the liquid droplet. Interestingly, we find and , which are the same scaling relations found for liquid droplets bouncing on inclined liquid layers28. The restitution coefficients deviate from these trends when t c /t w >1. In these cases, sphere vibrations interact with the water via the matryoshka cavity causing a significant decrease in bouncing efficiency. For both elastic skipping spheres and bouncing droplets, the restitution coefficients are always less than one as part of the initial translational kinetic energy goes into post-impact vibrations in the sphere or droplet28,32.

Figure 7: Restitution coefficients of elastic spheres skipping on water. Experiments reveal the dependence of the square of the normal (a) and tangential (b) restitution coefficients on (grey-filled symbols denote mode 1− transitional skips; marker shapes same as for Fig. 6; ρ*=1±0.08). We find that and scale with in the same way that for bouncing liquid drops28 they scale with . Thus, shear modulus G for skipping elastic solid spheres plays the role of the Laplace pressure σ/R d for bouncing liquid droplets. Because rebound is more efficient in the tangential direction, β o decreases (thus becoming more favourable) with each skip. In addition, both restitution coefficients increase with decreasing impact speed until the sphere enters the mode 1− skipping and transitional regimes. These two effects combine to enable multiple skip trajectories (Figs 6 and 8, Supplementary Movie 3). The demise of a multiple skip event starts when t c /t w >1 and a matryoshka cavity forms (below red lines); and decrease rapidly when this occurs. This interpretation is confirmed by a multiple skip experiment; the purple markers on each plot correspond to the multiple skip experimental data shown in Figs 6 and 8. Full size image