#Defining the Problem

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.


```
from IPython.display import SVG
SVG(filename='mass_spring_damper.svg')
```

Start by loading in the core functionality of both SymPy and Mechanics.


```
import sympy as sym
import sympy.physics.mechanics as me
```

We can make use of the pretty printing of our results by loading SymPy's printing extension.


```
from sympy.physics.vector import init_vprinting
init_vprinting()
```

We'll start by defining the variables we will need for this problem:
- $x$: distance of the particle from the ceiling
- $v$: speed of the particle
- $m$: mass of the particle
- $c$: damping coefficient of the damper
- $k$: stiffness of the spring
- $g$: acceleration due to gravity
- $t$: time


```
x, v = me.dynamicsymbols('x v')
```


```
m, c, k, g, t = sym.symbols('m c k g t')
```

Now, we define a Newtonian reference frame that represents the ceiling which the particle is attached to, $C$.


```
ceiling = me.ReferenceFrame('C')
```

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.


```
O = me.Point('O')
P = me.Point('P')
```

The velocity of point $O$ in the ceiling is zero.


```
O.set_vel(ceiling, 0)
```

Point $P$ can move downward in the $y$ direction and its velocity is $v$ in the downward direction.


```
P.set_pos(O, x * ceiling.x)
P.set_vel(ceiling, v * ceiling.x)
P.vel(ceiling)
```

There are three forces acting on the particle. Those due to the acceleration of gravity, the damper, and the spring.


```
damping = -c * P.vel(ceiling)
stiffness = -k * P.pos_from(O)
gravity = m * g * ceiling.x
forces = damping + stiffness + gravity
forces
```

Now we can use Newton's second law, $0=F-ma$, to form the equation of motion of the system.


```
zero = me.dot(forces - m * P.acc(ceiling), ceiling.x)
zero
```

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}$.


```
dv_by_dt = sym.solve(zero, v.diff(t))[0]
dx_by_dt = v
dv_by_dt, dx_by_dt
```

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 that we found manually above. First, define a particle that represents the mass attached to the damper and spring.


```
mass = me.Particle('mass', P, m)
```

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}$.


```
kane = me.KanesMethod(ceiling, q_ind=[x], u_ind=[v], kd_eqs=[v - x.diff(t)])
```

Now Kane's equations can be computed, and we can obtain $F_r$ and $F_r^*$.


```
fr, frstar = kane.kanes_equations([(P, forces)], [mass])
print fr
print frstar
```

The equations are also available in the form $M\frac{d}{dt}[q,u]^T=f(\dot{u}, u, q)$ and we can extract the mass matrix, $M$, and the forcing equations, $f$.


```
M = kane.mass_matrix_full
f = kane.forcing_full
M, f
```

Finally, we can form the first order differential equations of motion $\frac{d}{dt}[q,u]^T=M^{-1}f(\dot{u}, u, q)$, which is the same as previously found.


```
M.inv() * f
```

#Simulating the system#

Now that we have defined the mass-spring-damper system, we are going to simulate it.

PyDy's `System` is a wrapper that holds the Kanes object to integrate the equations of motion using numerical values of constants.


```
from pydy.system import System
```


```
sys = System(kane)    #Initializing System with the KanesMethod's object
```

Now, we specify the numerical values of the constants and the initial values of states in the form of a dict.


```
constants = {m:10, g:9.8, c:5, k:10}    
init_conditions = {x:0, v:0}

sys.constants = constants
sys.initial_conditions = init_conditions
```

We must generate a time vector over which the integration will be carried out.


```
import numpy as np
sys.times = np.linspace(0, 8, 100)
```


```
x0 = sys.integrate()
```

#Visualizing the System

PyDy has a native module `pydy.viz` which is used to visualize a System in a GUI.


```
from pydy.viz import *
from numpy import pi
```

For visualizing the system, we need to create shapes for the objects we wish to visualize, and map each of them 
to a VisualizationFrame, which holds the position and orientation of the object.


```
bob = Sphere(2, color="red", material="metal")
bob_vframe = VisualizationFrame(ceiling, P, bob)
```


```
ceil = Circle(radius=10, color="white", material="metal")
#Creating a new reference frame, to give correct orientation to ceiling.
ceil_reference_frame = ceiling.orientnew("ceil_new", 'Axis', [pi/2, ceiling.z])
ceil_vframe = VisualizationFrame(ceil_reference_frame, O, ceil)
```

Now we initialize a Scene. A Scene contains all the information required to visualize a System onto a canvas.
It takes a ReferenceFrame and Point as arguments.


```
scene = Scene(ceiling, O)
```

We provide the VisualizationFrames, which we want to visualize as a list to scene.


```
scene.visualization_frames = [bob_vframe, ceil_vframe]

```

The default camera of Scene has the z axis of the base frame pointing out of the screen, and the y axis pointing up. We want the x axis to point downwards, so we supply a new camera that will achieve this.


```
primary_camera = PerspectiveCamera(ceiling.orientnew('k','Axis', [pi/2, ceiling.z] ), O.locatenew('a',100*ceil_reference_frame.z))
scene.cameras = [primary_camera]
```

The `generate_visualization_json_system` method generates the required data for the animations, in the form of JSON files. These JSON files are needed before calling the display methods.


```
scene.generate_visualization_json_system(sys)
```

Now, we call the display method.


```
scene.display_ipython()
```


```

```

11fbe3c8a41d2a05948b0ea6cacb2bcd424561bc
ipynb
Jupyter Notebook
examples/mass_spring_damper/mass_spring_damper.ipynb
jellysheep/pydy
687480067a27bfd250c21a0318c049681b8c076b
BSD-3-Clause
2019-06-27T05:30:36.000Z
2019-06-27T05:30:36.000Z
examples/mass_spring_damper/mass_spring_damper.ipynb
Munyola/pydy
687480067a27bfd250c21a0318c049681b8c076b
BSD-3-Clause
examples/mass_spring_damper/mass_spring_damper.ipynb
Munyola/pydy
687480067a27bfd250c21a0318c049681b8c076b
BSD-3-Clause
2019-06-27T05:29:50.000Z
2019-06-27T05:29:50.000Z
Qwen/Qwen-72B
1. YES
2. YES