Mass Spring System:

In this example, we are going to study a simple mass-spring system. This example serves to demonstrate the visualization handling abilities of mechanics-visualization module which is still in development process

Dependencies:

For this system, following python modules are required:

1) sympy=0.7.2

2) numpy=1.6.2



In [49]:

    
%load_ext autoreload
%autoreload 2
from sympy.physics.mechanics import *
from sympy import symbols,cos,Eq,sqrt,sin
import viz                       #This module is included with the notebook
from IPython.display import *

That is all we are going to need to carry out this visualization.

Now, we will create symbols for coordinates, speed, force etc of the mass(bob).



In [50]:

    
q = dynamicsymbols('q')                  # Generalized coordinate
u = dynamicsymbols('u')                  # Generalized speed
f = dynamicsymbols('f')                  # Net Force applied to the box
m = symbols('m')                         # Mass of body
g, k, t, L = symbols('g k t L')               # Gravity,spring constant, time and length of unstretched spring

Now we create an inertial reference frame, assign an origin, and set origin velocity to zero.



In [51]:

    
I = ReferenceFrame('I')            # Inertial reference frame
O = Point('O')                     # Origin point
O.set_vel(I, 0)                    # Origin's velocity is zero

Now we define the bob as a Particle ..



In [52]:

    
P = Point('P')                      # Center Of Mass of the box
P.set_pos(O, q * I.y)               # Set the position of P
P.set_vel(I, u * I.y )              # Set the velocity of P
bob = Particle('bob',P,m)

We assign forces acting on the Particle. Coordinate axes are defined to be as : x-axis: running from left to right, y-axis: running from top to bottom, z-axis: running from behind the plane to above the plane.

Following the above notion, Gravity is taken to be in (+I.y) direction, and force due to spring(F=Kx), in -I.y direction ..



In [53]:

    
f = m*g* I.y - k*(q)*I.y    #Net force on bob.
force = [(P,f)]

Setting up equation of motion for the body:



In [54]:

    
diff_eq = [q.diff(t) - u]

Now, we call the Kane's Method to generate the equation of motion of the bob



In [55]:

    
kane = KanesMethod(I, q_ind=[q], u_ind=[u], kd_eqs=diff_eq) # Initialize the object
fr, frstar = kane.kanes_equations(force, [bob])

The Differential equation of the motion of bob is given by the expression: fr + frstar = 0, i.e. fr = frstar



In [56]:

    
fr, frstar









    Out[56]:





([-g*m + k*q(t)], [m*Derivative(u(t), t)])

Now, We need to solve this differential equation, for this we would be utilizing code generator sub-module in sympy mechanics, which is under development, and is not merged in the package.

First packing the parameters in a list.



In [57]:

    
parameters = [m, k, g, L]     #Mass, stiffness, gravity
parameter_vals = [1, 1, 9.8, 1]

Now we import the code generator module.



In [58]:

    
import code_gen
right_hand_side = code_gen.numeric_right_hand_side(kane, parameters)

The numeric_right_hand_side method takes the Kane's object and parameters in symbolic form, and returns the differential equation in a numeric form, which can be integrated using odeint() method in SciPy.

Now we integrate right_hand_side using odeint



In [65]:

    
from scipy.integrate import odeint
x0 = [1,0]        #Initial conditions, q=1,u=0.
t = [i*0.1 for i in range(0,10)]    #Taking 10 time intervals of 0.1 sec
numeric_vals = odeint(right_hand_side, x0, t, args=(parameter_vals,)) 
numeric_vals









    Out[65]:





array([[ 1.        ,  0.        ],
       [ 1.04396334,  0.87853405],
       [ 1.17541411,  1.7482901 ],
       [ 1.39303891,  2.6005778 ],
       [ 1.69466326,  3.42688139],
       [ 2.07727345,  4.21894471],
       [ 2.53704657,  4.96885373],
       [ 3.06938873,  5.66911561],
       [ 3.66898093,  6.31273356],
       [ 4.32983223,  6.89327675]])

The variable contains the positions and velocities, we need only positions for this 1D simulation, so we separate them, and pack them alongwith time.



In [63]:

    
Q = []
for entries in numeric_vals:
    Q.append(entries[0])
 
simulation_coords = [Q,t]     #Packing positions and time in a single array

Simulation:

For simulation and visualization, we are going to utilize, the Javascript rendering capabilities of IPython notebooks to render WebGL's. We are using Three.js Javascript library for 3D visualizations.



In [64]:

    
viz.visualize(simulation_coords)