Curious dynamics of a Golf Ball Bounce

Introduction

If you have ever tried looking up any analysis of a ball bounce you would have most likely found either simple problems of a rigid ball bouncing off a (rigid) table, or more complicated problems regarding tennis balls, footballs or basketballs. Even if the ball is allowed to deform during the bounce, one common denominator of those studies is the rigid surface which the ball impacts.

Enter the game of golf. The golf ball is often imagined as a rock-hard object, but it is not quite the case (I encourage you to watch [5]). The ball becomes quite elastic during its contact with the club, but remains approximately rigid as it bounces off the course. What is unusual, compared to the previous approaches, is that the turf itself becomes compliant. Such setting (a rigid ball against a compliant surface) has not been studied extensively, giving us a room to explore.

To understand why this problem would be of anybody’s interest let us firstly understand what qualifies as a golf ball. There are six rules in the official governance regarding the golf ball [6]:

  1. General which states that a ball must be a golf ball in a traditional sense;
  2. Weight which specifies the upper limit for ball’s mass;
  3. Size which specifies the lower limit for ball’s diameter;
  4. Spherical symmetry which requires the ball to behave symmetrically;
  5. Initial Velocity which specifies the upper limit for ball’s initial velocity under given launch conditions;
  6. Overall distance which specifies the upper limit for ball’s overall distance (including all bounces and the roll) under given launch conditions.

 

Rules 2-5 can easily be verified by simple measurements, rule 1 may be less precise but still hardly ambiguous. What may cause a headache is the overall distance rule. How would one go about measuring it? A slightest change of the course incline or a blow of wind could change the ruling either way…

This is why a precise model for the bounce would be of benefit to the manufacturers as well as the governing bodies of golf. In this project we are after not only a simple model for the most traditional or standard bounce conditions. We are seeking a model which will explain a wide range of initial conditions as well as a number of phenomena which can be observed on a golf course. During the impact a ball can either slide or roll, it can bounce forwards but it is not unusual for it to bounce backwards too, and finally it can also change the orientation of its spin — these are just some of the behaviours we wish to account for in our modelling and which will remain at the centre of our attention throughout this project.

Current literature

Although not much has been done, some results are available. A starting point is often a simple model for the bounce of a rigid ball against a rigid surface, such as the one described in [2]. This usually assumes a simple Coulomb force acting as a friction between the ball and the surface, and the lift off is determined by a coefficient of restitution r, meaning that the ball will lift off when the ratio of its lift off normal velocity to the initial normal speed is r. In [4] Penner used such formulation to explain a golf ball bounce — according to his hypothesis multiple points of contact due to a compliant surface can be simplified to a single force acting at an angle — see Figure 1. In this case the analysis of the bounce focuses on finding the angle of incline and applying the rigid bounce theory at such angle. Given limited data that Penner had it has been concluded that this angle will be inversely proportional to the speed and angle on the ball’s inbound trajectory (spoiler alert: we shall not discuss this in detail, but our data does not seem to support this claim).

 

penner

Figure 1: An illustration of Penner’s hypothesis — multiple forces acting on the surface in contact (left) are equivalent to a single force acting at an angle (right).

A larger and more in-depth analysis, similar to the one we carry out here, has been carried out by Haake in [3]. The study focused on determining the physical properties which characterise a golf turf. Haake proposed some simple bounce models to explain some phenomena, however these have been supported by only a limited set of the collected data.

Finally, a number of models and experiments have been conducted by Rod Cross, such as the ones presented in [1]. These, however, rarely focus on explaining the visco-elastic behaviour of the turf and explain the observed phenomena using the coefficient of restitutions (both in the normal and tangential directions) model. There is a larger experimental validation of these models, but explanations behind them can rarely be expanded to a wide range of initial conditions.

There is therefore a wide gap in current literature between the modelling and experimental validation. The more precise or physically valid methods are not well supported by experimental data, whereas the models with some experimental validation are often focusing on a narrow set of conditions.

This is the void which we are now attempting to patch up. As mentioned before, our models should be robust, but we are also hoping to run a series or experimental campaigns which will validate them. In other words — our model should be driven by both data and physics…

Simple approaches at the modelling of the compliant behaviour of the ground

One of the simplest way to model visco-elastic properties is by the use of simple springs and dampers. Unsurprisingly, this also is a route we took. We present two simple mechanisms based on the Kelvin-Voigt formulation (where a spring and a damper are connected in parallel) and the Maxwell mechanism (where the spring and the damper are connected in series). For now we shall also make a simple assumption that the ball spins in the plane of its motion.

In both models the normal elasto-plastic behaviour is modelled using the Kelvin-Voigt formulation. The analysis of our data (which will be described in more detail in the next parts) showed that such formulation agreed better with the observed behaviour. This is further confirmed in [3], where it is also deemed that Maxwell-like behaviour in the normal direction is unphysical for the golf ball and turf interaction.

Kelvin-Voigt model

In this simple model we set two uncoupled mechanisms, one acting in the normal and the other in the horizontal directions. Both of these mechanisms will be modelled by a spring and damper connected in parallel — see Figure 2. It is then a straight forward exercise in the Lagrangian mechanics to describe the motion of the ball in contact with such surface, yielding our equations of motion:

    \begin{equation*} \begin{aligned} \ddot{x} +\frac{2d_x}{\varepsilon_x} \dot{x} +\frac{1}{\varepsilon_x^2}x &= \lambda_T,\\ \ddot{y} +\frac{2d_y}{\varepsilon_y} \dot{y} +\frac{1}{\varepsilon_y^2}y &= -g,\\ \frac{2}{5}R^2\,\dot{\omega}&=R \lambda_T, \end{aligned} \end{equation*}

 

kv_kv

Figure 2: Elasto-plastic model based on Kelvin-Voigt formulation.

where (x,y) is the position of the centre of mass of the ball, \omega is its angular speed and R is the radius of the ball. Here, we have manipulated our parameters so that they are reduced to just two parameters in each direction: \varepsilon_{x,y} are parameters that can be thought of as stiffness coefficient. Small values of \varepsilon will imply a stiffer surface in the respective direction. On the other hand d_{x,y} can be thought of as damping ratio — larger values of d will imply larger energy loss during the impact. There is also a more general understanding of the behaviour of these systems that arises from such arrangement of the parameters: with d>1 the explicit solutions to these differential equations will take a form of exponential functions which tend to a horizontal asymptote, whereas with 0<d<1 the solution will be an oscillatory transcendental function. We will therefore refer to the former as an overdamped mechanism, and the latter as an underdamped one.

We also introduced a tangential reaction \lambda_T which regulates the behaviour of the ball in contact with the surface. In our analysis we have assumed this to be a simple Coulomb friction force, that is \left|\lambda_T \right| = \mu \lambda_N, when the ball is slipping and \left|\lambda_T \right| < \mu \lambda_N when the ball is rolling on the surface, where \lambda_N are the forces acting in the normal direction. The parameter \mu (together with previously introduced \varepsilon_{x,y} and d_{x,y}) will then characterise the surface.

This very simplistic formulation gives us a wide range of possibilities in terms of possible behaviours it can model. Introduction of the pushing force in the horizontal direction, which is the main change from a simple rigid bounce formulation, gives us a wider spectrum of what can happen to the ball. Without a doubt the ball can bounce backward and change its spin, but at the same time it can still be in contact with the surface while slipping only.

Mathematically, this is a very simple initial value ODE problem. In theory, an analytic result could be obtained, but does not introduce much insight into the work. On the other hand we know that the solutions are linearly dependent on the initial conditions, in particular in the normal direction — the switch condition in the horizontal direction caused by the Coulomb law will introduce piecewise linearity only. This should prove to be a useful property when estimating the parameters of the model from the data.
One will see that working with particularly rigid surface (i.e. where \varepsilon_y\ll 1) leads to a singular perturbation problem. Let us consider this case (in vertical direction only). We introduce a new time variable \tau = \varepsilon_y^{-1} t so that with the rescaled time the vertical equation of motion becomes

    \begin{equation*} y'' + 2d_y y' + y =-\varepsilon_y^2 g, \end{equation*}

where the prime now denotes the derivative with respect to \tau. With the initial conditions

    \begin{equation*} y(0) = 0,\qquad \dot{y}(0)= \dot{y}_0, \end{equation*}

the above solves to

    \begin{equation*} y = \frac{\dot{y}_0 \e^{-\tau d_y}}{\sqrt{1-d_y^2}} \sin \left[ \tau \sqrt{1-d^2}\right] +\varepsilon_y^2 g \left[\cos \left[ \tau \sqrt{1-d^2}\right] + \frac{d_y}{\sqrt{1-d_y^2}} \sin \left[ \tau \sqrt{1-d^2}\right]-1 \right], \end{equation*}

where we have assumed d_y<1 (this appears to be the natural choice of the parameter given the aforementioned physical meaning of it and the fact that we wish our ball to lift off). Let us now focus on the limit case where \varepsilon_y \rightarrow 0. We note that in this setting effects of the gravity are neglected, and this indeed agrees with the rigid bounce theory. Furthermore, by solving the lift off condition y''(\tau) =0, we obtain an explicit formulation of the time of lift off

    \begin{equation*}\label{eq:lift_off_time_kv} \tau_F = \frac{\arctan\left[\frac{2\,d_y\,\sqrt{1-d_y^2}}{2\,d_y^2-1}\right] + k\pi}{\sqrt{1-d_y^2}}, \end{equation*}

where

    \begin{equation*} k=\begin{cases} 1& \mbox{ if }\quad 2d^2-1<0,\\ 0& \mbox{ if }\quad 2d^2-1>0, \end{cases} \end{equation*}

which assures that the `correct’ zero of our time function is used. Note now that the lift off time no longer depends on the initial velocity, which is a result of the linear dependence of the solution on the initial condition.

Finally, the above formulation allows us to find the coefficient of restitution
and this is

    \begin{equation*}\label{eq:limit_cor} r = -\frac{\dot{y}(\tau_F)}{\dot{y}_0} = \e^{-\tau_F d_y}, \end{equation*}

with \tau_F as previously defined. This indeed is a constant for a given surface, as the only parameter in the expression is d_y, which is the property of the ground. In other words, in the limit \varepsilon_y\rightarrow 0, implying a rigid ground, we obtain a constant coefficient of restitution, which agrees with the theory of the rigid bounce.

Kelvin-Voigt & Maxwell model

 

kv_m

Figure 3: Elasto-plastic model based on Kelvin-Voigt and Maxwell formulations.

A similar, but ever so slightly different way of describing the same situation is by placing the spring and the damper in the horizontal direction to be placed in series (see Figure 3).  With such change we increase the plasticity of the ground in the horizontal direction — the damper responds only to a push force. In this setting we have a force from the movement of the ball pushing it in the positive horizontal direction, but nothing pushing it backwards. This is a change when compared with the previous model, where the push was exerted by the spring connected to the damper in parallel.

 

The energy methods again allow us for the formulation of the equations of motion

    \begin{equation*} \begin{aligned} \dddot{x} +\frac{2d_x}{\varepsilon_x} \ddot{x} +\frac{1}{\varepsilon_x^2}\dot{x} &= \dot{\lambda_T} +\frac{2d_x}{\varepsilon_x} \lambda_T,\\ \ddot{y} +\frac{2d_y}{\varepsilon_y} \dot{y} +\frac{1}{\varepsilon_y^2}y &= -g,\\ \frac{2}{5}R^2\,\ddot{\theta}&=R \lambda_T. \end{aligned} \end{equation*}

The good news is that the equation of motions do not differ that much from our previously studied Kelvin-Voigt system. We still have the parameters \varepsilon and d which correspond to the stiffness and damping ratio. With the addition of the Coulomb friction law governing the behaviour of \lambda_T we still have a piecewise-smooth system of linear ODEs. In the limit \varepsilon \rightarrow 0 we recover the rigid bounce theory in the vertical direction (when solving the singular perturbation problem). In terms of the analysis of the equations — these are not any more difficult that what we have studied before nor do they introduce a particular insight to the behaviour of the ball.

What remains for us to do with these models is to see how well do they fit the data…

Experimental data

In October 2019 we carried out our first experimental campaign. In this simple setting golf balls were fired from a modified baseball canon and bounced off a plank of wood with a layer of artificial grass on top of it. We selected 55 initial conditions, ranging over different speeds, spins and incline angles. Each choice of initial condition has been repeated 6 times and the bounces were recorded using a high-speed camera (at 10,000 frames per second).

The resulting 330 recordings were analysed using image analysis methods to identify the velocity of the ball and its spin before and after the bounce. As the ball was fired from the launcher by two vertically placed wheels the assumption that the ball spins in the plane of its motion appeared to hold in the experimental setting.

Rigid bounce theory

A natural place to start would be the simple theory. We therefore attempt a fitting of two parameters that define a rigid surface in the rigid bounce theory — these parameters are the coefficient of restitution r and the Coulomb’s coefficient of friction \mu. The problem with the theory is that the model is a piecewise linear one, where the switch place depends on the parameters themselves. We thus implemented an iterative optimisation algorithm, in which the switch condition is updated after each iteration and new estimates for r and \mu are made. We determine the goodness of fit by implementing the least squares method

    \begin{equation*} f=\sum_{i} \left|\bm{P}_F^{i} - \bm{p}_F^{i}\right|^2, \end{equation*}

where \bm{P}_F = [ \dot{X}_F, \dot{Y}_F, \Omega_F], that is a vector of measured velocities (horizontal, vertical and angular) at lift off, and \bm{p}_F is the vector of corresponding velocities as predicted by the model. One should also mention that an appropriate scaling (taking radius of the ball to be of unit length) is implemented, so that the sizes of each of the vector components are comparable.

Unsurprisingly, we were able to arrive at some optimal values of r and \mu, but the estimate gave a whooping 50\% average relative error value! This should not be too unexpected though, one is unlikely going to be able to model an artificial grass surface as a rigid one!

Simple elasto-plastic models

Let us therefore return to our presented visco-elastic models, which were based on simple spring-and-damper mechanisms. Each of these included five parameters, but the models included the Coulomb friction law, which dictates a non-smooth switch between the slipping and rolling behaviours. For this precise reason, running an optimisation regime for the parameter fitting is unstable and highly sensitive to initial conditions. A structured approach is one of the methods that is currently being explored in this project, but we have instead estimated the parameters using a Monte Carlo approach. For each of the models 5,000 samples of parameters were selected at random from pre-defined domains. These domains have been established based on our physical assumptions, relating to the average time the bounce took, the distance the ball travelled and certain other behaviours it exhibited. For the best estimates the algorithm then looked in the vicinity of the given approximation to find a local minimum of the objective function.

With such conditions, best parameter estimates yielded an average error of 28.5 \% for the fit into the Kelvin-Voigt model and 29.4\% for the Kelvin-Voigt & Maxwell model. Although this is a significant improvement from the rigid bounce model, the result is far from a desired one. We also noticed certain inconsistencies amongst the best estimates — in each of the models the set of 10 best estimates included both overdamped (d_x>1) and the underdamped (d_x<1) cases in the horizontal direction. Given that these present significantly different behaviours of the mechanism itself this is an inconsistency that is yet to be explored.

Linear fits

So far we have been exploring the goodness of previously discussed models, with the models being our starting point. Let us now turn this around and let’s try a reverse approach, where we shall try to explore the data and discuss the possibilities arising from there.

The analysis up to this point, both technical and that of available data, suggests a significant change in the behaviour depending on the slipping vs. rolling scenario. And so the first and simplest approach would be to divide the data into two sets. In one set we shall keep all the data where the lift off happens with a positive velocity of the supposed contact point (therefore imitating a possible slipping scenario) and the second set will be the one where the said velocity negative (therefore suggesting the ball entered the rolling motion). One can easily fit a linear transformation into these sets, so that \bm{p}_F = A \bm{p}_0 — that is the lift off velocities are a simple linear transform of the initial values. Note how we omit any offset here — this is assuring a vague physical validity (a ball with no velocities before bounce will not suddenly lift off with some velocities after the `bounce’). Such fit (performed separately for the previously discussed two sets) gives an average error of 12\% — the lowest so far! The problem however is that the transforms are discontinuous between the two sets.

And so, let us now turn to a linear fit, with two constraints imposed:

  • a ball with zero velocity will not lift off (origin is mapped to the origin);
  • the resulting linear transform will be piecewise-smooth and continuous across the switch condition.

That is, we are looking for a transform of the form

    \begin{equation*} \bm{p}_F = \left\{ \begin{matrix} A \bm{p}_0 + \bm{d}_1 & \mbox{if } \bm{w} \cdot \bm{p}_0 + v >0\\ C \bm{p}_0 + \bm{d}_2 & \mbox{if } \bm{w} \cdot \bm{p}_0 +v <0\\ \end{matrix} \right. \end{equation*}

where A, C \in \mathbb{R}^{3\times 3}, \, \bm{d}_1, \bm{d}_2, \bm{w} \in \mathbb{R}^3,\, v\in \mathbb{R}, such that either \bm{d}_1 = \bm{0} or \bm{d}_2 = \bm{0} and A \bm{p}_0 + \bm{d}_1 = C \bm{p}_0 + \bm{d}_2 when \bm{w} \cdot \bm{p}_0 +v =0.

The constraints allow us to reform the problem, so that the map we look for takes the form

    \begin{equation*}\label{eq:linear_black_box} \bm{p}_F = \begin{cases} A\,\bm{p}_0 & \mbox{if} \quad \boldsymbol{w} \cdot \bm{p}_0 +v \gtrless 0,\\ (A -\bm{k} \otimes \bm{w})\,\bm{p}_0 -v \bm{k} \otimes \bm{w}\,\bm{i}& \mbox{if} \quad \bm{w} \cdot\bm{p}_0 +v \lessgtr 0, \end{cases} \end{equation*}

where we have assumed the first element of the vector \bm{w} to be w_1 = 1 for simplicity, and \bm{i} = [1,0,0]. This means that we now have a linear transformation with 15 parameters that need to be fitted within the optimisation regime. The problem is, once again, coming from the fact that the switch condition includes the parameters themselves, which in terms of parameter fitting creates steep gradient and thus badly fitted results.

We approach this problem by smoothing out the switch condition with a smooth sigmoid function and another iterative regime, which increases the steepens of the switch condition at each step and re-estimates the parameters.

The resulting estimation gives us the average error of 15\% — this is worse than a simple linear fit, but at the same time we have gained physical validity. The question now remains about all the other physical principles underlying such linear fit.

So where is it all heading?

And so one cannot finish the consideration of the topic now without a general wonder — where is it all leading us? The work is clearly incomplete — we still lack a model that would explain all the phenomena and match the data. But our intuition and the experimental results don’t give a clear answer as to where we will go next.

Physical intuition pushes us more and more into the world of non-linearity. A ball that is fully underground will behave very differently to the one that simply slides on the surface. We expect that the forces acting on the ball will depend on the surface area of the ball in contact with the ground. There are some models that we are hoping to explore soon (namely the ones previously proposed by Haake in [3] as well as the Hertzian contact models), however this type of behaviour unavoidably leads to non-linear analysis. Is that a bad thing? Well, it is given our data! Linear fits currently outperform any attempts to fit parameters into the models. If we try focusing on simply exploiting the data currently available — nothing should push us into the non-linear domain. Furthermore, our linear data has not yet been fully explored — we are still lacking a good explanation for the fits that we achieve…

But we must remember one thing — our data has been collected in a rather simplified setting. At the end of the day, the game of golf hardly ever sees a ball bouncing off a slab of wood with a layer of artificial grass on top of it. And although such experiential setting was highly beneficial for the early stages of the development of our models, we are now ready to collect and analyse the data on the bounce off a real golf turf. The experimental setting for such campaign is a lot more elaborate (once a piece of turf has been hit with a ball it is damaged and hence the setting needs to be moved) and hence the collection of such data will be a longer process. But we are hoping that it will shine some light on the direction in which the modelling can progress…

Acknowledgements

This work is partially funded and developed in cooperation with The R&A Ltd. We also acknowledge further funding from EPSRC.

References

[1] R. Cross. Backward bounce of a spinning ball. European Journal of Physics, 39(4):045007, 2018.
[2] C. B. Daish. The physics of ball games. Hodder and Stoughton, London, 1981.
[3] S. J. Haake. An apparatus for measuring the physical properties of golf turf and their application in the field. Ph.D. Thesis, The University of Aston in Birmingham, 1989.
[4] A. R. Penner. The run of a golf ball. Canadian Journal of Physics, 80(8):931{940, 2002.
[5] SmarterEveryDay. How Hard Can You Hit a Golf Ball? https://www.youtube.com/watch?v=JT0wx27J9xs, 2019. [Online; accessed 22-February-2021].
[6] The R & A and USGA. The Equipment Rules. R & A Rules Limited and The United States Golf Association, 2000.


This article was written by Stanisław Biber, an Engineering Mathematics PhD candidate at the University of Bristol.  For any questions relating to this work, please contact s.biber@bristol.ac.uk.