Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] 3 link pendulum #10

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

[WIP] 3 link pendulum #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
06398c4
Added a 3 link pendulum problem, with massless links. and
tarzzz Jun 18, 2013
11a0059
Added README
tarzzz Jun 18, 2013
d5f2d44
Update three_link_pendulum.py
tarzzz Jun 18, 2013
a61e602
Modified script for two dimensional joints
tarzzz Jun 19, 2013
e69d5cf
Merge branch '3-link-pendulum' of https://github.com/tarzzz/pydy_exam…
tarzzz Jun 19, 2013
8d379c2
modified Inertia and Velocity symbols. Changed code to acc to pep8 st…
tarzzz Jun 23, 2013
b7571e7
Updated according to PEP8 style guide
tarzzz Jun 23, 2013
944ab76
modified README, added a figure
Jul 11, 2013
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions three_link_pendulum/README
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
This example is the symbolic equations of motions derivation of a three
link conical pendulum.
The links of the pendulum are also considered to have mass, and are represented
as RigidBodies. The the end points of the links are represented by Particles.
The pendulum has six degrees of freedom.

To run this program type in terminal:

python three_link_pendulum.py

Binary file added three_link_pendulum/three_link_pendulum.gif
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
161 changes: 161 additions & 0 deletions three_link_pendulum/three_link_pendulum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
#For this one we assume RigidBody links
from sympy import symbols, sympify
from sympy.physics.mechanics import *

#Number of links = 3
N_links = 3

#Number of masses = 3
N_bobs = 3

#Defining Dynamic Symbols ................

#Generalized coordinates(angular) ...

alpha = dynamicsymbols('alpha1 alpha2 alpha3')
beta = dynamicsymbols('beta1 beta2 beta3')

#Generalized speeds(angular) ...
vel_alpha = dynamicsymbols('vel_alpha1 vel_alpha2 vel_alpha3')
vel_beta = dynamicsymbols('vel_beta1 vel_beta2 vel_beta3')

#Mass of each bob:
m = symbols('m:' + str(N_bobs))


#Length and mass of each link ..
l = symbols('l:' + str(N_links))
M = symbols('M:' + str(N_links))

#For storing Inertia for each link :
Ixx = symbols('Ixx:' + str(N_links))
Iyy = symbols('Iyy:' + str(N_links))
Izz = symbols('Izz:' + str(N_links))
Ixy = symbols('Ixy:' + str(N_links))
Iyz = symbols('Iyz:' + str(N_links))
Ixz = symbols('Ixz:' + str(N_links))

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A general inertia matrix has 9 entries: three moments of inertia and 6 products of inertia. You are only introducing the moments of inertia about the x and y axes. Where are the other entries? Is the link assumed symmetric? Depending on the symmetry you can reduce the inertia down to the three moments of inertia but since you have no Z inertia at all this seems problematic. I'd prefer you to add all of the inertia tensor entries here and then if you assume symmetry about various planes you can set the appropriate products of inertia equal to zero when you get to the numerical simulation.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only 6 entries will be unique, inertia matrices are always symmetric (i.e., Ixy == Iyx, Ixz == Izx, Iyz == Izy). If you create a 3x3 with 9 unique entries and somehow you accidentally make it non-symmetric, it won't be representative of an inertia matrix and is probably a mistake. The inertia() function in mechanics may be useful here.

#gravity and time ....
g, t = symbols('g t')


#Now defining an Inertial ReferenceFrame First ....

I = ReferenceFrame('I')

#And some other frames ...

A = I.orientnew('A', 'Body', [alpha[0], beta[0], 0], 'ZXY')

B = A.orientnew('B', 'Body', [alpha[1], beta[1], 0], 'ZXY')

C = B.orientnew('C', 'Body', [alpha[2], beta[2], 0], 'ZXY')

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The .orientnew method is a little more compact here.

A = ReferenceFrame('A')
B = A.orientnew('B', 'Body', (alpha[0], beta[0], 0), 'ZXY')

Also please read through the PEP8 style guide if you haven't: http://www.python.org/dev/peps/pep-0008/

You should have spaces after your commas.


#Setting angular velocities of new frames ...
A.set_ang_vel(I, vel_alpha[0] * I.z + vel_beta[0] * I.x)
B.set_ang_vel(A, vel_alpha[1] * A.z + vel_beta[1] * A.x)
C.set_ang_vel(B, vel_alpha[2] * B.z + vel_beta[2] * B.x)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you had defined generalized speed variables you could insert them here instead of the derivatives.

u0 = DynamicsSymbols('u0')
u1 = DynamicsSymbols('u1')
A.set_ang_vel(I, u0 * I.z + u1 * I.x)

Also did you set the angular velocity of I equal to zero. I think you need to, but not sure.



# An Origin point, with velocity = 0
O = Point('O')
O.set_vel(I, 0)

#Three more points, for masses ..
P1 = O.locatenew('P1', l[0] * A.y)
P2 = O.locatenew('P2', l[1] * B.y)
P3 = O.locatenew('P3', l[2] * C.y)

#Setting velocities of points with v2pt theory ...
P1.v2pt_theory(O, I, A)
P2.v2pt_theory(P1, I, B)
P3.v2pt_theory(P2, I, C)
points = [P1, P2, P3]

Pa1 = Particle('Pa1', points[0], m[0])
Pa2 = Particle('Pa2', points[1], m[1])
Pa3 = Particle('Pa3', points[2], m[2])
particles = [Pa1, Pa2, Pa3]



#defining points for links(RigidBodies)
#Assuming CoM as l/2 ...
P_link1 = O.locatenew('P_link1', l[0]/2 * A.y)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pep8 would be

0.locatenew('P_link1', l[0] / 2 * A.y)

P_link2 = O.locatenew('P_link1', l[1]/2 * B.y)
P_link3 = O.locatenew('P_link1', l[2]/2 * C.y)
#setting velocities of these points with v2pt theory ...
P_link1.v2pt_theory(O, I, A)
P_link2.v2pt_theory(P_link1, I, B)
P_link3.v2pt_theory(P_link2, I, C)

points_rigid_body = [P_link1, P_link2, P_link3]


#defining inertia tensors for links

inertia_link1 = inertia(A, Ixx[0], Iyy[0], Izz[0], ixy = Ixy[0], iyz = Iyz[0], izx = Ixz[0])
inertia_link2 = inertia(B, Ixx[1], Iyy[1], Izz[1], ixy = Ixy[1], iyz = Iyz[1], izx = Ixz[1])
inertia_link3 = inertia(C, Ixx[2], Iyy[2], Izz[2], ixy = Ixy[2], iyz = Iyz[2], izx = Ixz[2])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not write a loop to do these three lines?


#Defining links as Rigid bodies ...

link1 = RigidBody('link1', P_link1, A, M[0], (inertia_link1, O))
link2 = RigidBody('link2', P_link2, B, M[1], (inertia_link2, P1))
link3 = RigidBody('link3', P_link3, C, M[2], (inertia_link3, P2))
links = [link1, link2, link3]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't a loop be nicer here? If you make the right loops you will be able to specify N at the top of the script and everything will run.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do feel that for an illustrative example, loops could slow down the comprehension of a budding user.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fine too. But the example is half loops and half not, currently. It should be either or. The nice part about having loops is that you can pick any N at the beginning of the problem and it should derive your system.



#Applying forces on all particles , and adding all forces in a list..
forces = []
for particle in particles:

mass = particle.get_mass()
point = particle.get_point()
forces.append((point, -mass * g * I.y) )
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make a comment somwhere describing which direction you are defining to be "down" (i.e. aligned with gravity).


#Applying forces on rigidbodies ..
for link in links:
mass = link.get_mass()
point = link.get_masscenter()
forces.append((point, -mass * g * I.y) )
kinetic_differentials = []
for i in range(0,N_bobs):
kinetic_differentials.append(vel_alpha[i] - alpha[i].diff(t))
kinetic_differentials.append(vel_beta[i] - beta[i].diff(t))

#Adding particles and links in the same system ...
total_system = []
for particle in particles:
total_system.append(particle)

for link in links:
total_system.append(link)

q = []
for angle in alpha:
q.append(angle)
for angle in beta:
q.append(angle)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for angle in alpha + beta:
q.append(angle)

print q
u = []

for vel in vel_alpha:
u.append(vel)
for vel in vel_beta:
u.append(vel)

print u
kane = KanesMethod(I, q_ind = q, u_ind = u, kd_eqs = kinetic_differentials)
fr, frstar = kane.kanes_equations(forces, total_system)

print fr