{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 6. Sinusoidal Forcing Of A Spring-Mass-Damper System\n", "\n", "\n", "This notebook introduces external forcing of a vibratory system, where the\n", "external force is modeled as a sinusoidal input to the bottom of a bus driver's\n", "seat.\n", "\n", "After the completion of this assignment students will be able to:\n", "\n", "- excite a system with a sinusoidal input\n", "- understand the difference in transient and steady state solutions\n", "- use autocorrelation to determine period\n", "- relate the frequency response to the time series\n", "- create a frequency response plot\n", "- define resonance and determine the parameters that cause resonance\n", "\n", "# Introduction\n", "\n", "We have only been dealing with the free response of a system given some initial coordinate (position, orientation) or speed value (velocity, angular velocity). The behavior observed is the \"natural\" behavior of the system, i.e. the motion that it would exhibit if no external loads act on the system. For the simple single degree of freedom systems we have dealt with so far it is possible to apply a time varying external force or torque to the system to influence the time trajectory of the single coordinate. For example let's investigate the simplest system in vibrations, a particle that can slide laterally but is fixed to a wall via a linear spring and linear damper and has a sinusoidal force applied to it as it moves:\n", "\n", "![](fig/06/Mass-Spring-Damper.svg)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from resonance.linear_systems import MassSpringDamperSystem" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sys = MassSpringDamperSystem()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mass, $m$, stiffness of the spring, $k$, and the damping coefficient of the dashpot, $c$, can be set on this model." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.constants" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "# Undamped Forced Motion\n", "\n", "Each system has a function called `sinusoidal_forcing_response()` that works in a similar fashion to `free_response()`. Two new pieces of information are required: the forcing amplitude, $F_o$, and the forcing frequency, $\\omega$ as the first two arguments. Set $F_o=.0$ N and $\\omega=2\\pi$ rad/s and generate a trajectory for 20 seconds." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "solution2": "hidden" }, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(1.0, 2 * np.pi, 20.0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that there is no longer a simple sinusoidal oscillation. The trajactory seems that it may still be *periodic* but the motion is now some combination of the natural motion and the motion due to the forcing term. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj.plot(subplots=True);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "traj = sys.sinusoidal_forcing_response(0.0, 0.0, 20.0)\n", "lines = ax.plot(traj.index, traj.position)\n", "ax.set_xlabel('Time [s]')\n", "ax.set_ylabel('Position [m]')\n", "ax.set_ylim((-0.5, 0.5))\n", "\n", "def adjust_forcing(amplitude=0.0, frequency=0.0):\n", " traj = sys.sinusoidal_forcing_response(amplitude, frequency, 20.0)\n", " lines[0].set_data(traj.index, traj.position)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ipywidgets import interact\n", "\n", "interact(adjust_forcing, amplitude=(0.0, 10.0, 0.2), frequency=(0.0, 10 * np.pi, 0.3));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Forcing at and around the natural frequency\n", "\n", "You may have notice some particularly interesting motion in a particular regime of frequencies. It turns out that if you force the system at a frequency near the natural frequency of the system. Simulate the free response of the system and find the natural frequency in radians per second. Investigate the forced trajectory with a 1 Newton forcing amplitude and a forcing frequency slightly less than, slightly more than, and exactly at the natural frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "sys.coordinates['position'] = 0.1 # m\n", "traj = sys.free_response(5.0)\n", "traj.plot(subplots=True);" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Estimate the natural frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "from resonance.functions import estimate_period\n", "2 * np.pi / estimate_period(traj.index, traj.position)" ] }, { "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": [ "# Just lower than the natural frequency" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "sys.coordinates['position'] = 0.0 # m\n", "\n", "traj = sys.sinusoidal_forcing_response(1.0, 9.8, 40.0)\n", "\n", "traj.plot(subplots=True);" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "# Just higher than the natural frequency" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(1.0, 10.2, 40.0)\n", "traj.plot(subplots=True);" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "This is called *beating*. Beating signals seem to have two primary frequencies, one that is much slower than the other. There is a rapid frequency that is close to the forced frequency and then another that can be estimated with:\n", "\n", "$$\\omega_{beat} = |\\omega_n - \\omega|$$\n", "\n", "**Exercise**\n", "\n", "Calculate the beat frequency and beat period from the above equation and see if it matches the trajectory in the above plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "w_beat = np.abs(9.9959766250584341 - 10.2)\n", "print(w_beat)\n", "\n", "print(2 * np.pi / w_beat)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write you answer here" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "# Exactly at the natural frequency\n", "\n", "See what happens at exactly the natural frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(1.0, 9.9959766250584341, 100.0)\n", "\n", "traj.plot(subplots=True);" ] }, { "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": [ "# Damping\n", "\n", "There is also interesting vibrational behavior when both damping and forcing are introduced to the system. Set $m=100$ kg, $c=100$ kg/s, and $k=910$ N/m and simulate the system's free response for 20 seconds to remind yourself of the behavior. Use an initial position of 0.001 m and initial velocity of 0.02 m/s. What type of vibration is this? (over-, under-, critically-damped, etc)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "solution2": "hidden" }, "outputs": [], "source": [ "sys.constants['mass'] = 100 # kg\n", "sys.constants['damping'] = 100 # kg/s\n", "sys.constants['stiffness'] = 910 # N/m\n", "\n", "sys.coordinates['position'] = 0.001 # m\n", "sys.speeds['velocity'] = 0.02 # m/s" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "solution2": "hidden" }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj = sys.free_response(20)\n", "traj.plot(subplots=True);" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "# Natural Frequency\n", "\n", "The exact value of a linear system's natural frequency can be calculate with:\n", "\n", "$$\\omega_n = \\sqrt{\\frac{k}{m}}$$\n", "\n", "where $m$ is the system mass and $k$ is the system's stiffness. Calculate this value using the equation above." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "np.sqrt(sys.constants['stiffness'] / sys.constants['mass'])" ] }, { "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": [ "# Forcing a damped system below it's natural frequency\n", "\n", "Subject the system to $F(t) = 10 \\cos(1t)$ with the same initial conditions as the previous case. This frequency is well below the natural frequency. Simulate for 5, 10, 20, and 50 seconds." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(10.0, 1.0, 5)\n", "traj.plot(subplots=True);" ] }, { "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": [ "# Just below the natural frequency\n", "\n", "Subject the system to $F(t) = 10 \\cos(3t)$ with the same initial conditions as the previous case. This frequency is well below the natural frequency. Simulate for 5, 10, 20, and 50 seconds." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(10.0, 3.0, 5)\n", "traj.plot(subplots=True);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Higher than natural frequency\n", "\n", "Subject the system to $F(t) = 10 \\cos(10t)$ with the same initial conditions as the previous case. This frequency is well below the natural frequency. Simulate for 5, 10, 20, and 50 seconds." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj = sys.sinusoidal_forcing_response(10.0, 10.0, 5)\n", "traj.plot(subplots=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Excersise**\n", "\n", "What do you observe about how the vibration behavior changes with respect to the forcing amplitude and frequency? How does the amplitude and frequency of the position trajectory change with respect to the forcing function? Create a interactive plot where you can adjust the forcing amplitude and frequency and observe the change in position over a 30 second simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Transient and Steady State Behavior\n", "\n", "When you increase the simulation time long enough in all of the above forced responses, you should see that there is more unpredictable behavior at the beginning of the response and that in the later portion of the response the system gets into a very predictable sinusoidal motion. What is the frequency of the motion after about 10 seconds of simulation? Do you recognize this frequency?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The underdamped regime of motion for this system occurs with a damping coefficient value as about $60.33 < c < 603.32$ kg/s." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy.optimize import curve_fit" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def cosine_func(times, amp, freq, phase_angle):\n", " return amp * np.cos(freq * times - phase_angle)" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Using the curve fit tool and a fitting function that looks like $X\\cos(\\omega t - \\theta)$ to fit the steady state behavior of the position trajectory ($t > 10$s or so, use `traj[m:]` to get just the last `m` seconds). With the damping ratio set at 100 kg/s make a plot of the position amplitude as a function of forcing frequences $1.0 < \\omega < 10.0$ rad/s. Use a forcing amplitude of 10 N. Plot a vertical line using `ax.axvline()`.\n", "\n", "```python\n", "Fo = 10.0 # N\n", "\n", "frequencies = np.linspace(1.0, 10.0, num=100)\n", " \n", "sys.constants['damping'] = 100 # kg/s\n", " \n", "amplitudes = []\n", " \n", "for omega in frequencies:\n", " # write the code for the loop here\n", "\n", "fig, ax = plt.subplots(1, 1, sharex=True)\n", "ax.set_xlabel('$\\omega$ [rad/s]')\n", "ax.set_ylabel('Steady state amplitude [m]') \n", "\n", "# write the plot commands here\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "solution2": "hidden" }, "outputs": [], "source": [ "Fo = 10.0 # N\n", "\n", "frequencies = np.linspace(1.0, 10.0, num=100)\n", " \n", "sys.constants['damping'] = 100 # kg/s\n", " \n", "amplitudes = []\n", " \n", "for omega in frequencies:\n", " traj = sys.sinusoidal_forcing_response(Fo, omega, 20.0)\n", " popt, pcov = curve_fit(cosine_func,\n", " traj[10:].index, traj[10:].position,\n", " p0=(Fo, omega, 0.05))\n", " amplitudes.append(abs(popt[0]))\n", "\n", "amplitudes = np.array(amplitudes)\n", "\n", "fig, ax = plt.subplots(1, 1, sharex=True)\n", "ax.set_xlabel('$\\omega$ [rad/s]')\n", "ax.set_ylabel('Steady state amplitude [m]') \n", "\n", "ax.axvline(np.sqrt(sys.constants['stiffness'] / sys.constants['mass']), color='black')\n", "ax.plot(frequencies, amplitudes);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write you answer here" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "**Exercise**\n", "\n", "Now let's plot a new line on the plot for different values of damping coefficients, instead of a single value. You can use a *nested loop* to cycle through frequencies for each damping ratio. Use 5 equally spaced values of $60 < c < 600$ and the same parameters as above. Add a legend to show which color line corresponds to what damping coefficient.\n", "\n", "```python\n", "Fo = 10.0 # N\n", "\n", "dampings = np.linspace(60.0, 600.0, num=5)\n", "frequencies = np.linspace(1.0, 10.0, num=100)\n", "results = []\n", "\n", "for c in dampings:\n", " \n", " # set the damping\n", " \n", " amplitudes = []\n", " \n", " for omega in frequencies:\n", " # write this portion\n", " \n", " amplitudes = np.array(amplitudes)\n", " \n", " # store the amplitudes in the results list\n", "\n", "fig, ax = plt.subplots(1, 1, sharex=True)\n", "ax.set_xlabel('$\\omega$ [rad/s]')\n", "ax.set_ylabel('Steady state amplitude [m]') \n", "\n", "ax.axvline(np.sqrt(sys.constants['stiffness'] / sys.constants['mass']), color='black')\n", "\n", "for amps in results:\n", " ax.plot(frequencies, amps)\n", "\n", "# add legend\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "solution2": "hidden" }, "outputs": [], "source": [ "Fo = 10.0 # N\n", "\n", "dampings = np.linspace(60.0, 600.0, num=5)\n", "frequencies = np.linspace(1.0, 10.0, num=100)\n", "results = []\n", "\n", "for c in dampings:\n", " \n", " sys.constants['damping'] = c\n", " \n", " amplitudes = []\n", " \n", " for omega in frequencies:\n", " traj = sys.sinusoidal_forcing_response(Fo, omega, 20.0)\n", " popt, pcov = curve_fit(cosine_func,\n", " traj[10:].index, traj[10:].position,\n", " p0=(Fo, omega, 0.05))\n", " amplitudes.append(np.abs(popt[0]))\n", " \n", " amplitudes = np.array(amplitudes)\n", " \n", " results.append(amplitudes)\n", "\n", "fig, ax = plt.subplots(1, 1, sharex=True)\n", "ax.set_xlabel('$\\omega$ [rad/s]')\n", "ax.set_ylabel('Steady state amplitude [m]') \n", "\n", "ax.axvline(np.sqrt(sys.constants['stiffness'] / sys.constants['mass']), color='black')\n", "\n", "for amps in results:\n", " ax.plot(frequencies, amps)\n", "\n", "ax.legend(dampings);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The frequency response\n", "\n", "The plot you created above is called a *frequency response plot* it shows how the magnitude of the steady state response to a sinusoidal forcing. The X axis is the forcing frequency and the Y axis is the amplitude of the position oscillation. The frequency response plot also typically includes a plot of the position's steady state phase shift angle as a function of the forcing frequency too. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sys.constants['damping'] = 60\n", "sys.frequency_response_plot(10.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exerice**\n", "\n", "Create an interactive plot using the `sys.frequency_response_plot()` function with a slider for the damping coefficient that goes from $60 < c < 600$ with an appropriate interval spacing. You can access the lines in the plot like so:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "axes = sys.frequency_response_plot(10.0)\n", "amp_line = axes[0].lines[1]\n", "phase_line = axes[1].lines[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can get the frequencies in the plot from:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "frequencies = amp_line.get_xdata()" ] }, { "cell_type": "markdown", "metadata": { "solution2": "hidden", "solution2_first": true }, "source": [ "Write a plot update function that looks like:\n", "\n", "```python\n", "def update(damping):\n", " # code that modifies the plot\n", "```\n", "\n", "You can use `sys.frequency_response()` to get the amplitude and phase values." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "solution2": "hidden" }, "outputs": [], "source": [ "def update(damping):\n", " sys.constants['damping'] = damping\n", " amp, phase = sys.frequency_response(frequencies, 10.0)\n", " amp_line.set_ydata(amp)\n", " phase_line.set_ydata(np.rad2deg(phase))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# write answer here" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interact(update, damping=(60, 600))" ] } ], "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 }