(C) 2004 Roberto Lot

Step 1 - Define the physical constituents of the systems: points, bodies, etc.

Define the pin joint A

> A := origin(ground): show(A);

Define the body frame

> BF := rotate('Z',theta(t)) * translate(L,0,0);

Define the body

> pendulum := make_BODY(BF,m,Ix,Iy,Iz): show(pendulum);

Define acceleration due to gravity

> _gravity := make_VECTOR(ground,g,0,0): show(_gravity);

Reaction force (body coordinates)

> RF := make_FORCE(BF,RA,TA,0, A, pendulum): show(RF);

Define the externally applied torque

> AT := make_TORQUE(ground,0,0,MZ, pendulum): show(AT);

Step 2 - Once the definition of the system is complete you can easily derive the equation of motion

mbs := { pendulum, RF, AT}; short-hand summary of the system model (body, force, torque)

Netwon's equations

> eqnsN := newton_equations(mbs,BF): show(eqnsN);

Euler's equations with respect to center of mass, in ground frame

> eqnsE := euler_equations(mbs, CoM(pendulum), ground): show(eqnsE);

Euler's equations with respect to the point A, in moving frame

> eqnsE := euler_equations(mbs,A,BF): show(eqnsE);

Equation(s) of motion for the system

> eqn := Zcomp(eqnsE);

Step 3 - Numerical intergration and solution for specified values

Integration of the equation of motion

> data := [m=4, L=0.2, g=9.81, Iz=0.12, MZ=.2+2.1*t];

> ode := subs(data,eqn);

Define the time domain range

> t0:=1; tend:=3; nt:=100;

> tvec := array([seq(1.*i/nt*(tend-t0)+t0 , i=0..nt)]): creates a sequence of numerical integration steps

Define the initial conditions

> ic := theta(t0)=0.2, D(theta)(t0)=0.0;

Integrate the equation(s) of motion

> motion := dsolve( { ode , ic }, type=numeric, output=tvec):

> outvar := evalm(motion[1,1]);

> solz := evalm(motion[2,1]):

Plot the solution

plot([ [[solz[n,1],solz[n,2]] $n=1..nt+1], [[solz[n,1],solz[n,3]] $n=1..nt+1] ] , color=[red,green]);

Calculation of the reaction forces

> RA = solve(Xcomp(eqnsN),RA);

> TA = solve(Ycomp(eqnsN),TA);