sec4.2.mws

Classical Mechanics with Maple

Section 4.2: Systems of Plane Rigid Bodies

Dr. Harald Kammerer
maple@jademountain.de

Initialisation

> restart;

> libname:="C:/mylib/m6dynlib","C:/mylib/m6dynfig",libname:

> with(linalg):with(plots):with(plottools):with(dynamics);with(figures_chapter_4);

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

[LAGRANGE, MI, gravitational_energy, impact, kineti...

[Fig_4_1, Fig_4_10, Fig_4_11, Fig_4_12, Fig_4_13, F...
[Fig_4_1, Fig_4_10, Fig_4_11, Fig_4_12, Fig_4_13, F...

4.2 Plane Rigid Bodies

In this section we derive some geometrical properties of the apportionment of the mass of rigid bodies. These properties don't depend on the motion of the body but they may affect it. Here we consider a special kind of body. We consider bodies with a symmetry plane with respect to their mass distribution, i.e., they are cylindrical (or prismatic). We will call these planar rigid bodies . In this special case it is enough to consider only the symmetry plane.

4.2.1 Kinematics of Planar Motion of Rigid Bodies

Velocity

Now we consider a rigid body that is moving with respect to a fixed (x,y,z) -reference system with origin O . We talk about planar motions when all points of the considered body are moving in parallel planes. Usually we define the z -axis orthogonal to this plane. Additionally we define a local body - fixed (X,Y,Z) -system, again with the Z -axis orthogonal to the plane of motion. So the Z -axis and the z -axis are always parallel. The origin of the moving system is the body-fixed point P . In Fig. 8 we see the situation.

> display(Fig_4_5(),scaling=constrained,axes=none,title="Figure 8");

[Maple Plot]

The position of the rigid body is defined by the vector rP between O and P and the angle phi between the global x -axis and the local X -axis. (For clarity this angle is not shown in the figure above, but it should be clear what we mean).

The velocity of the body is given by the vector

vP = diff(rP,t)

of the point P with respect to the inertial system and the vector of the angular velocity

omega = diff(phi,t)*eZ

As in section 2, when we described the motion in relative coordinates, we consider now an additional point Q . But now the point Q is fixed on the body and so it is fixed in the local coordinate system. Before we continue we repeat some details about the formulation for motions in relative coordinates. As in the prior section, we define the unit vectors in the direction of the x- and the y- axis as

> ex := vector(3,[1, 0, 0])

ex := vector([1, 0, 0])

> ey := vector(3,[0, 1, 0])

ey := vector([0, 1, 0])

The unit vectors in the direction of the X - and the Y- axis of the moving system are defined by

> eX := vector(3,[eXx, eXy, 0])

eX := vector([eXx, eXy, 0])

> eY := vector(3,[eYx, eYy, 0])

eY := vector([eYx, eYy, 0])

As before we desribe the derivative with respect to time by appending " _d " to the variable and saying for example: " r dot".

The current position and the motion of the considered point Q is described by the vector rQ , which is affected by the current position of the moving coordinate system, described by the vector rP , and the relative position of the point Q with respect to the origin P of the moving system, described by the vector r (rho).

> rP := xP*ex+yP*ey

> evalm(rP);

vector([xP, yP, 0])

> rho := Xrho*eX+Yrho*eY

> evalm(rho);

vector([Xrho*eXx+Yrho*eYx, Xrho*eXy+Yrho*eYy, 0])

Xrho*eXx is the x-component of the projection of the vector rho on the X-axis of the moving coordinate system, Yrho*eYx is the x-component of the projection of the vector rho on the Y-axis of the moving coordinate system. Accordingly Xrho*eXy is the y-component of the projection of the vector rho on the X-axis of the moving coordinate system, Yrho*eYy is the y-component of the projection of the vector rho on the Y-axis of the moving coordinate system.

Again as in section 2 we make the substitutions

> subs1 := Xrho*eXx = XQx

> subs2 := Yrho*eYx = YQx

> subs3 := Xrho*eXy = XQy

> subs4 := Yrho*eYy = YQy

For the vector rQ we get

> rQ := rP+rho

rQ := xP*ex+yP*ey+Xrho*eX+Yrho*eY

This yields by use of subs1 to subs4 (we use the appended letter s to distinguish the expression with and without substutution)

> rQs := subs({subs1, subs2, subs3, subs4},evalm(rQ))...

rQs := vector([xP+XQx+YQx, yP+XQy+YQy, 0])

The rotation of the (X,Y,Z) -system is described by the vector

> Omega := vector(3,[0, 0, omega])

Omega := vector([0, 0, omega])

For the velocity of Q we go three componenets in section 2. At first there was the velocity of the point P

> vP := vector(3,[xP_d, yP_d, 0])

vP := vector([xP_d, yP_d, 0])

Next there is the changing of the position of Q inside the relative system. In contrast to section 2 there is in the current situation no relative motion of Q and so we get

> vrel := vector(3,[0, 0, 0])

vrel := vector([0, 0, 0])

At last we have the part of the velocity which is caused by the rotation of the moving coordinate system

> vrot := crossprod(Omega,rho)

vrot := vector([-omega*(Xrho*eXy+Yrho*eYy), omega*(...

or with use of the substitutions above

> vrots := subs({subs1, subs2, subs3, subs4},evalm(vr...

vrots := vector([-omega*(XQy+YQy), omega*(XQx+YQx),...

For the point Q we get the total velocity

> vQ := vP+vrel+vrots

vQ := vP+vrel+vrots

> evalm(vQ)

vector([xP_d-omega*(XQy+YQy), yP_d+omega*(XQx+YQx),...

or written as components

> vx := evalm(vQ)[1]

vx := xP_d-omega*(XQy+YQy)

> vy := evalm(vQ)[2]

vy := yP_d+omega*(XQx+YQx)

> vz := evalm(vQ)[3]

vz := 0

We see that this consideration is a special case of the description of motions in relative coordinates.

Next we consider the same situation as before. But now we choose another point inside the rigid body as the origin of the local coordinate system. As shown in Fig. 9 we define an additional point inside the rigid body. The question is now: Is the motion of the system dependent on the chosen reference points?

> display(Fig_4_6(),scaling=constrained,axes=none,title="Figure 9");

[Maple Plot]

The position of with respect to the fixed coordinate system is given by the vector rP´ , the position of Q with respect to is given by the vector `rho´` . At last the position of with respect to P is described by the vector rPP´ .

With the relations above we want to formulate the velocity of point Q when we describe the situation by use of point instead of point P . At first we have to make some preparatory work, which is the same as we done above but with new variables.

The unit vectors of the local coordinate system in point are defined by

> `eX´` := vector(3,[`eXx´`, `eXy´`, 0])

`eX´` := vector([`eXx´`, `eXy´`, 0])

> `eY´` := vector(3,[`eYx´`, `eYy´`, 0])

`eY´` := vector([`eYx´`, `eYy´`, 0])

The current position of the point Q is affected by the current position of the moving coordinate system, described by the vector rP´ , and the relative position of the point Q with respect to the origin P´ of the moving system, described by the vector rho´ .

> `rP´` := `xP´`*ex+`yP´`*ey

> evalm(rP´);

vector([`xP´`, `yP´`, 0])

> `rho´` := `Xrho´`*`eX´`+`Yrho´`*`eY´`

> evalm(rho´);

vector([`Xrho´`*`eXx´`+`Yrho´`*`eYx´`, `Xrho´`*`eXy...

The position of with respect to P can be written as

> rPP´:=rho-rho´:

> evalm(rPP´);

vector([Xrho*eXx+Yrho*eYx-`Xrho´`*`eXx´`-`Yrho´`*`e...

Again we make some substitutions

> `subs1´` := `Xrho´`*`eXx´` = `XQx´`

> `subs2´` := `Yrho´`*`eYx´` = `YQx´`

> `subs3´` := `Xrho´`*`eXy´` = `XQy´`

> `subs4´` := `Yrho´`*`eYy´` = `YQy´`

And for the vector rQ´ we get

> `rQ´` := `rP´`+`rho´`

`rQ´` := `xP´`*ex+`yP´`*ey+`Xrho´`*`eX´`+`Yrho´`*`e...

This yields by using of subs1 to subs4 (we use the appended letter s to distinguish the expression with and without substutution)

> `rQs´` := subs({`subs1´`, `subs2´`, `subs3´`, `subs...

`rQs´` := vector([`xP´`+`XQx´`+`YQx´`, `yP´`+`XQy´`...

The rotation of the ( X´,Y´,Z´)- system with respect to the fixed system is described by the vector

> `Omega´` := vector(3,[0, 0, `omega´`])

`Omega´` := vector([0, 0, `omega´`])

The velocity of the point is

> `vP´` := vector(3,[`xP´_d`, `yP´_d`, 0])

`vP´` := vector([`xP´_d`, `yP´_d`, 0])

The component of the velocity which is caused by the rotation of the moving ( X´,Y´,Z´)- coordinate system is

> `vrot´` := crossprod(`Omega´`,`rho´`)

`vrot´` := vector([-`omega´`*(`Xrho´`*`eXy´`+`Yrho´...

or with use of the substitutions above

> `vrots´` := subs({`subs1´`, `subs2´`, `subs3´`, `su...

`vrots´` := vector([-`omega´`*(`XQy´`+`YQy´`), `ome...

The total velocity of Q is

> `vQ´` := `vP´`+`vrots´`

`vQ´` := `vP´`+`vrots´`

> evalm(`vQ´`)

vector([`xP´_d`-`omega´`*(`XQy´`+`YQy´`), `yP´_d`+`...

or written as components

> `vx´` := evalm(`vQ´`)[1]

`vx´` := `xP´_d`-`omega´`*(`XQy´`+`YQy´`)

> `vy´` := evalm(`vQ´`)[2]

`vy´` := `yP´_d`+`omega´`*(`XQx´`+`YQx´`)

> `vz´` := evalm(`vQ´`)[3]

`vz´` := 0

With respect to P we get the velocity of

> vP´:=vP+crossprod(Omega´,rPP´);

`vP´` := vP+vector([-`omega´`*(Xrho*eXy+Yrho*eYy-`X...

This yields for the velocity of Q

> evalm(vQ´);

vector([xP_d-`omega´`*(Xrho*eXy+Yrho*eYy-`Xrho´`*`e...
vector([xP_d-`omega´`*(Xrho*eXy+Yrho*eYy-`Xrho´`*`e...

We compare this result with the velocity of Q by use of the description with respect to point P

> evalm(vQ);

vector([xP_d-omega*(XQy+YQy), yP_d+omega*(XQx+YQx),...

Of course all components of the velocity must be the same, independent of the reference point, so we have the equations

> eq1:=evalm(vQ)[1]=evalm(vQ´)[1];

eq1 := xP_d-omega*(XQy+YQy) = xP_d-`omega´`*(Xrho*e...

> eq2:=evalm(vQ)[2]=evalm(vQ´)[2];

eq2 := yP_d+omega*(XQx+YQx) = yP_d+`omega´`*(Xrho*e...

for the x and the y component. The solution of both equation yields with use of the substitutions above

> solve({subs({`subs1´`, `subs2´`, `subs3´`, `subs4´`,subs1, subs2, subs3, subs4},eq1)},{omega´});

{`omega´` = omega}

> solve({subs({`subs1´`, `subs2´`, `subs3´`, `subs4´`,subs1, subs2, subs3, subs4},eq2)},{omega´});

{`omega´` = omega}

That means that the rotational velocity of a rigid body is independent of the position of the reference points .

Next we consider the question: Is there a point M which ha s zero velocity?

When we know the velocity of point P we know that the velocity of point M is defined by

v[M] = v[P]+crossprod(Omega,rho[M])

When we postulate that the velocity of point M is v[M] = 0 we can write

v[P] = -crossprod(Omega,rho[M])

When we supposed that v[P] and Omega are known then rho[M] is defined by this relation. rho[M] is the vector from point P to M . In case we know the velocity of a second point we get point M from the intersection of both orthogonal lines to the velocity vectors v[P] and v[`P´`] as shown in Fig. 10.

> display(Fig_4_7(),scaling=constrained,axes=none,title="Figure 10");

[Maple Plot]

Otherwise we get the velocity of every point by the relation

v[P] = crossprod(Omega,R)

where R is the vector from M to the point P. Every general motion can momentarily be considered as a rotation around the momentary center of rotation . The word "momentary" shows that the center of rotation in general changes its position during the motion. This fact is shown on the basis of the following example.

Example: Man on Ladder

Let's consider the following situation: A man climbs a ladder. When he reaches the top point, the one end of the ladder slips off. Don't worry, nothing happens to the man, as you can see in Fig. 11, where the situation is animated. (Click again the picture and press "play" in the menu bar.)

> display(Fig_4_8(),insequence=true,scaling=constrained,axes=none,title="Figure 11");

[Maple Plot]

The motion of the left part of the ladder is clear, we have a rotation around the left support of the ladder. This point is fixed.

The motion of the right part of the ladder is a combination of a translation and a rotation. As described above we can formulate every momentary situation during the motion as a rotation around a momentary center of rotation. We get this momentary center of rotation M as follows.

Because the motion of the left part of the ladder is known as a rotation around the left support we know the motion of the joint between both parts of the ladder. The velocity of this joint points orthogonal to the vector from the left support to the joint. Because this joint also belongs to the right part of the ladder, the momentary center of the rotation of this part of the ladder must lie on the line orthogonal to the velocity vector of the joint.

Second, we know that the right support of the ladder slips on the ground, which in this example should be the horizontal. The velocity vector of the right support points in the horizontal direction, too. The momentary center of rotation lies again on the orthogonal line to the velocity vector, which means in this case on a vertical line at the right support.

Now we know the position of the momentary center of rotation of the right part of the ladder: it lies on the intersection of both lines. In Fig. 12 these two lines are drawn dotted in red. The momentary center of rotation is marked by the red point. When you click on the figure and press "play" on the menu bar you see the trajectory of the position of the momentary center of rotation. It is a circular arc around the left support of the ladder with the radius the double the length of the ladder.

> display(Fig_4_9(),insequence=true,scaling=constrained,axes=none,title="Figure 12");

[Maple Plot]

Finally it should be said that the motion of the man is not a rotation but a translational motion on a circular arc. For completeness it should also be said that a translational motion can be considered as a rotation with infinite distance between the momentary center of rotation and the body.

Acceleration

Next we calculate the acceleration of an arbitrary fixed point Q of a rigid body with respect to the inertial system. We know about the velocity of Q

vQ := vP+crossprod(Omega,rho)

Notice that vQ, vP, aQ, aP and Omega are vectors.

Differentiating with respect to time yields

aQ := diff(vQ,t)

aQ := aP+crossprod(diff(Omega,t),rho)+crossprod(Ome...

As shown above it is diff(rho,t) = crossprod(Omega,rho) and so

aQ := aP+crossprod(diff(Omega,t),rho)+crossprod(Ome...

or because

Omega = diff(phi,t)*eZ

and so

aQ := crossprod(Omega,crossprod(Omega,rho)) = -diff...

we can write

aQ := aP+crossprod(diff(Omega,t),rho)-diff(phi,t)^2...

Remember: eZ is the unit vector in Z -direction.

This formulations are not exactly as they are to be used in MAPLE, but for better presentation we make some simplifications.

The total acceleration of Q so has three components: first the velocity of point P , the radial acceleration -diff(phi,t)^2*rho and the tangential acceleration crossprod(diff(Omega,t),rho) . The last two parts result from the rotation of point Q around point P . We say diff(Omega,t) = diff(phi,`$`(t,2))*eZ is the angular acceleration of the body with respect to the inertial system.

As the relations for the velocity this results are analogous to the relations for the relative motion in section 2, but without any relative motion between point Q and the local body-fixed (X,Y,Z) -system.

4.2.1 The Center of Mass

Here it is assumed that the center of gravity is known from the statics. The center of mass is defined equivalently, but we constitute the formulation with respect to infinite mass particles of the body instead of infinite partial areas.

Consider the planar rigid body shown in Fig. 13. We choose an arbitrary (X,Y,Z)- coordinate system with origin P and the Z -axis perpendicular to the symmetry plane of the body.

> display(Fig_4_10(),scaling=constrained,axes=none,title="Figure 13");

[Maple Plot]

The center of mass is defined by

rho[S] := int(rho,m)/int(1,m)

Where int(1,m) = m is the total mass of the body. " rho " is written "rho" in the figure. The integrals in the equation are taken over the total mass of the body. In case that P = S it is rho[S] = 0 and so int(rho,m) = 0 .

When the mass of the body is uniform over the area the center of mass is the same point as the center of gravity.

4.2.2 Moment of Inertia

The moment of inertia of a rigid body is defined by

J[P] := int(rho*rho,m) = int(X^2+Y^2,m)

with the integral over the total mass of the body.

In the special case of prismatic bodies with constant height h and density gamma we get

dm = gamma*h*dA

with dA the cross section of the element dm in the (X,Y) -plane.

So we get the relation

J[P] = gamma*h*I[P,polar]

with the polar second moment of area I[P,polar] what is known from the subject " Strength of materials ".

Theorem of Steiner

Let J[S] be the known moment of inertia with respect to the center of mass S . We are now asking for the moment of inertia of the body with respect to an arbitrary point P . The situation is shown in Fig. 14.

> display(Fig_4_11(),scaling=constrained,axes=none,title="Figure 14");

[Maple Plot]

From the figure we get the relations

X = x-a

Y = y-b

a^2+b^2 = r^2

The integrals are over the entire surface. From the definition above we know

J[P] = int(X^2+Y^2,m)

J[P] := int(x^2+y^2,m)+(a^2+b^2)*int(1,m)-2*a*int(x...

The last two parts 2*a*int(x,m) and 2*b*int(y,m) are zero because x and y are the coordinates with respect to the center of mass. Then we get the theorem of Steiner

J[P] = J[S]+m*r^2

Notice that the moment of inertia with respect to the center of mass has always the smallest value in comparison with the moment of inertia with respect to any other point.

In the following examples we deduce the moments of inertia of some special bodies. The results of these examples and the moments of inertia for some more plane rigid bodies are included in the dynamics package. See MI for more information.

Example: Slender Rod

We consider the rod shown in Fig. 15. The dimensions of the rod orthogonal to the X-axis shall be negligible.

> display(Fig_4_12(),scaling=constrained,axes=none,title="Figure 15");

[Maple Plot]

We are at first interested in the moment of inertia of the rod with respect to point P . Therefore we consider a small mass element

dm = m/L*dX

Here m is the total mass of the rod and L is its length. Substitute this into

J[P] := int(X^2,m)

yields

> J[P] := int(X^2*m/L,X = 0 .. L)

J[P] := 1/3*m*L^2

The distance between the points P and S is L/2 . For the moment of inertia with respect to the center of mass we get from the theorem of Steiner

> J[S] := J[P]-m*(L/2)^2

J[S] := 1/12*m*L^2

Because we consider only planar systems, the depth of the rod has no influence on the result. So this relation is also valid for thin plates.

Example: Semi-Circular Disk

Next we consider the moment of inertia of a semi-circular disk with radius R , thickness t and mass m as shown in Fig. 16.

> display(Fig_4_13(),scaling=constrained,axes=none,title="Figure 16");

[Maple Plot]

The distance between the center of mass and the center of the circle is given by

> s := 4*R/3/Pi

For the calculations consider Fig. 17.

> display(Fig_4_14(),scaling=constrained,axes=none,title="Figure 17");

[Maple Plot]

The moment of inertia with respect to point M is defined by

J[M] := int(x^2+y^2,m) = int(r^2,m)

The thickness of the disk has to be constant, the material homogenous, so we can write for the infinite mass element which is marked in Fig. 17

dm := t*rho*r*dr*dphi

with the infinitesimal angle dphi (=d f ) .

The density of the material is

> rho := m/t/(R^2*Pi/2)

The integral in the definition of the moment of inertia is now divided into two integrals:

> J[M] := int(int(r^2*t*rho*r,r = 0 .. R),phi = 0 .. ...

J[M] := 1/2*m*R^2

We get the moment of inertia of the disk with respect to the center of mass by use of the theorem of Steiner

> J[S] := J[M]+m*s^2

J[S] := 1/2*m*R^2+16/9*m*R^2/Pi^2

> collect(J[S],{m,R})=evalf(J[S]);

(1/2+16/9/Pi^2)*R^2*m = .6801265486*m*R^2

4.2.3 The Angular Momentum of Rigid Bodies

We continue now the considerations form subsection 4.1.2 (The Angular Momentum for Systems of Mass Particles). One of the final results was the axiom of the angular momentum for rigid bodies

Lo = int(crossprod(r,rd),m)

with rd = diff(r,t) the derivative of r with respect to time.

To calculate the angular momentum we need the position and the velocity of every mass particle of the rigid body. We know from above that the motion of a rigid body is known when the velocity v[P] of an arbitrary point P and the angular velocity Omega is given. For an arbitrary mass particle then

r := r[P]+rho

diff(r,t) := v = v[P]+crossprod(Omega,rho)

is valid. See Fig. 18 for the details.

> display(Fig_4_15(),scaling=constrained,axes=none,title="Figure 18");

[Maple Plot]

The only value which denotes the mass particle is the vector rho , which specifies the position of the mass particle with respect to point P . So we can write

Lo = int(crossprod(r[P]+rho,v[P]+crossprod(Omega,rh...

Lo = crossprod(r[P],v[P])*int(1,m)+crossprod(r[P],c...

Further we know

int(1,m) = m

int(rho,m) = m*rho[S]

with rho[S] the vector from point P to the center of mass S . The use of the local unit vectors eX , eY and eZ yields with rho = X*eX+Y*eY and Omega = omega*eZ (notice: Omega is a vector while omega is a scalar value)

int(crossprod(rho,crossprod(Omega,rho)),m) = int(X^...

int(X^2+Y^2,m)*omega = J[P]*Omega

Now we get for the angular momentum

Lo = crossprod(r[P],v[P])*m+crossprod(r[P],crosspro...

For the first two terms we can write

crossprod(r[P],v[P]+crossprod(Omega,rho[S]))*m = cr...

with v[S] the velocity of the center of mass.

The angular momentum is now

Lo = (crossprod(r[P],v[S])+crossprod(rho[S],v[P]))*...

For the axiom of the angular momentum we need the derivative with respect to time

diff(Lo,t) = (crossprod(r[P],a[S])+crossprod(v[P],v...

Because

crossprod(v[P],v[P]) = 0

and

crossprod(v[P],crossprod(Omega,rho[S])) = -crosspro...

we get

diff(Lo,t) = (crossprod(r[P],a[S])+crossprod(rho[S]...

M[O,res] = M[P,res]+crossprod(r[P],F[res])

Together with the first axiom of the center of gravity we get

M[O,res] = M[P,res]+crossprod(r[P],a[S]*m)

Finally we get for the axiom of the angular momentum

J[P]*diff(Omega,t)+crossprod(rho[S],a[P])*m = M[P,r...

4.2.4 Establish the Equations of Motion

Next we consider two special cases.

Rotation of a rigid body around a fixed axis

In this case the body fixed point P is coincidently with the fixed point O . So we have a[P] = 0 . Further we write omega = diff(phi,t) and get

J[P]*diff(phi,`$`(t,2)) = M[z,P,res]

General plane motion

In this case we choose the center of mass S for the point P , then rho[S] = 0 is valid. Then we get the scalar equation

J[S]*diff(phi,`$`(t,2)) = M[z,S,res]

We name this equation the second axiom of the center of gravity . The first axiom describes in case of a general motion the translation of the center of mass, the second axiom the rotation around the center of mass.

In both special cases we see that the form of the equations is analogous to Newton's law. In case of translations the mass of the body specifies the inertia, in case of rotations the moment of inertia J[P] or J[S] . When we write the equations analogous to section 3.1 we get

M[z,P,res]-J[P]*diff(phi,`$`(t,2)) = 0

and

M[z,S,res]-J[S]*diff(phi,`$`(t,2)) = 0

This means that we can consider the moment caused by the rotational inertia of the body as an active moment and handle the situation similar to the statical exercises of the momentum balance at rigid bodies.

In the special case of rotations around a fixed axis we need only to consider the momentum balance with respect to the rotation axis. In case of general plane motions we need to use the momentum balance and the force balance. When we use cartesian coordinates we get the three scalar conditions for the balance of the body, called the

equations of motion

F[x,res]-m*diff(x,`$`(t,2)) = 0

F[y,res]-m*diff(y,`$`(t,2)) = 0

M[z,S,res]-J[S]*diff(phi,`$`(t,2)) = 0

These three equations are the same as we know from the statics only with the additional part to account for the influence of the inertia.

Example: Cuboid And Cylinder on a Slope

> unassign('J'):

We consider the situation shown in Fig. 19. There are two rigid bodies, a cuboid and a cylinder, both on a slope with gradient alpha . Both bodies have the same mass m . The radius of the cylinder is R , the length of the edges of the cuboid is 2*R . In the initial position at t = 0 the distance between the bodies is d . At last we assume that there is always adherence between the cylinder and the ground, and there is no friction between the cuboid and the ground ( mu = 0 ).

The question is now: Will the distance between the bodies increase or decrease? And when it decreases: After what time do the bodies have contact?

> display(Fig_4_16()[1],scaling=constrained,axes=none,title="Figure 19");

[Maple Plot]

First we consider the motion of the cylinder. We choose the coordinate system as shown in Fig. 20. In this example we have a general plane motion. There is no rotation around a fixed axis but the momentary pole of the rotation is changing during the motion.

> display(Fig_4_16()[2],scaling=constrained,axes=none,title="Figure 20");

[Maple Plot]

We consider now the three equations of the equilibrium for the forces in x - and y -direction and the moment around the center of mass S . Therefore we cut the cylinder free and consider all active and passive forces, as shown in Fig. 21.

> display(Fig_4_16()[3],scaling=constrained,axes=none,title="Figure 21");

[Maple Plot]

In Fig. 21 we have the inertial forces

> Tx := m*diff(x(t),`$`(t,2))

> Ty := m*diff(y(t),`$`(t,2))

We assume that the acceleration is measured in the same direction as the displacements. So the inertial forces are assumed to point against the coordinate axis when they are positive.

The inertia moment (Do not confuse with the moment of inertia) is

> Mz := J[S]*diff(phi(t),`$`(t,2))

Additionally we have the weight force

> G := m*g

with the acceleration of gravity g and the contact forces N and H . So we get from Fig. 21 the three equations for the equilibrium

> EQ_Fx := Tx-G*sin(alpha)+H = 0

EQ_Fx := m*diff(x(t),`$`(t,2))-m*g*sin(alpha)+H = 0...

> EQ_Fy := Ty+G*cos(alpha)-N = 0

EQ_Fy := m*diff(y(t),`$`(t,2))+m*g*cos(alpha)-N = 0...

> EQ_MS := Mz-H*R = 0

EQ_MS := J[S]*diff(phi(t),`$`(t,2))-H*R = 0

Further we have some additive geometrical conditions. We know that the distance between the center of mass S and the ground is R. This means

> y(t):=0;

y(t) := 0

and so the acceleration in the y -direction is

> diff(y(t),`$`(t,2))

0

From the condition that the cylinder rolls without sliding, we get

> phi(t) := x(t)/R

and so

> diff(phi(t),`$`(t,2))

diff(x(t),`$`(t,2))/R

This yields for the equilibrium equations

> EQ_Fx;

m*diff(x(t),`$`(t,2))-m*g*sin(alpha)+H = 0

> EQ_Fy;

m*g*cos(alpha)-N = 0

> EQ_MS;

J[S]*diff(x(t),`$`(t,2))/R-H*R = 0

From the equation EQ_Fy we get the information that the normal force N is independent from the motion. It has the same value as in the case of statical considerations.We get for N

> N:=solve(EQ_Fy,N);

N := m*g*cos(alpha)

From the equation EQ_MS we get the tangential contact force H

> H:=solve(EQ_MS,H);

H := J[S]*diff(x(t),`$`(t,2))/R^2

This force is dependent on the motion, as we see in the term diff(x(t),`$`(t,2)) in the relation.

From the function MI from the dynamics package we get the moment of inertia

> J[S] := MI('solidcylinder',m,R);

J[S] := 1/2*m*R^2

At last we get from the equation EQ_Fx the acceleration of the cylinder

> a_cylinder(t):=solve(EQ_Fx,diff(x(t),t$2));

a_cylinder(t) := 2/3*g*sin(alpha)

Next we consider the situation of the cuboid. Therefore we reset some variables first.

> unassign('N','Tx','Ty','G','EQ_Fx','EQ_Fy','y(t)'):

The situation of the cuboid is shown in Fig. 22.

> display(Fig_4_16()[4],scaling=constrained,axes=none,title="Figure 22");

[Maple Plot]

Again we consider first the equations of the equilibrium, whereby we don't need to regard any balance of the moments, because there is no rotation. We consider the cut free body in Fig. 23.

> display(Fig_4_16()[5],scaling=constrained,axes=none,title="Figure 23");

[Maple Plot]

Because the mass of the cuboid is again m and we use the same coordinates to describe the situation as before we get the same inertia forces in Fig. 23 as for the cylinder

> Tx := m*diff(x(t),`$`(t,2))

> Ty := m*diff(y(t),`$`(t,2))

And the weight force is again

> G := m*g

We get now the equilibrium equations

> EQ_Fx := Tx-G*sin(alpha) = 0

EQ_Fx := m*diff(x(t),`$`(t,2))-m*g*sin(alpha) = 0

> EQ_Fy := Ty+G*cos(alpha)-N = 0

EQ_Fy := m*diff(y(t),`$`(t,2))+m*g*cos(alpha)-N = 0...

Further we have same geometrical condition for the motion in the y -directuion as before, that is that the distance between the center of mass S and the ground is R. This means

> y(t):=0;

y(t) := 0

and so the acceleration is y -direction is

> diff(y(t),`$`(t,2))

0

This yields for the equilibrium equations

> EQ_Fx;

m*diff(x(t),`$`(t,2))-m*g*sin(alpha) = 0

> EQ_Fy;

m*g*cos(alpha)-N = 0

From the equation EQ_Fy we get the normal force N

> N:=solve(EQ_Fy,N);

N := m*g*cos(alpha)

This is the same result as that for the motion of the cylinder.

And from the equation EQ_Fx we get the acceleration of the cuboid

> a_cuboid(t):=solve(EQ_Fx,diff(x(t),t$2));

a_cuboid(t) := g*sin(alpha)

Comparison with the acceleration of the cylinder

> a_cylinder(t);

2/3*g*sin(alpha)

shows that the accelerations are different. The rolling cylinder is slower than the cuboid. So we need to answer the question, when do the bodies make contact? . To answer this question we have to calculate the displacement as a function of time for both bodies. To get this we integrate the acceleration twice.

> v_cylinder(t):=int(a_cylinder(t),t)+C1;

v_cylinder(t) := 2/3*g*sin(alpha)*t+C1

> s_cylinder(t):=int(v_cylinder(t),t)+C2;

s_cylinder(t) := 1/3*g*sin(alpha)*t^2+C1*t+C2

> v_cuboid(t):=int(a_cuboid(t),t)+D1;

v_cuboid(t) := g*sin(alpha)*t+D1

> s_cuboid(t):=int(v_cuboid(t),t)+D2;

s_cuboid(t) := 1/2*g*sin(alpha)*t^2+D1*t+D2

Here we use the coordinate s_cylinder and s_cuboid to describe the displacement and assume that both coordinates start in both bodies' center of mass at the beginning. From the integration we get two constants for the motion of each body. To compute these we need four inital conditions. We assume that both bodies are at rest at the beginning of the time period. That means

> IC1:=subs(t=0,v_cylinder(t))=0;

IC1 := C1 = 0

> IC2:=subs(t=0,s_cylinder(t))=0;

IC2 := C2 = 0

> IC3:=subs(t=0,v_cuboid(t))=0;

IC3 := D1 = 0

> IC4:=subs(t=0,s_cuboid(t))=0;

IC4 := D2 = 0

> sol:=solve({IC1,IC2,IC3,IC4},{C1,C2,D1,D2});

sol := {C1 = 0, D1 = 0, D2 = 0, C2 = 0}

So we get for the displacement of the cylinder

> s_cyl(t):=subs(sol,s_cylinder(t));

s_cyl(t) := 1/3*g*sin(alpha)*t^2

and for the cuboid

> s_cub(t):=subs(sol,s_cuboid(t));

s_cub(t) := 1/2*g*sin(alpha)*t^2

We must compute, after what time is the displacement of the cuboid by the value d (initial distance between both bodies) larger than the displacement of the cylinder? This means

> Cond:=s_cub(t)=s_cyl(t)+d;

Cond := 1/2*g*sin(alpha)*t^2 = 1/3*g*sin(alpha)*t^2...

And the solution of this equation yields

> T:=solve(Cond,t);

T := 1/g/sin(alpha)*6^(1/2)*(g*sin(alpha)*d)^(1/2),...

We know that time has to be positive, so we need only consider the positive result.

At last we consider the special case of alpha = 30 °, that means

> alpha := Pi/180*30;

alpha := 1/6*Pi

and we assume for the initial distance

> d:=10:

For the acceleration of gravity we know

> g:=9.81:

and so we get the value

> Time:=evalf(abs(T[1]));

Time := 3.497487084

We see that the motion of both bodies and in the end the time of the contact don't depend on the mass of the bodies.

4.2.5 Balance of Energy, Conservation of Energy

> unassign('J'):

In section 3.2 we had considered the balance of energy and the conservation of energy of mass particles. Now we continue this for rigid bodies. Fig. 24 shows an arbitrary rigid body with the notation which is used now.

> display(Fig_4_17(),scaling=constrained,axes=none,title="Figure 24");

[Maple Plot]

We describe the general motion of the body by the translation of the center of mass and the rotation of the body around its center of mass. The velocity of the center of mass is described by the variable vV[S] (vS in the figure). To avoid confusion we use in this section the notation vV to describe the velocity vectors and v to describe the scalar value of the velocity. The rotation is described by the vector Omega = omega*eZ with the unit vector eZ along the Z-axis.

For an arbitrary mass element of the rigid body (in the figure yellow colered) we get the velocity vector

vV = vV[S]+crossprod(Omega,rho)

or with crossprod(Omega,rho) = vV[Rot] we can write

vV = vV[S]+vV[Rot]

Considering the mass element as a mass particle yields the kinetic energy (see subsection 3.2)

> T := kinetic_energy(v,dm)

T := 1/2*v*dm^2

For the total rigid body we get the kinetic energy

T := int(v^2,m)/2

The square of the velocity is given by

v^2 = dotprod(vV,vV)

v^2 = dotprod(vV[S],vV[S])+2*dotprod(vV[S],crosspro...

We use the relations

Omega = omega*eZ

rho = X*eX+Y*eY

with the unit vectors of the body fixed coordinate system eX , eY and eZ and thus

crossprod(Omega,rho) = omega*(X*eX+Y*eY)

dotprod(crossprod(Omega,rho),crossprod(Omega,rho)) ...

So we get from the description with use of vectors to the scalar notation. Additionally we use the known relations

int(1,m) = m, int(X,m) = 0, int(Y,m) = 0, int(X^2+y...

Thus we get for the kinetic energy

T := 1/2*m*v[S]^2+1/2*J[S]*omega^2

The kinetic energy of a rigid body is seperated into two parts:

The kinetic energy of the translation with velocity v[S]

T[Trans] := 1/2*m*v[S]^2

and the kinetic energy of the rotational motion of the body around its center of mass with the rotational velocity omega

T[Rot] := 1*J[S]*omega^2/2

We can use the function kinetic_energy from the package dynamics to calculate both parts, the translational and the rotational kinetic energy, because the format of both expressions is the same.

> T[Trans] := kinetic_energy(m,v[S])

T[Trans] := 1/2*m*v[S]^2

> T[Rot] := kinetic_energy(J[S],omega)

T[Rot] := 1/2*J[S]*omega^2

In the case that we consider a motion of a rigid body in form of a rotation around a fixed axis through point P it is convenient to describe the motion with respect to this axis. With the vector rho[P] from S to P we know from 4.2.1 the velocity of the center of mass

V[S] = crossprod(Omega,rho[P])

and thus

v[S]^2 = omega^2*rho[P]^2

For the kinetic energy we get

T = (J[S]+m*rho[P]^2)*omega^2/2

and according the axiom of Steiner this means

> T = kinetic_energy(J[P],omega)

T = 1/2*J[P]*omega^2

As we see we can use again our function from the dynamics -package for the calculation.

Next we consider the balance of energy of the rigid body. Here we regard the rigid body as a system of mass particles dm . We know that the passive forces do no work (see 3.2). Further we know from 4.1.1 that the inner forces are equal by pairs and point against each other. At last there is no relative motion between the mass particles, by definition of a rigid body. The balance of energy for the rigid body is given by the summation of the work of all n bodies acting forces on the one side of the equation, and the integration of the kinetic energy over the total mass of the body on the other side

sum(int(F[i],r[i] = r[i1] .. r[i2]),i = 1 .. n) = i...

The indices 1 and 2 denote here two different states.

For the right side of the equation we can write according to the results above

T[2,Trans]+T[2,Rot]-T[1,Trans]-T[1,Rot]

The left side of the equation is the work of all active forces on the path from 1 to 2 . It will be rearranged as follows:

With

dr[i] = vV[i]*dt, vV[i] = vV[S]+crossprod(Omega,rho...

we get

sum(int(F[i]*vV[S],t = t[1] .. t[2])+int(dotprod(F[...

(exactly the term in the first integral on the right side F[i]*dr[S] should be written as dotproduct)

Under consideration of

sum(F[i],i = 1 .. n) = F[Res]

sum(crossprod(rho[i],F[i]),i = 1 .. n) = M[Res]

Omega*dt = eZ*d*phi

dotprod(M[Res],eZ) = scalar(M[Res])

Notice: The last relation means in words: The dotproduct of the resultant moment vector and the unit vector in Z-direction is equal to the scalar value of the resultant moment.

With these relations we get

int(F[Res],r[S] = r[S1] .. r[S2])+int(M[Res],phi = ...

Here we have taken in that the active forces are the combination of conservative and non-conservative forces (see subsection 3.2).

The energy balance has now the same format as in the consideration of mass particles

T[2]+V[2] = T[1]+V[1]+W[1 .. 2]

But now we have to see that the kinetic energy is composed of the two parts

T = T[Trans]+T[Rot]

Beside this we have the work

W[1 .. 2] = int(F[Res],r[S] = r[S1] .. r[S2])+int(M...

which includes the work of non-conservative moments.

The conservation of energy results when the work W[1 .. 2] vanishes.

Example: Connected Cuboid And Cylinder

> unassign('H','R','N','L','s','g','alpha','W(s)','omega','E','v'):

We consider a system with a cuboid and a cylinder, each with mass m and connected by a rope. The cylinder lies on a slope with angle alpha . Gravity causes the system to move. Have a look at Fig. 25, where the initial situation is shown. Here the velocity of all parts of the system is zero. Assume that the first part of the surface where the cuboid lies has no friction (green color) and on the second part there is friction with coefficient mu between ground and cuboid (brown color). Between the cylinder and the ground always exist adherence. Additionally we assume that the rope never sags and the displacement of the cuboid is the same as the displacement of the cylinder. So we use only one coordinate, s , to describe the momentary situation.

> display(Fig_4_18()[1],scaling=constrained,axes=none,title="Figure 25");

[Maple Plot]

It is important to see what happens while the cuboid passes the transition from the part without friction to the part with friction, That means the situation between s = L1 and s = L1+2*R . These situations are shown in the Fig.'s 26 and 27.

> display(Fig_4_18()[2],scaling=constrained,axes=none,title="Figure 26");

[Maple Plot]

> display(Fig_4_18()[3],scaling=constrained,axes=none,title="Figure 27");

[Maple Plot]

Here we think it should be a good approximation of reality to assume that the friction force increases linearly from zero up to the maximum value. Further we assume that the friction force is not dependent on the velocity but solely on the law for the friction force FR , given by

FR := mu*N

with the normal force

N := m*g

We can write the friction force as a function of the position of the cuboid by

> FR(s) := piecewise(s < L1,0,L1 <= s and s < L1+2*R,...

FR(s) := PIECEWISE([0, s < L1],[1/2*mu*m*g/R*(s-L1)...

We will see later in this example by use of some concrete values that this is correct.

The end of the motion is reached when the friction force has dissipated all energy or when the cylinder runs against the wall at the end of the slope as shown in Fig. 28. That means that s = L3 is reached.

> display(Fig_4_18()[4],scaling=constrained,axes=none,title="Figure 28");

[Maple Plot]

To formulate the expressions for the energy of the system in the different states we need to define a reference horizontal. The altitude of the cuboid doesn't change during the motion, so we don't need to consider it. It seems to be useful to define the position o