Bringing constraint forces into evidence
Iceman_5
I'm interested in the forces exerted by the joints between the solids.
Luke
points/frames to have velocities/angular velocites which they cannot
actually possess, without introducing any generalized extra
generalized coordinates. Then you need to include the applied forces
at the point of interest (or torque for the reference frame of
interest). These extra generalized speeds result in extra partial
velocities and partial angular velocities, which lead to extra motion
equations which can be solved for the constraint forces.
There is a section in chapter 4 (I think) of Kane's 1985 Dynamics book
that describes this in detail. I would be happy to add this to the
pendulum or the rolling disc or torus example. I will work on this
over the weekend and update that example. Or if you have time and
figure it before then, email me a patch and I can review it and make
it part of the example.
~Luke
Iceman_5
commented in my post with the trigonometric expressions and I don't
know how to solve it. I post here my code, maybe that could help. It
could be a good example if all the problems are solved finally. It's a
little bit complicated due to the number of reactions, but I think the
code is well organized and someone can understand it.
from sympy import *
from pydy import *
# Create a Newtonian reference frame
N = NewtonianReferenceFrame('N')
#declare parameters
Inertia_parametersn = 'Ixc Iyc Izc Jxyc Jxzc Jyzc Ixd Iyd Izd Jyzd Ixe
Iye Ize '
Mass_parametersn = 'mc md me '
Geom_parametersn = 'R l4 l6 lgd betta l1 h1 '
auxiliar_paramsn = 's_b c_b '
Forces_extn = 'Fe1 Fe2 '
Forces_lign = 'FSC_x FSC_y FSC_z FCD1_r FCD1_y FCD1_t FCD2_r FCD2_y
FCD2_t FSD1_r FSD2_r '
Torques_lign = 'MSC_x MSC_z MCD1_x MCD2_x MSD1_x MSD2_x'
paramsn = Inertia_parametersn + Mass_parametersn + Geom_parametersn +
auxiliar_paramsn + Forces_extn + Forces_lign \
+ Torques_lign
params = N.declare_parameters(paramsn)
(Ixc, Iyc, Izc, Jxyc, Jxzc, Jyzc, Ixd, Iyd, Izd, Jyzd, Ixe, Iye, Ize)
= params[0:13]
(mc, md, me) = params[13:16]
(R, l4, l6, lgd, betta, l1, h1) = params[16:23]
(s_b, c_b) = params[23:25]
(Fe1, Fe2) = params[25:27]
(FSC_x, FSC_y, FSC_z, FCD1_r, FCD1_y, FCD1_t, FCD2_r, FCD2_y, FCD2_t,
FSD1_r, FSD2_r) = params[27:38]
(MSC_x, MSC_z, MCD1_x, MCD2_x, MSD1_x, MSD2_x) = params[38:44]
#declare generalized speeds and coordinates
q, qd = N.declare_coords('q', 1)
u, ud = N.declare_speeds('u', 14)
#----------------------------------------------------------------------
# Orient the three rigid bodies, assign inertias:
# We define the solid crankshaft
C = N.rotate("C", 2, q[0], I=(Ixc, Iyc, Izc, Jxyc, Jyzc, Jxzc))
C.G = N.O.locate('CG', 0, C, mc, force = FSC_x*N[1] + FSC_y*N[2] +
FSC_z*N[3])
# We define an auxiliary reference frame in the crankshaft
Caux = N.rotate("Caux", 2, q[0]+betta)
# We define the solid disc 1
D1 = N.rotate('D1', 2, -q[0], I=(Ixd, Iyd, Izd, 0, Jyzd, 0))
D1.G= N.O.locate('D1G', (-l6/2-l4+lgd)*C[2] + (R/2)*C[3], D1, md)
#mass center of the disc 1
# We define the solid disc 2
D2aux=N.rotate('D2aux', 3, pi)
D2 = D2aux.rotate('D2', 2, q[0]+betta, I=(Ixd, Iyd, Izd, 0, Jyzd, 0))
D2.G= N.O.locate('D2G', (l6/2+l4-lgd)*Caux[2] + (R/2)*Caux[3], D2, md)
#-----------------------------------------------------------------------
#locate the important points and apply forces
D1.P1 = D1.G.locate('P1', -lgd*D1[2] + (R/2)*D1[3], D1, force = -Fe1*N
[3])
D2.P2 = D2.G.locate('P2', -lgd*D2[2] + (R/2)*D2[3], D2, force = -Fe2*N
[3])
D1.S = D1.G.locate('CSD1', -lgd*D1[2] + (R/2)*C[3], D1, force = -
FSD1_r*C[3])
D2.S = D2.G.locate('CSD2', -lgd*D2[2] + (R/2)*Caux[3], D2, force = -
FSD2_r*Caux[3])
#these are the contact points, so each time they are fixed to the
discs but they aren't always the same point of the disc
D1.C = D1.G.locate('RD1C', (l4-lgd)*D1[2], D1, force = -FCD1_r*C[3] -
FCD1_t*C[1] - FCD1_y*C[2])
D2.C = D2.G.locate('RD2C', (l4-lgd)*D2[2], D2, force = -FCD2_r*Caux[3]
- FCD2_t*Caux[1] - FCD2_y*Caux[2])
C.D1 = C.G.locate('RCD1', -l6/2*C[2] + R/2*C[3], C, force = FCD1_r*C
[3] + FCD1_t*C[1] + FCD1_y*C[2])
C.D2 = C.G.locate('RCD2', l6/2*Caux[2] + R/2*Caux[3], C, force =
FCD2_r*Caux[3] + FCD2_t*Caux[1] + FCD2_y*Caux[2])
#apply torques
C.apply_torque(MSC_x*N[1] + MSC_z*N[3])
C.apply_torque(MCD1_x*C[1], other = D1)
C.apply_torque(MCD2_x*Caux[1], other = D2)
D1.apply_torque(MSD1_x*C[1])
D1.apply_torque(MSD2_x*Caux[1])
#---------------------------------------------------------------------------
## Set velocities and angular velocities
C.abs_ang_vel = Vector(u[0]*C[2] + u[4]*N[1] + u[5]*N[3])
D1.abs_ang_vel = Vector(-u[0]*D1[2] + u[8]*C[1] + u[9]*C[3])
D2.abs_ang_vel = Vector(u[0]*D2[2] + u[12]*Caux[1] + u[13]*Caux[3])
C.G.abs_vel = Vector(u[1]*N[1] + u[2]*N[2] + u[3]*N[3])
D1.G.abs_vel = Vector(R/2*u[0]*C[1] + u[6]*C[3] + u[7]*C[2])
D2.G.abs_vel = Vector(R/2*u[0]*Caux[1]+ u[10]*Caux[3] - u[11]*Caux[2])
D1.P1.abs_vel = D1.G.abs_vel + dt(D1.P1.rel(N.O), N).subs({qd[0]: u
[0]})
D2.P2.abs_vel = D2.G.abs_vel + dt(D2.P2.rel(N.O), N).subs({qd[0]: u
[0]})
D1.S.abs_vel = D1.G.abs_vel + dt(D1.S.rel(N.O), N).subs({qd[0]: u[0]})
D2.S.abs_vel = D2.G.abs_vel + dt(D2.S.rel(N.O), N).subs({qd[0]: u[0]})
### Define generalized speeds
u_defs = N.define_speeds([Eq(u[0], dot(C.ang_vel(), C[2]))])
print u_defs
## Set accelerations and angular accelerations
C.abs_ang_acc = dt(C.abs_ang_vel, N)
D1.abs_ang_acc = dt(D1.abs_ang_vel, N)
D2.abs_ang_acc = dt(D2.abs_ang_vel, N)
C.G.abs_acc = dt(C.G.abs_vel, N)
D1.G.abs_acc = dt(D1.G.abs_vel, N)
D2.G.abs_acc = dt(D2.G.abs_vel, N)
D1.P1.abs_acc = dt(D1.P1.abs_vel, N)
D2.P2.abs_acc = dt(D2.P2.abs_vel, N)
D1.S.abs_acc = dt(D1.S.abs_vel, N)
D2.S.abs_acc = dt(D2.S.abs_vel, N)
kanes_eqns = N.form_kanes_equations()
If you try to execute the code you will get the error with the
trigonometric expressions that can't be simplified.