Plane motion of a rigid body

    Scilab is also well adapted for solving the equations of motion of a body of a system of rigid bodies. These equations are systems of differential equations, linear or non linear, of the second order with respect to time.

    As an example we consider the motion of a rigid body on an inclined plane (figure hereafter). This problem is analysed in Chapter 23 (Section 23.2) of the text-book Mechanics of Rigid Bodies de l'ouvrage Mécanique des Solides Rigides (2006), Part V

 
 
 
   

  The equations of motion can be written  (27.30) under the form:                            

 

 

where x and y are the coordinates of the mass center G, ψ is the rotation angle of the body (S), ft and fr the damping coefficients  in translation and rotation, respectively.

   We report hereafter an example of the command file of Scilab and the results obtained for the trajectories:

    -- of the mass center (red curve);

    -- of a point attached to the body (blue curve), with coordinates (4m, 0, 0) relatively to the system attached to the body;

for different values of the inclination of the plane and of the damping coefficients.

    The differential equations are solved by using a Runge-Kutta process of order  4.

    For an extended analysis, please report to Chapter 27 of the text-book Mechanics of Rigid Bodies, Part VI.

 

Fichier.sce

//Plane motion of a rigid body

ta=5 // final time
alpha=20;//inclination of the plane
alpha=alpha*%pi/180;
ft=0.5
fr=0.8
x0=0;y0=0;psi0=0//initial coordinates and angle
vx0=0;vy0=5;vpsi0=360;//initial velocities
g=9.807;
n=200;//step number for the calculation
t0=0;
pas=(ta-t0)/n;
t(1)=t0;x(1)=x0;y(1)=y0;psi(1)=psi0*%pi/180;vx(1)=vx0;vy(1)=vy0;vpsi(1)=vpsi0*%pi/180;
for i=1:n
    ti=t(i);xi=x(i);yi=y(i);psii=psi(i);vxi=vx(i);vyi=vy(i);vpsii=vpsi(i);
    Fvxi=g*sin(alpha)-ft*vxi;
    Fvyi=-ft*vyi;
    Fvpsii=-fr*vpsii;
    phi1x=pas*vxi;phi1y=pas*vyi;phi1psi=pas*vpsii;
    phi1vx=pas*Fvxi;phi1vy=pas*Fvyi;phi1vpsi=pas*Fvpsii;
    ti=ti+pas/2;
    xi=x(i)+phi1x/2;yi=y(i)+phi1y/2;psii=psi(i)+phi1psi/2;
    vxi=vx(i)+phi1vx/2;vyi=vy(i)+phi1vy/2;vpsii=vpsi(i)+phi1vpsi/2;
    Fvxi=g*sin(alpha)-ft*vxi;
    Fvyi=-ft*vyi;
    Fvpsii=-fr*vpsii;
    phi2x=pas*vxi;phi2y=pas*vyi;phi2psi=pas*vpsii;
    phi2vx=pas*Fvxi;phi2vy=pas*Fvyi;phi2vpsi=pas*Fvpsii;
    xi=x(i)+phi2x/2;yi=y(i)+phi2y/2;psii=psi(i)+phi2psi/2;
    vxi=vx(i)+phi2vx/2;vyi=vy(i)+phi2vy/2;vpsii=vpsi(i)+phi2vpsi/2;
    Fvxi=g*sin(alpha)-ft*vxi;
    Fvyi=-ft*vyi;
    Fvpsii=-fr*vpsii;
    phi3x=pas*vxi;phi3y=pas*vyi;phi3psi=pas*vpsii;
    phi3vx=pas*Fvxi;phi3vy=pas*Fvyi;phi3vpsi=pas*Fvpsii;
    ti=ti+pas/2;
    xi=x(i)+phi3x;yi=y(i)+phi3y;psii=psi(i)+phi3psi;
    vxi=vx(i)+phi3vx;vyi=vy(i)+phi3vy;vpsii=vpsi(i)+phi3vpsi;
    Fvxi=g*sin(alpha)-ft*vxi;
    Fvyi=-ft*vyi;
    Fvpsii=-fr*vpsii;
    phi4x=pas*vxi;phi4y=pas*vyi;phi4psi=pas*vpsii;
    phi4vx=pas*Fvxi;phi4vy=pas*Fvyi;phi4vpsi=pas*Fvpsii;

    t(i+1)=ti;
    x(i+1)=x(i)+(phi1x+2*phi2x+2*phi3x+phi4x)/6;
    y(i+1)=y(i)+(phi1y+2*phi2y+2*phi3y+phi4y)/6;
    psi(i+1)=psi(i)+(phi1psi+2*phi2psi+2*phi3psi+phi4psi)/6;
    vx(i+1)=vx(i)+(phi1vx+2*phi2vx+2*phi3vx+phi4vx)/6;
    vy(i+1)=vy(i)+(phi1vy+2*phi2vy+2*phi3vy+phi4vy)/6;
    vpsi(i+1)=vpsi(i)+(phi1vpsi+2*phi2vpsi+2*phi3vpsi+phi4vpsi)/6;
end
//motion of the point M(xs,ys)
xs=4;ys=0;
for i=1:n
    ci=cos(psi(i));si=sin(psi(i));
    xm(i)=x(i)+xs*ci-ys*si;
    ym(i)=y(i)+xs*si+ys*ci;
end
clf

plot(xm,ym,"b",x,y,"r")
xtitle("Inclinaison = 20°  fr = 0.5 /s  ft = 0.8 /s","m","m");

 

Trajectories obtained