{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 4. Clock Pendulum with Air Drag Damping\n", "\n", "This notebook introduces the third fundamental characteristic of vibration: energy dissipation through damping. A compound pendulum system representing a clock pendulum is implemented that allows students to vary the damping parameters and visualize the three regimes of linear damping.\n", "\n", "After the completion of this assignment students will be able to:\n", "\n", "- understand the concept of damped natural frequency and its relationship to\n", " mass/inertia, stiffness, and damping\n", "- state the three fundamental characteristics that make a system vibrate\n", "- compute the free response of a linear system with viscous-damping in all\n", " three damping regimes\n", "- identify critically damped, underdamped, and overdamped behavior\n", "- determine whether a linear system is over/under/critically damped given its\n", " dynamic properties\n", "- understand the difference between underdamping, overdamping, and crticial\n", " damping\n", " \n", "## Introduction\n", "\n", "Many clocks use a pendulum to keep time. Pendulum's have a very constant osicallation period and if designed with quality components take very little energy to run. There is a downside to pendulums though. Any friction in the pendulum's pivot joint or any air drag on the pendulum itself will cause the pendulum to slow stop oscillating and this energy dissaption will affect the period of oscillation. `resonance` includes a simple clock pendulum that represents a clock pendulum that looks like:\n", "\n", "![](fig/04/clock-pendulum.png)\n", "\n", "Import the pendulum as so:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from resonance.linear_systems import ClockPendulumSystem" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sys = ClockPendulumSystem()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check out its constants:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.constants" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the coordinates and speeds:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.coordinates" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.speeds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The system can be simulated as usual if the coordinates or speeds are set to some initial value." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sys.coordinates['angle'] = np.deg2rad(5.0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "traj = sys.free_response(5.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And a plot can be shown:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "ax.plot(traj.index, traj.angle)\n", "ax.set_ylabel('Angle [rad]')\n", "ax.set_xlabel('Time [s]');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above simulation shows that we get a sinusoid oscillation with a period of about 1 second, which is good." ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Creative an interactive plot of the angle trajectory with a slider for the `viscous_damping` coefficient. The slider shoudl go from 0.0 to 5.0 with step of 0.1. The code should follow this pattern, as before:\n", "\n", "```python\n", "fig, ax = plt.subplots(1, 1)\n", "\n", "sim_line = ax.plot(traj.index, traj.angle)[0]\n", "ax.set_ylim((-sys.coordinates['angle'], sys.coordinates['angle']))\n", "\n", "ax.set_ylabel('Angle [rad]')\n", "ax.set_xlabel('Time [s]')\n", "\n", "def plot_trajectory(viscous_damping=0.0):\n", " # fill out this function so that the plot will update with a slider.\n", "\n", "plot_trajectory()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "sim_line = ax.plot(traj.index, traj.angle)[0]\n", "ax.set_ylim((-sys.coordinates['angle'], sys.coordinates['angle']))\n", "\n", "ax.set_ylabel('Angle [rad]')\n", "ax.set_xlabel('Time [s]')\n", "\n", "def plot_trajectory(viscous_damping=0.0):\n", " sys.constants['viscous_damping'] = viscous_damping\n", " traj = sys.free_response(5.0)\n", " sim_line.set_data(traj.index, traj.angle)\n", " fig.canvas.draw()\n", "\n", "plot_trajectory()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from ipywidgets import interact" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "widget = interact(plot_trajectory, damping=(0.0, 5.0, 0.1));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Visualize the Motion" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sys.constants['viscous_damping'] = 0.2\n", "frames_per_second = 30\n", "traj = sys.free_response(5.0, sample_rate=frames_per_second)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.animate_configuration(interval=1000 / frames_per_second) # sample time in milliseconds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**\n", "\n", "Try out some of the values for the viscous damping that had interesting trajectories and see what the animation looks like." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Oscillation Period and Viscous Damping\n", "\n", "You may have noticed that the period seems to change with different viscous damping values. It is worth investigating this more thoroughly." ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Use your function for estimating the period of a trajectory in a loop to collect period estimates for 30 values of viscous damping from 0.0 to 5.0. The code for the loop should be structured like:\n", "\n", "```python\n", "viscous_damping_vals = np.linspace(0.0, 5.0, num=30)\n", "periods = []\n", "for c in viscous_damping_vals:\n", " sys.constants['viscous_damping'] = c\n", " traj = sys.free_response(5.0)\n", " periods.append(estimate_period(traj.index, traj.angle))\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "from resonance.functions import estimate_period\n", "\n", "viscous_damping_vals = np.linspace(0.0, 5.0, num=30)\n", "periods = []\n", "for c in viscous_damping_vals:\n", " sys.constants['viscous_damping'] = c\n", " traj = sys.free_response(5.0)\n", " periods.append(estimate_period(traj.index, traj.angle))\n", " \n", "fig, ax = plt.subplots(1, 1)\n", "ax.set_xlabel('Damping Coefficient [Ns/m]')\n", "ax.set_ylabel('Period [s]')\n", "ax.plot(viscous_damping_vals, periods);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Have a look at the `periods` list and see if anything is unusual. Use the same loop as above but investigate viscoud damping values around the value that causes issues and see if you can determine how high the viscous damping value can be for a valid result." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "viscous_damping_vals = np.linspace(1.5, 1.7, num=50)\n", "periods = []\n", "for c in viscous_damping_vals:\n", " sys.constants['viscous_damping'] = c\n", " traj = sys.free_response(5.0)\n", " periods.append(estimate_period(traj.index, traj.angle))\n", "\n", "viscous_damping_vals[np.isnan(periods)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At any given configuration dynamic systems potentially have both potential energy and kinetic energy. Recall that the potential energy of a mass that is lifted above ground in a gravitational field is:\n", "\n", "$$ PE = mgh$$\n", "\n", "where $m$ is the mass of the object, $g$ is the acceleration due to gravity, and $h$ is the height above the ground. Additionally, the translational and rotational kinetic energy of a planar rigid body can be expressed as:\n", "\n", "$$ KE = \\frac{m v^2}{2} + \\frac{I \\omega^2}{2} $$\n", "\n", "where $m$ is the mass of the body, $v$ is the magnitude of the linear velocity of the center of mass, $I$ is the centroidal moment of inertia of the rigid body, and $\\omega$ is the angular velocity of the rigidbody.\n", "\n", "For example the bob of the pendulum has:\n", "\n", "$$ PE_{bob} = m_{bob} g l (1 - \\cos{\\theta}) $$\n", "\n", "and\n", "\n", "$$ KE_{bob} = \\frac{m_{bob} (l \\dot{\\theta})^2}{2} + \\frac{m_{bob} r^2 \\dot{\\theta}^2}{4} $$\n", "\n", "You can add a measurement for this that looks like:\n", "\n", "```python\n", "def kinetic_energy(bob_mass, bob_radius, rod_length, bob_height, rod_mass, angle_vel):\n", " \n", " v_bob = rod_length * angle_vel\n", " I_bob = bob_mass * bob_radius**2 / 2.0\n", " KE_bob = bob_mass * v_bob**2 / 2.0 + I_bob * angle_vel**2 / 2.0 \n", " \n", " v_rod =\n", " I_rod =\n", " KE_rod = \n", " \n", " return KE_rod + KE_bob\n", "```" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Add a measurement to the system that outputs the kinetic energy in the trajectory data frame from `free_response`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "def kinetic_energy(bob_mass, bob_radius, rod_length, bob_height, rod_mass, angle_vel):\n", " v_rod_com = rod_length / 2.0 * angle_vel\n", " v_bob = rod_length * angle_vel\n", " I_rod = rod_mass * rod_length**2 / 12 # rod's centroidal moment of inertia\n", " I_bob = bob_mass * bob_radius**2 / 2.0\n", " KE_rod = rod_mass * v_rod_com**2 / 2.0 + I_rod * angle_vel**2 / 2.0\n", " KE_bob = bob_mass * v_bob**2 / 2.0 + I_bob * angle_vel**2 / 2.0\n", " return KE_rod + KE_bob\n", "\n", "sys.add_measurement('kinetic_energy', kinetic_energy)\n", "\n", "def potential_energy(bob_mass, rod_mass, rod_length, bob_height, acc_due_to_gravity, angle):\n", " PE_bob = bob_mass * acc_due_to_gravity * (rod_length - rod_length * np.cos(angle))\n", " PE_rod = rod_mass * acc_due_to_gravity * (rod_length / 2 - rod_length / 2 * np.cos(angle))\n", " return PE_bob + PE_rod\n", "\n", "sys.add_measurement('potential_energy', potential_energy)\n", "\n", "def total_energy(potential_energy, kinetic_energy):\n", " return potential_energy + kinetic_energy\n", "\n", "sys.add_measurement('total_energy', total_energy)\n", "\n", "sys.measurements['kinetic_energy'], sys.measurements['potential_energy'], sys.measurements['total_energy']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Investigate what happens to the kinetic energy if the viscous damping value is increased." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "sys.constants['viscous_damping'] = 0.2\n", "traj = sys.free_response(5.0)\n", "plt.figure()\n", "traj.kinetic_energy.plot();" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Log Decrement\n", "\n", "The time constant of an exponential decaying oscillation, $\\tau$, is defined as:\n", "\n", "$$ x(t) = e^{t/\\tau} $$\n", "\n", "$\\tau$ is the time it takes for the function to decay to $1 - 1/e$ which is about 63.2% of the starting value.\n", "\n", "Another useful value that is related is the log decrement $\\delta$, which is the ratio of the period of osciallation to the time constant. It is defined as:\n", "\n", "$$\\delta = ln \\frac{x(t)}{x(t+T)} = T / \\tau $$" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Using a value of 0.2 for the viscous damping and 2 degrees as an initial angle, create a free response and determine the log decrement of decayed osciallation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "sys.constants['viscous_damping'] = 0.2\n", "sys.coordinates['angle'] = np.deg2rad(2.0)\n", "sys.speeds['angle_vel'] = 0.0\n", "traj = sys.free_response(10.0, sample_rate=200)\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "ax.plot(traj.index, np.rad2deg(traj.angle))\n", "ax.set_xlabel('Time [s]')\n", "ax.set_ylabel('Angle [deg]');\n", "\n", "t = 0.75\n", "T = estimate_period(traj.index, traj.angle)\n", "delta = np.log(traj.angle[t] / traj.angle[round(t + T, 3)])\n", "delta" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Underdamped\n", "\n", "When the viscous damping value is relatively low, a very nice decayed oscillation is present. This is called an underdamped oscillation because there is vibration, but it still dissipates." ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Create a single plot that shows the free response at these different viscoud damping values: 0.0, 0.08, 0.2, 1.6. Use an initial angle of 5 degrees and inital angular velocity of 0 deg/s. Include a legend so that it is clear which lines represent which values. Use a loop to reduce the amount of typing needed." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "for c in [0.0, 0.08, 0.2, 1.6]:\n", " sys.constants['viscous_damping'] = c\n", " traj = sys.free_response(5.0)\n", " ax.plot(traj.index, traj.angle, label='c = {:0.3f} [Ns/m]'.format(c))\n", " \n", "ax.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Critically Damped\n", "\n", "Above in the plot you created of the period of oscillation versus the viscous damping value, you should have discovered that when the value is at 1.6780754836568228 (or somewhere close) that there are no longer any oscillations. This boundary between oscillation and not oscillating wrt to the viscous damping value is called \"critically damped\" motion." ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Make a single plot of angle trajectories with the viscous damping value set to 1.6780754836568228 and an initial angle of 0.1 degrees. Plot three lines with the initial angular velocity as -10 degs, 0 degs, and 10 degrees." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "sys.constants['viscous_damping'] = 1.6780754836568228\n", "sys.coordinates['angle'] = np.deg2rad(0.1)\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "\n", "for v0 in [-np.deg2rad(10), 0, np.deg2rad(10)]:\n", " sys.speeds['angle_vel'] = v0\n", " traj = sys.free_response(5.0)\n", " ax.plot(traj.index, traj.angle, label='v0 = {:0.3f} [rad/s]'.format(v0))\n", " \n", "ax.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Overdamped\n", "\n", "Finally, if the viscous damping value is greater than the critical damping value, the motion is called over damped. There is no oscillation and with very high values of damping the system will rapidly decay." ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Create a single plot with the viscous damping value at 2.0 and these three sets of initial conditions:\n", "\n", "- $\\theta_0$ = 1 deg, $\\dot{\\theta}_0$ = 0 deg/s\n", "- $\\theta_0$ = 0 deg, $\\dot{\\theta}_0$ = 10 deg/s\n", "- $\\theta_0$ = -1 deg, $\\dot{\\theta}_0$ = 0 deg/s" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "sys.constants['viscous_damping'] = 2.0\n", "\n", "initial = ((np.deg2rad(1), 0.0),\n", " (0.0, np.deg2rad(10)),\n", " (-np.deg2rad(1), 0.0))\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "ax.set_xlabel('Time [s]')\n", "ax.set_ylabel('Angle [deg]')\n", "\n", "for x0, v0 in initial:\n", " sys.coordinates['angle'] = x0\n", " sys.speeds['angle_vel'] = v0\n", " traj = sys.free_response(5.0)\n", " lab_temp = 'x0 = {:0.1f} [deg], v0 = {:0.1f} [deg/s]'\n", " ax.plot(traj.index, np.rad2deg(traj.angle),\n", " label=lab_temp.format(np.rad2deg(x0), np.rad2deg(v0)))\n", " \n", "ax.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }