*Nicola Bertin – nicola.bertin.5@studenti.unipd.it – Degree in Mechanical Engineering*

**Introduction**

The self-stability of bicycle and two-wheeled vehicles in general, is the results of many effects and parameters. In particular the main effects that made a bicycle self-stable are: position of masses (and inertia), trail and gyroscopic effects, taken together.

At “Delft University of Technology” a special bicycle named “TMS (Two mass skate) bike” was studied, a bicycle that is different from a standard everyday bike, in order to understand how these effect work together.

**Objectives**

The aim of this project is to recreate and study in Adams View environment the behavior of the real TMS bicycle.

The analyses which are carried out are:

– determination of cornering stiffness of wheels

– bahavior at different speed

– frequency analysis

– yaw rate and lean rate impulse behavior

**Procedure**

A paper that discusses previous work was downloaded, and the important parameters such masses and geometries and inertia values were extrapolated.

The entire 3d model was downloaded following the links in the paper http://bicycle.tudelft.nl/stablebicycle/: it was in Solidworks file format and was created and used from the same university for their analyses. The model was exported in Parasolid file format and imported in Adams.

The contact wheels-ground was handled with a point-to-plane formulation, with equations that simulate in a basic manner the behavior of a rigid tyre (in this case made of metal).

After that, various analyses were carried out to understand the parameters of such kind of tyre, and if the equations were sufficient to describe the behavior of the bike.

Finally, with other analyses, the experimental data were compared with those of the analyzed model.

**Geometry**

In the paper there are the dimensions, geometric properties and inertia parameters of the structure:

The 3d model was downloaded following the link in http://bicycle.tudelft.nl/stablebicycle/, the assembly was in Solidworks file format.

A reference frame named **“Absolute”** was placed in the point of contact of the rear wheel with the ground, and its axes were aligned with the same reported in the paper:

General dimensions were checked, like the wheel base.

In the model there are also the materials of the components and related density, so the weight of the whole system was checked:

###### Mass properties of GB7Sv9c

Configuration: Default

Coordinate system: Absolute

###### Mass = 8.83710 kilograms

###### Volume = 0.00170 cubic meters

###### Surface area = 0.91822 square meters

###### Center of mass: ( meters )

X = 0.56700

Y = 0.00000

Z = -0.39356

###### Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )

Taken at the center of mass.

Ix = (0.82606, 0.00001, -0.56358) Px = 0.16632

Iy = (0.56358, 0.00033, 0.82606) Py = 2.60452

Iz = (0.00020, -1.00000, 0.00026) Pz = 2.75102

###### Moments of inertia: ( kilograms * square meters )

###### Taken at the center of mass and aligned with the output coordinate system.

Lxx = 0.94075 Lxy = 0.00006 Lxz = -1.13511

Lyx = 0.00006 Lyy = 2.75102 Lyz = 0.00002

Lzx = -1.13511 Lzy = 0.00002 Lzz = 1.83009

###### Moments of inertia: ( kilograms * square meters )

Taken at the output coordinate system.

Ixx = 2.30956 Ixy = 0.00006 Ixz = -3.10713

Iyx = 0.00006 Iyy = 6.96088 Iyz = 0.00002

Izx = -3.10713 Izy = 0.00002 Izz = 4.67114

The ball bearings was semplified to limit the amount of work to do in Adams with the setup of constraints between the parts:

The density of torus component is set to maintain the same original total weight of the bearing.

The rear and front parameters are reported. However, in the 3d model the front frame subassembly doesn’t have the shaft and spacers, because these three component are in the rear frame subassembly.

It is believed these component more correctly belong to the front frame because they are fixed with it when it moves (rotate).

After “moving” these parts from the rear to the front frame subassembly, the weight and inertia tensor are congruent with those reported in the paper.

Rear frame:

###### Mass properties of GB7Sv9c

###### Configuration: Default

Coordinate system: Absolute

###### Mass = 6.42539 kilograms

###### Volume = 0.00126 cubic meters

###### Surface area = 0.73922 square meters

###### Center of mass: ( meters )

X = 0.50440

Y = 0.00000

Z = -0.42786

###### Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )

Taken at the center of mass.

Ix = (0.81896, 0.00001, -0.57385) Px = 0.04349

Iy = (0.57385, 0.00129, 0.81896) Py = 2.56756

Iz = (0.00075, -1.00000, 0.00105) Pz = 2.59280

###### Moments of inertia: ( kilograms * square meters )

Taken at the center of mass and aligned with the output coordinate system.

Lxx = 0.87468 Lxy = 0.00005 Lxz = -1.18621

Lyx = 0.00005 Lyy = 2.59280 Lyz = 0.00000

Lzx = -1.18621 Lzy = 0.00000 Lzz = 1.73638

###### Moments of inertia: ( kilograms * square meters )

Taken at the output coordinate system.

Ixx = 2.05093 Ixy = 0.00005 Ixz = -2.57288

Iyx = 0.00005 Iyy = 5.40378 Iyz = 0.00000

Izx = -2.57288 Izy = 0.00000 Izz = 3.37110

Front frame:

###### Mass properties of GB7Sv9c

Configuration: Default

Coordinate system: Absolute

###### Mass = 2.41171 kilograms

###### Volume = 0.00044 cubic meters

###### Surface area = 0.17900 square meters

###### Center of mass: ( meters )

X = 0.73380

Y = 0.00000

Z = -0.30219

###### Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )

Taken at the center of mass.

Ix = (-0.01492, -0.00043, -0.99989) Px = 0.00142

Iy = (0.03717, 0.99931, -0.00099) Py = 0.03825

Iz = (0.99920, -0.03718, -0.01490) Pz = 0.03839

###### Moments of inertia: ( kilograms * square meters )

Taken at the center of mass and aligned with the output coordinate system.

Lxx = 0.03838 Lxy = 0.00001 Lxz = 0.00055

Lyx = 0.00001 Lyy = 0.03825 Lyz = 0.00002

Lzx = 0.00055 Lzy = 0.00002 Lzz = 0.00143

###### Moments of inertia: ( kilograms * square meters )

Taken at the output coordinate system.

Ixx = 0.25862 Ixy = 0.00001 Ixz = -0.53424

Iyx = 0.00001 Iyy = 1.55710 Iyz = 0.00002

Izx = -0.53424 Izy = 0.00002 Izz = 1.30004

Wheels:

###### Mass properties of Aluwielcpl3

Configuration: Default

Coordinate system: Absolute

###### Mass = 0.17406 kilograms

###### Volume = 0.00006 cubic meters

###### Surface area = 0.02555 square meters

###### Center of mass: ( meters )

X = 0.00000

Y = 0.00000

Z = 0.00000

###### Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )

Taken at the center of mass.

Ix = (0.00000, 0.00000, 1.00000) Px = 0.00009

Iy = (1.00000, 0.00000, 0.00000) Py = 0.00009

Iz = (0.00000, 1.00000, 0.00000) Pz = 0.00018

###### Moments of inertia: ( kilograms * square meters )

Taken at the center of mass and aligned with the output coordinate system.

Lxx = 0.00009 Lxy = 0.00000 Lxz = 0.00000

Lyx = 0.00000 Lyy = 0.00018 Lyz = 0.00000

Lzx = 0.00000 Lzy = 0.00000 Lzz = 0.00009

###### Moments of inertia: ( kilograms * square meters )

Taken at the output coordinate system.

Ixx = 0.00009 Ixy = 0.00000 Ixz = 0.00000

Iyx = 0.00000 Iyy = 0.00018 Iyz = 0.00000

Izx = 0.00000 Izy = 0.00000 Izz = 0.00009

Counterwheels:

###### Mass properties of Alucounterwiel4

Configuration: Default

Coordinate system: Absolute

###### Mass = 0.16664 kilograms

###### Volume = 0.00006 cubic meters

###### Surface area = 0.02948 square meters

###### Center of mass: ( meters )

X = 0.00000

Y = 0.00000

Z = 0.00000

###### Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )

Taken at the center of mass.

Ix = (0.00000, 0.00000, -1.00000) Px = 0.00009

Iy = (-1.00000, 0.00000, 0.00000) Py = 0.00009

Iz = (0.00000, 1.00000, 0.00000) Pz = 0.00016

###### Moments of inertia: ( kilograms * square meters )

Taken at the center of mass and aligned with the output coordinate system.

Lxx = 0.00009 Lxy = 0.00000 Lxz = 0.00000

Lyx = 0.00000 Lyy = 0.00016 Lyz = 0.00000

Lzx = 0.00000 Lzy = 0.00000 Lzz = 0.00009

###### Moments of inertia: ( kilograms * square meters )

Taken at the output coordinate system.

Ixx = 0.00009 Ixy = 0.00000 Ixz = 0.00000

Iyx = 0.00000 Iyy = 0.00016 Iyz = 0.00000

Izx = 0.00000 Izy = 0.00000 Izz = 0.00009

The Solidworks’ convention for the inertia tensor is the following:

The entire assembly was exported in **Parasolid** file format with ascii encoding (*.x_t), with this format the parts were imported separately in Adams, and not in a only part (solids within) just like if other format was used instead (e.g.: step *.stp).

More than this, the Solidworks material density of each part was automatically included in the exported file, so in Adams the weight, positions of center of mass and inertia tensor was the same.

Setting in Adams the same reference frame **“Absolute”** of the paper, the following is obtained:

As can be seen, the center of mass of the structures in Y direction is not coincident with this axes, but decentralized even though a little bit.

The gravity is direct along the Z axis, for a standard Adams value of **-9.80665.
**

The units are:

The parameters for the rear frame only, reference frame **“absolute”**:

Rear frame with respect of its center of mass **“rear_cm”**, aligned with **“absolute”**:

The parameters for the front frame only, reference frame **“absolute”**:

Front frame with respect of its center of mass **“front_cm”**, aligned with **“absolute”**:

Rotating wheels (rear or front is the same), with respect of its center of mass **“front_whl_cm”**, aligned with **“absolute”**:

Counterotating wheels (rear or front is the same), with respect of its center of mass **“front_ctrwhl_cm”**, aligned with **“absolute”**:

The paper’s inertia moment along Y seem to have an error in the exponent of the power of 10, because the mantissa value is the same, only the exponent change. The mass is the also the same.

The Adams’ convention for the inertia tensor is the following:

Solidworks and Adams convention are the same.

**Constraint**

**“Fixed”**between the parts: even though in the reality there is the possibility of regulating the middle part of the structure, in Adams all the parts are fixed each other:

- On all 4 wheel, rotating and counterotating, a
**“revolute joint”**(between the support frame, front or rear, and the wheel) is set up for simplicity and functionality (in reality the wheel’s shaft is supported by two bearings):

**“Gear”**coupling with transmission ratio equal to one is set up between rotating and counterotating wheels:

In this way solid contact formulation is avoided, make the simulation more demanding.

The described cancellation of gyroscopic effect is kept anyway with this technique.

It is considered negligible the reduced mass to add to the total structure mass because of the inertia of the wheels, as well as the low value of this last parameter, also because this will alter the behavior of the structure in other directions in addition to the longitudinal. The structure is also launched with an initial speed and **is not considered the deceleration** related to the rolling resistance coefficient described in the paper, in this way less than other friction phenomena (which are not set) the structure maintains its launch velocity.

- Rotating wheels are set up with an angular speed during the simulation such that the peripheral speed corresponds to the longitudinal velocity of the point of contact, i.e.: hypothesis of pure rolling.

- A
**“revolute joint”**is set up on the hinge that connects the rear frame with the front frame, for the same reason of simplicity in the treatment is set with a single constraint although in reality such a hinge is formed by two bearings. The coulping is located between the support structure of the bearings and the shaft on which are mounted:

The joint is located in the middle of the shaft and the axis of rotation is aligned with the real hinge, and because the geometry of the assembly is not a simplified one, this ensures that the trail of

4 mm on the front wheel is respected.

The model is verified with Adams model check, with positive results:

The system has **7 degrees of freedom** because: 6 are those of the rear frame free to move in space, 1 of the front frame, after that the wheels have constraint of pure rolling (motion that removes the only rotational DOF that have wheels and counterwheels).Finally though there are the wheels contact with the ground, they don’t reduce the DOF because the Adams formulation to handle this type of contact is to impose the normal forces to the contact surface.

**Final modeling of wheel-ground contact**

The **point on plane** formulation was chosen, so two points are created on the bottom of two rotating wheels (the ones who touch the ground).

These points are not fixed with the wheels, otherwise they will rotate with them continuously changing the position, but they are fixed with the wheels’ supporting structure.

In this way the points remain in the same position (with respect of the support) like real contact point do.

To these points a unilateral constraint point on plane is applied, it prevent the vertical translation against the ground/plane with the action of a normal force.

No friction is set up:

Two forces are applied on each point, one transverse and one longitudinal with the wheel’s support (rear and front):

These forces follow the point, they change their magnitude but their direction remains fixed in space.

These two forces are the components of the lateral force of friction due to lateral slip of the wheel.

Their equations are:

**Determination of cornering stiffness Ks**

With the simplified equation governing the friction on the wheels, it is determined by trials the value of coefficient Ks that give the same behavior of the real structure.

The papers shows that the **critical speed** which split the zone of stability from the one of instability is **2.3 m/s**:

Various analyses are carried out launching the structure at this velocity, varying on each the coefficient from **1 to 30**, trying to find the stability.

To assess whether the structure is stable or not the maximum roll reached is measured after a reasonably long simulation time, **60 seconds**, allowing the structure the possibility of falling.

After **2 sec** from the start of the simulation a **impulse** is applied to the upper weight of structure, with a value of **0.5 N** for **0.05 s** in the direction Y, to speed up the process for oscillation of the bike to avoid natural divergence from unstable vertical position, because it would take much more time.

The value of the pulse is considered little so it doesn’t disturb the system too much.

The impulse is sum of two step function, with rising and falling time of **1/10** of the total duration of the impulse:

Small step size is used to better catch the impulse defined above. **Step size = 0.0025
**The values of Ks for each analysis are:

**1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0**

The maximum angle of roll reached during the simulation is measured, if it goes

**over 90°**a sensor stop the simulation for that value of Ks, to avoid the continuing of simulation because the structure is considered fallen:

Between 5 and 10 there is the value of Ks with which the bicycle becomes stable, for values greater than 10 the structure is stable at velocity of 2.3 m/s.

For the next analyses the **Ks is set up to 10** and a check is performed at 2.2 m/s where the bicycle must be unstable.

With 2.2 m/s launch and following impulse it is obtained:

After the impulse at t = 2 s the structure react with unstability followed by a fall, **Ks = 10** can be considered a good aproximation of (minimum) cornering stiffness of the wheels with the ground.

**Behavior at different speeds**

Analyses are carried out at various velocities and the maximum angle of roll is measured.

The speeds are:

**0.1, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 2.2, 2.3, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0
**It is obtained:

After 2.3 m/s the bicycle is stable.

Zooming in the stability region it is observed:

The values are very close to zero, which means the bicycle has not reached angles of roll so high that made her falling down, but it can be seen looking at the diagram in the time domain (always with vel. = 2.3 – 10 m/s) as this value is reached immediately after the impulse, and also as the maximum value reached changes with the speed, and even as the behavior between oscillating or damped changes with the speed:

Front rotation (steering axis) chart:

The trend follows the one shown in the paper chart where it can be seen at different speeds the value of real part of eigenvalue, which is the on who returns the damping capacity of bicycle at a stress:

Notice how the maximum damping is around 3-4 m/s as returned by the graph of simulations, for which the oscillations are no longer present to speed around 3 m/s.

It can be noticed in front rotation chart that at speeds greater than 4.0 m/s there is another oscillation in the first time after the impulse.

Measuring with FFT the value for 10.0 m/s velocity, in order to compare this with the value returned from the eigenmodes analysis, it is obtained:

With a peak around 1.84 Hz.

**Frequency determination**

The oscillation frequencies of the bicycle roll and front frame rotation at different speeds are determined, applying a **FFT** to the charts where the bicycle is unstable and calculating the **eigenmodes** where the bicycle is stable.

For both behavior the impulse is not given, in the stable region in order to have bicycle still and stable as possible for the calculating of eigenvalues, and in the unstable region because otherwise the bicycle would fall shortly after. In this region the structure is left free to swing due to the asymmetry in the vertical longitudinal plane and to achieve **one complete oscillation** after a time of **1 s** for contacts adjustment. For example at 1.5 m/s the roll has the following trend:

With a peak at 0.53 Hz.

It is obtained:

n.a. = no one (1) complete oscillation after 1 s from the simulation beginning (when t > 1 s).

For the speeds in which the bicycle is stable the eigenvalues will be determined after a very long time and for different moments of time. The assumption is that, if the bike is in steady state conditions, the same eigenvalues should be obtained.

The procedure is as follows. The structure is launched at different speeds. It is then simulated for 1200 s with step size = 0.01, and 250 eigenvalues are found every 1 s after the end (hence the final time is 1450 s). However, the eigenvalues have considerable variability from each other (with 1 s gap).

A graph is plotted which displays the value of eigenvalues in the complex plane.

It can be noticed that there are different number of eigenvalues that form curves and lines, these eigenvalues in the plane can be collected in different groups to calculate the average.

Completing the table for the speed greater than 2.3 m/s:

It can be noticed that the real part have a minimum around 3-4 m/s like the chart from the paper.

The imaginary part is constant from 3 to 7 m/s and it seem to be shifted to the right, with respect of the paper’s chart. Moreover the maximum absolute values reached for both real and imaginary part are not very similar to the ones in the paper’s chart.

**Yaw-lean rate impulse behavior**

Finally in graphical form the yaw-rate and lean-rate is plotted at different speed equal and greater than 2.3 m/s, where the damping is present.

The magnitude of impulse is adjusted in order to obtain similar values of yaw and lean rate reported in the paper, in particular 7 N.

The curves needs to be mirrored in the X axis, because maybe in the experimental test the impulse was given in the opposite direction from the one given in the simulation:

At 3.5 m/s:

It can be noticed that the trends of the analysised yaw-lean rate are more similar to the measured yaw-lean rate because the presence of oscillation, that are not present in the simulated yaw-lean rate.

**Conclusions**

It was possible to recreate the behavior of the TMS bicycle at different speeds, by using wheel-ground friction through simple equations. The effect of friction and the rolling resistance have been neglected. Trends are similar, however oscillation frequencies differ.

**References**

[1]** **J. D. G. Kooijman, J. P. Meijaard, Jim M. Papadopoulos, Andy Ruina,* A. L. Schwab, “A Bicycle Can Be Self-Stable Without Gyroscopic or Caster Effects”, Science 332, 339 (2011)