{
"metadata": {
"name": "mass_spring_damper"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we will derive the equations of motion for the classic mass, spring, damper system under the influence of gravity. The following figure gives a pictorial description of the problem."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from IPython.display import SVG\n",
"SVG(filename='mass_spring_damper.svg')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "pyout",
"prompt_number": 1,
"svg": [
""
],
"text": [
""
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Start by loading in the core functionality of both SymPy and Mechanics."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import sympy as sym\n",
"import sympy.physics.mechanics as me"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can make use of the pretty printing of our results by loading the SymPy printing extension."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%load_ext sympy.interactive.ipythonprinting"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll start by defining the variables will need for this problem. $x$ is the distance of the particle from the ceiling, $v$ is the speed of the particle, $m$ is the mass of the particle, $c$ is the damping coefficient of the damper, $k$ is the stiffness of the spring, $g$ is the acceleration due to gravity, and $t$ represents time."
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"x, v = me.dynamicsymbols('x v')"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"m, c, k, g, t = sym.symbols('m c k g t')"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we define a Newtonian reference frame that represents the ceiling which the particle is attached to, $C$."
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"ceiling = me.ReferenceFrame('C')"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will need two points, one to represent the original position of the particle which stays fixed in the ceiling frame, $o$, and the second one, $p$ which is aligned with the particle as it moves."
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"o = me.Point('o')\n",
"p = me.Point('p')"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The velocity of point $o$ in the ceiling is zero."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"o.set_vel(ceiling, 0)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Point $p$ can move downward in the $x$ direction and it's velocity is $v$ in the downward direction."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"p.set_pos(o, x * ceiling.x)\n",
"p.set_vel(ceiling, v * ceiling.x)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are three forces acting on the particle. Those due to the acceleration of gravity, the damper, and the spring."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"damping = -c * p.vel(ceiling)\n",
"stiffness = -k * p.pos_from(o)\n",
"gravity = m * g * ceiling.x\n",
"forces = damping + stiffness + gravity\n",
"forces"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "pyout",
"prompt_number": 10,
"text": [
"(-c\u22c5v + g\u22c5m - k\u22c5x)*\u001b[94m\u001b[1mc_x\u001b[0;0m\u001b[0;0m"
]
}
],
"prompt_number": 10
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use Newton's second law, $0=F-ma$, to form the equation of motion of the system."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"zero = me.dot(forces - m * p.acc(ceiling), ceiling.x)\n",
"zero"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$- c \\operatorname{v}{\\left (t \\right )} + g m - k \\operatorname{x}{\\left (t \\right )} - m \\frac{\\partial}{\\partial t} \\operatorname{v}{\\left (t \\right )}$$"
],
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAOMAAAAfCAYAAAAV3nouAAAABHNCSVQICAgIfAhkiAAABx9JREFU\neJzt3HuwVVUdwPHP5QqBJN0UcgSLwjErR5hMM0nhoviICqeaYpoyMTJzpobsIb0cmCYrEk2nh9pk\n3l5IrymnEgdqoqwmxkqimp5i00uHdDSyqEDoj9/anX3PPcd77rmbcy6X9Z25c1hr7732b//Wb6/1\n+/3W2vTKZA5d5uBKnI7f4h/dFSeTOXR5AqZhCX6DSd0VJ5PJwFewtJsCTOjmzTOZMcSvcFw3BRir\nL+NMTB7mnOl4YgdkOdhYgB9iP75acdvjqV8m43pswudwiurc1HGjp2l4U4P64/BXzE7lHrwXEzsk\n18HERPwbb62wzfHWL1/DgJB3IvbiwgraHVd6Wo3HN6h/szCw8ogzB1e0eZ+5OKzNa8c688XM+NwK\n2+xUv3SCBXhEzE5wtNDX0ypou209jTU3dYYQ9pEGx87EVvFABTvwDO25F2/BsW1cdzDQj3/ipxW1\n18l+6QSnC5kfSOWX4k7cO8p2R6WnsfYynifinUacge81qP+5eNBMjX78SLheVTDe+uUBMVjBLLxb\nNTP5qPRUtZt2kohTdmE3/oP3i9HgeKzEU0XAvCFdsxwfxAnCvVpbam8ZLhaB7pOwCHfg6/hYOmcb\nTsO3K36WRkzCu8Sz/A2/xjnYjPuwGPNwEY7Ey4X7Mx/rsBGX46j0PJPE8+2pUMaJ6X5Xp/LjhOvU\nm2R6fZIJnoPX4FHRL6/DpegTRroa9+h+v1ygWt1+MT3DerHWuFwMXs1oxXb/rvt6+j+vEAHqvFRe\nIHY0XCBm4BvFy78S20vXvVYo9Xki+ze1QduXihd7SoNjc/HRNuQdEMptlSnYgi+U6t4pZp9luC7V\n3YXvi0GpJ9Wtwk5cU7pnr9DPRSOUeziKePEsYZBXCUO4Isl6VDpvjjCIwjsaELtQ5uP52KeWAOpk\nv9QzSXd126rtMko9VeWmPhufxTvws1TXJzr3JzhXjFx7cX6qL/iUmv8+NQldzyLREbsbHNutM7HJ\nWjwLl5Tq7hY63CeMpEcY+f3COIoZaI8YzdfjD6nu0fR3dMVy9qf73S3cr3XCWKeJUf3BdN7bRH/t\nS+WpeEi4WX/EteIFLY51q18W6q5uW7Vdxoj9bhYd2Nvk+CzhPh0rlPTiuuMfF8q+VW3kLnO/GOEb\ncYpwhUfKgNZnxieLTr+urv59+CWOEaPeXGEk9bHSrYbGEnPSuS9sVeAW2SSMcp2YEZtRnzn8s+Y6\n7mS/1NNt3bZqu4UsbeupHDOehE+UGh6ObXiDyCCdjRuSsI34S/q9UMSTG0vHesSIvB+/FwumD5aO\nnyhGuO82aXtmuq4Zn1Zzncs8RaT+/9vg2Aoxoxe8TOhqY915C0VQfl8qnyVGuq115/Xjlrq680Us\nXf9c7fYDtXjxI8JDWYubNI6HypnDE4TRfafJParol3afq0rdtkOrtsuBsd8RcVoS5uIWzv2WobtC\nFquNNotEgqHMG8WsVKzd9ImZquADYiQcKQNanxlvFM/YV6qbLDr8laW62wwNxJ+Zrj23rn6zmi5m\na+5VjIQiXjw5lZeJNHuxttXX6CJcJtyrw0t15a1hneyXZoxWt/vb/CsYznYZpZ6qiBl3pt9Gn5/M\nxgtK5Vn4Xd05L8Lt6d93CnekzJlilCzWblaqpex7hFuwY8RSj4yHxfM9XKpbKDKVxYg3QSStttRd\nu0h0xg9KdUem+s+n8irNvYqR0C9G722pvEfEMUek8rXpdwo+JGYrIiO8Hf9K5Ql4e6ndbvdLFbrt\nafOvYDjbZZR6quJlvFfEKQsaCLE6CViwXUzLBSvEiFcY4l4RrJ9cOqdXLTA/VbgqheuyVATuB5ov\ni1mjiMFmC7f8HpFBJpJYfRobzF1q61rEjNwr9HaqiDuroF/or0jKPJR+d4kMabHOtUS8bCeKRefj\nDU48vEe49wXd7pcqdTtHDEprDLbF4RjOdum+nhCKugU3i/Wt65Ow9S/7TDHV34AP4yVN2luh5lrN\nEyPfNSI9X7Q5w2AXcaQMGNnSxiX4hphR1uDHIptWsBS/MDQztkWktsv04kvC/V2juqz2VryqVJ6A\nT4pM91Wl+0wX/XW1eJ6p4uW7ScSb5zRpvxP90ogqdftY3zA+XezOaUSrtkv39HTQMmBkL2OZw0Xi\nZ3lFsmS6Q/03jDcbvHTVUcbadrhOssvgfYLNmGHoR6dLhL+/qWqhMh2l/hvGxXKfjmk2iFmwyDQe\nI2LFy7smUaZd6r9hvEMkeJYI9/xPYiPEGd0Qbrx+QlQlt4kY40qxjjdbpKvr1xwzY58NIiN+nrD9\n3SLreruICSeKvaaZTOYAMtw3jJ8R67KZTOYAs8rgTQOXGfxJ0w6xdDVB4y1tB5xDOYGTObR4rG8Y\np4u8wE682vD/f00mkxkFR4jEzXp8U2ROCw4TGzuWi33WmUwmk8lkMplMJpPJZDKZzFD+B5R/CE33\nm63mAAAAAElFTkSuQmCC\n",
"prompt_number": 11,
"text": [
" d \n",
"-c\u22c5v(t) + g\u22c5m - k\u22c5x(t) - m\u22c5\u2500\u2500(v(t))\n",
" dt "
]
}
],
"prompt_number": 11
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can then form the first order equations of motion by solving for $\\frac{dv}{dt}$ and introducing the kinematical differential equation, $v=\\frac{dx}{dt}$."
]
},
{
"cell_type": "code",
"collapsed": true,
"input": [
"dvdt = sym.solve(zero, v.diff(t))[0]\n",
"dxdt = v\n",
"dvdt, dxdt"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$\\begin{pmatrix}\\frac{- c \\operatorname{v}{\\left (t \\right )} + g m - k \\operatorname{x}{\\left (t \\right )}}{m}, & \\operatorname{v}{\\left (t \\right )}\\end{pmatrix}$$"
],
"output_type": "pyout",
"prompt_number": 12,
"text": [
"\u239b-c\u22c5v(t) + g\u22c5m - k\u22c5x(t) \u239e\n",
"\u239c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500, v(t)\u239f\n",
"\u239d m \u23a0"
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Forming the equations of motion can also be done with the automated methods available in the Mechanics package: `LagrangesMethod` and `KanesMethod`. Here we will make use of Kane's method to find the same equations of motion which we found manually above. First define a particle which represents the mass attached to the damper and spring."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"mass = me.Particle('mass', p, m)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 13
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can construct a `KanesMethod` object by passing in the generalized coordinate, $x$, the generalized speed, $v$, and the kinematical differential equation which relates the two, $0=v-\\frac{dx}{dt}$."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"kane = me.KanesMethod(ceiling, q_ind=[x], u_ind=[v], kd_eqs=[v - x.diff(t)])"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 14
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now Kane's equations can be computed which returns $F_r$ and $F_r^*$."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"kane.kanes_equations([(p, forces)], [mass])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$\\begin{pmatrix}\\left[\\begin{smallmatrix}c \\operatorname{v}{\\left (t \\right )} - g m + k \\operatorname{x}{\\left (t \\right )}\\end{smallmatrix}\\right], & \\left[\\begin{smallmatrix}m \\frac{\\partial}{\\partial t} \\operatorname{v}{\\left (t \\right )}\\end{smallmatrix}\\right]\\end{pmatrix}$$"
],
"output_type": "pyout",
"prompt_number": 15,
"text": [
"\u239b \u23a1 d \u23a4\u239e\n",
"\u239c[c\u22c5v(t) - g\u22c5m + k\u22c5x(t)], \u23a2m\u22c5\u2500\u2500(v(t))\u23a5\u239f\n",
"\u239d \u23a3 dt \u23a6\u23a0"
]
}
],
"prompt_number": 15
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equations are also available in the form $M\\frac{dq,u}{dt}=f(\\dot{u}, u, q)$ and we can extract the mass matrix, $M$, and the forcing equations, $f$."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M = kane.mass_matrix_full\n",
"f = kane.forcing_full\n",
"M, f"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$\\begin{pmatrix}\\left[\\begin{smallmatrix}1 & 0\\\\0 & m\\end{smallmatrix}\\right], & \\left[\\begin{smallmatrix}\\operatorname{v}{\\left (t \\right )}\\\\- c \\operatorname{v}{\\left (t \\right )} + g m - k \\operatorname{x}{\\left (t \\right )}\\end{smallmatrix}\\right]\\end{pmatrix}$$"
],
"output_type": "pyout",
"prompt_number": 16,
"text": [
"\u239b\u23a11 0\u23a4, \u23a1 v(t) \u23a4\u239e\n",
"\u239c\u23a2 \u23a5 \u23a2 \u23a5\u239f\n",
"\u239d\u23a30 m\u23a6 \u23a3-c\u22c5v(t) + g\u22c5m - k\u22c5x(t)\u23a6\u23a0"
]
}
],
"prompt_number": 16
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we can form the first order differential equations of motion $\\frac{dq,u}{dt}=M^{-1}f(\\dot{u}, u, q)$, which is the same as previously found."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"M.inv() * f"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$\\left[\\begin{smallmatrix}\\operatorname{v}{\\left (t \\right )}\\\\\\frac{- c \\operatorname{v}{\\left (t \\right )} + g m - k \\operatorname{x}{\\left (t \\right )}}{m}\\end{smallmatrix}\\right]$$"
],
"output_type": "pyout",
"prompt_number": 17,
"text": [
"\u23a1 v(t) \u23a4\n",
"\u23a2 \u23a5\n",
"\u23a2-c\u22c5v(t) + g\u22c5m - k\u22c5x(t)\u23a5\n",
"\u23a2\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23a5\n",
"\u23a3 m \u23a6"
]
}
],
"prompt_number": 17
}
],
"metadata": {}
}
]
}