This notebook simulates particles that are coupled together by springs that obey Hooke's law $F=-k x$. It allows the user to place particles with different initial conditions and to couple together particles of their choosing with springs. The notebook will then produce a simulation of the desired situation. Gravity can be toggled on and off and so can collisions of particles. A particle may also be fixed at its initial position to act as an anchor for other particles.
The cell below imports all of the required libraries and sets up variables that are used throughout the notebook. Some of which may be changed including the time step of the simulation, the total simulation run time, default particle radius, the acceleration due to gravity and whether springs should be visible or not.
#NAME: Coupled Points
#DESCRIPTION: Interactive 3D simulation of systems of particles and springs.
from math import sqrt
from vpython import *
from __future__ import division, print_function
from ipywidgets import widgets
#should not be changed
points=[]
springs=[]
GRAVITY=False
COLLISIONS=False
flag=0 #start 1, end -1, or 0 if not attached to spring
# can be changed
DELTAt=0.001
tMAX=100.0
R=0.4
g=-9.81
HELICES=True #whether springs should be visible or not
This cell defines several functions that are used in various parts of this notebook. The first 5 functions are used to manipulate Python list objects as if they were 3D vectors, allowing us to subtract, add, scalar multiply, find the magnitude of, take a dot product and also to find a unit vector parallel to the vectors that are inputted. The conVec function converts the Python lists in to VPython vectors so they can be used by the VPython library.
The fl and i functions are used to convert strings in to floats and ints respectively with the added bonus of handling empty strings. These are used to convert the text inputted via the GUI into particle and spring variables.
The allowedRange function is used to determine if the particles the user wishes to attach a spring to actually exist.
def sub(e1,e2):
return [e1[0]-e2[0],e1[1]-e2[1],e1[2]-e2[2]]
def add(e1,e2):
return [e1[0]+e2[0],e1[1]+e2[1],e1[2]+e2[2]]
def sMul(a,e1):
return ([a*e1[0],a*e1[1],a*e1[2]])
def dot(e1,e2):
return (e1[0]*e2[0]+e1[1]*e2[1]+e1[2]*e2[2])
def NORM(e1):
return sqrt(e1[0]**2+e1[1]**2+e1[2]**2)
def unit(e1):
try:
return sMul(1.0/NORM(e1),e1)
except ZeroDivisionError:
return [0.0,0.0,0.0]
def conVec(e1):
vec=vector(0.0,0.0,0.0)
vec.x=e1[0]
vec.y=e1[1]
vec.z=e1[2]
return vec
def fl(string,m=0,r=0):
try:
return float(string)
except ValueError:
if(m):
return 1.0
else:
if(r):
return R
else:
return 0.0
def INT(string):
try:
return int(string)
except ValueError:
return 0
def allowedRange(a):
if(a>=0 and a<len(points)):
return True
else:
return False
The class definition below defines the various properties of a spherical particle and also contains functions used to update the velocity and position of the particle due to the various forces acting on it. While computing these updates we assume that the force remains constant over a small time step defined in the first cell as DELTAt allowing us to use the constant acceleration equations
$$\\ \textbf{v}=\textbf{v}_{0}+\frac{\textbf{F}}{m} \delta t \\$$ $$\ \textbf{r}=\textbf{r}{0}+\textbf{v}{0} \delta t +\frac{\textbf{F}}{2m} \delta t ^{2} \$$
to update the velocity and position of the particle.
class point:
def __init__(self,r0,v0=[0.0,0.0,0.0],m0=1.0,fixed=0,RADIUS=R):
self.r=r0
self.v=v0
self.radius=RADIUS
if m0==0.0:
self.m=0.0001
else:
self.m=m0
self.F=[0.0,0.0,0.0]
self.FIXED=fixed
if(self.FIXED):
self.v=[0.0,0.0,0.0]
def vUpdate(self):
self.v=add(self.v,sMul(DELTAt/self.m,self.F))
def rUpdate(self):
self.r = add(self.r,add(sMul(DELTAt,self.v),sMul(0.5*(DELTAt**2)/self.m,self.F)))
Below is a class definition for the spring object. It has properties such as unstretched length, spring constant and the indexes (for the points array defined in the first cell) of the particles the spring is attached to. It also contains functions to compute the length of the spring, the magnitude of the force supplied by the spring and an update function which computes the updated lengths,force magnitudes and directions after each time step. The springs are assumed to be hookean and so obey Hookes law
$$\ \textbf{F}=k (L-L_{0}) \hat{\textbf{r}}.\$$
(Where $\hat{\textbf{r}}$ is the direction of the force, which is also computed by the update function).
class spring:
def length(self):
return NORM(sub(points[self.start].r,points[self.end].r))
def force(self):
return self.k*(self.L-self.L0)
def __init__(self,point1,point2,l0,K):
self.start=point1
self.end=point2
self.L0=l0
self.L=self.length()
self.k=K
self.F=self.force()
self.rs=unit(sub(points[self.end].r,points[self.start].r))
self.re=unit(sub(points[self.start].r,points[self.end].r))
def update(self):
self.L=self.length()
self.F=self.force()
self.rs=unit(sub(points[self.end].r,points[self.start].r))
self.re=unit(sub(points[self.start].r,points[self.end].r))
Below we have two functions which handle collisions between particles.
$$\ \$$
def impulse(i,j):
if(i==j):
pass
#else:
# n=unit(sub(points[i].r,points[j].r))
# component=sMul(2.0*dot(n,sub(points[i].v,points[j].v))/(points[i].m*((1.0/points[i].m)+(1.0/points[j].m))),n)
# points[i].v=add(points[i].v,component)
# points[i].r=add(points[i].r,sMul(5.0*DELTAt,points[i].v))
def collisions(i):
for j in range(0,len(points)):
if(NORM(sub(points[i].r,points[j].r))<=(points[i].radius+points[j].radius)):
impulse(i,j)
break
This cell contains a function which is responsible for computing the net force acting on each particle, including gravity if the toggle is activated and then calls the position and velocity update functions attached to each particle. It also calls the collision functions if the collisions toggle is activated to determine if any particles have collided and if they have, it then updates each particle's momentum accordingly.
def updatePoints():
for i in range(0,len(points)):
points[i].F=[0.0,0.0,0.0]
for j in range(0,len(springs)):
if(springs[j].start==i):
flag=1
elif(springs[j].end==i):
flag=-1
else:
flag=0
if(points[i].FIXED):
points[i].F=[0.0,0.0,0.0]
else:
if(flag==1):
points[i].F=add(points[i].F,sMul(springs[j].F,springs[j].rs))
elif(flag==-1):
points[i].F=add(points[i].F,sMul(springs[j].F,springs[j].re))
else:
pass
if (GRAVITY and (not points[i].FIXED)):
points[i].F=add(points[i].F,sMul(g*points[i].m,[0,1,0]))
if(COLLISIONS):
collisions(i)
points[i].vUpdate()
points[i].rUpdate()
for i in range(0,len(springs)):
springs[i].update()
The class below is responsible for defining the 3D scene containing the particles and springs. It also contains the main loop which is where time is incremented and the particles positions are computed and plotted after each time step.
class vis3D:
def __init__(self):
self.BALL=[]
self.scene = canvas()
for i in range(0,len(points)):
self.BALL.append(sphere(pos=conVec(points[i].r), radius=points[i].radius))
if(points[i].FIXED):
self.BALL[i].color=color.red
else:
self.BALL[i].color=color.green
if HELICES:
self.HELIX=[]
for i in range(0,len(springs)):
self.HELIX.append(helix(pos=conVec(points[springs[i].start].r),axis=conVec(sMul(springs[i].L,springs[i].rs)), radius=0.4, color=color.blue))
def update(self):
updatePoints()
for i in range(0,len(points)):
self.BALL[i].pos=conVec(points[i].r)
if HELICES:
for i in range(0,len(springs)):
self.HELIX[i].pos=conVec(points[springs[i].start].r)
self.HELIX[i].axis=conVec(sMul(springs[i].L,springs[i].rs))
def run(self):
t=0.0
self.scene.background=color.white
display(self.scene)
while t<tMAX:
rate(int(1.0/DELTAt))
self.update()
t+=DELTAt
The final class is an object that defines the user interface of this notebook, it also has functions that create springs and particles when certain buttons are pressed.
class gui:
def __init__(self):
self.pos = [widgets.Text(description='X',width=100),widgets.Text(description='Y',width=100 ),widgets.Text(description='Z',width=100)]
self.POS = widgets.HBox(children=self.pos)
self.vel=[widgets.Text(description='Vx',width=100),widgets.Text(description='Vy',width=100 ),widgets.Text(description='Vz',width=100)]
self.VEL=widgets.HBox(children=self.vel)
self.misc=[widgets.Text(description='Mass',width=100),widgets.Text(description='Radius',width=100),widgets.widget_bool.Checkbox(description='Fixed',width=100)]
self.MISC=widgets.HBox(children=self.misc)
self.create=widgets.Button(description="Create Point",width=100)
self.NEXT = widgets.Button(description="Next",width=100)
self.sprAtt = [widgets.Text(description='Start',width=100),widgets.Text(description='End',width=100 )]
self.SPRATT = widgets.HBox(children=self.sprAtt)
self.sprProp = [widgets.Text(description='L0',width=100),widgets.Text(description='K',width=100 )]
self.SPRPROP = widgets.HBox(children=self.sprProp)
self.createSpr=widgets.Button(description="Create Spring",width=100)
self.grav=[widgets.widget_bool.Checkbox(description='Gravity',width=100),widgets.widget_bool.Checkbox(description='Collisions',width=100)]
self.GRAV=widgets.HBox(children=self.grav)
self.START=widgets.Button(description="Start",width=100)
self.create.on_click(self.addPoint)
self.NEXT.on_click(self.nxt)
self.createSpr.on_click(self.addSpring)
self.START.on_click(self.start)
def display(self):
display(self.POS,self.VEL,self.MISC,self.create,self.NEXT)
def addPoint(self,b):
points.append(point([fl(self.pos[0].value),fl(self.pos[1].value),fl(self.pos[2].value)],[fl(self.vel[0].value),fl(self.vel[1].value),fl(self.vel[2].value)],fl(self.misc[0].value,m=1),INT(self.misc[2].value),fl(self.misc[1].value,r=1)))
if(points[len(points)-1].FIXED):
print("Fixed Particle " +str(int(len(points)))+" Created. r0="+str(points[len(points)-1].r)+"m. v0="+str(points[len(points)-1].v)+"m/s. mass="+str(points[len(points)-1].m)+"kg. radius="+str(points[len(points)-1].radius)+"m.")
else:
print("Movable Particle " +str(int(len(points)))+" Created. r0="+str(points[len(points)-1].r)+"m. v0="+str(points[len(points)-1].v)+"m/s. mass="+str(points[len(points)-1].m)+"kg. radius="+str(points[len(points)-1].radius)+"m.")
def nxt(self,b):
display(self.SPRATT,self.SPRPROP,self.GRAV,self.createSpr,self.START)
#make plot of point location numbered
def addSpring(self,b):
if(not(allowedRange(INT(self.sprAtt[0].value)-1) and allowedRange(INT(self.sprAtt[1].value)-1))):
print("Couldn't Create Spring")
else:
springs.append(spring(INT(self.sprAtt[0].value)-1,INT(self.sprAtt[1].value)-1,fl(self.sprProp[0].value,m=1),fl(self.sprProp[1].value,m=1)))
print("Spring Created Between Particles " +str(INT(self.sprAtt[0].value))+" and " + str(INT(self.sprAtt[1].value))+". L0="+str(springs[len(springs)-1].L0)+"m. K="+str(springs[len(springs)-1].k)+"Nm.")
def start(self,b):
if self.grav[0].value:
global GRAVITY
GRAVITY = True
if self.grav[1].value:
global COLLISIONS
COLLISIONS = False#True
self.visual=vis3D()
self.visual.run()
Select cell in the menu bar above and then run all to start the GUI where particles and their connections (springs) can be defined.
GUI=gui()
GUI.display()