2.1. Experimental set-up

The experimental set-up consists of a clear, translucent acrylic container and high-speed monochrome cameras. The cross-section of the container is large, 200 $\times$ 200 $\mathrm{mm}^2$ , 20 times larger than the diameter of the largest sphere utilized in the present study. The depth of the container is 400 mm. By creating a suction pressure at the end of a needle tip, the spheres are firmly held and are released by interrupting this suction pressure. The free-falling sphere impacts the liquid with no rotation, as confirmed by the images captured prior to sphere impact. All the experiments are performed in a closed room with an ambient temperature of 25 $^\circ$ C.

Figure 1. Two orthogonally and synchronized cameras record the motion of a rigid sphere after impacting the air–water interface. Panel (a) illustrates schematically the experimental arrangement and the coordinate system used for the imaging data; (b) a sample trajectory of a sphere with a density ratio ( $\rho ^*$ ) of 2.16 and diameter 10 mm is viewed by the two cameras in the $z$ – $y$ and $z$ – $x$ planes. Image sequences of 10 mm spheres with an initial Reynolds numbers of 15 700 are shown for the density ratios (c) $\rho ^* = 2.16$ , (d) $\rho ^* = 3.26$ , (e) $\rho ^* = 6.08$ and (f) $\rho ^* = 7.92$ . The time step between consecutive images is 60 ms for (c), 40 ms for (d), 25 ms for (e) and 20 ms for (f). The black dashed line represents a pure vertical trajectory. The circular data points in the last images indicate the position of the sphere at equal time intervals.

Experiments were conducted using two cameras (Phantom VEO E-340-L with Nikon microlenses of 28 mm, spatial resolution of 192 $\unicode{x03BC}$ m pixel^–1 and 24–85 mm, spatial resolution of 199 $\unicode{x03BC}$ m/pixel) placed orthogonally to one another. Backlighting consisted of a 30 W monochromatic light source and a diffuser. The two cameras were synchronized and recorded the sphere motion at 1000 f.p.s., with a resolution of 1152 $\times$ 1100 pixels (100 mm $\times$ 160 mm in the object plane) and an exposure time of 400 $\mu s$ . The exposure time leads to a maximum relative motion blur of 28 % of the sphere diameter for the smallest sphere (4 mm) with the highest initial velocity (2.8 m s^–1). However, the motion blur for the 10 mm sphere is only 11 % upon impact, and for all spheres, this motion blur decreases rapidly after impact since the velocity immediately goes through a strong deceleration phase. Furthermore, the edge detection routine remains the same at all time steps, so the error through motion blur for the relative motion is neglected. The data are collected while the sphere descends until it traverses out of the field of view of either of the cameras. A pictorial view of the experimental set-up is shown in figure 1(a). Figure 1(b) shows a sample trajectory of a 10 mm sphere with density ratio of $\rho ^*=$ 2.16, as observed with the two cameras.

To further illustrate the raw data with which we are working, we present figure 1(c–f), showing example trajectory traces from a single camera for the 10 mm diameter spheres, all with the same initial Reynolds number ( $Re_i$ defined as $v_{{i}}D/\nu$ ), but for four different density ratios ( $\rho ^*$ ). It is apparent from these visualizations that the lighter spheres (figure 1 c,d) exhibit higher lateral displacements away from the vertical, whereas the heavier spheres (figure 1 e,f) have a more ballistic-like trajectory, as expected. The associated lift forces causing these lateral displacements will be quantitatively derived below.

The diameter ( $D$ ) of the spheres are 4 mm, 6 mm, and 10 mm, as measured using a vernier caliper, and their mass ( $m$ ) is determined using an electronic weight balance (Ohaus), yielding $\rho _{{s}}$ . The impact velocity ( $v_{{i}} \approx \sqrt { 2gh }$ ) is derived from the initial release height ( $h$ ) of the sphere, measured from its centre to the air–water interface. Three impact velocities of 1.4, 1.98, and 2.8 $m/s$ are considered in the present study. The Reynolds number ( $Re_i$ = $v_{{i}}D/\nu$ ) upon impact ranges from 6300 to 31 500. Given that the critical Reynolds number for a sphere in steady flow lies well above 10 $^5$ , we assume that the boundary layers on all spheres throughout all phases of the trajectory remain laminar.

For these definitions, $\nu$ is the kinematic viscosity (0.89 $ \mathrm{mm^2\, s^{-1}}$ ) of the fluid, and $g$ is the gravitational acceleration. Throughout the following discussion, length scales are rendered dimensionless using the sphere diameter $D$ , velocities with the impact velocity $v_{{i}}$ and time scales using $D/v_{{i}}$ . Dimensionless quantities are designated with the superscript ‘ $*$ ’, and unit vectors are written in boldface font.

This range of Reynolds number was chosen because sphere impact at lower initial Reynolds numbers showed no significant difference in trajectory behaviour to those conducted at higher Re, and higher Reynolds number impacts led to an entry air cavity being formed behind the penetrating sphere. The entry phenomenon was outside the scope of the present study, but has been investigated at length in connection with the impact and penetration of superhydrophic spheres (Speirs et al. Reference Speirs, Mansoor, Belden and Truscott2019). The density ratio was not chosen beyond 7.92, since at this value, no significant change was observed from the value 6.08. At these high density ratios the sphere exhibits a nearly vertical trajectory, which is seen in figure 1(e,f).

2.2 Methodology

For the presentation of results and the subsequent analysis invoking a force model, it is convenient to work in a natural coordinate system, i.e. in a coordinate system in which the unit vectors of the accompanying triad of the pathline are used as basis vectors. This coordinate system is pictured in figure 2(a). The unit vector tangential to the pathline is given as

(2.1) \begin{equation} \mathbf{t}=\frac {\boldsymbol{v}_{{s}}}{| \boldsymbol{v}_{{s}} |}, \end{equation}

where $\boldsymbol{v}_{{s}}$ is the velocity vector of the sphere along the pathline $s$ . Here $\sigma$ is the coordinate in the direction $\mathbf{t}$ , $n$ is the coordinate in the direction of the principal normal vector $\mathbf{n}_\sigma = R\mathrm{d}\mathbf{t}/\mathrm{d}\sigma$ and $b$ the coordinate in the direction of the binormal unit vector $\mathbf{b}_\sigma =\mathbf{t} \times \mathbf{n}_\sigma$ . Here $R$ is the radius of curvature of the pathline in the plane spanned by the normal vectors $\mathbf{t}$ and $\mathbf{n}_\sigma$ . The velocity vector $\boldsymbol{v}_{{s}}$ is understood to also represent the slip velocity in the equation of motion, since the pool is quiescent upon sphere impact.

Figure 2. (a) Definition of the natural coordinate system based on the sphere pathline. (b) Definition of the plane of curvature and the total lift magnitude $F_{{L}}$ .

2.3. Image processing

The video images are processed using ImageJ software to subtract the background and to create a binarized image of the sphere inside the liquid. Subsequently, an in-house MATLAB code is employed to determine the position of the sphere. The detection of the bottom most point of the sphere is accomplished by utilizing an edge detection technique. The air–water interface is established as the reference point for spatial coordinates ( $z=0$ ). Similarly, the time instant a sphere makes first contact with the air–water surface is established as the reference point for time ( $t=0$ ). This time is taken as the first frame in which contact with the liquid has been made.

The accumulated dimensionless path length that the sphere covers is denoted $s^*=s/D$ , whereby $s=0$ at the air–water interface. The dimensionless lateral displacement at each time step is expressed as $r^*=\sqrt {x^2 + y^2}/D$ . The displacement data are smoothed using an in-house MATLAB function ’loess’, a method that involves using linear regression in a locally weighted scatter plot. Subsequently, velocity magnitude $|\boldsymbol{v}_{{s}} |$ is computed by differentiating the smoothed displacement data after first fitting the data with a quintic spline function. The dimensionless velocity magnitude is given as $v_{{s}}^*=| \boldsymbol{v}_{{s}}|/v_{{i}}$ . Acceleration data was obtained using a differentiation of the smoothed velocity data. Time is made dimensionless using the diameter of the sphere and the impact velocity, i.e. $t^*=t v_{{i}}/D$ .

Given the $x$ , $y$ and $z$ coordinates of the sphere as a function of time, the curvature of the pathline can be computed as

(2.2) \begin{equation} \kappa = \frac {1}{R} = \frac {\sqrt {(z''y'-y''z')^2 + (x''z'-z''x')^2 + (y''x'-x''y')^2}}{(x'^2 + y'^2 + z'^2 )^{3/2}}, \end{equation}

where the primes indicate first and second differentiation with respect to time. These differentials are computed numerically using second-order central differences. The angular frequency along the curved pathline is then given as $\omega = | \boldsymbol{v}_{{s}}|/R$ , which is a necessary quantity in computing the centrifugal force acting on the sphere, which will be related to the lift force, examined in the following section.

2.4. Force model

In formulating the equation of motion for the sphere, the BBO equation (Zhu & Fan Reference Zhu and Fan1998) is used, which includes body forces ( $F_{{G}}$ – weight, $F_{{B}}$ – buoyancy), apparent forces ( $F_{{H}}$ – Basset or history term, $F_{{A}}$ – added mass), hydrodynamic forces ( $F_{{D}}$ – viscous and pressure forces combined as drag) and inertial forces ( $F_{{I}}$ ). We will neglect the Saffman lift force (Saffman Reference Saffman1965), applicable only in sheared flow, and any rotational-lift force (Rubinow & Keller Reference Rubinow and Keller1961), applicable only with rotation/spin of the sphere.

The scalar momentum equation expressed along the direction of motion/pathline can be written as (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011)

(2.3) \begin{align} \underbrace {\frac {1}{6} \rho _{{s}} \pi D^3 \frac {\mathrm{d}v_{{s}}}{\mathrm{d}t}}_{{\textit{Inertial}\,\textit{force}}} & = \underbrace {\frac {1}{6} \pi D^3 (\rho _{{s}} - \rho _{{l}}) g \cos {\chi }}_{{\textit{Body}\,\textit{forces}}} - \underbrace {\frac {1}{8} C_{{D}} \rho _{{l}} \pi D^2 v_{s} ^2}_{{\textit{Drag}\,\textit{force}}} - \underbrace {\frac {1}{6} C_{{A}} \rho _l \pi D^3 \frac {\mathrm{d}v_{{s}}}{\mathrm{d}t}}_{{\textit{Added}\,\textit{mass}}} \nonumber \\ & \quad - \underbrace {\frac {3}{2} D ^2 \sqrt {\pi \mu \rho _{{l}}} \int _{0}^{t} \frac {1}{\sqrt {t-\zeta }} \frac {\mathrm{d}v_{{s}}}{\mathrm{d}\zeta }\mathrm{d}\zeta }_{{\textit{Boussinesq}{-}\textit{Basset}\,\textit{term}}}, \end{align}

where $C_{{D}}$ is the coefficient of drag, and $C_{{A}}$ is the coefficient of the added mass force. Note that no hydrodynamic lift force has been included in (2.3). This force will be introduced and discussed below.

The added mass coefficient $C_{{A}}$ is usually taken as 0.5 (Guo Reference Guo2011) and expresses the kinetic energy imparted into the surrounding fluid through acceleration/deceleration of the sphere. The Boussinesq–Basset term captures the viscous force change due to boundary-layer development on an accelerating or decelerating submerged body. These viscous forces are not expected to be significant relative to the inertial and pressure forces involved over large portions of the trajectory; hence, the Boussinesq–Basset term will be neglected in the following analysis (Nouri et al. Reference Nouri, Ganji and Hatami2014). This is not to say that the transient boundary layer development does not play a central role in determining the trajectory and speed of the sphere, but it is expected to be more through the separation and wake behaviour due to the state of the boundary layer; hence, through the resulting pressure distribution around the sphere. These forces would make themselves apparent in the above equation as variations in the drag force.

The initial Reynolds numbers encountered in this study all lie below approximately $3 \times 10^4$ and throughout most of the sphere trajectory the values are much lower. For very similar Reynolds numbers, Truscott et al. (Reference Truscott, Epps and Techet2012) used a drag coefficient of 0.5. They viewed the unsteady added mass as part of the pressure force acting on the sphere, i.e. in the net hydrodynamic force. In the present study, the experimental data allows all of the remaining terms in (2.3) to be evaluated; thus, the value of $C_{{D}}$ can and will be computed at each time step.

The terminal velocity of a sphere can be determined by establishing an equilibrium in which the net force acting on the sphere is reduced to zero. The expression for the dimensionless terminal velocity ( $v_{{t}}^*$ ) can be obtained by equating the buoyancy force with the drag force plus the weight:

(2.4) \begin{align} v_{{t}}^* = \frac {1}{v_{{i}}}\sqrt { \frac {4 (\rho ^* - 1)gD }{3C_{{D}}} }. \end{align}

In this equation, a drag coefficient must be prescribed, and a value of $C_{{D}}=0.5$ has been used, assuming at this stage no significant acceleration or deceleration.

If the pathline has non-zero curvature ( $\kappa$ ), this implies a force ( $\overrightarrow {F}_{{\def\negativespace{} L}}$ ) acting in the $-\mathbf{n}_\sigma$ direction, i.e. perpendicular to $\mathbf{b}_\sigma$ , as depicted in figure 2(b). Since lift force is defined as the force acting perpendicular to the direction of motion, the total lift force magnitude, $||\overrightarrow {F}_{{\def\negativespace{} L}}||=F_{{L}}$ , can be computed as the sum of the buoyancy ( $F_{{B}}$ ), gravity ( $F_{{G}}$ ) and the hydrodynamic (lift) force ( $F_{\textit{HL}}$ ) acting in the $-\mathbf{n}_\sigma$ direction. Thus, to compute the hydrodynamic lift force from the sphere trajectory, the direction of the unit normal $\mathbf{n}_\sigma$ must be determined. Recognizing that $\mathbf{n}_\sigma =\mathbf{b}_\sigma \times \mathbf{t}$ , it is necessary to first compute the unit normal $\mathbf{b}_\sigma$ from the trajectory data. This is done using three points along the trajectory, pictured in figure 2(b) as points PQR, whereby point P is the position at which the local force is to be estimated. Point R represents the position just before point P, and point Q is the position following point P along the trajectory. The vector $ \overrightarrow {\mathbf{PQ}}$ is derived by computing the forward difference between the coordinates of point P and point Q. Similarly, the vector $ \overrightarrow {\mathbf{PR}}$ is determined from the points P and R. This procedure is applied at every point along the trajectory of the sphere from the air–water interface.

Knowing the $x,y,z$ coordinates of the points P, Q and R, the unit normal to the subscribed triangle is given by $\mathbf{b}_\sigma = \overrightarrow {\mathbf{PR}} \times \overrightarrow {\mathbf{PQ}}/||\overrightarrow {\mathbf{PR}}\times \overrightarrow {\mathbf{PQ}}||$ . The trajectory unit normal $\mathbf{t}$ , can be computed using forward differencing around point P allowing $-\mathbf{n}_\sigma =-\mathbf{b}_\sigma \times \mathbf{t}$ to be computed. The angle $\psi$ is the angle between $-\mathbf{n}_\sigma$ and the direction of gravity (see figure 2 a). Only the components of buoyancy and gravity along the direction $-\mathbf{n}_\sigma$ can contribute to the total magnitude of the lift force $F_{{L}}$ , which itself must be equal to the centripetal force; hence, $(F_{{G}}-F_{{B}})\cos {\psi } + F_{\textit{HL}}=F_{{L}}=m\omega ^2R$ , and with the angular frequency given by $\omega = | \boldsymbol{v}_{{s}} |/R$ , we have

(2.5) \begin{equation} (F_{{G}}-F_{{B}})\cos {\psi } + F_{\textit{HL}}=\frac {m | \boldsymbol{v}_{{s}} |^2}{R}. \end{equation}

Note that this equation is a scalar equation, since through the $\cos \psi$ factor, the gravitational and buoyancy forces have been projected onto the $-\mathbf{n}_\sigma$ vector. Moreover, although the total lift force defined in this manner will always be positive and directed towards the origin of the local curvature radius, as shown in figure 2(b), the hydrodynamic lift force could become negative. A negative hydrodynamic lift force would act to reduce the local curvature of the sphere trajectory, i.e. straighten the trajectory. If now the buoyancy and gravity forces are known, then the hydrodynamic portion of the lift force $F_{\textit{HL}}$ can be computed and plotted as a function of time or displacement of the sphere. The hydrodynamic lift force is made dimensionless in the form of a lift coefficient, $C_{{L}}$ ,

(2.6) \begin{equation} C_{{L}}=\frac {F_{\textit{HL}}}{\frac {1}{2}\rho _{{l}}| \boldsymbol{v}_{{s}} |^2 A}, \end{equation}

where $A$ is the projected area in the direction of motion, in this case $\pi D^2/4$ . In some of the results presented below, an abrupt jump in lift coefficient was observed. An explanation of this result is given in the Appendix.

Before presenting the results of the experiments, the origin of a lift force arising from an asymmetric wake and the interrelation between drag and lift from such wakes will be phenomenologically discussed. The initial Reynolds number lies in the approximate range $6300 \lt Re_i \lt 31\,500$ ; however, the sphere decelerates during its trajectory to values of approximately 10 % of the impact velocity, thus, the total encountered Reynolds number range is approximately $630 \lt Re_i \lt 31\,500$ . This range lies in the Newton regime of drag coefficient, in which $C_{{D}}$ takes an almost constant value over all Reynolds numbers and which is significantly below transitional Reynolds numbers for even roughened spheres. On the other hand, the data on which this statement is based come from experiments in which the flow is steady, i.e. no acceleration or deceleration. Thus, accepted drag coefficients for steady flow may not necessarily be applicable over the entire sphere trajectories of the present experiments.

Although the flow and wake may be statistically symmetric around the sphere in this Reynolds number range, the instantaneous wake can be highly asymmetric. Thus, even though the time averaged lift coefficients in this Reynolds number range may be zero, instantaneous fluctuations of lift can occur and have been demonstrated through numerical simulations (Constantinescu & Squires Reference Constantinescu and Squires2004; Yun et al. Reference Yun, Kim and Choi2006). The transition of the boundary layer to a turbulent state after separation from the sphere can be irregular, causing it to temporarily reattach to the surface of the sphere and separate farther downstream (Hadžić et al. Reference Hadžić, Bakić, Perić, Šajn and Kosel2002). Both Taneda (Reference Taneda1978) and Hadžić et al. (Reference Hadžić, Bakić, Perić, Šajn and Kosel2002) observe a progressive wave motion around the sphere for $10^4 \lt {Re}\lt 3.8 \times 10^5$ by which the separation points rotate around the sphere randomly. Achenbach (Reference Achenbach1974) also observed this for a very similar Reynolds number range. Such an irregular boundary-layer separation will also affect the drag, since the drag arises from the integration of the pressure around the sphere and is thus highly correlated with the location of flow separation: a fluctuating separation will yield a fluctuating drag force. However, an asymmetric separation and vortex shedding will not only influence the drag, but the asymmetry will also mean that the resultant force will no longer be aligned with the flow direction; this force will therefore, have both drag and lift components. For a free moving body this results in a change of trajectory and a new orientation of the drag and lift forces in a lab-fixed coordinate system.

In summary, our force model is cast in a natural coordinate system and all changes in the motion speed of the sphere are attributed to a change in drag plus the body forces along the direction of motion. All direction changes are attributed to a hydrodynamic lift force plus the body forces acting in the direction of local pathline curvature. A positive lift force increases the pathline curvature, a negative lift force decreases the pathline curvature.

3.1. Kinematics of sphere motion

Before discussing the results of the parametric variations conducted in this study, a sample set of data will be examined to illustrate the kinematics of the sphere motion, as shown in figure 3. This figure shows a three-dimensional rendition of the trajectory (see figure 3 a), dimensionless pathline distance ( $s^*$ ) (figure 3 b), dimensionless instantaneous velocity magnitude ( $v_s^*$ ) (figure 3 c) and the dimensionless lateral deviation $r^*$ (figure 3 d) for a 10 mm sphere ( $\rho ^* = 2.16$ ) impacting the pool with a Reynolds number of 31 500. The dimensionless velocity magnitude ( $v_s^*$ ) of the sphere in a three-dimensional rendition is depicted using a colour bar in figure 3(a). Two dashed lines have been added to figure 3(b–d), the first denoting the dimensionless time at which the sphere first deviates more than 10 % of its diameter from the vertical trajectory ( $r^*\gt 0.1$ ), and the second denoting the time at which the sphere attains a minimum dimensionless velocity ( $v_s^*$ ). These timelines divide the total time span into three phases, designated here as follows.

1. (i) Submersion phase, exhibiting a nearly vertical trajectory and extending up to the time at which the lateral displacement remains below 10 % of the diameter of the sphere.

2. (ii) Deceleration phase, during which the sphere velocity continues to decrease to a minimum value.

3. (iii) Settling phase, the remaining time during which the sphere velocity tends towards its nominal terminal velocity in the vertical direction.

Notably, the end of the deceleration phase (minimum velocity) coincides with the time at which the lateral deviation of the pathline ( $r^*$ ) begins to exhibit strong variations between repetitions of the same experiment. This behaviour is illustrated in figure 4, in which typical lateral deviation curves are shown for repeated experiments using the 4 mm, 6 mm and 10 mm diameter spheres with $\rho ^*=2.16$ . Any changes in $r^*$ must necessarily be associated with a lift force, as discussed above.

In figure 4 the variation of the angles $\psi$ (figure 4 d–f) and $\chi$ (figure 4 g–i) are also shown over dimensionless time. These results will be further discussed in § 3.2; however, their physical interpretation is briefly described here. The angle $\chi$ is simply the angle between the instantaneous pathline and the vertical. It starts with the value zero and ends at zero if the sphere has reached its terminal velocity downward. The angle $\psi$ is the angle between the plane of pathline curvature and gravity; hence, this angle expresses to what extent the body forces contribute to instantaneous lift. This angle must start at 90 $^\circ$ and would also end at 90 $^\circ$ if the sphere was moving vertically.

Figure 3. Motion kinematics of a 10 mm sphere ( $\rho ^*=2.16$ ) entering the pool with an initial Reynolds number of 31 500. (a) Visualization of three-dimensional trajectory; (b) dimensionless pathline distance over dimensionless time; (c) instantaneous dimensionless velocity over dimensionless time; (d) dimensionless lateral deviation over dimensionless time. The horizontal dashed line in this graph represents the terminal velocity computed according to (2.4). The dimensionless velocity ( $v_s^*$ ) of the sphere is represented by the colour bar in (a). Data points are spaced equally in time with time step 6 ms.

Figure 4. Examples of measured dimensionless lateral deviations ( $r^*$ ), angle $\psi$ and angle $\chi$ when repeating an experiment: (a,d,g) 4 mm sphere at $Re_i = 8900$ ; (b,e,h) 6 mm sphere at $Re_i = 18\,900$ ; (c,f,i) 10 mm sphere at $Re_i = 22\,300$ . All spheres have $\rho ^*=2.16$ . The time step between the two successive points in (a), (b) and (c) is 6 ms. For clarity, angle $\psi$ is plotted as lines in (d)–(f), and angle $\chi$ is plotted as lines in (g)–(i).

The interpretation of $r^*$ is further illustrated and discussed in the dimensionless trajectory plots shown figure 5. Here the dimensionless lateral deviation $r^*$ is plotted against the dimensionless depth $z^*$ for several density ratios $\rho ^*$ and for the three sphere diameters. Note that the axes in these graphs vary since, in dimensionless terms, the observation volume of the acrylic container depends on the sphere diameter. It is apparent that not all of the spheres reach this final state within the available observation volume of the acrylic container. On the other hand, this is not considered a limitation, since, as will be shown in the next section, the fluctuating lift force has decreased significantly before the terminal velocity is reached.

Figure 5. Sample dimensionless lateral deviation ( $r^*$ ) versus corresponding dimensionless depth ( $z^*$ ) for spheres of various density ratios. (a) A 4 mm sphere at $Re_i = 6300$ ; (b) 6 mm sphere at $Re_i = 9400$ ; (c) 10 mm sphere at $Re_i =15\,700$ . The time step between trajectory points is 10 ms.

The dimensionless time boundaries between the three phases of motion – submersion, deceleration, settling – are remarkably constant over initial Reynolds number, sphere diameter and density ratio, as shown in figure 6 for the dimensionless time at which the submersion phase and deceleration phase ends. Truscott et al. (Reference Truscott, Epps and Techet2012) also measured the penetration trajectory of hydrophilic spheres with various density ratios, all impacting with the same Reynolds number. Their measured dimensionless lateral displacement is in excellent agreement with the end of the submersion phase shown in figure 6. Only one sphere trajectory in the work of Truscott et al. (Reference Truscott, Epps and Techet2012) reached a clear termination of the deceleration phase (minimum velocity in their figure 6 b), and this occurred at a dimensionless time of $t^*=18$ , also in good agreement with our value.

Figure 6. Dimensionless time denoting end of submersion phase and end of deceleration phase of sphere motion as a function of initial Reynolds number and density ratio. The error bars express one standard deviation computed from four repetitions of the same experiment. For densities $\rho ^* = 6.08$ and $7.92$ of the 10 mm sphere, the termination of the deceleration phase is not captured, due to the sphere moving beyond the field of view. For the density ratio 7.92, the termination of the deceleration phase for the 4 mm and 6 mm sphere was also not clear, since the velocity did not exhibit a distinct minimum.

The similarity among the many experiments is further underlined in figure 7, where the dimensionless magnitude of velocity and acceleration of spheres with a density ratio of $\rho ^* = 2.16$ over time and for all sphere diameters and initial Reynolds numbers are shown. Although the previous results in figure 4(a–c) and 5 indicate large variations in specific trajectories, in dimensionless terms, both the velocity and acceleration exhibit remarkable uniformity, confirming that the dimensionless scaling is appropriate. Furthermore, this uniformity is evident also for the other sphere density ratios. As such, these universal curves represent a predictive tool for estimating the deceleration of penetrating spheres within the parameter limits described above.

Figure 7. Magnitude of dimensionless (a) velocity ( $v_s^*$ ) and (b) acceleration ( $a_s^*$ ) for various sphere diameters and impact $Re_i$ ( $\rho ^* = 2.16$ ). The dashed blue line represents the end of the submersion phase, while the magenta dashed line illustrates the end of the deceleration phase.

Figure 8. Penetration of a 10 mm sphere of density ratio $\rho ^*=2.16$ at $Re_i=22\,300$ . (a) Three-dimensional rendition of pathline in dimensionless time; (b) dimensionless lateral distance $r^*$ for four repetitions of the same experiment; (c) drag coefficient as a function of dimensionless time. (d) Lift coefficient as a function of dimensionless time. Panels (e) and (f) represent angle $\psi$ and angle $\chi$ as a function of dimensionless time, respectively. The graphs (c)–(f) correspond to the red curve shown in (b). The colour bar depicted in (a) illustrates the variation in dimensionless velocity ( ${v_s}^*$ ) of the sphere. The time step between the two successive points plotted in (b) is 6 ms.

3.2. Dynamics of sphere motion

In discussing the dynamics of sphere motion, we again begin by examining an example set of data, in this case pertaining to the 10 mm sphere impacting at a Reynolds number of 22 300, and shown in figure 8. In this figure and in subsequent figures, for clarity, the graphs are plotted only with lines and no symbols. For each of the experimental conditions, four repetitions were performed, and the entire digitized trajectories ( $x$ , $y$ , $z$ coordinates) over time are available and documented in Billa et al. (Reference Billa, Josyula, Tropea and Mahapatra2024). Although there exists a high degree of randomness in the actual trajectories for each experimental repetition, the strong commonalities alluded to above will now be discussed in terms of the three phases, beginning with the submersion phase.

3.2.2. Deceleration phase

We now move to the deceleration phase, marked by the first deviations from the initial vertical trajectory (see figures 8 a and 8 b). The trajectory begins this phase with an abrupt curvature, which is reflected in the rapidly changing angles $\psi$ (figure 8 e) and $\chi$ (figure 8 f), evidently a result of a rapid increase in both drag and lift. Since the first trajectory curvature arises from the vertical state, the lift at this stage can only be attributed to hydrodynamic lift, since the body forces are acting vertically. Assuming the wall shear stress contributes little to the overall drag or lift, the integrated pressure over the sphere with a skewed wake would result in a drag and lift that would no longer be aligned with the motion axis of the sphere. The large values of drag and lift in this phase are, therefore, clearly related to asymmetric wake effects since these would result in a skewed base pressure area on the rear of the sphere.

Whereas the drag coefficient rises to a maximum value and then decreases again towards the end of the deceleration phase (figure 8 c), the lift coefficient exhibits an abrupt jump from positive to negative values (figure 8 d). This is very typical of all other data sets and can be explained by examining changes in the angle $\psi$ . As discussed in the Appendix, when the trajectory goes through a ’projected’ inflection point, this angle can exhibit sharp jumps in magnitude, since the orientation of curvature, i.e. $\mathbf{n}_\sigma$ , will change direction. Such a change is seen in figure 8(e), corresponding to a sharp drop in lift coefficient. We conclude that this arises due to a reorientation of the wake, such that the asymmetry changes orientation on the sphere. During this change, the wake is momentarily symmetric, leading to a short period of zero hydrodynamic lift. However, whereas the wake orientation is changing, the wake area apparently remains larger, resulting in a persistently large drag coefficient, i.e. the wake base pressure still acts over an area undiminished in magnitude. Once through the point of changing the sign of $\mathbf{n}_\sigma$ , the lift force is again high and in the direction of the new unit vector $-\mathbf{n}_\sigma$ , i.e. positive.

What is particularly noteworthy is that the lateral deviation $r^*$ exhibits a constant slope over the entire deceleration phase, almost identical in all repetitions of experiments at the same initial conditions. This infers a sustained hydrodynamic lift force in a constant direction. This deduction must be explained. The net gravity and buoyancy force acts downward, and would act to make the trajectory vertical. If the sphere is not moving vertically downwards, then these body forces contribute to lift through angle $\psi$ . Whereas a vertical trajectory results in a constant value of $r^*$ , this is not observed in the deceleration phase. Therefore, over this entire deceleration phase, there must be a sustained and constant hydrodynamic lift force in both magnitude and direction, counteracting the ever-present body force acting downwards. This suggests that the asymmetric vortex shedding, which is understood to be the origin of the lift force, once established, remains approximately constant in its orientation on the sphere throughout this deceleration phase.

Figure 9. Lateral deviations ( $r^*$ ) and the ratio of the vertical component of the hydrodynamic lift force to the total body forces ( $\phi$ ) are shown for various initial Reynolds numbers. Results are shown for spheres of diameter $D =4$ mm in (a) and (d), $D = 6$ mm in (b) and (e) and $D = 10$ mm in (c) and (f). In (a (b) and (c), the time interval between successive points is 4 ms. For clarity, $\phi$ is plotted as continuous lines in (d, (e) and (f).

In this regard, the behaviour of the three different sized spheres shown in figure 4 at the end of the deceleration phase and in the settling phase is distinctly different. Whereas the trajectories of the 10 mm sphere (see figure 4 c) tend towards a vertical motion ( $r^*$ constant) following the deceleration phase, the 6 mm sphere (see figure 4 b) exhibits larger fluctuations in trajectory inclination, albeit fluctuating around an approximate constant value. The 4 mm sphere (see figure 4 a) continues on an inclined trajectory throughout the entire field of view. In the decelerating phase, the body forces are competing with the hydrodynamic lift force to achieve a vertical trajectory, as described above. Having computed the hydrodynamic lift force coefficient at all times, the ratio of the vertical component of the hydrodynamic lift force to the sum of the body forces can be given as

(3.1) \begin{equation} \phi =\frac {F_{\textit{HL}}\cos {(180-\psi )}}{F_{{G}}-F_{{B}}} = \frac {3}{4g}\frac {C_{{L}} v_{{s}}^2 \cos {(180-\psi )}}{D(\rho ^* -1)}, \end{equation}

and is shown for the three sphere diameters 4, 6 and 10 mm, and for varying initial Reynolds numbers in figure 9. From this diagram, this value does tend towards unity throughout the deceleration phase. From (3.1) it can be seen that the body forces will tend to dominate for larger spheres and from the diagram, the values for the 10 mm sphere appear to decrease earlier in dimensionless time, suggesting that the body forces begin to dominate the hydrodynamic lift forces earlier for the heavier sphere. This is intuitively correct, since the sphere is decelerating and the hydrodynamic lift force scales with velocity squared and diameter squared, whereas the body forces scale with the cube of the diameter and remain constant. What remains unclear is the exact reason that such a stable hydrodynamic lift force, i.e. an asymmetric wake, is maintained throughout the deceleration phase, and for the 4 mm sphere even throughout the settling phase. Further insight would require a time resolved measurement of the velocity field around the sphere throughout this phase.

Phenomenologically, these observations are in good agreement with experiments from Truscott et al. (Reference Truscott, Epps and Techet2012), where also all hydrophilic spheres initially exhibited a sharp deviation from the vertical, accompanied by rapid deceleration. Similarly, they attribute the trajectory change from the vertical to asymmetrical vortex shedding, as observed and explained also by Horowitz & Williamson (Reference Horowitz and Williamson2010). The exact direction of trajectory change was non-repeatable, as in the present experiments. What the present experiments reveal, is that while the direction of trajectory change is non-repeatable, the growth of lateral displacement remains constant in dimensionless time until a minimum velocity is reached, as indicated in figure 7(a).

The deceleration phase ends when the absolute velocity magnitude reaches a minimum, and again, this occurs for all investigated cases at an approximately constant value of $t^* \approx$ 20–21, as seen in figure 6. Typically, this minimum velocity is close to and sometimes slightly less than the terminal velocity, given in (2.4). For instance, in figure 3 for the 10 mm sphere, a slight acceleration is seen following the deceleration phase. This is likely due to a relaxation of the asymmetric wake into a symmetric, steady state condition, which is then maintained throughout the final settling phase. Note that Truscott et al. (Reference Truscott, Epps and Techet2012) also observe for lighter spheres a velocity slightly less than the terminal velocity, followed by a light acceleration (their figure 5). They call this ‘underdamped behaviour’, but do not elaborate on its physical origins.

Assuming the end of the deceleration phase corresponds to conditions at which inertia no longer contributes significantly to the force balance (2.3), it can be concluded that the time to dissipate the initial kinetic energy scales remarkably well with $D/U$ , independent of initial Reynolds number and density ratio. Retrospectively, this is not surprising, since no other length or velocities scales are involved in the problem.

It is now consequential to plot the drag coefficient evolution as a function of deceleration. This is done for the two spheres $\rho ^*=2.16$ and 3.26 in figure 10 for all instantaneous values within the deceleration phase, i.e. excluding the submersion and settling phases, assuming these phases with residual entry effects or with no inertial effects may exhibit a different behaviour. Although some irregularity and scatter are observed in this data, all graphs show a very similar pattern: first, increasing $C_{{D}}$ with increasing deceleration followed by decreasing $C_{{D}}$ with decreasing deceleration. The following is postulated from this data. Under initial deceleration at the beginning of the deceleration phase, the drag coefficient increases due to the earlier boundary-layer separation, as described above, thus increasing the deceleration. This is then a self-augmenting and perpetuating process, which continues until some limit is reached, close to the maximum drag coefficient. Thereafter, any perturbation causing the boundary layer to attach later will decrease the wake area, the drag coefficient and the deceleration; again a self-augmenting process, until the sphere reaches its minimum velocity at the end of the deceleration phase. As seen in this figure, the limiting or maximum drag coefficient is larger for the lighter sphere ( $\rho ^*=2.16$ ), presumably since inertia is more dominant for the heavier sphere ( $\rho ^*=3.26$ ). This then explains the rise and fall of the drag coefficient seen in the deceleration phase, for instance, shown exemplary in figure 8. For the spheres, $\rho ^*=6.08$ and 7.92 such curves are not obtained simply because the inertial force overwhelms any small variation of deceleration.

Figure 10. Dependence of instantaneous drag coefficient on instantaneous deceleration for times within the deceleration phase: (a,b,c) $\rho ^*=2.16$ ; (d,e,f) $\rho ^*=3.26$ . The arrows indicate an increasing time in the approximate interval $t^*=10{-}20$ . The legends in (a), (b) and (c) apply to (d), (e) and (f), respectively.

This self-augmenting behaviour is inherent to any instability and gives rise to hysteresis. This is evident in the curves of figure 10, where for the same value of deceleration, two values of the drag coefficient can be obtained, depending on whether the deceleration is increasing or decreasing. This is, therefore, a history effect arising from the boundary-layer separation adjusting to varying values of the sphere velocity; however, this history effect is more influential than the Boussinesq–Basset term, which only considers viscous drag and not pressure drag, as in this case.

To explain the above results, we now consider the flow around the sphere under conditions of strong acceleration/deceleration, particularly the behaviour of the boundary layer near its separation point. To observe the flow patterns around the sphere, a uniform coating of red dye was applied to the sphere before it was released into the quiescent liquid with a $Re_i$ of 15 700. Figure 11(a) shows the resulting flow visualization around the sphere. Under steady flow conditions (which is not the case in the deceleration phase), separation for a laminar boundary layer would be expected near an angle of 80 $^\circ$ –84 $^\circ$ . The expected velocity profile of the boundary layer in the vicinity of separation is shown schematically in figure 11(c), whereby the solid lines represent the boundary layer velocity profile at some instant in time. If the sphere experiences a strong acceleration, the outer flow would effectively increase in velocity, and over a short period of time, the separation point would move rearwards, decreasing the wake area. This is indicated in the figure with a dashed line marked ‘accelerating sphere’. With a deceleration of the sphere, the opposite trend would be expected, i.e. the separation would occur earlier, resulting in a larger wake area over which the base pressure would be exerted. In the figure, this is indicated with the dashed velocity profile labelled ‘decelerating sphere’. This would result in a larger base pressure area; hence, a higher drag, consistent with the observations of the present study throughout the deceleration phase.

Figure 11. Schematic interpretation of the origins of lift and drag through wake asymmetry. (a) Visualization of a free-falling sphere ( $\rho ^* = 2.16$ , $D = 10$ mm and $Re_i = 15700$ ), showing a change in trajectory and indicating an asymmetry of the left-hand (P $_{{L}}$ ) and right-hand (P $_{{R}}$ ) points of separation. (b) Schematic (and exaggerated) visualization of asymmetric wake area arising from asymmetric separation points, indicating that relative to the flow direction, a lift force component arises. (c) Schematic (and exaggerated) interpretation of how a deceleration or acceleration of the sphere might influence the boundary-layer velocity profile in the vicinity of a separation point.

It is apparent from the data in figure 10, that the drag coefficient of a decelerating sphere at some instantaneous Reynolds number can be significantly higher than for a sphere experiencing a steady flow at the same Reynolds number, usually of the order 0.5. Although there exist numerous studies of sphere drag under decelerating conditions, most of these are devoted to low Reynolds number flows (Temkin & Mehta Reference Temkin and Mehta1982; Liu et al. Reference Liu, Liu, Qu, Lin, Hu and Agarwal2018; Velazquez & Barrero-Gil Reference Velazquez and Barrero-Gill2024). Some studies of decelerating bluff bodies at similar Reynolds numbers have indicated that the drag coefficient can differ significantly from values for steady flow (Potvin et al. Reference Potvin, Peek and Brocato2003), while others remain rather inconclusive about how drag changes for decelerating bodies (Marchildon & Gauvin Reference Marchildon and Gauvin1979). Thus, at present it is difficult to find corroborating data to our experimental results.

Similarly, little quantitative literature exists on unsteady lift forces in this Reynolds number range, as confirmed in the exhaustive survey of spherical particles in unbounded flows given by Shi & Rzehak (Reference Shi and Rzehak2019). Nevertheless, it is clear that bifurcations and symmetry breaking in the wake of spheres occur already at relatively low Reynolds numbers and result in a random change of direction (Fabre et al. Reference Fabre, Auguste and Magnaudet2008). The magnitude of these lift forces in the present Reynolds number range has, to the authors’ knowledge, not been previously measured.

This transient state of the boundary layer during sphere deceleration is analogous to Stokes first problem in which the similarity variable $\eta =y/2\sqrt {\nu t}$ is used to describe the invoked velocity profile above a plate, suddenly accelerated to some finite velocity (Spurk & Aksel Reference Spurk and Aksel2007). Interpreting $y$ as the boundary-layer thickness in the present case, the time necessary for the diffusion of a strong change of outer flow velocity into the same dimensionless state near the sphere surface would scale with $y^2$ . If, for instance, one boundary layer were only half as thick as a second boundary layer, a change of outer velocity would be felt in $0.5^2$ , i.e. a quarter of the time. Given that the boundary-layer thickness over the sphere will scale approximately with Re $^{1/2}$ , the spheres with higher initial Reynolds numbers would exhibit a lower boundary-layer thickness. Hence, this should lead to an earlier and stronger influence on the separation point, resulting in a higher drag. This is consistently seen in the data of figure 10.

The spheres with a higher density ratio will have higher inertial forces for the same Reynolds number; hence, any wake variations will have a relatively lower effect on the overall sphere motion, resulting in less modulation of the computed drag coefficient using (2.3). The results shown in figure 10 confirm that spheres with higher density ratios exhibit lower maximum drag coefficients.

The settling phase is prominently characterized with widely varying values of lateral deviation $r^*$ , as shown in figure 4, but also evident in figure 8(b). During this phase, the velocity is very low, and inertia no longer plays a dominant role. Thus, the trajectory is now highly susceptible to even small lift forces arising from wake asymmetries. Thus, the trajectory curvatures become more frequent and pronounced. This leads to variations in drag coefficient, but now only seen as smaller variations near values typical of steady-state flow, i.e. $C_{{D}} \approx 0.5$ . Similarly, the lift coefficient reduces in magnitude, approaching values near zero. These variations, both in drag and lift, are still attributed to a fluctuating and non-symmetric wake structure, but its consequence is now much less significant.

Referring back to the fluid mechanic interpretation of the deceleration phase and in particular to the yet unexplained stability of the wake asymmetry orientation throughout this phase, it is evident that this stability no longer exists in the settling phase. Not only the dimensionless lateral displacement exhibits large fluctuations, but also the orientation angles $\psi$ and $\chi$ are continually changing. This is seen in figure 8, but is also evident in figure 4(d–i), in which similar results are shown for several other spheres and impact conditions. Recalling the physical meaning of the angles $\psi$ and $\chi$ given above in § 3.1, it is evident from these measurements that in this phase, after which the velocity has reached a minimum value, the vortex shedding and instantaneous wake and associated base pressure area exhibit much less stability in their orientation with respect to the flow direction. The continual variations of $\psi$ and $\chi$ , while small in magnitude, suggest a random change in pathline orientation, i.e. a random change in lift force orientation. This would appear to be a consequence of random asymmetry directions of the wake and/or vortex shedding. Given the low velocities and correspondingly low inertial forces, even small changes in wake orientation will have a larger effect on the trajectory in this settling phase.