> restart: libname := libname, "C:/MBSymba": with(MBSymba_r6):
> PDEtools[declare]([theta(t)],prime=t,quiet):
Definition of the mechanical systems
Define the body frame
> TP := rotate('Z',theta(t)) * translate(0,-L,0);
Define the body
> pendulum := make_BODY(TP,m,0,0,IG): show(pendulum);
Define acceleration due to gravity
> _gravity := make_VECTOR(ground,0,-g,0): show(_gravity);
Define the pin joint A
> A := origin(ground): show(A);
Reaction force (body coordinates)
> make_VECTOR(TP,TA,RA,0):
RF := make_FORCE(%, A, pendulum): show(RF);
Define the damping torque
> make_VECTOR(ground,0,0,-c*diff(theta(t),t)):
DT := make_TORQUE(%, pendulum): show(DT);
Newton's equations
> Describe(newton_equations);
> eqnsN := newton_equations({ pendulum, RF, DT}): show(eqnsN);
# Compute the Newton's equations of the multibody system <MBS>, which may
# include bodies and forces.
newton_equations( MBS::{set, list}, verbose::boolean := false, $ ) :: VECTOR
> eqnsN := project(eqnsN,TP): show(eqnsN);
Euler's equations with respect to the fixed point A, in moving frame
> Describe(euler_equations);
# Euler's equations of the multibody system <MBS> with respect the (optional)
# point <AA>.
euler_equations( MBS::{set, OBJECT, list}, AA::POINT := NULL,
verbose::boolean := false, $ )
> eqnsE := euler_equations({ pendulum, RF, DT},CoM(pendulum)): show(eqnsE);
> eqnsE := euler_equations({ pendulum, RF, DT},A): show(eqnsE);
Equation(s) of motion for the system
> eqn := comp_Z(eqnsE);
Numerical integration
Integration of the equation of motion
> data := [m=4, L=0.2, g=9.81, IG=0.12, c=0.5];
Define the initial conditions
> ic := theta(0)=2.8, D(theta)(0)=0.0;
Integrate the equation(s) of motion
> motion := dsolve( { subs(data,eqn) , ic }, type=numeric, range=0..5):
> plots[odeplot](motion, [[t,theta(t)],[t,D(theta)(t)]], color=[red,blue],numpoints=200);
Reaction forces
Netwon's equations
> Describe(euler_equations);
eqnsN := newton_equations({ pendulum, RF, DT}): show(eqnsN);
# Euler's equations of the multibody system <MBS> with respect the (optional)
# point <AA>.
euler_equations( MBS::{set, OBJECT, list}, AA::POINT := NULL,
verbose::boolean := false, $ )
Calculation of the reaction forces
> RA = solve(comp_Y(eqnsN),RA);
> TA = solve(comp_X(eqnsN),TA);
Lagrange's equations
> KE := simplify(kinetic_energy(pendulum));
> PE := gravitational_energy(pendulum);
> eqnL := lagrange(KE-PE,theta(t),t) - generalized_force({DT},theta(t));
sintetically
> eqnL := lagrange_equations({pendulum,DT},[theta(t)]);
>
>