A mathematical model for simulating cycling: applied to track cycling

A review of existing mathematical models for velodrome cycling suggests that cyclists and cycling coaches could benefit from an improved simulation tool. A continuous mathematical model for cycling has been developed that includes calculated slip and steering angles and, therefore, allows for resulting variation in rolling resistance. The model focuses on aspects that are particular, but not unique, to velodrome cycling but could be used for any cycling event. Validation of the model is provided by power meter, wheel speed and timing data obtained from two different studies and eight different athletes. The model is shown to predict the lap by lap performance of six elite female athletes to an average accuracy of 0.36% and the finishing times of two elite athletes competing in a 3-km individual pursuit track cycling event to an average accuracy of 0.20%. Possible reasons for these errors are presented. The impact of speed on steering input is discussed as an example application of the model.


Motivation
Cycling is a sport which lends itself to performance analysis. The relative ease of data collection means that competitive teams carry out much analysis and spend significant resources determining optimal choices of equipment, athlete or strategy. As power meters and other measurements have become more accurate, the desire for accurate mathematical models has grown. A review of existing predictive cycling models has revealed scope for improvement in the modelling of tyre forces in track cycling in particular.

Literature review
Since the release of Schoberer Rad Meßtechnik's (Jülich, Germany) first power meter in the early 1990s, mathematical models of cycling have become increasingly sophisticated.
Olds et al. [1] presented a comprehensive equation for power demand accounting for aerodynamic drag (including head winds), rolling resistance and equipment specifications. Olds et al. [2] improved this model, including the impact of drafting, crosswinds and the kinetic energy of limbs and wheels. This revised model was validated by testing 41 athletes over a 26-km time trial. Using measured power as an input the model predicted finishing time to an accuracy of 5%. Martin et al. [3] developed a model more complicated in its calculation of aerodynamic drag, with wheel drag varying with velocity, but otherwise less sophisticated. This model had an accuracy of 3% when validated for six test subjects. Basset et al. [4] derived a model to compare the several different world hour record attempts in the 1990s. This model included some limiting assumptions: i.e. equal groundspeed and air speed, and frontal surface area, a constant fraction of body surface area.
It was not until Martin et al. [5] and Lukes et al. [6] who presented models for velodrome cycling that cornering on a banked track was addressed. Martin et al. [5] considered the track as circular with constant radius and banking angle. In contrast, Lukes et al. [6] modelled the velodrome as two straights and two corners, albeit with no transition between the two. A similar approach was taken by Caddy et al. [7] in an investigation into the impact of cyclist posture on event performance. All only approximate the true shape of a velodrome. Lukes et al. [8] refined their approximation by splitting the velodrome geometry into eight sections rather than four. This improved model also included tyre scrubbing effects. It could predict finishing times with an accuracy of 2%. In a separate study investigating the aerodynamics of track cycling, Underwood [9,10] created the most accurate track cycling model to date. Using measured power and field-derived values for drag area, C d A, the model predicted elite individual pursuit finishing times with an average error of 0.42%.

Paper overview
Although four of the above models [4,5,[8][9][10] concern velodrome cycling, they all neglect some aerodynamic aspects and model the bends simplistically. To address these limitations, this paper describes a more complete mathematical model for all cycling events. It proposes a method of predicting the slip and steer angles necessary to navigate turns (allowing for banking), as well as allowing for rotating bodies in the system. The model is derived, a method of implementation is explained and then the results of two different validation studies are presented. The impact of speed on steering input is discussed as an example application of the model.

Model principles
The focus of this model is velodrome cycling; however, the equations derived are generic and could be applied to other cycle sports. The model is a significant extension of the velodrome-specific model presented in a previous study [7].
The model is a quasi-steady-state analysis that assumes instantaneous equilibrium of the cyclist at each time step but allows for changes in speed and configuration between time steps. At each time step, the cyclist is assumed to be following a path of known local curvature at constant speed, i.e. with all accelerations other than centripetal neglected. The rate of work done against dissipative forces is calculated based on this instantaneous equilibrium. Any difference between the cyclist's input mechanical work and that dissipated is attributed to changes of gravitational potential and/or kinetic energy (any change in the latter implying acceleration).
Forces are resolved both tangential and perpendicular to the cyclist's direction of motion at every time step. The lean angle, tyre slip angles and steering angle are calculated. The model assumes that the heading angle of the cyclist, χ, and the steering input, δ, are small and will, therefore, only be considered in determining the tyre slip angles (which are of a similar magnitude).

Governing equation
The governing equation for this model is an energy balance: Over a time period, δt, the available mechanical work is the product of the cyclist's input power, P in , and the efficiency, η, of the bicycle transmission. Drag forces dissipate much of this energy, E diss . The remaining power results in changes in the total kinetic, ΔT, and/or potential, ΔV, energies. Figure 1 shows forces on the cyclist viewed in a direction parallel to the ground surface and perpendicular to the direction of motion (see also Fig. 2, a view along the direction of motion). β is the banking angle of the ground surface; λ is the angle of vertical inclination of the cyclist's direction of motion and F λ is the consequent component of the  cyclist's weight acting in their direction of motion. F d is the aerodynamic drag force and ζ is the angle from the direction of motion through which it acts. F R and F N are rolling resistance and normal contact forces, respectively, where subscripts 1 or 2 refer to the front or rear tyre. F F is the propulsive force acting at the contact patch of the rear tyre. F x is the horizontal component of the centripetal, F c , and drag, F d , forces (Fig. 2). F y is the sum of the forces acting on the cyclist at an angle λ to the vertical and perpendicular to the direction of motion (Fig. 2). Both F x and F y have components acting in the plane of Fig. 1 (as shown). Ω is the angular velocity of the cyclist and bicycle. Figure 2 shows the forces acting on a cyclist that determine the angle of lean, θ. This view is of a plane normal to the direction of motion, i.e. at angle λ to the vertical. Within this plane, the weight of the cyclist and bike has component F w .

Cyclist dynamics
As the wheels navigate a bend with instantaneous radius of curvature R w , the cyclist's centre of gravity moves on a path with instantaneous radius R CG . The corresponding centripetal force, F c , acts in a direction perpendicular to both the direction of motion and the axis of rotation and thus at an angle, κ, to the horizontal. Taking moments about the wheel contact point, and using an iterative formula for θ can be derived: where m is the total mass of the cyclist and the bicycle, g is the gravitational acceleration, v CG is the velocity of the cyclist/bicycle centre of gravity and h CG the distance between the centre of gravity and the wheel/ground contact line.
The roll angle of the cyclist, φ, is given by

Aerodynamic drag
Aerodynamic drag can exceed 90% of a cyclist's resistance to motion [11,12]. In this model, the aerodynamic drag force is calculated via using the drag area, C d A, air density, ρ, velocity of the centre of drag, v d , and local air velocity, v air . v d/air is the velocity of the centre of drag relative to the air and is found using The drag force acts in the same direction as v d/air , at an angle ζ to the direction of motion and parallel to the ground surface (it is assumed that v air is in a plane parallel to the ground surface).

Slip, camber and steer angles
To determine the necessary steering input, the tyre loading and slip and camber angles must be calculated (Fig. 3a). The tyre loading depends on the overall equilibrium of the cyclist (see Figs. 1, 2, 3b). The resultant reaction force, P, acting through the tyres and its components can be derived by where F N and F S are the total normal contact force and side force, respectively. F N and F S can be separated into forces acting through front and rear tyres by taking moments about the front tyre contact patch, resolving in two directions and using the lengths a and b (see Fig. 3b).
Thus it can be found that Having determined the side force required at each tyre, it is possible to calculate the slip and camber angles by the use of which follows from Eqs. (14)- (18). The ratio of side to normal forces is also given by where i and i are slip and camber angles of the wheel in question, C i the cornering stiffness (/rad) and C i the camber stiffness (/rad) of the tyre.
Camber and slip angles for front and rear tyres are defined by [13] where ε is the steering rake angle. u n and u t are the normal and tangential components of the bicycle velocity with respect to the bicycle frame, see Fig. 3a and where χ is the heading angle. The relationship between ground steer angle, δ′ and steer angle, δ, is given by [13] sin 1 = sin + sin cos ,  29), it is now possible to obtain a function of the bicycle and track geometry, tyre coefficients and roll and heading angles that can be solved iteratively, by substituting Eq. (28) for the steer angle, δ:

Rolling resistance
The total rolling resistance is given by where C rr1 and C rr2 are coefficients of rolling resistance for front and rear tyre, respectively. These coefficients depend (28) tan � = cos cos − sin sin . (31) (32)

Potential energy
Work done against gravity is determined by changes in potential energy, V, which is equal to where the overall height of the centre of gravity, z, depends on both the varying height of the path of the wheels, h w , and the lean angle of the cyclist, θ.

Kinetic energy
The kinetic energy of the system is given by where m w and I w are the mass and moment of inertia of one bicycle wheel, v cw and ω w are the translational and angular velocities of the wheel, respectively, m c and I c are the mass and moment of inertia of the cyclist/bicycle, and v c and ω c are the translational and angular velocities of the cyclist/ bicycle. Rotational kinetic energy of the limbs, pedals and cranks is neglected; Olds [2] showed them to account for only 0.07% of total kinetic energy, compared to 2% for the wheels. To make use of Eq. (34), the different velocities must be determined from the geometry of the cyclist and bicycle (Fig. 4).
By assuming no longitudinal slip for the tyre/ground contact, the angular velocity of the wheel relative to the cyclist, ω w/c , can be found by where r is the outer radius of the tyre and the components are in the first, second and third directions (Fig. 2). The cyclist and bicycle frames rotate as a rigid body, with an angular velocity of The components of the wheels' and cyclist's angular velocities in the first, second and third directions are given by Due to the quasi-steady-state approach adopted for this model, the component of angular velocity in the third direction (i.e. dθ/dt, the rate of change of lean angle) is assumed to be zero at each instant.
Moments of inertia about the centre of gravity in the first, second and third directions (Fig. 2) are assumed to be principal moments of inertia for both wheel and cyclist; thus: With all terms in Eq. (34) determined, the total kinetic energy, T, can be defined in terms of v CG by where and

Numerical solution and implementation
A straightforward implementation is forward integration of the acceleration, a CG , over fixed time increments, δt. Forces and configuration are calculated assuming . 4 Showing velocities of the cyclist and bicycle in a plane parallel to the bicycle frame instantaneous steady-state cornering using the method above; both are assumed constant throughout the time step. a CG is calculated, and also assumed constant throughout the time step. This approach is a computationally efficient approximation that should be sufficiently accurate if δt remains small (i.e. less than 2 s).
The rate at which energy is lost to dissipative forces is the product of the force magnitudes and the corresponding velocities and can, therefore, be calculated by Differentiating Eq. (33) means that the rate of change of potential energy can be determined by Note that the quasi-steady-state approximation assumes that the angular velocity of roll, dθ/dt, is approximately zero at each time step; therefore, in the simplest implementation dθ/dt must be estimated via linear extrapolation from previous time steps.
If the input power P in of the cyclist is known, the power P T associated with a change in kinetic energy can be calculated by modifying Eq. (1) to If K is assumed to be approximately constant over a time step then P T can be calculated by The acceleration of the cyclist's centre of gravity is given by The velocity, displacement, configuration and forces acting on the cyclist at the beginning of the next time step can then be determined.
To avoid the discontinuity that arises from v CG equalling zero at the start of the initial time step, Eq. (48) can be modified to where G is the gear ratio and Q the starting torque.

Method
The model was implemented in Matlab (Mathworks, Cambridge, UK) and validated by two different methods. The investigation has been approved by the Cambridge University Engineering Department Research Ethics Committee. First, the tool was used to predict the lap times of six different elite female athletes cycling at approximately constant speed. These athletes each took part in three sub maximal efforts at the Manchester velodrome on three separate days. Throughout each effort, the athletes were asked to maintain a specified speed and to follow the 250-m datum line as closely as possible. The specified speed was varied for each session and athlete. In total, the six athletes completed 174 laps at speeds between 42 and 51 km/h. The recorded power data and measured athlete characteristics were then used to predict the athlete's performance throughout the 174 laps. All of the participants gave informed consent for their data to be used in this investigation.
Second, the tool was used to predict the finishing time of two elite female cyclists competing in the 3 km Individual Pursuit (3KIP) event at the 2017 UEC European Track Championships (ETC2017) in Berlin from the input power recorded for each cyclist during the same event.

Athlete power
Input power of the athletes was recorded using a power meter (Schoberer Rad Meßtechnik, Germany) which had been calibrated according to the manufacturer's instructions.

Atmospheric conditions
Air density, ρ, was calculated from local atmospheric conditions at the time using Teten's formulation [15]. Gravitational acceleration, g, was determined for the velodrome locations [16].

Track geometry and cyclist trajectory
Track geometry (banking angles, radii, inclinations) was required for two different velodromes. The Manchester (UK) velodrome geometry was found from a survey of the 250-m datum line using a TC403L total station (Leica Geosystems, Heerbrugg, Switzerland). The geometry of the Berlin velodrome was determined from a combination of expert knowledge and information given by the track designers: Schuermann Architects (Muenster, Germany). It was assumed that the wheels of the cyclists exactly followed the datum line. In contrast to other studies [4,5,7,8] the altitude, h w , of the datum line was allowed to vary.

Drag area and aerodynamic drag
Bulk airflow is caused by cyclists circling a velodrome. Throughout the validation process, this airflow was assumed to remain constant in both magnitude and direction. The magnitude of the airflow was the average of that measured during the validation session using an anemometer. The direction of the airflow was assumed tangential to the cyclist's motion, i.e. ζ = 0.
Due to lack of access to a wind tunnel, the cyclist's drag area C d A was derived from field testing. Each athlete underwent aerodynamic testing with identical equipment and maintaining the same position as that used in each effort but at an earlier date. The protocol outlined by Fitton et al. [17] was used in each instance. C d A was assumed constant throughout each effort.
Note that the centre of drag (the point through which the aerodynamic drag force acts), is assumed the same as the centre of gravity. The cyclists' positions, and the fact that their bodies typically account for 70% of their total aerodynamic drag [11], means that the centre of drag is very close to the centre of gravity.

Mass and inertia
In the first part of this validation, the mass of the cyclists and their equipment was measured before and after each session and the average used in the simulation. In the second part mass was measured only once, as close to the event as practical.
The model requires the centre of gravity location and moment of inertia of the cyclist. These inputs were determined by measuring an average-sized elite cyclist, who did not take part in either study but was part of the same team, from a high-definition photo and then modelling each limb, the torso and the head as separate ellipsoids. The same cyclist was weighed and each ellipsoid assigned a proportion of the cyclist's total mass typical for an average human [18]. Using two more photos of the same athlete in their cycling position, three-dimensional coordinates were assigned to each limb. With the limb dimensions, estimated masses and positions, it was possible to determine the centre of mass and the inertia of the cyclist and bicycle (without wheels), I C . The values for other athletes were determined by scaling for their relative physical characteristics. Wheels were assumed to be uniform discs to calculate inertia.

Efficiency of the bicycle
Sources of inefficiency on a bicycle include drivetrain, frame flexibility and wheel bearings. For this validation, a fixed mechanical efficiency, η, of 98% [19] has been assumed.

Tyre properties
Tyre properties C α , C γ and C rr were determined in a previous study [14]. C α and C γ were found to be constant and independent of tyre normal force, whereas C rr was found to be a function of both the loading and orientation (α and γ) of the tyre. C rr has, therefore, been modelled to increase non-linearly with both α and γ, and most significantly with the former.

Results
For part one of this validation, lap times were predicted with an average accuracy of 0.36%. The maximum and minimum errors of the model were 0.98% and 0.001%, respectively. The standard deviation of the errors was 0.22%. The prediction was less than the actual lap time in 86% of cases; a probable explanation is the cyclists' imperfect handling ability. Analysis of the measured wheel speed data revealed that the elite cyclists travelled 0.7% further than the track length.
In the second validation study, finishing times of the two events were predicted with an average accuracy of 0.20% and individual split times with an average absolute accuracy of 0.24% (Table 1). Figure 5 compares simulated and recorded wheel speed data for Athlete A. Again the simulation under predicted the finishing time of both athletes. An additional contributing factor may be the assumption of constant C d A throughout the event. In a 3KIP C d A may be greater at the start, as the athlete pedals out of the saddle, and at the end of the event, where the tiring cyclist may not maintain their position.
The error associated with the simulation's prediction of finishing times and, importantly, split times is lower than in any previous comparable study. The model described by Lukes et al. [8] was capable of predicting finishing time in a 4KIP to within 2% and individual split times (0-1 km, 1-2 km, 2-3 km, 3-4 km) with a slightly larger error than that. Underwood's [9,10] proposed model was able to predict finishing times for elite athletes competing in the 3KIP and 4KIP to within 0.42%. When investigating Underwood's model's prediction of split times (0-1 km, 1-2 km, 2-3 km), however, the available data suggest higher errors of approximately 2.5%.
Significant factors contributing to the accuracy of the model presented here include the consideration of rotational kinetic energy and varying tyre forces and the care taken in measuring inputs, particularly coefficients of aerodynamic drag and rolling resistance.

Example application
One novel aspect of this model is the capability to predict tyre slip angles and the necessary steering input, δ, to navigate a particular trajectory. Using the geometry of the Manchester velodrome datum line, the impact of speed on δ at the bend apex has been predicted for Athlete A (Fig. 6).
At low speeds, the model predicts a low δ despite the low lean angle, θ (Fig. 7). As speed increases δ also increases to a peak of ~ 1.7° at ~ 50 km/h. This approximately coincides with the speed at which roll angle, φ, equals zero. As speed and φ further increase δ is predicted to decrease.

Conclusions
A mathematical model of simulating cycling has been developed. The model includes aspects of particular relevance to velodrome cycling. Via two different validation studies, the accuracy of the model has been shown to surpass previous comparable models: errors in predicted lap times are consistently less than 0.36%. A key advantage of the model is the calculation of steer and tyre slip angles; this enables the rolling resistance to be predicted more accurately. This makes it possible, for example, to comment on the impact of handling ability and tyre choice on event performance.  Fig. 6 The impact of speed on the necessary steering input at the apex of the bend in the Manchester Velodrome Fig. 7 The impact of speed on the cyclist's lean and roll angle at the apex of the bend in the Manchester Velodrome