Trifilar pendulum

by Roberto Lot

> restart: mylib :="C:/LIBS": libname := mylib,libname: with(MBSymba_r5):

Warning, Non-Linear Modeling option has been choosen

> PDEtools[declare]([x(t),y(t),z(t), ux(t),vx(t),wx(t), uy(t),vy(t),wy(t), uz(t),vz(t),wz(t)] , prime=t, quiet):

pendulum description

suspended plate

frame TN, fixed to the plate

> TN := < < ux(t),uy(t),uz(t), 0> |
< vx(t),vy(t),vz(t), 0> |
< wx(t),wy(t),wz(t), 0> |
< x(t), y(t), z(t), 1> >;

center of mass

> G := make_POINT(TN, l/4,0,0):# show(G);

plate inertia

> plate := make_BODY(G,m,Ix,Iy,Iz,Cyz,Cxz,Cxy): show(plate);

gravity force

> gravity_force := make_FORCE (ground, 0,0,-m*g, G, plate):
show(gravity_force);

plate vertex definition

plate vertex

> A := make_POINT(TN,l*sqrt(3)/2, 0, 0): show(A,'synthetic');
B := make_POINT(TN, 0, l/2, 0): show(B,'synthetic');
C := make_POINT(TN, 0, -l/2, 0): show(C,'synthetic');

check vertex disctance

> AB := make_VECTOR(A,B): BC := make_VECTOR(B,C): CA := make_VECTOR(C,A):
dot_prod(AB,AB), dot_prod(BC,BC), dot_prod(CA,CA);

ground vertex definition

> A0 := make_POINT(ground,l*sqrt(3)/2, 0, h): show(A0,'synthetic');
B0 := make_POINT(ground, 0, l/2, h): show(B0,'synthetic');
C0 := make_POINT(ground, 0, -l/2, h): show(C0,'synthetic');

Coordinates and Constraints

dependent coordinates

> q := [x(t),y(t),z(t), ux(t),vx(t),wx(t), uy(t),vy(t),wy(t), uz(t),vz(t),wz(t)];
nq := nops(q);

since 12 coordinates are used for describing a 3 dof system, 9 independent constraints equations must be found

rigid body constraints

these constraints are obtainde from the orthonormal property of the rotation submatrix

> evalm(inv_frame(TN) &* TN - 1);

> constraints := convert(%,set) minus {0}:
for i from 1 to nops(constraints) do i, constraints[i]; od;

ground-plate vertex distance constraints

> A0A := make_VECTOR(A0,A): c1 := dot_prod(A0A,A0A)-h^2;

> B0B := make_VECTOR(B0,B): c2 := dot_prod(B0B,B0B)-h^2;

> C0C := make_VECTOR(C0,C): c3 := dot_prod(C0C,C0C)-h^2;

update constraints

> constraints := convert( constraints union { c1,c2,c3} ,list):

Non-Linear Lagrange's equations

kinetic energy

> KE := kinetic_energy(plate);

generalized forces

> Qforce := map(simplify,generalized_force(gravity_force,q));

Lagrange's equations

> LL := KE-GE + sum('lambda||k(t)'*constraints['k'],'k'=1..nops(constraints)): lagrangian

> leqns := [seq(lagrange(LL,q[i],t) - Qforce[i] , i=1..nq)]:

> for i from 1 to nq do i,q[i]; leqns[i]; od;

The result is a set on n=12 Lagrange's equations (2nd order differential equations) plus f=9 constraint equations (algebraic equations).
The unknowns are n=12 coordinates plus f=9 Lagrange's multipliers

Baumgarte's method

The Differential-Algebraic set of equations is being transformed into a pure Differential set of equations

> baumgarte := convert( evalm(omega^2*constraints + 2*zeta*omega*diff(constraints,t) + diff(constraints,t,t)) ,list):
for i from 1 to nops(baumgarte) do i; baumgarte[i]; od;