Mechanical Systems.mws
Analysis and Simulation of Simple Dynamic Systems
by Forhad Ahmed
Columbus Alternative High School
2632 McGuffey Road
Columbus, Ohio 43211
U.S.A
email
: thesoccerstar2k3@hotmail.com
This worksheet demonstrates the use of Maple for studying dynamic systems.
Background and Program Code
Usage
Dynamic Systems
Dynamic Systems (as they are used in
this worksheet) are mechanical systems comprising of masses and constraints.
Constraints are a set of rules that have to be followed as the masses
are allowed to move freely. For example, a very simple dynamic system
might be that of a pendulum. The constraint of a pendulum is that one
mass has to stay at a fixed distance from a fixed mass (rigidrod constraint).
Examples of more complex system are a double pendulum and a triple pendulum.
A Simple Example:
Study of dynamic systems is done by analyzing the differential equations that
govern every mass in the system. For a very simple example, consider
only one point object of mass 2 [kg] located at the origin, and the only force
acting on it is a 2 [N] force 60 degrees above the xaxis.
This is indeed a very simple system and Newton's second law can be used easily
to find the differential equations of motion. The solutions to these
equations will be the functions that represent the motion of the mass object.
The solution to the preceding problem
can be found easily using analytic methods, but in order to do so, initial
conditions for the system must be defined. These are the values of functions
at initial times needed to successfully integrate the equations and find a
particular solution. The initial conditions needed for each mass in
a system are its initial x and y positions (x0, y0) and its initial velocities
(vx0 , vy0) In this case, the initial position is the origin (0,0)
and the initial velocity is also (vx0 , vy0) = (0,0) because the object is
not given an initial velocity.
When analyzing more complex systems however, the differential equations are
not so simple. The more masses, constraints, vibrations, pin joints
etc. are present in the system, the more complex the equations are. Solving
them by hand is nearly impossible, but computers have the ability to find
numerical answers given the proper instructions and the equations.
Program Design
This entire project is based on a program
code written in Maple 7. Maple enables solving of differential equations
(DE's) in either symbolic or numeric form. For most simple DE's, Maple
is able to find analytic solutions but for complex DE's, Maple uses its numerical
algorithms to find a numeric answer. My objective as a programmer is
to set up the differential equations that govern a dynamic system that the
user defines through Maple. Defining can be thought of as describing
to the computer a system of masses, constraints and initial conditions, and
the job of the computer, then, is to simulate and analyze the motion of the
system.
Program Diagram:
Program Usage
The Maple program is made up of modules
or procedures that perform specific tasks when called upon. The user
defines a dynamic system using the following commands in the proper format
indicated below. Also shown here are different commands that are used
to utilize this program.
1) (restart/initialize): Initialize
( ):
2) (define gravity and
viscosity): GravityViscosity
(g, ):
3) (define a mass object):
Mass (n,
x0, y0, vx0, vy0, m):
4) (define a Hooke's spring ): Spring
(A1, A2, k, L):
5) (add horizontal forceto a mass object):
AddForceX (n, Fx):
6) (add
vertical force: AddForce Y(n,
Fy):
to
a mass object)
7)
(anchor a mass object so: Fixed
(n):
that it remains fixed)
8) (add constraint that allows a :
SlideX (n):
mass object to
only slide
horizontally )
9) (add constraint that allows a :
SlideY (n):
vertically )
10) (cause a mass object to vibrate:
VibrateX (n,
A,
)
horizontally)
11) (cause a mass object to vibrate:
VibrateY (n,
A,
)
vertically)
12) (create differential equations
CreateEquation s(
)
for the specified
dynamic system. Issue this command
after all masses, constraints
etc. have been defined)
13) solve the DE's that have been: SolveEquations
( )
created and store the solution in
a list format)
14) (animate the simulation): RunSimulation
(tf, dt )
15) (visualize graphs of position: Visualize
(v1, v2, tf)
Initializations
For running this program, the following packages are required:
> 
restart:
with(plots):
with(plottools):
with(ListTools);
with(DEtools):

Warning, the name changecoords has been redefined
Warning, the name arrow has been redefined
Warning, the assigned name Group now has a global binding
Warning, the name translate has been redefined
Program Code
Execute the following code in order for this program to work correctly
> 
# 'Initialize all variables and set all flags to 0' module #;
Initialize:=proc()
global i,F,ic,m,noo,nos,nof,nslx,nsly,nvibx,nviby,picture;
for i from 1 to 250 do:
F[i]:=0:
ic[i]:=0:
m[i]:=0:
noo:=0:
nos:=0:
nof:=0:
nslx:=0:
nsly:=0:
nvibx:=0:
nviby:=0:
picture:=0:
end do:
end proc:

> 
# module for defining gravity and air resistance #
GravityViscosity:=proc(GG,muu)

> 
# module for defining mass objects and their initial conditions #
Mass:=proc(n,x0,y0,vx0,vy0,m0)
global F,ic,m,noo:

> 
F[n]:=array([mu*diff(x[n](t),t),m0*gmu*diff(y[n](t),t)]):

> 
ic[n]:={x[n](0)=x0,y[n](0)=y0,D(x[n])(0)=vx0,D(y[n])(0)=vy0}:
m[n]:=m0;
noo:=noo+1

> 
# module for defining springs and their properties #
Spring:=proc(n1,n2,k,L)

> 
global d,F,ic,temp1,temp2,temp3,temp4,springarray,nos:
d:=L;
nos:=nos+1;
springarray[nos]:=array([n1,n2,k,L]);

> 
temp1:=(F[n1][1]):
temp1:=temp1+(x[n1](t)x[n2](t))/sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)*(k)*(sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)d)mu*diff(x[n1](t),t):
(F[n1][1]):=temp1:

> 
temp2:=(F[n1][2]):
temp2:=temp2+(y[n1](t)y[n2](t))/sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)*(k)*(sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)d)mu*diff(y[n1](t),t):
(F[n1][2]):=temp2:

> 
temp3:=(F[n2][1]):
temp3:=temp3+(x[n2](t)x[n1](t))/sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)*(k)*(sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)d)mu*diff(x[n2](t),t):
(F[n2][1]):=temp3:

> 
temp4:=(F[n2][2]):
temp4:=temp4+(y[n2](t)y[n1](t))/sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)*(k)*(sqrt((x[n1](t)x[n2](t))^2+(y[n1](t)y[n2](t))^2)d)mu*diff(y[n2](t),t):
(F[n2][2]):=temp4

> 
# module for adding a horizontal force to a mass object #
AddForceX:=proc(n,fx)
global F;
F[n][1]:=F[n][1]+fx;
end proc:

> 
# module for adding a vertical force to a mass object #
AddForceY:=proc(n,fy)
global F:
F[n][2]:=F[n][2]+fy;
end proc:

> 
# module for anchoring mass objects (fixing them in place) #
Fixed:=proc(n0)
global F,pin,nof;
nof:=nof+1;
pin[nof]:=n0:
F[n0][1]:=0:
F[n0][2]:=0:
end proc:

> 
# module for allowing mass objects to slide horizontally only #
SlideX:=proc(n0)
global F,nslx,sliderx;
nslx:=nslx+1;
sliderx[nslx]:=n0;
F[n0][2]:=0;
end proc:

> 
# module for allowing mass objects to slide vertically only #
SlideY:=proc(n0)
global F,nsly,slidery;
nsly:=nsly+1;
slidery[nsly]:=n0;
F[n0][1]:=0;
end proc:

> 
# module for causing a mass object to vibrate horizontally #
VibrateX:=proc(n,amp,om)
global F,nvibx,arrx;
nvibx:=nvibx+1;
arrx[nvibx]:=n;
F[n][1]:=m[n]*amp*om^2*(cos(om*t))/2;
F[n][2]:=0;
end proc:
# module for causing a mass object to vibrate vertically #
VibrateY:=proc(n,amp,om)
global F,nviby,arry;
nviby:=nviby+1;
arry[nviby]:=n;
F[n][2]:=m[n]*amp*om^2*cos(om*t)/2;
F[n][1]:=0;
end proc:

> 
# module for creating the system differential equations #
CreateEquations:=proc()
global eq,i,h,k,H,u;
for i from 1 to noo do;
eq[i]:={m[i]*diff(x[i](t),t,t)=F[i][1],m[i]*diff(y[i](t),t,t)=F[i][2]}:
end do;
H:=[seq([eq[i][1],eq[i][2],ic[i][1],ic[i][2],ic[i][3],ic[i][4]],i=1..noo)]:
u:= Flatten(H);
end proc:

> 
# module for solving the differential equations for the system #
SolveEquations:=proc()
global Sys,eq,sol,i,X,Y,u;
Flatten(H);
u:=Flatten(H):
Sys:=u:
sol:=dsolve(Sys,numeric,output=listprocedure):
for i from 1 to noo do;
X[i]:=eval(x[i](t),sol);
Y[i]:=eval(y[i](t),sol);
end do:
end proc:

> 
# module for animating the dynamic system #
RunSimulation:=proc(tf,delt)
global massobject,springobject1,h,k,picture,sliderobjectx,sliderobjecty,vibobjectx,vibobjecty,pinobject:
local x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,y5,y6,y7,y8,y9,cc,s,i,j,middlenumber,size,a1,a2,dis,z,jj:
cc:=1:
for i from 0 by delt to tf do:

> 
for j from 1 to noo do:

> 
middlenumber:=trunc((4*noo+1)/2)+1:

> 
massobject[j]:=disk([h[j],k[j]],.07,color=blue):

> 
end do:
for z from 1 to nos do:
a1:=springarray[z][1];
a2:=springarray[z][2];
size:=springarray[z][4];
x1:=h[a1];
y1:=k[a1];
x2:=h[a2];
y2:=k[a2];
dis:=sqrt((x2x1)^2(y2y1)^2);
x3:= (x1 + x2)/2: y3 := (y1 + y2)/2:
x4:= (x1 + x3)/2: y4 := (y1 + y3)/2:
x5:= (x3 + x2)/2: y5 := (y3 + y2)/2:
x6:= (x4 + x1)/2: y6 := (y4 + y1)/2:
x7:= (x3 + x4)/2: y7 := (y4 + y3)/2:
x8:= (x3 + x5)/2: y8 := (y3 + y5)/2:
x9:= (x5 + x2)/2: y9 := (y5 + y2)/2:
springobject1[z]:={point([x3,y3],color=red),point([x4,y4],color=red),point([x5,y5],color=red),point([x6,y6],color=red),point([x7,y7],color=red),point([x8,y8],color=red),point([x9,y9],color=red)};
end do:
for j from 1 to nslx do;
jj:=sliderx[j]; sliderobjectx[j]:=line([h[jj]1/2,k[jj]],[h[jj]+1/2,k[jj]],linestyle=7)
end do:
for j from 1 to nsly do:
jj:=slidery[j];
sliderobjecty[j]:=line([h[jj],k[jj]+1/2],[h[jj],k[jj]1/2],linestyle=7)
end do:
for j from 1 to nvibx do:
jj:=arrx[j];
vibobjectx[j]:=
{line([h[jj]+.15,k[jj]],[h[jj].15,k[jj]],color=red),
line([h[jj]+.15,k[jj]],[h[jj]+.10,k[jj]+.05],color=red),
line([h[jj]+.15,k[jj]],[h[jj]+.10,k[jj].05],color=red),
line([h[jj].15,k[jj]],[h[jj].10,k[jj]+.05],color=red),
line([h[jj].15,k[jj]],[h[jj].10,k[jj].05],color=red)}
end do:
for j from 1 to nviby do:
jj:=arry[j];
vibobjecty[j]:=
{line([h[jj],k[jj]+.15],[h[jj],k[jj].15],color=red),
line([h[jj],k[jj]+.15],[h[jj]+.05,k[jj]+.10],color=red),
line([h[jj],k[jj]+.15],[h[jj].05,k[jj]+.10],color=red),
line([h[jj],k[jj].15],[h[jj]+.05,k[jj].10],color=red),
line([h[jj],k[jj].15],[h[jj].05,k[jj].10],color=red)}
end do:
for j from 1 to nof do;
jj:=pin[j];
pinobject[j]:={line([h[jj].1,k[jj]],[h[jj]+.1,k[jj]],thickness=2,color=red),line([h[jj],k[jj]+.1],[h[jj],k[jj].1],thickness=2,color=red)};
end do;

> 
picture[cc]:=display(seq(massobject[q],q=1..noo),seq(springobject1[w],w=1..nos),seq(pinobject[w],w=1..nof),seq(sliderobjectx[w],w=1..nslx),seq(sliderobjecty[w],w=1..nsly),seq(vibobjectx[w],w=1..nvibx),seq(vibobjecty[w],w=1..nviby),scaling=constrained):

> 
cc:=cc+1;
print("Solving system of equations at time (t)"=i);

> 
end do:
display(seq(picture[q],q=1..cc1),scaling=constrained,axes=boxed,insequence=true):

> 
# module for visualizing a mass object's position vs. time and more #
Visualize:=proc(a,b,tf)
local temp1,temp2,temp3;
temp3:=Flatten([seq([x[i](t),y[i](t)],i=1..noo)]);
temp1:=Flatten([seq([eq[i][1],eq[i][2]],i=1..noo)]);
temp2:=[Flatten([seq([ic[i][1],ic[i][2],ic[i][3],ic[i][4]],i=1..noo)])];
DEplot(temp1,temp3,t=0..tf,temp2,scene=[a,b],stepsize=.005,thickness=1,linecolour=red);
display(%,axesfont=[HELVETICA,7],labelfont=[HELVETICA,7],scaling=constrained);
end proc:

Case Study 1 (Damped Oscillation)
This program is used for studying simple mechanical systems, especially systems that involve vibrations. Consider a simple system consisting of two mass objects of mass 1 [kg]; one fixed and one free (with massnumbers 1 and 2 respectively) , connected by a spring of stiffness 50 [N/m] and length one. The free mass is given an initial velocity downwards of 0.5 [m/s]. Viscosity (airresistance) is present and the coefficient of dampening is
= 0.3. Gravity is present and the gravitational constant is the usual 9.8 [m/
]. How will this system behave?
In the above picture of the model, the masses are blue circles and the spring is the set of dots connecting the masses. Here is the Maple input for defining this system (the masses are located at (0,0) and (0,1); the red cross indicates that the mass is anchored or fixed):
> 
GravityViscosity(9.8,0.3):

> 
Mass(1,0,0,0,0,1):
Mass(2,0,1,0,.5,1):
Spring(1,2,50,1):

After entering these commands, Maple has created the differential equations that govern the motion of the system. Issuing the
RunSimulation
(finaltime, increment) command animates the system from t = 0 to t = finaltime with the given increment. Let's run this simulation upto time = 6 seconds with an increment of .1
This program also allows you to plot the position of a mass object versus time. Suppose you want to graph the yposition of the bottom mass (mass # 2) versus time upto time = 6 seconds. This following command is used:
> 
Visualize(t,y[2](t),6);

As predicted, the airresistance causes the vibration of mass2 to decrease with time.
Case Study 2 (Resonance)
Consider a similar model to the previous one; instead of the top mass being fixed, it is forced to vibrate vertically with an amplitude of .2 meters and a frequency of 6. Also, the length of the spring is increased to 3 meters. There is no air resistance. Here is a visual representation of the model. The red arrows indicate that the mass is vibrating:
The top mass (vibrating) is mass 1 and the bottom is mass 2. Their coordinates are (0,0) and (0,3) respectively. Here is the Maple input for defining this model:
> 
GravityViscosity(9.8,0):

> 
Mass(1,0,0,0,0,1):
Mass(2,0,3,0,0,1):
Spring(1,2,50,3):

Animate the model until time = 4 seconds with .1 second increments:
Graph the yposition of the bottom mass (mass 2) versus time upto time = 20 seconds:
> 
Visualize(t,y[2](t),20);

It can be seen that the amplitude of the vibration of the bottom mass increases, then drecreases and increases etc. within the 20 second period. However, if the vibration of the top mass matched the natural freqency of the system, the bottom mass will continue vibrate with increasing amplitude (resonance). The natural frequency of a springmass system is
so in this case that frequency would be
or about 7.07. Let's try this model again only changing the vibrational frequency of the top mass to match the natural frequency so resonance can occur.
Maple Input:
> 
GravityViscosity(9.8,0):

> 
Mass(1,0,0,0,0,1):
Mass(2,0,3,0,0,1):
Spring(1,2,50,3):

Graph the yposition of mass 2 versus time upto time = 8 seconds:
> 
Visualize(t,y[2](t),8);

The bottom mass object's vibrational amplitude increses with time: it is resonating.
Case Study 3 (Beadonwire)
Consider a mass object with a mass of 0.1 [kg] located at the origin and it is only allowed to slide horizontally as if it were a bead on a frictionless wire. It is connected to a fixed mass (located at (1,0)) and a free mass (located at (1, 1)) by springs of stiffness 1 and 200 [N/m] respectively. One of the two connections (with the stiffer spring) forms a pendulum in which the free mass is allowed to swing freely. Our goal is to model this scenario. Here is a picture of the model:
The masses are labeled for clarity. The horizontal dotted line through mass1 indicates that it is only allowed to slide horizontally. Mass1 is connected to mass2 and mass3. The first connection (between 1 and 2) is made with a spring of stiffness 5 [N/m] and the second connection (between 1 and 3) is made with a spring of stiffnes 200 [N/m] (stiff springs are used as approximations of rigid rods for modeling a pendulm). Ther gravitational constant is 9.8 and there is no air resistance (
= 0).
Maple Input:
> 
GravityViscosity(9.8,0):

> 
Mass(1,0,0,0,0,.1):
Mass(2,1,0,0,0,.1):
Mass(3,1,1,0,0,.1):
Spring(1,2,5,1):
Spring(1,3,200,1.414):

Let's run this simulation until time = 3 seconds with .1 second increments:
Here is a graph of the horizontal vibration of mass1 versus time upto time = 8 seconds:
> 
Visualize(t,x[1](t),8);

Here is a graph of the horizontal vibration of mass3 (swinging mass) versus its vertical vibration upto time = 3 seconds:
> 
Visualize(x[3](t),y[3](t),3);

Case Study 4 (Vibrating Tower)
Consider a system of 8 masses connected by springs. They are aligned vertically (resembling a tower) and the lowest of the masses are forced to vibrate horizontally with a given amplitude and frequency:
All mass objects (in blue) have a massnumber (from 1 to 8) and a mass of 0.1 [kg] and the springs have stiffness 400. The bottom two masses are 1 & 2; the top 2 masses are 7 & 8. The bottom two masses (1 & 2) are forced to vibrate at an amplitude of .2 meters and a frequency of 10. The goal is to see how much the tower vibrates and how the vibration can be minimized. This program allows you to model this dynamic system and analyze its motion. Here is the Maple input:
> 
GravityViscosity(10,0):

> 
Mass(1,0,0,0,0,.1):
Mass(2,1,0,0,0,.1):
Mass(3,0,2,0,0,.1):
Mass(4,1,2,0,0,.1):
Mass(5,0,4,0,0,.1):
Mass(6,1,4,0,0,.1):
Mass(7,0,6,0,0,.1):
Mass(8,1,6,0,0,.1):
Spring(1,3,400,2):
Spring(2,4,400,2):
Spring(1,4,400,sqrt(5)):
Spring(2,3,400,sqrt(5)):
Spring(4,3,400,1):
Spring(4,6,400,2):
Spring(3,5,400,2):
Spring(3,6,400,sqrt(5)):
Spring(4,5,400,sqrt(5)):
Spring(5,6,400,1):
Spring(5,7,400,2):
Spring(7,8,400,1):
Spring(8,6,400,2):
Spring(5,8,400,sqrt(5)):
Spring(6,7,400,sqrt(5)):

> 
VibrateX(1,.2,10):
VibrateX(2,.2,10):

Let's run this simulation upto time = 3 seconds with increments of .1 and see the motion (this may take a LONG time):
This program also allows you to plot the position of a mass object versus time. Suppose you want to graph the xposition of the top right mass (mass # 8) versus time upto time = 2 seconds. This following command is used:
> 
Visualize(t,x[8](t),2);

The preceding is a graph of the xposition of the top right mass object versus time upto 2 seconds. It can be seen that witin that time period, the maximum displacement of the mass object was about .4 meters. Now you can try different values for stiffness of the springs and try to see which ones cause the tower to vibrate less. There is a value of stiffness of the springs that, if given the right frequency of vibration of the two bottom masses, will cause the tower to resonate and become unstable! Its up to you to experiment.
Practice Problem
The preceding examples are only a small part of what this program can do. This worksheet may be considered as a sum of many worksheets that are present in this website because it allows for the study of not just one case but of infinately many cases. You can define a system as simple or as complex as you like and this program will model it for you! The beauty lies within its versatility.
Practice Problem
: Use this worksheet to solve the following problem:
Consider the above model of a bridge. Number the masses (in blue) in any way you like and connect them as they are connected in the picture by springs of different stiffnesses. All mass objects have a mass of 0.1 [kg]. The two masses that have gray arrowheads pointing to them are fixed masses. Define all the masses and all the springs. Set the gravitational constant to 9.8 and the coefficient of airresistance to 0.1. Try different values for the stiffness of the springs and notice which ones lead to a stable bridge and which ones cause the bridge to collapse. (For example, weak springs in the middle of the bridge will probably cause the structure to collapse). Have patience using this program because animating the preceding model will take time.
I am a high school student and would be honored to hear criticisms of this worksheet from professors, mathematicians and scientists. Also, if you have any questions or comments about this program, please feel free to email me. I will be glad to know and help you.
Also, I know Maple has brought about Maplets as a tool for integration of the Maple platform with Java applets. I do not own a copy of Maplets but I know enough about it to say that it would be a great tool to combine with the program in this worksheet. The user can define masses, springs and other constraints through the Maplets Graphical Interface as input and the program can display the animation of graphs of motion as output. If anyone is interested in making such a Maplet and wants to use this program as the basis of their project, feel free to contact me. I think it would be a great application for studying mechanics!
Disclaimer:
While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.