Trifilar Pendulum in Natural Coordinates

frame TN, fixed to the plate

> TN := matrix(4,4,[[ux,vx,wx(t),xA(t)], [uy,vy,wy(t),yA(t)],[uz,vz,wz(t),zA(t)],[0,0,0,1]]);

plate vertex definition

plate vertex

> A := origin(TN):
B := make_POINT(TN, l*sqrt(3)/2, -l/2, 0):
C := make_POINT(TN, l*sqrt(3)/2, +l/2, 0):
show(A,'synthetic'), show(B,'synthetic'), 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);

equivalence of definitions

> BB := make_VECTOR( make_POINT(ground,xB(t),yB(t),zB(t)), B):

> CC := make_VECTOR( make_POINT(ground,xC(t),yC(t),zC(t)), C):
show(BB); show(CC);

> solz := solve ( {Comps(BB), Comps(CC)} , {ux, vx, uy, vy, uz, vz});

> TN := subs(solz,evalm(TN));

suspended plate

center of mass

> #G := make_POINT(TN, l/4,0,0):
G := make_POINT(ground, (xA(t)+xB(t)+xC(t))/3, (yA(t)+yB(t)+yC(t))/3, (zA(t)+zB(t)+zC(t))/3):
show(G);

plate inertia

> plate := make_BODY(G,m,Jxx,Jyy,Jzz,Jyz,Jxz,Jxy): show(plate,'synthetic');

gravity force

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

ground vertex definition

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

test

> A0B0 := make_VECTOR(A0,B0): B0C0 := make_VECTOR(B0,C0): C0A0 := make_VECTOR(C0,A0):
dot_prod(A0B0,A0B0), dot_prod(B0C0,B0C0), dot_prod(C0A0,C0A0);

Coordinates and Constraints

dependent coordinates

> q := [xA(t),yA(t),zA(t), xB(t),yB(t),zB(t), xC(t),yC(t),zC(t), wx(t),wy(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];
#simplify(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;

The Lagrange's multipliers are algebraic unknowns, but they may be easily transformed into a differential form

> SL := [seq(lambda||i(t)= gamma||i(t)+tau*diff(gamma||i(t),t), i=1..nops(constraints))];

assembly variables and equations

> all_vars := q union [seq(gamma||i(t), i=1..nops(constraints))];
all_eqns := subs(SL,leqns) union baumgarte:

> #for i from 1 to nops(all_vars) do i,all_vars[i]; all_eqns[i]; od;

ODEs numerical integrations

This problem is to big for beeing solved in Maple, it should be converted into a Fortran program (or MatLab or something else).

The problem complexity is given by both by the number of equations and the number of operations required for their evaluation

change variables

> SX := [seq(op([diff(all_vars[i],t,t)=X2[1],diff(all_vars[i],t)=X1[1],all_vars[i]=X[1]]),i=1..nops(all_vars))]:

evaluate complexity

> nops(all_eqns);

> codegen[cost](subs(SX,evalf(all_eqns)));