Estimating a Bicycle Wheel's Radial Inertia

This notebook introduces the relationship of vibrational motion to the inertia of a system and how one might estimate the inertia given measured vibrational data. The determination of the radial and axial moments of inertia of a bicycle wheel are used as a motivating example.

After the completion of this assignment students will be able to:

  • import measurement data into Python
  • visualize the vibrational measurements
  • use interactive plots to cause a simulation to match measurements
  • understand the concept of natural frequency and its relationship to mass/inertia and stiffness
  • state two of the three fundamental characteristics that govern vibration (mass/inertia and stiffness)

Inertia and Vibration

One of the fundamental properties that affects how systems vibrate is the inertia of a system. Inertia is colloquially defined as the tendency to do nothing or to remain unchanged. A vibrating system is changing with respect to time, thus we may infer that inertia will try to prevent the system from changing. There are two specific types of inertia that we will address in this class: mass which is a resistance to linear motion and moment of inertia which is a resistenace to rotational motion. The moment of inertia is more specifically a descriptor for the distribution of mass in a system.

Most people are familiar with an ordinary bicycle wheel, so we will look into the inertial characteristics of bicycle wheels and how this relates to vibration. A bicycle wheel is mostly symetric about a plane that passes through the cross sectional center line of the circular rim in addition to being symmetric about any radial plane normal to the wheel plane. Recalling from dynamics, this means that the moment of inertia will have three prinicipal moments of inertia, two of which are equal due to the nature of the symmetry. The inertia about the axle of the wheel is the rotational inertia (resists the wheel rolling) and the inertia about any radial axis resists motions like turning the handlebars, etc. We demonstrated in the previous lesson that the inertia of the book affected the period of oscillation, and thus if a bicycle wheel is vibrated its inertia should affect the vibration in some way also.

Torsional Pendulum

Vibrations of a mass or inertia occur when there is a force that acts on the system which displaces it from its equilibrium. A common way to measure inertia is to use a torsional spring with a known linear force/displacement ratio, the spring constant or stiffness, to resist motion and thus create this force. If we attach a torsional spring along a radial line of the bicycle wheel, through its center of mass and glue the other end to the ceiling it is possible to create a torsional pendulum. If you twist the wheel about the axis of the torsional spring it will start vibrating.

The video belows shows such a setup. A bicycle wheel is attached to the ceiling by a slender steel rod. The rod's axis passes through the center of mass of the wheel. When the bicycle wheel is given an initial angular displacement from the equilibrium and then released, it vibrates.

In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('uuc1_RTYZLw', width=640, height=480)
Out[1]:

A free body diagram can be sketched of this system:

image

Exercise

This system is slightly different than the book and cup system in the previous lesson. There are two fundamental properties of the system that make this it vibrate. What are these two things and how is it differen than the book on the cup?

Write your answer here

Data

During the experiment shown in the above video I recorded some important constants about the system and measured the angular velocity of the system about the rod's axis with a angular velocity gyro (which also functions internally based on vibrations). A gyro is a micro electromechanical device that outputs a voltage proportional to the angular velocity it experiences. The data is stored in a plain text comma separated value (CSV) file. The Unix head command can show the first set of lines and can be invoke in Jupyter using the prepended !.

In [2]:
!head bicycle-wheel-radial-inertia-rate-gyro-measurement.csv

time,angular_velocity
0.0,3.236790061330015
0.0020000000000000018,3.2304788385296246
0.003999999999999997,3.2195776355287498
0.005999999999999998,3.2092501800752573
0.008,3.1926115018876944
0.009999999999999995,3.1880215217073844
0.011999999999999997,3.1765465712743266
0.013999999999999999,3.1490066903384144
0.016,3.123188052093832

The data from the measurement can be loaded into a Panda's DataFrame with Panda's read_csv() function. The time index is in seconds and the angular velocity is in radians per second.

In [3]:
import pandas as pd
In [4]:
radial_gyro_meas = pd.read_csv('bicycle-wheel-radial-inertia-rate-gyro-measurement.csv', index_col='time')

This file has 1001 rows:

In [5]:
len(radial_gyro_meas)
Out[5]:
1001

The head() and tail() functions associated with a DataFrame (identical to the Unix head and tail) can be used to inspect the beginning and end of the data:

In [6]:
radial_gyro_meas.head()
Out[6]:
angular_velocity
time
0.000 3.236790
0.002 3.230479
0.004 3.219578
0.006 3.209250
0.008 3.192612
In [7]:
radial_gyro_meas.tail()
Out[7]:
angular_velocity
time
1.992 -0.271675
1.994 -0.207989
1.996 -0.144303
1.998 -0.085207
2.000 -0.030127

This can then be plotted to see what the measurement looks like. I use the '.' style to show the individual measurements more clearly. You can zoom in to see the individual measurements. Also, remember to enable the plotting display with:

In [8]:
%matplotlib inline
In [9]:
radial_gyro_meas.plot(style='.');

Exercise

Use the zoom feature and make a best estimate of the period of oscillation, i.e. the time for one cycle to occur (time between two peaks).

In [10]:
# write your answer here

There is a convenience function estimate_period() that attempts a rough estimate that is helpful. Does this give a similar answer to your visual estimate?

In [11]:
from resonance.functions import estimate_period
In [12]:
estimate_period?
In [13]:
estimate_period(radial_gyro_meas.index, radial_gyro_meas.angular_velocity)
Out[13]:
0.42100000000000004

Simulating the system

There is a system in the resonance package that represents a basic torsional pendulum that can be used to simulate this system:

In [14]:
from resonance.linear_systems import TorsionalPendulumSystem
In [15]:
sys = TorsionalPendulumSystem()

Note that the constants are not yet set:

In [16]:
sys.constants
Out[16]:
{'rotational_inertia': 0.0,
 'torsional_damping': 0.0,
 'torsional_stiffness': 0.0}

Exercise

Calculate the stiffness of the circular cross section rod based on your knowledge of torsional elastic mechanics. Here are some of the torsion rod's geometry and material properties:

  • Length of the torsion rod, $l$ : 1.05 m
  • Diameter of the torsion rod, $d$ : 6.35 mm
  • Modulus of Rigidity of steel, $G$ : 77 GPa

Set the system's torsional_stiffness to the value you calculate.

In [17]:
import numpy as np
l = 1.05  # m
d = 0.00635  # m
G = 77.0E9  # Pa
J = np.pi * d**4 / 32
sys.constants['torsional_stiffness'] = G * J / l
sys.constants
Out[17]:
{'rotational_inertia': 0.0,
 'torsional_damping': 0.0,
 'torsional_stiffness': 11.705668520051944}
In [18]:
# write your answer here

Exercise

As a starting estimate of the inertia of the bicycle wheel, use the axial moment of inertia of an infinitely thin hoop of mass, $m$. See Wikipedia's List of moments of inertia for the equation. The bicycle wheel used in the data collection had these values:

  • Outer radius of the bicycle wheel, $R$ : 0.336 m
  • Mass of the bicycle wheel, $m$ : 1.55 kg

Set the rotational inertia value in the system to the value you calculate.

In [19]:
mass = 1.55  # kg
radius = 0.336  # m
radial_inertia = mass * radius**2 / 4.0 # kg m**2
sys.constants['rotational_inertia'] = radial_inertia
sys.constants
Out[19]:
{'rotational_inertia': 0.04374720000000001,
 'torsional_damping': 0.0,
 'torsional_stiffness': 11.705668520051944}
In [ ]:

Exercise

Use the above constants and simulate the system with an initial angular velocity that matches the initial measured angular velocity in the measurements and store the result in a variable named trajectory. You can set the initial angular velocity by accessing the System.speeds dictionary and assigning a value. One way to select a single value from a Pandas DataFrame is to use the DataFrame.loc[] notation. This will allow you to select a row by it's index (a time value). Calculate the period of oscillation of this simulation and see if it is similar to the period you estimated from the data above. If the period is different, why might this be so?

In [20]:
sys.speeds['torsion_angle_vel'] = radial_gyro_meas.loc[0.0]['angular_velocity']

duration = radial_gyro_meas.index[-1] - radial_gyro_meas.index[0]

trajectory = sys.free_response(duration)

trajectory.head()
Out[20]:
torsion_angle torsion_angle_acc torsion_angle_vel
time
0.00 0.000000 -0.000000 3.236790
0.01 0.032224 -8.622277 3.193582
0.02 0.063587 -17.014358 3.065113
0.03 0.093253 -24.952190 2.854811
0.04 0.120429 -32.223852 2.568291
In [21]:
trajectory.plot(subplots=True)
Out[21]:
array([<matplotlib.axes._subplots.AxesSubplot object at 0x7f84b2d3fd30>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7f84b2cdab50>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7f84b2ce79d0>],
      dtype=object)

Interactive Plots

You should see that the computational system does not quite predict the actual period of the oscillation. One useful tool for interactively changing the constants of the system is a Jupyter "widget". Widgets let you create sliders, drop downs, value entry boxes, etc. to interactively change the results of the simulation. Below is an example of how you can create a widget that updates a plot interactively. In this case we'd like to adjust the inertia value until the period of oscillation of the simulation matches the period of the actual data.

The first step is to create a plot of the data, as we did above, but make sure to assign the axis that command creates to a variable ax so that we can add more information to the plot.

The second step is to define a function which has an input that you want to adjust interactively, in our case the radial inertia. This function should compute the new simulation trajectory, set the values for the simulation line with the new data, and finally redraw the figure.

In [22]:
import matplotlib.pyplot as plt
In [23]:
fig, ax = plt.subplots(1, 1)

# create initial plot
meas_lines = ax.plot(radial_gyro_meas.index, radial_gyro_meas.angular_velocity, '.')

# add a line to the plot
sim_lines = ax.plot(trajectory.index, trajectory.torsion_angle_vel)

ax.set_ylabel('Angular Rate [rad/s]')
ax.set_xlabel('Time [s]')
ax.legend(['Measured', 'Simulation'], loc=1)

# now create a function that has an input that you'd like to varying interactively
# the function should compute the updated trajectory with this new value and change the data on the line
def plot_trajectory_comparison(rotational_inertia=0.04):
    # set the new inertia value
    sys.constants['rotational_inertia'] = rotational_inertia
    # simulate the system with new value
    trajectory = sys.free_response(duration)
    # set the x and y data of the simulation line to new data
    sim_lines[0].set_data(trajectory.index, trajectory.torsion_angle_vel)
    # redraw the figure with the updated line
    fig.canvas.draw()

# call the function to make the initial plot
plot_trajectory_comparison()

Now the fun part. You can use the interact function from the ipywidgets package to make the plot interactive. You pass in the function that does the plot updating and then a range for the slider for the axial inertia as (smallest_value, largest_value, step_size). This will create a slider that updates the plot as you change the value.

In [24]:
from ipywidgets import interact
In [25]:
widget = interact(plot_trajectory_comparison, rotational_inertia=(0.01, 0.2, 0.0005));

Here is a way to get to the value that the widget is set too (it's a bit obtuse, but we'll show different methods later):

In [26]:
widget.widget.children
Out[26]:
(FloatSlider(value=0.04, description='rotational_inertia', max=0.2, min=0.01, step=0.0005),
 Output())
In [27]:
widget.widget.children[0].value
Out[27]:
0.04

How does inertia relate to the period?

By now it should be pretty clear that there is a relationship between the inertia of the system and the period of oscillation. It would be nice to plot the period versus the change in inertia to try and determine what the relationship is. Say we want to check a range of inertia values from X to Y, we can create those values with:

In [28]:
import numpy as np
In [29]:
inertias = np.linspace(0.001, 0.2, num=50)

Python For Loops

Instead of typing the simulation code out for each of the 50 inertia values, we can use a loop to iterate through each value and save some typing. Python loops are constructed with the for command and, in the simplest case, the in command. Here is an example that iterates through the inertia values, computes what the hoop radius would be if the mass was fixed, and prints the value of the radius on each iteration:

In [30]:
for inertia in inertias:
    radius = np.sqrt(2 * inertia / mass)
    print(radius)

0.03592106040535498
0.08081220356417686
0.10849378742191074
0.13042695876774693
0.14916930393903552
0.16580642599703171
0.18092003915034385
0.19486495382806718
0.20787650154519038
0.22012026158198236
0.2317179734573327
0.24276224765689153
0.25332547947586864
0.2634655351576967
0.27322953332351235
0.28265644790883004
0.29177895264855885
0.3006247609112442
0.30921762003028797
0.31757806316156084
0.32572398724334245
0.3336711038282367
0.34143329538023487
0.3490229001920354
0.35645094265931226
0.3637273211995012
0.37086096296485277
0.37785995225145347
0.38473163787206727
0.3914827235574724
0.3981193445541539
0.4046471329102455
0.41107127342680455
0.41739655185565144
0.4236273966178773
0.429767915076839
0.43582192520999435
0.44179298337343326
0.44768440873259313
0.45349930483574535
0.4592405787283579
0.46491095794249665
0.47051300564305626
0.4760491341694841
0.48152161717597713
0.48693260054346443
0.4922841122119156
0.49757807106073987
0.5028162949475479
0.508000508000762

It is also often useful to store the computed values in a list or array as the loop iterates so you can use the values in future computations. Here is a way to store the values in a Python list and then convert the list to a NumPy array.

In [31]:
radii = []  # create an empty list
for inertia in inertias:
    radius = np.sqrt(2 * inertia / mass)
    radii.append(radius)  # add the radius to the end of the list
radii = np.array(radii)  # convert the list to a NumPy array
radii
Out[31]:
array([0.03592106, 0.0808122 , 0.10849379, 0.13042696, 0.1491693 ,
       0.16580643, 0.18092004, 0.19486495, 0.2078765 , 0.22012026,
       0.23171797, 0.24276225, 0.25332548, 0.26346554, 0.27322953,
       0.28265645, 0.29177895, 0.30062476, 0.30921762, 0.31757806,
       0.32572399, 0.3336711 , 0.3414333 , 0.3490229 , 0.35645094,
       0.36372732, 0.37086096, 0.37785995, 0.38473164, 0.39148272,
       0.39811934, 0.40464713, 0.41107127, 0.41739655, 0.4236274 ,
       0.42976792, 0.43582193, 0.44179298, 0.44768441, 0.4534993 ,
       0.45924058, 0.46491096, 0.47051301, 0.47604913, 0.48152162,
       0.4869326 , 0.49228411, 0.49757807, 0.50281629, 0.50800051])

If you want to work with the NumPy array up front you can create an empty NumPy array and then fill it by using the correct index value. This index value can be exposed using the enumerate() function. For example:

In [32]:
radii = np.zeros_like(inertias)
radii
Out[32]:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [33]:
for i, inertia in enumerate(inertias):
    radius = np.sqrt(2 * inertia / mass)
    radii[i] = radius
radii
Out[33]:
array([0.03592106, 0.0808122 , 0.10849379, 0.13042696, 0.1491693 ,
       0.16580643, 0.18092004, 0.19486495, 0.2078765 , 0.22012026,
       0.23171797, 0.24276225, 0.25332548, 0.26346554, 0.27322953,
       0.28265645, 0.29177895, 0.30062476, 0.30921762, 0.31757806,
       0.32572399, 0.3336711 , 0.3414333 , 0.3490229 , 0.35645094,
       0.36372732, 0.37086096, 0.37785995, 0.38473164, 0.39148272,
       0.39811934, 0.40464713, 0.41107127, 0.41739655, 0.4236274 ,
       0.42976792, 0.43582193, 0.44179298, 0.44768441, 0.4534993 ,
       0.45924058, 0.46491096, 0.47051301, 0.47604913, 0.48152162,
       0.4869326 , 0.49228411, 0.49757807, 0.50281629, 0.50800051])

It is also worth noting that NumPy provides vectorized functions and that loops are not needed for the above example, for example:

In [34]:
radii = np.sqrt(2 * inertias / mass)
radii
Out[34]:
array([0.03592106, 0.0808122 , 0.10849379, 0.13042696, 0.1491693 ,
       0.16580643, 0.18092004, 0.19486495, 0.2078765 , 0.22012026,
       0.23171797, 0.24276225, 0.25332548, 0.26346554, 0.27322953,
       0.28265645, 0.29177895, 0.30062476, 0.30921762, 0.31757806,
       0.32572399, 0.3336711 , 0.3414333 , 0.3490229 , 0.35645094,
       0.36372732, 0.37086096, 0.37785995, 0.38473164, 0.39148272,
       0.39811934, 0.40464713, 0.41107127, 0.41739655, 0.4236274 ,
       0.42976792, 0.43582193, 0.44179298, 0.44768441, 0.4534993 ,
       0.45924058, 0.46491096, 0.47051301, 0.47604913, 0.48152162,
       0.4869326 , 0.49228411, 0.49757807, 0.50281629, 0.50800051])

It is a good idea to use vectorized functions instead of loops when working with arrays if you can, because they are much faster:

In [35]:
%%timeit -n 1000

radii = []
for inertia in inertias:
    radius = np.sqrt(2 * inertia / mass)
    radii.append(radius)
radii = np.array(radii)

92.4 µs ± 3.33 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [36]:
%%timeit -n 1000

np.sqrt(2.0 * inertias / mass)

2.45 µs ± 87.3 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Exercise

Use a loop to construct a list of periods for different inertia values. After you have both arrays, plot the inertias on the $X$ axis and the frequencies on the $Y$ axis. Is there any functional relationship that describes the relationship between the variables?

In [37]:
periods = []
for inertia in inertias:
    sys.constants['rotational_inertia'] = inertia
    trajectory = sys.free_response(duration)
    periods.append(estimate_period(trajectory.index, trajectory.torsion_angle))

fig, ax = plt.subplots(1, 1)
ax.plot(inertias, periods)
ax.set_xlabel('Inertia [$kg \cdot m^2$]')
ax.set_ylabel('Period [s]');
In [ ]:

Trying plotting the inertia versus the square root of the inertia to see if there are any similarities.

In [38]:
plt.figure()
plt.plot(inertias, np.sqrt(inertias))
plt.xlabel('$I$')
plt.ylabel('$\sqrt{I}$');

It turns out that the period of oscillation, $T$, is proportional to the square root of the oscillating inertia.

$$ T \propto \sqrt{I} $$

and more precisely the proportionality coefficient ends up being:

$$ \frac{2\pi}{\sqrt{k}} $$

we will discuss why this is the case when we get to modeling in later notebooks.

In [39]:
plt.figure()
x = np.linspace(0, 0.2, 100)
plt.plot(inertias, 2 * np.pi / np.sqrt(sys.constants['torsional_stiffness']) * np.sqrt(inertias))
plt.plot(inertias, periods)
plt.xlabel('$I$')
plt.ylabel(r'$\frac{2 \pi}{\sqrt{k}} \sqrt{I}$');

Frequency and Period

We've been talking about the period of oscillation so far, which is defined as the time between a repeating configuration, i.e. seconds per cycle of oscillation. The inverse of period is frequency: the number of cycles per second. The unit Hertz is a special unit for cycles per second.

Exercise

Given a period of 1.35 seconds, (1) what is the frequency in Hertz? and (2) what is the is frequency in radians per second? (3) create two functions, one that coverts period to frequency in radians per second called period2freq and one that does the opposite called freq2period.

In [40]:
T = 1.35  # s
f = 1.0 / 1.35  # H
print('Frequency: {} H'.format(f))

omega = 2 * np.pi * f  # rad/s
print('Frequency: {} rad/s'.format(omega))

def period2freq(period):
    return 1.0 / period * 2.0 * np.pi

def freq2period(freq):
    return 1.0 / freq * 2.0 * np.pi

print('Period from frequency: {}'.format(period2freq(1.0)))

print('Frequency from period: {}'.format(freq2period(6.283185)))

Frequency: 0.7407407407407407 H
Frequency: 4.654211338651545 rad/s
Period from frequency: 6.283185307179586
Frequency from period: 1.000000048889152
In [ ]:

The fundamental frequency of osciallation of the generalized coordinate of a single degree of freedom system is called the natural frequency. The natural frequency is named because it is the frequency a system tends to oscillate at when no driving forces are applied.

Decaying Oscillations

Exercise

The system also has a constant called torsional_damping. This is constant represents the proportion of a force that is linearly related to the angular velocity of the oscillation of the wheel. The faster it rotates the more force is applied to slow it down. This will cause more damping in the simulation the higher the value is. It is a resonable representation of the damping from air resistance on the wheel. Make an interactive plot like you did above with two adjustable constants: torsional_damping and rotational_inertia. Adjust the parameters until the simulation matches. Does the simulation match the best fit from above?

In [41]:
ax = radial_gyro_meas.plot(style='.')
line = ax.plot(trajectory.index, trajectory.torsion_angle_vel)
ax.set_ylabel('Angular Rate [rad/s]')
ax.set_xlabel('Time [s]')
ax.legend(['Measured', 'Simulation'], loc=1)

def plot_trajectory_comparison(rotational_inertia=0.044, torsional_damping=0.0):
   sys.constants['rotational_inertia'] = rotational_inertia  # set the new inertia value
   sys.constants['torsional_damping'] = torsional_damping
   traj = sys.free_response(duration)  # simulate the system with new value
   line[0].set_data(traj.index, traj.torsion_angle_vel)  # set the x and y data of the simulation line to new data
   plt.gcf().canvas.draw()  # redraw the figure with the updated line

# call the function to make the initial plot
plot_trajectory_comparison()
In [42]:
# write solution here
In [43]:
widget = interact(plot_trajectory_comparison,
                  rotational_inertia=(0.01, 0.2, 0.001),
                  torsional_damping=(0.0, 0.05, 0.001));

<Figure size 432x288 with 0 Axes>
In [44]:
widget.widget.children[0].value
Out[44]:
0.044
In [45]:
widget.widget.children[1].value
Out[45]:
0.0