Finite element simulation of pole vaulting

Pole vaulting is one of the most spectacular disciplines in athletics. The evolution of world record heights is strongly influenced by the development of advanced poles and subsequent materials. Employing advanced, load-adjusted composites has resulted in a steady increase of the world record height. This study provides a framework for finite element simulations of pole vaulting with focus on the initial and boundary conditions as well as finite element choices. The influence of the pole bending stiffness on the achievable height is systematically simulated. Higher effective bending stiffness leads to higher pole vaulting heights. However, if a certain stiffness is exceeded, the vaulter will not be able to bend the pole enough which leads to failed attempts.


Introduction
Pole vaulting is an astonishing sport where athletes use poles as an auxiliary equipment to clear the bar at heights of 5-6 m and more. Risen from a purpose driven technique, e.g., overcoming castle walls during siege, crossing irrigation ditches [1] and transport [2], pole vaulting has become one of the most spectacular and high-tech disciplines at the Olympic games. It is a dynamic, elegant and fascinating sport.
To clear the bar, pole vaulters need to transfer the kinetic energy of the approach into potential energy via the pole. Two factors determine the success of a vault: (1) physical abilities of the athlete (speed, strength and springiness) and (2) properties of the pole (low weight and ability to store and return energy in the most efficient way).
In times where the physical abilities of the athletes reach their limits, the pole and, thus, the material it is made of, becomes a key factor. Figure 1 depicts the development of the world record in pole vaulting over time.
In the early 20th century, bamboo poles were utilized which are light but also very stiff and, thus, not able to transfer the kinetic energy from the approach in an efficient manner. In the 1940s, aluminum and steel led to new world records-however, these materials show high stiffness as well and therefore did not significantly improve the vaulting height. In the beginning of the 1960s, lightweight fiber-reinforced plastics (with lower bending stiffness) became increasingly popular. They allow large elastic deflections. As a consequence, they are able to transfer more energy resulting in a more pronounced catapult effect. In the following years, the world record steadily rose due to improved fiber arrangements and vaulting techniques. In the 1990s, a new plateau of nearly 6.15 m was reached. This indicates that human physiology and current materials may have reached a limit and enhancements of the pole's material are necessary to further improve the world record. The aim of this work is to provide a finite element (FE) model allowing for simulation aided pole design.

Model specifications
FE simulations of the pole vaulting process were conducted with the commercial FE code ABAQUS 6.14. The initial boundary value problem of the dynamic equation of motion was solved with an implicit time integration method under consideration of finite deformations. Exemplary FE models of the pole and the vaulter were obtained and generated in ABAQUS/CAE. To capture the kinematics of the pole vault procedure, a setting was employed that couples pole and vaulter allowing for a relative motion.

Finite element model of the vaulter
A comprehensive overview of the mechanics of pole vaulting can be found in Frère et al. [4] describing the process in four steps: run-up, take-off, pole bending and pole straightening. A crucial part in the jumping process simulation is the description of the vaulter. Muscle work of the vaulter increases the performance [5,6]. In addition, the moment exerted on the pole by the vaulter influences the vaulting performance [7]-elite vaulters bend the pole such that its effective length 1 relative to the original length is reduced by ca. 30% [8]. As stated by Ekevad and Lundberg [9], a modeling approach representing the vaulter by a point mass with a fixed position relative to the pole is not sufficient. The complex motion in combination with muscle work needs to be considered.
Hubbard [10] applied a model of the vaulter consisting of three rigid segments representing different parts of the body. Ekevad and Lundberg [11] extended this approach to six segments in a 2D setting, sequentially connected by joints: forearm, upper arm, head, trunk, thigh and shank (cf. Fig. 2).
The two arms and legs were represented by one equivalent segment. A summary of the properties of the segments is shown in Table 1. In this paper, we extend the approach of Ekevad and Lundberg [11] by additionally accounting for a joint at the neck.
An exemplary motion of a vaulter was obtained by analyzing videos of a German elite pole vaulter 2 . Specific frames of the video were taken to represent different phases of the vault process. With the aid of CAD software several measurements were taken for each phase to describe the position of the vaulter's body segments. In accordance with the planar motion of the segments model, a side view was considered in the videos 3 . The motion was described relative to the upper pole tip. Thus, it can be used in simulations with poles of different lengths. To do so, three coordinate systems were introduced as shown in Fig. 2. First, the global coordinate system x-y was fixed and located at the lower tip of the pole. Second, a relative coordinate system x 1 -y 1 was established at the upper tip of the pole. In addition, a relative coordinate system x 2 -y 2 was utilized with an origin coinciding with the previous one. This coordinate system rotates to maintain the x 2 -axis tangential to the pole. The measured quantities are: -Angle u between ground and a line connecting the tips of the pole, -Height h p of upper pole tip, -Angle h 0 of the rotated coordinate system x 2 -y 2 , -Angle h 1 between the coordinate system x 2 -y 2 and segment A, and -Angles h 2 to h 6 between the segments A to F.
The angles h i between the body segments of the reference pole vault over time are plotted in Fig. 3. A linear interpolation was performed between the discrete values.
At a second stage, the body was reduced to a point mass (center of gravity) with m = 80 kg. Thus, the significant effort to control the individual segments of the vaulter with muscle torques is avoided allowing more computational time on material modeling. The position of the mass center x m in the coordinate system x 2 -y 2 was calculated from the measured angles in each phase by where x i is the position of the mass center of the ith segment of the body and m i is the mass of the ith segment. The point mass moves relative to the pole in coordinate system x 2 -y 2 applying ABAQUS' connector elements (cf. Sect. 2.3). The different phases of the vault and the position of the mass center (blue square) are illustrated in Fig. 4. We employed two variables during the pole vault to trigger the relative motion of the point mass: the pole angle u from take-off to the instance of maximum pole bending, and the relative height h m;rel of the point mass from maximum pole bending till pole release 4 . This is due to motion of the point mass relative to the pole that cannot be implemented as time based. It depends on properties of the pole such as the pole's length and stiffness-these properties, however, are to be varied during the simulations.

Finite element model of the pole
Modern, elite vaulting poles are manufactured out of lightweight materials consisting mainly of glass fiber-reinforced plastics. The fiberglass pole is rolled on a mandrel and, subsequently, wrapped to stiffen the pole. The perfect   4 Ekevad and Lundberg [11] control the vaulter motion by the angle u. This leads to instabilities from the instance of maximum pole bending to pole release where the vaulter movement changes, based on the almost constant angle u. Ekevad and Lundberg [11] tried to overcome this drawback by calculating an average rate of the change of the angle u in the first half and applying this for the second half. Their approach is non-physical: the stretching of the pole is not coupled to the relative motion of the vaulter.
pole differs for athletes, depending on their physical properties, abilities and vaulting technique-and plays an essential role in the vaulter's performance. As a first step, an isotropic hyperelastic material model of Neo-Hookean type was applied with the strain-energy function where l is the shear modulus, I 1 ¼ trðbÞ is the first invariant ofb,b ¼ J À2=3 b is the distortional component of the left Cauchy-Green strain tensor, j is the bulk modulus and J is the determinant of the deformation gradient F. Neo-Hookean hyperelasticity renders linear elasticity as long as it is valid but also allows for nonlinear behavior which occurs due to the large deflection of the pole. Therefore, it is the more general choice. The material Fig. 4 Positions of the vaulter's body segments and mass center (blue square) during different phases of the pole vault.
Coordinate system x 1 -y 1 is chosen for the presentation (color figure online) density q and the Poisson's ratio m were calculated with the rule of mixtures (following [13]) for a fiber volume fraction of 50% as in Davis and Kukureka [14]. The modulus of elasticity was obtained via a flex test simulation 5 . This kind of test is used in the pole industry to classify the stiffness of a pole. The pole is pin-supported at both ends and loaded with a mass of 22.7 kg in the center. The deflection in cm gives the flex number (a relative stiffness number) of the pole. The assumed geometry of the equipment, namely the pole, is shown in Fig. 5. Table 2 summarizes the chosen properties. In the study on the pole stiffness (see Sect. 3), the elastic modulus is varied in the neighborhood of the previously calculated value of the flex test.
We discretized the pole with continuum elements as well as beam elements and compared the outcome. Continuum elements can represent arbitrary structures and bodies and can resolve complex three-dimensional stress states. Further, they allow for incorporating user defined material laws which might be desired when accounting for microstructure, for example. Structural elements such as beams rely on assumptions like slenderness simplifying the model to one predominant dimension. This leads to a reduced number of degrees of freedom and therefore, higher computational efficiency. ABAQUS' Timoshenko beam elements capture large axial strains [15].

Coupling of pole and vaulter
The point mass needs to be connected to the pole such that it can move relatively to the upper pole tip according to the vaulter's motion in different phases of the vault. Moreover, the connection fails once the pole is stretched and the vaulter would release it to clear the bar. An IAAF requirement is that the vaulter must hold the pole in the grip area (no higher than 0.1524 m from the top of any pole or no lower than 0.4572 m from the top of the pole).
ABAQUS' connector elements allow constraints involving relative motion of the connected parts via Lagrange multipliers as additional solution variables and also failure of the connection [15]. Through the connector element, a relative position of the point mass was given as a constraint on the upper tip of the pole. Moreover, the connector elements were employed as a sensor to measure the position ½x p ; y p T of the upper pole tip and the height h m of the mass center. A user subroutine UAMP was used in each increment of the calculation to call the quantities measured by the sensors and enforce relative displacement components [16].
In the UAMP routine, the pole angle u and the relative height h m;rel of the mass center were calculated via As the sensor measurements are from the beginning of the time increment, the corresponding values of u and h m;rel at the end of the time increment were extrapolated. Subsequently, the relative position of the point mass in coordinate system x 2 -y 2 in the next time step was determined as:

Initial and boundary conditions of pole vaulting
Regarding the boundary conditions of the simulation, the pole was supported at the lower end fixing all translational degrees of freedom. This represents the contact to the planting box. Furthermore, a gravitational acceleration of g ¼ 9:81 m s À2 was assumed. The choice of the initial conditions is not trivial and has to be done carefully. Frequently reported in the literature, is the employment of only the take-off velocity v 0 on the vaulter (see e.g., [9,11,18]), while no initial velocity was applied to the pole. The initial velocity of the vaulter is split into two components: one in the horizontal direction due to the approach speed and one in the vertical direction due to the jump. The components are calculated with the take-off angle a. Table 3 summarizes exemplary values for initial conditions according to Linthorne [18].
Due to the 'hard contact' boundary condition at the lower end, such an initial condition results in artificial oscillations of the pole with large amplitudes. All modes of the pole, which is initially at rest, are excited by the impulse. To circumvent these artificial oscillations, we started the simulation at the instance right after the pole was planted and when it started to bend as the vaulter's mass compressed it. Then, in addition to the initial velocity of the point mass, also a velocity profile v 0;pole on the pole was necessary that induced the bending. This profile was obtained from a preliminary simulation.
In this preliminary simulation, only the take-off velocity v 0 of the point mass was applied as an initial condition and the density of the pole was artificially reduced by six orders of magnitude. Therefore, the oscillations still occurred in the beginning of the simulations but the amplitudes were small since the pole had nearly no inertia associated with it. In the following, the oscillations decayed quickly (approximately 0.01 s) and a homogeneous velocity profile developed in the pole (cf. Fig. 6).
The profile has velocity components in the horizontal and vertical directions and is not necessarily perpendicular to the pole. This profile was then extracted and employed as the initial condition for the simulation. The velocity profile's shape resembles the bending of the pole in its first mode. Thus, an initial velocity profile initiates the real deformation. In the pole vaulting simulation, the profile was approximated with a third order polynomial along the pole length L requiring: 1. Zero velocity at the lower end, 2. Velocity at the upper tip equals the preliminary simulation, 3. Maximum velocity node position matches the preliminary simulation, and 4. Maximum velocity value equals the preliminary simulation.
3 Results and discussion IAAF [19] demands the crossbar, that needs to be cleared, to be straight above the planting box or up to 0.8 m behind it. Thus, in the x-y coordinate system from Fig. 2, the vault height h max should be attained at an x-position in the range of 0.0-0.8 m. Therefore, pole vault attempts, where the x-position of the vaulter at the vault height h max is negative, are considered to be a failed trial. Figure 8 depicts the maximum heights h max of the point mass depending on the effective bending stiffness EI eff . Valid trials are indicated with a green square whereas failed trials are displayed with a red triangle. The general trend of a higher effective bending stiffness leading to an increased height is observable. However, if the pole is too stiff, the vaulter is not able to bend the pole in such a way that he is catapulted forward-leading to a failed approach. The x-position at the height is illustrated in addition to signify valid and failed trials. The maximum height that is achieved in a valid trial is h max ¼ 7:58 m with an effective bending stiffness of the pole of EI eff ¼ 2218 N m 2 , corresponding to an elastic modulus of E ¼ 36:52 GPa. Considering that only the vaulter's body, but not his/her center of gravity, needs to pass the bar, slightly higher heights are possible. External circumstances as well as the fact that an athlete is never a perfect machine will decrease the realistic value. These results correspond well with world records. However, since dissipative effects are neglected, the computed height is slightly larger. The kinetic energy dissipated in the form of heat through the vaulter's muscles significantly influences the maximum height h max [18]. In  addition, the vaulter dissipates energy with impact when the pole is planted. The impact forces increase with higher pole stiffness. Since the simulation starts after the pole planting, this dissipation is as well neglected which especially influences the simulation outcome with poles of higher stiffness. Moreover, energy is dissipated due to friction at the planting box and the viscoelastic properties of the pole [20]. Further no force and momentum restrictions were considered for the relative movement of the mass point, meaning limits of human muscular strength are neglected at this stage. Essentially, computational simulations can compare different poles (i.e., different materials, lengths) as well as support the development of new poles. The velocity of the point mass is displayed in Fig. 9. The velocity v x in the horizontal direction decreases in the beginning as the kinetic energy of the approach is transferred into strain energy by deforming the pole. Subsequently, the velocity v y in the vertical direction increases when the pole recovers to catapult the vaulter. After releasing the pole, the horizontal velocity v x is constant since no forces act on the point mass in that direction. Contrary, the upward vertical velocity v y of the vaulter decreases linearly under the effect of gravity. When the vertical velocity v y decreases to zero, the maximum height h max and accordingly the maximum potential energy is reached. The energy transformations correspond well to the experimentally based descriptions of Dillmann and Nelson [21]. The simulated velocities correspond well in trend and magnitude to the measured ones of Frère et al. [4] and Angulo-Kinzler et al. [22].

Finite element choices: continuum vs. beam elements
The aim of this contribution was to allow for more detailed material models taking into account the heterogeneity of fiber-reinforced plastics. If the pole's heterogeneous microstructure is to be taken into account, continuum elements 6 need to be applied. They are computationally more costly, but distributions of field variables can be identified precisely and they allow to account for a sophisticated microstructured material model at a later stage. If the underlying pole material is treated as a homogeneous continuum, the pole can be discretized with beam elements 7 . The use of beam elements is suitable since the pole has only one dominant dimension whereas the other two are relatively small. The advantage of this approach is its computational cheapness allowing fast model studies. Thus, beam elements were used for the simulations in the effective bending stiffness study performed in this work and the microstructure was neglected. For reasons of comparison, the performance of beam vs. continuum elements is illustrated exemplarily for one case in Fig. 9.

Conclusion
This study provides a framework for FE simulations of pole vaulting, suitable for investigations on the pole's material. The vaulter was accounted for as a point mass to The initial velocity profile is computed in a preliminary simulation 6 ABAQUS' C3D8R continuum elements (8 nodes with 3 displacement degrees of freedom per node and linear shape functions) were applied in this case [15]. These elements possess one Gaussian integration point. Reduced integration was used to avoid shear locking. For such elements, artificial hourglass deformation modes may occur. The occurrence of these modes was controlled via the artificial strain energy that has to be small compared to the strain or internal energy. 7 ABAQUS' B21 beam elements were employed which are located in 2D space and use linear shape functions [15]. They possess two nodes with three degrees of freedom each (two displacements and one rotation) and two integration points. Timoshenko beam formulation was utilized.
save computational effort in control that was necessary in other studies. The motion of the vaulter relative to the pole was implemented with an algorithm based on a reference vault of an elite pole vaulter. The focus was on the proper choice of boundary and initial conditions to avoid unphysical oscillations. An initial velocity profile of the pole was determined by running a pre-simulation with reduced pole weight and overcomes unnatural simulation outcomes.
Beam elements as well as continuum elements can be used for the discretization of the pole yielding similar results. Continuum elements give a more detailed insight in distributions of field variables such as stresses and strains. Furthermore, they provide the option to study and incorporate complex computationally expensive material models considering the micromechanical structure of the pole.
The simulation results resemble measurements of the pole vault process in trend and magnitude. Additionally, we examined the influence of the effective bending stiffness of the pole on the achievable height in pole vaulting using an isotropic Neo-Hookean (continuum elements)/linear elastic (beam elements) material model and varying the elastic modulus. We observed that up to a certain stiffness, higher effective bending stiffness leads to higher pole vaulting heights. However, overly stiff poles would not bend enough to catapult the vaulter forward leading to failed vault attempts. This naturally occurring phenomenon is well reflected in our model.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.