3. Measuring a Bicycle Wheel's 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
- use curve fitting to estimate the period of oscillation
- 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)
- discover the center of percussion of a compound pendulum
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.
from IPython.display import YouTubeVideo
YouTubeVideo('uuc1_RTYZLw', width=640, height=480)
A free body diagram can be sketched of this system:
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 !.
!head bicycle-wheel-radial-inertia-rate-gyro-measurement.csv
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.
import pandas as pd
radial_gyro_meas = pd.read_csv('bicycle-wheel-radial-inertia-rate-gyro-measurement.csv', index_col='time')
This file has 1001 rows:
len(radial_gyro_meas)
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:
radial_gyro_meas.head()
radial_gyro_meas.tail()
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:
%matplotlib inline
radial_gyro_meas.plot(style='.');
Exercise
Use your period estimation function from the previous lesson to try to get an estimate of the period of this oscillation.
from resonance.functions import estimate_period
estimate_period(radial_gyro_meas.index, radial_gyro_meas.angular_velocity)
# write your answer here
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:
from resonance.linear_systems import TorsionalPendulumSystem
sys = TorsionalPendulumSystem()
Note that the constants are not yet set:
sys.constants
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.
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
# 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.
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
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?
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)
estimate_period(trajectory.index, trajectory.torsion_angle)
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.
import matplotlib.pyplot as plt
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.
from ipywidgets import interact
widget = interact(plot_trajectory_comparison, rotational_inertia=(0.01, 0.2, 0.001));
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):
widget.widget.children
widget.widget.children[0].value
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:
import numpy as np
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:
for inertia in inertias:
radius = np.sqrt(2 * inertia / mass)
print(radius)
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.
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
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:
radii = np.zeros_like(inertias)
radii
for i, inertia in enumerate(inertias):
radius = np.sqrt(2 * inertia / mass)
radii[i] = radius
radii
It is also worth noting that NumPy provides vectorized functions and that loops are not needed for the above example, for example:
radii = np.sqrt(2 * inertias / mass)
radii
It is a good idea to use vectorized functions instead of loops when working with arrays if you can, because they are much faster:
%%timeit -n 1000
radii = []
for inertia in inertias:
radius = np.sqrt(2 * inertia / mass)
radii.append(radius)
radii = np.array(radii)
%%timeit -n 1000
np.sqrt(2.0 * inertias / mass)
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?
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]');
Trying plotting the inertia versus the square root of the inertia to see if there are any similarities.
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.
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.
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)))
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.
Curve Fitting¶
Above, we interactively adjusted the inertia such that the simulation's period matched the measured data. Another way to determine the period of oscillation more accurately is to try and find the best fitting curve through the data points. The Python package scipy provides a very conveinent function called curve_fit. This function implements a least squares method that finds an optimal fit based on parameterized function provided by the user. Let's see how it works.
from scipy.optimize import curve_fit
First define a Python function where the first input argument are the "x" values (an array) and the remaining arguments are the parameters that can be adjusted to find the optimal fit. In this case, a simple cosine function with an adjustable amplitude and frequency is created.
def cos_func(times, amplitude, frequency):
return amplitude * np.cos(frequency * times)
To check that it works, provide some values:
cos_func(np.linspace(0, 10), 5.0, 5.0)
Now, provide this function to curve_fit along with the measure data (x and y) and an initial guess for the amplitude and frequency. A good initial guess is important, as the optimal solution can't always be found from an arbitrary initial guess. The function curve_fit returns two items. The first is the optimal values of the two parametes and the second is the covariance matrix that gives an idea of how certain the value of the parameters are. We will just work with the first value for now.
popt, pcov = curve_fit(cos_func, # our function
radial_gyro_meas.index, # measured x values
radial_gyro_meas.angular_velocity, # measured y values
p0=(3.0, period2freq(0.44))) # the initial guess for the two parameters
Now we see the optimal values for the amplitude and frequency:
popt
It is useful to plot the optimal function versus the measured data:
fig, ax = plt.subplots(1, 1)
ax.plot(radial_gyro_meas.index, radial_gyro_meas.angular_velocity, '.')
ax.plot(radial_gyro_meas.index, cos_func(radial_gyro_meas.index, popt[0], popt[1]));
That looks like a pretty nice fit! The period of this fitted function can now be computed:
freq2period(popt[1])
Exercise
This system, just like the book on a cup system has an amplitude that decreases with respect to time. This decrease is often referred to as a decaying amplitude. We just fit the mathematical function:
$$ \dot{\theta}(t) = A \cos{\omega t} $$where $A$ is the amplitude of the angular rate $\dot{\theta}$ and $\omega$ is the frequency of oscillation in radians per second. One way to account for the decay in amplitude oscillation is to introduce a multiplicative factor of $e^{-\lambda t}$ into the equation. This will cause exponential decay at the rate $\lambda$ (units are 1 / s). The mathematical equation will look like:
$$ \dot{\theta}(t) = A e^{-\lambda t} \cos{\omega t} $$Recall that $e^{-\lambda t}$ looks like:
plt.figure()
t = np.linspace(0, 2, num=100)
plt.plot(t, np.exp(-1.5 * t));
Use decaying oscillation mathematical function to create a curve fitting function and find the values of $A$, $\lambda$, and $\omega$ that best fit the data. Calculate the period of oscillation and compare it to the period from the purely sinusoidal fit from above. Is there any difference in the period of oscillation?
def decaying_sinusoid(t, a, lam, w):
return a * np.exp(lam * t) * np.cos(w * t)
popt, pcov = curve_fit(decaying_sinusoid,
radial_gyro_meas.index,
radial_gyro_meas.angular_velocity,
p0=(3.0, -0.0002, period2freq(0.44)))
fig, ax = plt.subplots(1, 1)
ax.plot(radial_gyro_meas.index, radial_gyro_meas, '.')
ax.plot(radial_gyro_meas.index, decaying_sinusoid(radial_gyro_meas.index, popt[0], popt[1], popt[2]));
freq2period(popt[2])
# write solution here
The natural frequency of a decayed oscillation is slightly different than that of a non-decaying oscillation.
This function does an excellent job at fitting the measurements!
Exercise
The system also has a constant called viscous_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: viscous_damping and radial_inertia. Adjust the parameters until the simulation matches. Does the simulation match the best fit from above?
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()
# write solutin here
widget = interact(plot_trajectory_comparison,
rotational_inertia=(0.01, 0.2, 0.001),
torsional_damping=(0.0, 0.05, 0.001));
widget.widget.children[0].value
widget.widget.children[1].value
Compound Pendulum¶
At this point we have figured out the moment of inertia about any radial axis that passes through the center of the wheel. We would also like to find the moment of inertia about the wheel's axle, the spin moment of inertia. This moment of interia affects how fast the bicycle wheel can be acelerated and decelerated. We could potentally hang the bicycle wheel on the torsion rod such that the axle's axis aligns with the torsion bar's axis. But there is a simpler way that only requires a fulcrum point and gravity to act as our (constant force) "spring". The video below shows a children's bicycle wheel hanging by the inner diameter of the rim on a small circular rod and a angular velocity gyro is attached to the wheel with it's measurement axis align with the wheel's axle. This arrangement is called a compound pendulum. A compound pendulum is any single degree of freedom pendulum in whch the swinging portion can be considered a rigid body.
YouTubeVideo('D2tSoGqhtx0', width=640, height=480)
The figure below shows a free body diagram of a compound pendulum. In our case the wheel is the rigid body and the revolute joint is the fulcrum.
<img src="bicycle-wheel-axial-inertia-measurement-fbd.png" >
The data from the measurement of the same full sized bicycle wheel from the previous analysis can be loaded with:
axial_gyro_meas = pd.read_csv('bicycle-wheel-axial-inertia-rate-gyro-measurement.csv',
index_col='time')
axial_gyro_meas.head()
axial_gyro_meas.plot(style='.');
Note that this measurement has much less damping than the torsional pendulum. Once again our goal is to determine what the moment of the inertia of the bicycle wheel is. We can use a CompoundPendulumSystem to simulate this motion:
from resonance.linear_systems import CompoundPendulumSystem
cpend_sys = CompoundPendulumSystem()
The mass and radius of the wheel is the same as above. You can estimate the rotational moment of inertia about the axle as:
$$ I = m r^2 $$But the inertia about the joint has to be computed using the parallel axis theorem.
Exercise
Using the parallel axis thereom, compute an estimate of the inertia about the joint if the distance from the axle to the joint is 0.296 meters. Store both results in variables named inertia_about_axle and inertia_about_joint.
mass = 1.55 # kg
radius = 0.336 # m
inertia_about_axle = mass * radius**2
print(inertia_about_axle)
inertia_about_joint = inertia_about_axle + mass * 0.296**2
print(inertia_about_joint)
# write your solution here
Now simulate the pendulum with an initial angular velocity that is the same as measurement and see how well they match.
cpend_sys.constants
cpend_sys.constants['acc_due_to_gravity'] = 9.81 # m/**2
cpend_sys.constants['inertia_about_joint'] = inertia_about_joint # kg m**2
cpend_sys.constants['joint_to_mass_center'] = 0.296 # m
cpend_sys.constants['pendulum_mass'] = mass # kg
cpend_sys.speeds['angle_vel'] = axial_gyro_meas.loc[0.0]['angular_velocity']
trajectory = cpend_sys.free_response(5)
fig, ax = plt.subplots(1, 1)
ax.plot(axial_gyro_meas.index, axial_gyro_meas, '.')
line = ax.plot(trajectory.index, trajectory['angle_vel'])
def plot(inertia=0.5):
cpend_sys.constants['inertia_about_joint'] = inertia
traj = cpend_sys.free_response(5)
line[0].set_data(traj.index, traj['angle_vel'])
interact(plot, inertia=(0.0, 1.0, 0.001));
Exercise
Use your decaying_sinusoid function to find a best fit for this data.
popt, pcov = curve_fit(decaying_sinusoid,
axial_gyro_meas.index,
axial_gyro_meas.angular_velocity,
p0=(0.6, -0.0002, 1.0 / 0.44 * np.pi * 2))
ax = axial_gyro_meas.plot(style='.')
ax.plot(axial_gyro_meas.index, decaying_sinusoid(axial_gyro_meas.index, *popt));
freq2period(popt[2])
# write you answer here
Exercise
Finally, using the period of oscillation, $T$, from the data above, you can compute the inertia about the joint using the CompoundPendulumSystem or this relationship:
Store the result of this calculation in a variable called actual_inertia_about_the_joint.
# write your solution here
actual_inertia_about_joint = freq2period(popt[2])**2 / 4/ np.pi**2 * mass * 9.81 * 0.296
actual_inertia_about_joint
Exercise
Use the parallel axis thereom to find the moment of inertia about the axle (the location of the center of mass) and store it in a variable called acutal_inertia_about_axle.
actual_inertia_about_axle = actual_inertia_about_joint - mass * 0.296**2
actual_inertia_about_axle
# write your answer here
Center of Percussion¶
The center of percussion of a compound pendulum is described as the location on the pendulum that corrsponds to the length of a simple pendulum of the same mass which has the same frequency of oscillation. The following video demonstrates the principle:
YouTubeVideo('Dw3UpKQVhVY', width=640, height=480)
Using a SimplePendulumSystem we can simulate it alongside the compound pendulum to determine the center of percussion.
from resonance.linear_systems import SimplePendulumSystem
spend_sys = SimplePendulumSystem()
spend_sys.constants
spend_sys.constants['acc_due_to_gravity'] = 9.81
spend_sys.constants['pendulum_mass'] = 1.55
spend_sys.constants['pendulum_length'] = 0.3
spend_sys.speeds['angle_vel'] = axial_gyro_meas.loc[0.0]['angular_velocity']
spend_traj = spend_sys.free_response(5)
cpend_sys.constants['inertia_about_joint'] = actual_inertia_about_joint
cpend_traj = cpend_sys.free_response(5)
fig, ax = plt.subplots(1, 1)
ax.plot(cpend_traj.index, cpend_traj.angle)
line = ax.plot(spend_traj.index, spend_traj.angle)
def plot(pend_length=0.3):
spend_sys.constants['pendulum_length'] = pend_length
spend_traj = spend_sys.free_response(5)
line[0].set_data(spend_traj.index, spend_traj.angle)
plot()
interact(plot, pend_length=(0, 1.0, 0.01));
More on center of percussion¶
The baseball bat: https://physics.csuchico.edu/baseball/DrBaseball/SweetSpot/