4

Dynamic Analysis. Mass Matrices and External Forces

The formulation of the inertia and external forces appearing at any of the elements of a multibody system, in terms of the dependent coordinates that describe their position, velocity, and acceleration, is of fundamental importance for the solution of the dynamic analysis. The external and inertia forces of a body subjected to an acceleration field may be expressed in a diversity of ways. The most general way is finding the resultant force and torque about a specific point of the element. However, there are many ways of doing this. All of the ways are based on representing the inertia and external forces by means of equivalent force systems (the same resultant force and torque about any point). This will depend on the type of representation (coordinates) used for the multibody system. We deal in this chapter with the representation of the inertia and external forces generated in the elements of planar and three-dimensional multibody systems that are characterized by natural coordinates. The formulation of the simpler planar element is the start, which will serve as an introduction to the more complicated three-dimensional development. The reader needs a minimum background in analytical dynamics for the understanding of this and subsequent chapters. For this reason, a background on this topic is provided in the first section, that can certainly be skipped by those with a sufficient knowledge.

4.1 Background on Analytical Dynamics

This section is intended to provide the reader with a basic background on some fundamental principles of analytical dynamics that are important for the understanding of the rest of this chapter and subsequent ones. The interested reader who wants to attain a deeper knowledge on this topic is referred to other works dedicated to the study of classical mechanics, such as Goldstein (1980), and Bastero and Casellas (1976).

120

4.1 Background on Analytical Dynamics

121

y x q dq L

dy dx
Figure 4.1. A virtual displacement on a single pendulum.

4.1.1 Principle of Virtual Displacements

The principle of virtual displacements is a powerful principle which is suitable for the dynamic analysis of connected rigid and flexible multibody systems. Prior to the definition of the principle, we need to introduce some concepts. A virtual displacement is defined as an infinitesimal imaginary change of configuration of a system at a stationary time that is consistent with its boundary and constraint conditions. If the configuration is represented by the position vector q, the vector of generalized virtual displacement is customarily denoted by d q. Figure 4.1 shows a single pendulum of length L with a mass m at its tip. If the system is characterized by the independent coordinate q, a virtual displacement dq consistent with the boundary condition is simply an imaginary infinitesimal rotation at the hinge. If the system is characterized by the dependent Cartesian coordinates x and y of its tip, the virtual displacements dx and dy are not independent but interrelated through the constraint condition:

f x2 + y2 L2 = 0
The relationship between dx and dy can be obtained by imposing the condition that a virtual variation of the constraint be zero:

df = f q dq = 2x

2y

dx = 0 dy

(4.1)

Virtual quantities operate the same way as variations. Without entering into the details of calculus of variations (Reddy (1984)), a virtual quantity or variation d, acts like a differential operator but with respect to the dependent variables only, since the time variable is considered fixed. In other words, the multibody

122

4. Dynamic Analysis. Mass Matrices and External Forces

system is considered at a stationary position. As an example, if f is a function of q, q, and t, its variation is

df = f q d q + f q d q
whereas its differential is df = f q dq + f q dq + ft dt

(4.2)

(4.3)

The laws of variation of sums, products, ratios, and so forth are the same as those of differentiation. In addition, the variational operator can be interchanged with the differential and integral operators. Virtual work dW is defined as the work done by all the forces acting on a system, including the inertia forces, that undergoes a virtual displacement and can be expressed as
n

dW = Qi dqi
1

(4.4)

Each of the generalized forces Qi represents the virtual work done when d qi=1 and d q j =0 for ji. When the system is characterized by n independent coordinates, the principle of virtual displacements can be defined as
n

dW = Qi dqi = 0
1

(4.5)

meaning that the virtual work of all the forces acting on the system, including the inertia forces, must be zero. In multibody dynamics the forces Qi do not include the reaction forces, because these do not produce any virtual work. The reason being that reaction forces are couples of internal forces that act along the vector connecting their position coordinates but with opposite signs, thus canceling the corresponding virtual work.

4.1.2 Hamilton's Principle

Consider a system characterized by a set of n independent coordinates qi. Let L=TV be the system Lagrangian, where T and V are the kinetic and potential energy, respectively, and Wnc is the work done by the non-conservative forces. Hamilton's principle (Hamilton (1834)) establishes that the motion of the system from time t1 to time t2, at which the motion is specified, is such that the integral action A=
t2 t1

L dt +

t2 t1

Wnc dt

(4.6)

has a stationary value for the correct path of the motion. This means that the variation of the action A has to vanish:

4.1 Background on Analytical Dynamics

t2 t2

123

dA =
t1

dL dt +
t1

dWnc dt = 0

(4.7)

where the property that the variation is interchangeable with the integral operator has been used. In most cases, the representation of multibody systems is done by means of dependent coordinates that are interrelated through the constraint conditions. Let us assume that the system is characterized by a vector q of n dependent coordinates that satisfy m constraint conditions fk(q, t)=0, which we will assume are of the holonomic type. Hamilton's principle can still be generalized for these cases by means of the Lagrange multipliers technique. Accordingly, the action A is augmented with an additional term:
t2 t2 t2 m

A=
t1

L dt +
t1

W nc dt
t1

( f k l k) dt k=1

(4.8)

where lk are the Lagrange multipliers affected by a minus sign for convenience of the formulation. The stationary condition dA=0 now leads to
t2 t2 t2

dA =
t1

dL dt +
t1

dWnc dt
t1

k = 1 i=1

fk (dqi q
i

lk) dt = 0

(4.9)

The summation in the last term of (4.9) can also be expressed in matrix form as (dqTFT ql ), where l is the vector of the Lagrange multipliers and F q is the Jacobian matrix of the constraints.

4.1.3 Lagrange's Equations

The Lagrange's equations can be directly obtained from Hamilton's principle. . Knowing that T=T(q, q ) and V=V(q) their variations become

dT = T dqi + i=1 qi

T dqi
i=1

qi

T d qT T + d q T q q

(4.10a)

n dV = V dqi dqT V q i=1 qi

(4.10b) (4.10c)

dWnc = dqT Q e x

where Qex represents the external forces and those not coming from a potential. Note that both index and matrix notation have been used simultaneously. Continuing with matrix notation, the application of Hamilton's principle defined by equation (4.9) leads to

124

4. Dynamic Analysis. Mass Matrices and External Forces

t2

t1

T T dqT T V + Q ex F T dt = 0 q l + dq q q q

(4.11)

The last term can be integrated by parts to yield

t2 T dq T dt = dqT T q q t2 t2

t1 t1

t1

dqT d T dt dt q

(4.12)

The first term on the RHS vanishes because the motion is specified at the two ends t1 and t2. Thus the variations will be zero: dq(t1)=dq(t2)=0. The substitution of (4.12) into (4.11) yields
t2

t1

d q T d L L + F T q l Q e x dt = 0 dt q q

(4.13)

Although the coordinates q are not independent, the expression between parenthesis can always be made zero through the selection of the m Lagrange multipliers l . According to the fundamental lemma of the calculus of variations (Reddy (1984)), equation (4.13) leads to d L L + F T q l = Qex dt q q (4.14)

which along with the m constraint equations F (q)=0 constitutes a set of (m+n) differential algebraic equations of motion. In the next example we apply these equations to a mechanical system.
Example 4.1 Use the Lagrange's equations to write the equations of motion of a mechanical sysT tem with kinetic energy T = 1/2 q M q, potential energy V = V (q), external forces Qex, and whose constraint conditions are F (q)=0. The partial derivatives of the kinetic and potential energies are L q then d L = M q q + M q q dt q and the application of equation (4.14) leads to M q q + F q l = Q ex M q q + T q V q It is worth pointing out that the terms M q q and Tq are quadratic in the velocities with coefficients that may depend on q (See Example 4.2). The terms that inT

=Mq q ;

L q

T q

V q

4.1 Background on Analytical Dynamics

125

m2 (x2,y2) m1 (x1,y1) q2 x

q1

Figure 4.2. Mass sliding along a single pendulum. volve q i2 are called centrifugal , and those that involve (q i q j) are called Coriolis terms. The term V q involves q but not its derivatives. Example 4.2 Figure 4.2 depicts a mass m 1 that slides along a massless rod which also has a mass m 2 attached at its tip. Find the equations of motion of this system subject to gravity using the two independent coordinates q1 and q2. The kinetic energy of this system is 2 2 2 2 T = 1 m1 x 1 + y 1 + 1 m2 x 2 + y 2 2 2 Knowing that x 1 = q1 cos q2, y 1 = q1 sin q2, x 2 = L cos q2, and y 2 = L sin q2 the kinetic energy in terms of q1 and q2 becomes q1 m1 0 T = 1 q1 q2 2 0 m 2 L 2 + m 1 q1 2 q2 Similarly, the potential energy is V = m 1 g q1 sin q2 + m 2 g L sin q2 and the application of equation (4.14) leads to m 1q1q2 q1 m 1 g sin q2 = 2 2 m 1q1 + m 2L g cos q2 0 m 2L + m 1q1 q2 2m 1q1q1q2 Since q 1 and q 2 are independent coordinates, note that there are no Lagrange multipliers involved in the equations. Note also that the mass matrix depends on the coordinate q 2 , and this leads to velocity-dependent nonlinear terms on the RHS of the equation. In addition, the gravity forces contain nonlinear terms that involve transcendental functions with q2 as argument. m1 0
2

..

126

4. Dynamic Analysis. Mass Matrices and External Forces

Example 4.3 Repeat Example 4.2 using the Cartesian (natural) coordinates of the two masses m1 and m2 . The kinetic energy of the system in Cartesian coordinates takes a simpler expression: T = 1 x 1y 1x 2y 2 2 m1 0 0 0 0 m1 0 0 0 0 m2 0 0 0 0 m2 x1 y1 x2 y2

Similarly the potential energy is V = m1 g y 1 + m2 g y 2 The constraint conditions to be satisfied by the Cartesian coordinates are
2 2 x2 + y2 L =0 x1 y2 x2 y1 = 0 Since the mass matrix is now constant, the equations of motion take this simpler form: 2

x1 0 0 y2 m y1 1g + 0 x 2 l1 = 0 x2 2x 2 y 1 l2 m 2g 2y 2 x 1 y2 Although the number of equations in Example 4.3 has increased compared to the results of Example 4.2, first the mass matrix and gravity forces are constant. Secondly, the degree of nonlinearity has decreased, since there are neither velocity-dependent terms nor transcendental functions in the RHS. In addition, the Jacobian matrix of the constraints is linear in q . m1 0 0 0 0 m1 0 0 0 0 m2 0 0 0 0 m2

This simple example already illustrates a general fact: the formulation of the equations of motion in independent coordinates leads to a minimum set of highly nonlinear and coupled ordinary differential equations. On the other hand, the formulation using natural coordinates, at the expense of increasing the number of equations, results in a simpler and less coupled set of equations with milder nonlinearities.

4.1.4 Virtual Power

Virtual power also constitutes a powerful principle that will be extensively used in the formulations of this chapter. A virtual velocity vector is defined as a set of imaginary velocities at a stationary time that is consistent with the homogeneous form of the velocity constraint conditions, that is, having no RHS term including the partial derivatives with respect to time. Following the example of Figure 4.1, let q* be a set of virtual velocities which must satisfy the velocity constraint conditions:

4.1 Background on Analytical Dynamics

* *

127

T f = fq q = 2x 2 y

x y

=0

(4.15)

Contrary to the virtual displacements, the virtual velocities need not be infinitesimal since equation (4.15) involves q and not d q. A virtual velocity (finite) is a virtual displacement (infinitesimal) divided by the infinitesimal scalar dt. The principle of virtual power is heavily used for the algorithms and formulations presented in this book; hence it is worthwhile to describe it in some detail. Virtual power can be applied with dependent or with independent coordinates. Both ways will be presented next. Dependent coordinates. If q* constitutes a set of n dependent virtual velocities, the principle of virtual power can be formulated as: W = Fi qi q
i=1 * n * *T

F=0

(4.16)

where F is the vector of all the forces that produce virtual power, including the inertia ones: F=M q Q (4.17)

Vector Q includes the external forces and the velocity-dependent inertia forces (centrifugal and Coriolis), but it does not include internal constraint forces, since they do not produce virtual power. Therefore, equation (4.17) leads to a set of equilibrium equations ( M q Q = 0) in which the internal constraint forces are missing. These forces should appear in the equilibrium equations. In order to find the equilibrium equations from (4.16) and (4.17), we need to add a set of forces T in the direction of the constraint violations (FT q l ), where the columns of F q (rows of F q ) give the direction of constraint forces and l is the vector of their unknown magnitudes. As the virtual velocity vector q* belongs to the nullspace of F q , the product (q*TFT q l ) is zero and can be added to equation (4.16) yielding: W* = q
*T

(M q Q + FT q l) = 0
*

(4.18)

Only nm elements of the virtual velocity vector q can be arbitrarily selected. It is always possible to find the m components l (Lagrange multipliers) in such a way that the parenthesis of (4.18) becomes zero. Consequently, the complete set of force equilibrium equations is: M q + FT q l =Q (4.19)

Equation (4.19) is analogous to the equations of motion of Example 4.1, that were obtained from the Lagrange's equations (4.14). Independent coordinates. The virtual power principle can be applied also with independent virtual velocities. From equations (4.16) and (4.17), it could not be

128

4. Dynamic Analysis. Mass Matrices and External Forces

concluded that ( M q Q = 0), because the multiplying virtual velocities q* were not independent; hence they could not be chosen arbitrarily. However, it is possible to use now the transformations between dependent and independent velocities and accelerations introduced in Chapter 3 as carried out in the following example.
Example 4.4 Starting from equation (4.16), obtain the equations of motion with independent coordinates, using the velocity and acceleration transformations defined in Section 3.5. The virtual velocity vector q* must satisfy the homogeneous version of the velocity constraint equations: Fq q* = 0 and according to equation (3.29) there is an (nxf) matrix R such that (i)

q* = R z* (ii) On the other hand, according to equation (3.32), the following relationship between dependent and independent accelerations can be established: q = R z + Sc (iii) where R is the same matrix that appears in the previous expression and (Sc ) is a term that depends on the actual velocities and can be computed easily. Introducing equation (iii) and the transpose of equation (ii) in equation (4.16), we get z * T RT M R z + M Sc Q = 0 (iv) Since the virtual velocities z * can be chosen arbitrarily, it is possible to conclude that the term that multiplies them in expression (iv) must be zero. Consequently, we arrive at the following set of equations: R T M R z = RT Q M Sc (v)

This is an important result that will be developed with more detail in Chapter 5. In Chapter 8, it will be used as the basis of very efficient dynamic formulations. The application of the virtual power method with independent coordinates is also referred to in the bibliography as Kanes method or Kanes equations (Kane and Levinson (1985), and Huston (1990)). Other authors (Schiehlen (1984)) refer to it as Jourdains principle.

4.1.5 Canonical Equations

The Lagrange's equations lead to a set of n second order differential equations in the coordinates q. Hamilton introduced a transformation that leads to a set of 2n first order differential equations and they are called the canonical or Hamilton's equations. The study of these equations is important because it gives one a further insight into the multibody problems. In addition, they provide an alternative to the acceleration-based formulations at the time of their numerical implementa-

4.1 Background on Analytical Dynamics

129

tion. Also, the canonical equations constitute the foundation for the study of quantum and relativistic mechanics. The canonical momenta is defined as p = L q (4.20)

where L is the Lagrangian and q a set of dependent coordinates that characterize the system. According to this new variable, the Lagrange's equations (4.14) take the following form: p = Q ex + L F T (4.21) q l q The Lagrangian is a function of q, q, and t, and, consequently, its differential becomes T dL = dq L + dqT L + L dt (4.22) q t q Using equations (4.20) and (4.21): L dt dL = dq Tp + dqT p Q ex + F T q l + t (4.23)

Knowing that dq T p = d(p T q) dp T q, equation (4.23) can be transformed into L dt d p T q L = dpT q + dqT Qe x F T q l p t (4.24)

The expression H = p T q L is called the Hamiltonian function. The RHS of (4.24) tells us that it is an explicit function of p, q, and t. Consequently, dH = dp T H + dqT H + dt H p q t (4.25)

Finally, identifying the terms on the RHS of equations (4.24) and (4.25), we arrive at the canonical equations: H = q (4.26) p H = Q F T l p ex q q (4.27)

In the case of mechanical systems, the Lagrangian L is defined in terms of q, q, and t. Rather than following a lengthy process to form the Hamiltonian as an explicit function of q, p, and t, and then differentiating as in (4.26) and (4.27), the canonical equations can be directly obtained from (4.20) and (4.27). Since the system kinetic energy is a quadratic function of the generalized velocities, equations (4.20) and (4.27) directly lead to the following set of equations in matrix form:

130

4. Dynamic Analysis. Mass Matrices and External Forces

p=Mq p = L q + Q ex F T q l

(4.28) (4.29)

where M is the mass matrix, L q = T q V q is the partial derivative of the Lagrangian with respect to the coordinates, Fq is the Jacobian matrix of the constraints, and Qex the vector of applied external forces. The combination of equations (4.28)-(4.29) and the constraints conditions constitutes a system of 2 n+m differential and algebraic equations (DAEs). Although we arrive at n more equations than with the acceleration-based formulation (4.14), p can be obtained explicitly by (4.29). When comparing equations (4.29) and (4.14), one may see that the term (M q) that appears in the Lagrange's equations is not present in their canonical counterparts.
Example 4.5 Repeat Example 4.2 using the canonical equations. The application of equations (4.28) and (4.29) along with the expressions obtained in Example 4.2 directly leads to the following equations: p1 = p2 m 1 q1 q2 m 1 g sin q2 + m 1 q1 + m 2 L g cos q2 0 m 0 p1 = 1 2 2 p2 0 m 2 L + m 1 q1 q1 q2
2

Note how the term related to (M q) that appears in Example 4.2 is not present now.

4.2 Inertia Forces. Mass Matrix

We study in this section the formation of the inertia forces that arise in each of the elements of a multibody system undergoing a given motion. The aim of this study is the formation of the mass matrix of the most common elements that appear in the analysis of multibody systems. The form of these mass matrices will undoubtedly depend on the type of coordinates chosen for the representation of the multibody system, and we will only use the setting provided by the natural coordinates. The inertia forces will be represented by means of equivalent forces that are congruent with the natural coordinates of the element. Of all the analytical methods exposed in Section 4.1, we will use the virtual power method, because it leads to a direct formulation of the inertia forces and avoids the differentiation process inherent in the Lagrange's equations. The Lagrange's equations can lead to rather involved computations in cases when the kinetic energy is position * dependent. As mentioned in Section 4.1.4, if q are the virtual velocities, the virtual power of the inertia forces can be written as the following dot product:

4.2 Inertia Forces. Mass Matrix

131

y P y r ri r

i x

Figure 4.3. Inertial and local system of coordinates in a planar element with natural coordinates.

W * = q*T QI

(4.30)

where QI are the inertia forces with respect to the natural coordinates, and W* is the scalar virtual power. The expressions of the inertia forces for the elements of the planar and threedimensional multibody systems will be developed below. The formulation of the planar element is simpler and will serve as introduction for the three-dimensional cases. In both the planar and spatial cases, the aim is to establish the inertia forces as a product of a matrix (mass matrix) times the acceleration vector, so that the virtual power expression becomes W * = q*T M q (4.31)

where the matrix M is a square, symmetric, and positive definite matrix which depends on the inertia characteristics of the element: mass, position of the center of gravity, and inertia tensor, and in some cases, on the position q of the element as well.

4.2.1 Mass Matrix of Planar Bodies

Consider the planar element shown in Figure 4.3, whose motion is completely defined by that of the basic points i and j. Consider also an inertial coordinates system (x, y) and a moving one (x , y), with its origin fixed at the basic point i and axis x going through point j. The location of a generic point P of the element is defined by a position vector r in the inertial system and r in the moving system, so that

132

4. Dynamic Analysis. Mass Matrices and External Forces

r = ri + A r

(4.32)

where A is the rotation matrix. Since the element is rigid, the local position vector r remains constant with the motion of the element. If the position of an element is characterized by the Cartesian (natural) coordinates of the points i and j, that of the point P will be defined as follows: r = r i + A r = r i + c1 (rj ri) + c2 n (4.33) where c1 and c2 are the components of the vector r in the basis formed by the orthogonal vectors (rjri) and n. The vector n follows the direction of y and has the same module as (rjri). Accordingly, the components of the vector r become x = x i + c1 (x j x i) c2 (yj yi) y = yi + c1 (yj yi) + c2 (x j x i) Equations (4.34) can be expressed in matrix form as follows: 1c1 r= x = y c2 c2 1c1 c1 c2 c2 c1 xi yi xj yj Cq (4.35) (4.34a) (4.34b)

where qT = {xi yi xj yj} is the vector containing the natural coordinates of the element. Note again that the matrix C is constant for a given point P. It does not change with the system's motion or with time. Consequently, r= C q (4.36) r= C q (4.37)

The coefficients c1 and c2 that define the matrix C can be simply expressed as a function of the coordinates of the points i and j in the local reference frame in the following way: r = c1 (r j ri) + c2 n (4.38) Note, however, that ri=0 in the local reference frame, therefore (4.38) becomes r = r j | n c1 X c (4.39) c2 where the vector c contains the coefficients c1 and c2, and the matrix X has as columns the components of the vectors rj and n . Therefore, the following form is taken: x j yj L 0 X= = ij (4.40) yj xj 0 Lij where Lij is the distance between points i and j. The matrix X can always be inverted unless points i and j are coincident. Equation (4.39) can be used to solve for c: c = X 1 r (4.41)

4.2 Inertia Forces. Mass Matrix

133

We now can formulate the virtual power of the inertia forces generated within the element. These can be obtained as the integral, extended to all the elements, of the virtual power of the inertia force of a differential mass located at the generic point P: W* = r
W

r*T r d W

(4.42)

where r is the mass density. Making use of equations (4.36) and (4.37), equation (4.42) becomes W* = r
W

q*T C T C q dW

(4.43)

Since the vectors q*T and q are independent of W , they can be moved out of the integral to yield W =q
* *T

r
W

C C dW q

(4.44)

and comparing (4.31) with (4.44) the mass matrix can be established as M =r
W

C T C dW

(4.45)

Performing the product CTC in equation (4.44), we obtain

2 2 (1c1) + c2 0 (1c1)c1 c2 c2 2 2 2 c2 (1c1)c1 c2 0 (1c1) + c2 2 (1c1)c1 c2 c2 2 c2 (1c1)c1 c2 2 2 c1 + c2 0 2

M=
W

dW

(4.46)

0 2 2 c1 + c2

Note that the equation (4.46) involves the following integrals:

r dW = m
W

(4.47a)

r c dW = X 1
W W

r r dW = m X 1 rG
I I X T = 12 x xy L Ixy Iy

(4.47b)

r c cTdW = X 1
W W

r r r T dW

(4.47c)

where m is the total mass of the element, rG represents the local coordinates of the center of gravity, and Ix , Iy and Ix y are the moments and products of inertia with respect to x and y respectively. The integral in (4.47b) is the static moment or first order area moment which is equal to the mass times the coordinates of the center of gravity. Similarly, the integral in (4.47c) yields the moments of inertia. Substituting the results of (4.47) into (4.46), we obtain the final expression for the mass matrix M:

134

4. Dynamic Analysis. Mass Matrices and External Forces

m 2mx G + Ii2 L ij L ij 0 M= m x G Ii 2 L ij L ij m yG L ij

0 m 2mx G + Ii2 L ij L ij m yG L ij m x G Ii 2 L ij L ij

m x G Ii m yG 2 L ij L ij L ij m yG m x G Ii 2 L ij L ij L ij Ii L ij 0
2

(4.48)

0 Ii L ij
2

where Ii is the polar moment of inertia with respect to the basic point i. The mass matrix defined by equation (4.48) is completely general for any planar element, since any planar element will have at least two basic points with which the mass matrix can be formulated. The mass matrix thus obtained is always constant and this constitutes an important fact.
Example 4.6 Derive the mass matrix of a single bar element of total mass m and length L ij that has the center of gravity located at the middle point. Choosing the points i and j as the end points, one finds that the different parameters that appear in (4.48) take the following values: 2 and the direct application of (4.48) yields xG = L ij ; yG = 0 ; Ii = m L ij 3
2

M=

m 0 m 0 3 6 0 m 0 m 3 6 m 0 m 0 6 3 0 m 0 m 6 3

4.2.2 Mass Matrix of Three-Dimensional Bodies

The determination of the mass matrix of the three-dimensional bodies (or elements) follows a method similar to that used in the planar case. Three-dimensional bodies not only require more complicated algebraic manipulations than the planar ones, but in some cases the mass matrix is not constant and, therefore, requires a special study. The analysis will start with the most general element defined by two points and two non-coplanar unit vectors, which happens to have a constant mass matrix. The mass matrix of the other elements including those

4.2 Inertia Forces. Mass Matrix

135

z z P j r u r i O x x y v y

Figure 4.4. Inertial and local system of coordinates in a spatial element with natural coordinates.

that are not constant, may be derived from this main one by coordinate transformations in which the virtual power method plays a key role. Element with two points and two non-coplanar unit vectors. Consider the element of Figure 4.4 defined by two basic points i and j and two unit vectors u and v. Similar to the planar case, consider an inertial reference frame (x,y,z) and a moving (or local) one (x ,y ,z) rigidly attached to the element that has its origin located at a point O. Again, P represents a generic point of the element, and its location is defined by the position vector r in the inertial frame and r in the local frame. Since the element is rigid, the local relative position vector (r ri) remains constant with the motion of the element. If the position of the three dimensional body is characterized by the Cartesian coordinates of the points i and j and the unit vectors u and v, the position of the point P relative to point i will be defined as follows: r r i = c1 (rj ri) + c2 u + c3 v (4.49)

where c1, c2, and c3 are the components of the vector (r ri) in the basis formed by vectors (rj ri), u, and v. Equation (4.49) may be represented in matrix form as ri x rj r = y = 1c1 I3 c1 I3 c2 I3 c3 I 3 = C qe (4.50) u z v

136

4. Dynamic Analysis. Mass Matrices and External Forces

where I3 is the identity matrix of order (3x3) and qe is the vector of the natural coordinates. As in the planar case, the matrix C is independent of the system's motion and therefore remains constant with time. Again, the following relations are obtained: r = C qe (4.51) r = C qe (4.52)

The coefficients c1, c2, and c3 that define C can be expressed as a function of the coordinates of the points i and j and unit vectors u and v in the local reference frame in the following way: r r i = c1 (r j ri) + c2 u + c3 v Equation (4.53) expressed in matrix form becomes c1 r ri = r j ri | u | v c2 = X c c3 (4.53)

(4.54)

where the vector c contains the coefficients c1, c2, and c3, and the matrix X has as columns the components of the vectors (r j ri), u , and v . As in the planar case, the matrix X can always be inverted (provided (r j ri), u , and v are noncoplanar). Equation (4.54) can be used to solve for c: c =X
1

(r ri)

(4.55)

We have gathered all the information necessary to formulate the virtual power of the inertia forces generated within the spatial element. The integral form of the virtual power becomes W* = r
W

r*T r dW

(4.56)

where r is the mass density. Making use of equations (4.51) and (4.52), equation (4.56) becomes W =r
V

*T

C C q dV = q

*T

r
V

C C dV q

(4.57)

where again the vectors q*T and q, which are independent of V, have been moved out of the integral. Comparing (4.31) with (4.57) we obtain the expression for the mass matrix M =r
V

CT C dV

(4.58)

The substitution of C into equation (4.57) yields

4.2 Inertia Forces. Mass Matrix

137

i
Figure 4.5. Element with three basic points and a unit vector.

M=
V

(1c1)2 I3 (1c1 ) c1 I 3 (1c1 ) c2 I3 (1c1) c3 I3 2 c1 c2 I3 c1 c3 I3 (1c1 ) c1 I3 c1 I3 dV 2 c2 I 3 c2 c3 I3 (1c1 ) c2 I3 c1 c2 I3 2 (1c1) c3 I 3 c1 c3 I3 c2 c3 I 3 c3 I3

(4.59)

The integration over the element of the product CTC involves the following integrals:

r dV = m
V 1 V 1 V 1

(4.60a)

r c dV = X
V

r ( r ri) dV = m X
T

(r G ri) m a
1 T

(4.60b)

r c cTdV = X
V

r ( r ri) (r ri) dV X

=X

Ji X

(4.60c)

where m is the total mass of the element, and rG represents the local coordinates of the element's center of gravity. The matrix Z can be formed from the matrix Ji which contains all the information about the moments and products of inertia of the element as in the planar case. Substituting the results of (4.60) into (4.49), the following expression for the mass matrix is obtained: m 2ma1 + z 11 I 3 ma1 z 11 I 3 ma1 z 11 I 3 z 11 I 3 ma2 z 21 I 3 z 21 I 3 ma3 z 31 I 3 z 31 I 3 ma2 z 12 I 3 z 12 I 3 z 22 I 3 z 32 I 3 ma3 z 13 I 3 z 13 I 3 z 23 I 3 z 33 I 3

M=

(4.61)

The mass matrix defined by equation (4.61) is a constant and symmetric matrix formed by diagonal sub-matrices of size (33). This matrix depends on a minimum set of ten different values, since the inertia characteristics of a three-

138

4. Dynamic Analysis. Mass Matrices and External Forces

k l

Figure 4.6. Element with four basic points.

j u

i
Figure 4.7. Element with two basic points and a unit vector.

dimensional solid depend on ten parameters: the mass, the coordinates of the center of gravity, and the six different elements of the inertia tensor. Other Elements with Constant Mass Matrix. Consider the element of Figure 4.5, determined by three non-aligned basic points and by one unit vector not contained in the plane determined by the three points. This element may actually be made similar to the one considered in the previous section (two points and two unit vectors), by simply defining the missing unit vector v at point j as the vector that goes from j to k. Thus v = r k rj / L jk (4.62)

Now we can use (4.61) to obtain the virtual velocity and acceleration vectors: ri rj u v =
I3 0 0 0 0 I3 0 0 0 0 1 I3 1 I3 L jk L jk 0 0 I3 0

ri rj rk u V

ri rj rk u (4.63)

where the transformation matrix V is a (1212) constant matrix. Equation (4.63) may be expressed as q = V qnew (4.64)

4.2 Inertia Forces. Mass Matrix

139

where q represents the components of the two points and two unit vectors, and qnew are those of the new element represented by three points and one unit vector. Similarly, the relation between the virtual velocities will be q * = V q* new (4.65)

The expression for the virtual power given by equation (4.31) can be used to obtain the mass matrix of the element at hand. This can be done by simply substituting (4.64) and (4.65) into (4.31) to yield
* T W * = q *T M q = q*T new V M V qnew

(4.66)

and, therefore, the new mass matrix becomes M new = V T M V (4.67)

Since V is constant the new mass matrix will also be constant. A similar situation arises with the type of element depicted in Figure 4.6, which can also be made equivalent to an element with two points and two noncoplanar unit vectors. In fact, the two unit vectors u and v located at points i and j may be defined as u = rk ri / L ik v = rl rj / L jl as: ri rj u = v I3 0 1 I3 L jk 0 0 I3 0 1 I3 L jk 0 0 1 I3 L jk 0 0 0 0 1 I3 L jk ri rj rk = V q rl (4.68a) (4.68b)

Equations (4.68) allow for the definition of a new transformation matrix V

(4.69)

The matrix V is again constant and, therefore, it also relates virtual velocities and accelerations as in (4.64) and (4.65). The new mass matrix is thus obtained by the same equation (4.66) with V defined by (4.69). Element with two points and one non-aligned unit vector. Figure 4.7 shows an element with two basic points and one unit vector not aligned with them. The coordinate transformation between the standard element consisting of two points and two non-coplanar unit vectors, and the new one is not as simple to establish as with the elements considered in the previous section. In the element of Figure 4.7, the vector v cannot be constructed non-coplanar with u by just a linear combination of the position vectors of points i, j, and u. A geometric transformation is needed in this case, and the method of constructing the vector v is to start from the following cross product of vectors:

140

4. Dynamic Analysis. Mass Matrices and External Forces

v = c r i rj u

(4.70)

where c is a constant that makes the vector v have a unit module. The differentiation of equation (4.70) with respect to time yields v = c r i rj u + c r i rj u (4.71)

A transformation matrix V for velocities can now be constructed from (4.71) as follows: ri I3 0 rj 0 I3 = 0 0 u c u c u v 0 0 I3 c rij ri rj = V q new u

(4.72)

where u and rij are skew-symmetric matrices of order (3x3) that correspond to the cross product of vectors obtained with u and (rirj), respectively. These matrices are written as follows: 0 uz u y u = uz 0 u x uy u x 0 0 zizj y i y j rij = zizj 0 xixj y i y j xixj 0 (4.73)

(4.74)

Equation (4.72) can be differentiated with respect to time so as to obtain the corresponding equation for the accelerations. However, the matrix V defined by equation (4.72) is no longer constant and, therefore, the differentiation process is different from those appearing in the previous section: ri I3 rj 0 = 0 u c u v or q = V q new + V q new The relation between the virtual velocities is q = V qnew
* *

0 0 I3 0 0 I3 c u c rij

rj ri u

I3 + 0 0

0 I3 0

0 0 I3

c u c u c rij

ri rj u

(4.75)

(4.76)

(4.77)

By substituting expressions (4.76) and (4.77) into the virtual power expression (4.31), we obtain W =q
* *T

M q = q new V M ( V q new + V q new)

*T

(4.78)

4.2 Inertia Forces. Mass Matrix

141

j P i s
Figure 4.8. Element with two basic points.

where M is the constant mass matrix corresponding to the element with two points and two unit vectors. This equation contains the new mass matrix Mnew= V TM V and velocity-dependent inertia terms defined by V TM V q new to be added to the external forces in the right-hand side of the equations of motion. Note the simplicity by which these forces are obtained using the virtual power method in which only the computation of V is required, as compared to the Lagrange's equations. The latter would have required the differentiation of the new mass matrix with respect to both time and q. Element with two basic points. The use of an element with two basic points only makes sense when it is used: first, to maintain a constant distance between two points; and second, when its dimensions, other than the length, are negligible. As in the planar case, the motion of an element of this type is fully characterized by the motion of the two basic points (See Figure 4.8). Any rotation around the element axis is disregarded, an assumption that is valid only if the moment of inertia about that axis is negligible. The position vector of a generic point P can be defined in terms of the variable s (Figure 4.8) as: r = ri L s + rj s L L which may be expressed in matrix form as: ri r = Ls s = N (s) q L L rj (4.79)

(4.80)

Since the matrix N is independent of the motion, the velocity and acceleration vectors will be given by r = N( s ) q r = N (s) q Consequently, the virtual power of the inertia forces is in this case:
L

(4.81) (4.82)

W* =
0

*T

r r A ds

(4.83)

142

4. Dynamic Analysis. Mass Matrices and External Forces

a)

b)

Figure 4.9. a) Element with two basic points and two coplanar unit vectors, b) Element with four basic points and four coplanar unit vectors.

where A is the cross-sectional area of the bar. By substituting the results of expressions (4.81) and (4.82) in equation (4.83), the following expression is obtained: W =
0 * L

*T

N r A N q ds = q

*T

L 0

N r A N ds q

(4.84)

Therefore the mass matrix is

L

M=
0 L

N T N r A ds =
L

I3
0

Ls L

r A ds

I3
0

Ls s r A ds L L
L

(4.85)

=
L

I3
0

Ls s r A ds L L

I3
0

s 2 r A ds L2

where I3 is the unit (3x3) matrix. Assuming a constant density and cross sectional area for the element, the following final expression for the mass matrix is obtained: m I3 m I3 3 6 M= (4.86) m I3 m I3 6 3 Mass matrix of other three-dimensional elements. All the three-dimensional elements that may appear in practice, belong either to one of the groups studied previously or contain a set of points and vectors whose mass matrix is known. For example, the body in Figure 4.9b contains four basic points and four unit vectors. Fortunately, one does not have to worry about calculating mass matrices of elements as complicated as this one. It will be sufficient to take any two points with which two non-coplanar unit vectors are associated and attribute the corresponding mass matrix to the subset of points and unit vectors given by

4.2 Inertia Forces. Mass Matrix

143

equation (4.61). The body's system of local coordinates will be located according to the selected points and unit vectors. The inertia properties of the element, such as center of gravity and inertia tensor, must be referred to those coordinates. The element shown in Figure 4.9a, contains two basic points and two unit vectors that are coplanar. For this reason, the mass matrix of equation (4.61), which assumes that the unit vectors are non-coplanar, is not applicable. Therefore, one of the two unit vectors must be selected and used in the virtual power equation (4.78) corresponding to an element with two points and only one unit vector. Irrespective of the geometry and number of basic points and unit vectors of any element of a multibody system, a subset of points and vectors can always be found whose mass matrix corresponds to one of those calculated in the previous sections.

4.2.3 Kinetic Energy of an Element

The formation of the mass matrices of different elements has been studied in the previous sections through the application of the principle of virtual power. Expressions for the kinetic energy of a body may be convenient, either because it is of direct interest to evaluate the energy, or because one wishes to formulate the equations of motion through another type of formulation, such as the method of Lagrange. The way of formulating the kinetic energy is rather similar to that used for the virtual power of the inertia forces. The kinetic energy is defined by: T=1 2
v

rT r r dV = 1 r 2

q T C T C q dV
v

(4.87)

Since the natural velocities are independent of V, they can be moved out of the integral to yield 1 T (4.88) q Mq 2 where q are the natural velocities of the element. The matrix M given by equation (4.61) corresponds to the element that has two non-coplanar points and two unit vectors. In the other cases, the mass matrix may be formed by the coordinate transformation matrix V TM V. The use of the Lagrange's equations for the formulation of the equation of motions with non-constant mass matrices will lead to the differentiation process explained in Example 4.1, which may become involved in those cases where the mass matrix is coordinate-dependent. The kinetic energy is then given by T= T= 1 T T( q ) q V M V( q ) q 2 (4.89)

144

4. Dynamic Analysis. Mass Matrices and External Forces

j v1 fP i P

u1

Figure 4.10. Concentrated force fP at point P.

4.3 External Forces

We deal in this section with the formulation of external forces, such as concentrated loads and moments that are applied to the elements of the multibody system, and those generated from specific types of elements, like translational and rotational springs. Although the forces generated by springs are not really external, it is customary in multibody dynamics to include them as applied external forces, as opposed to structural mechanics and the finite element formulation in which those forces are considered as internal forces acting by means of the stiffness matrix.

4.3.1 Concentrated Forces and Torques

This section with the analysis of concentrated forces and torques. The problem is the representation of those forces and torques in terms of the natural coordinate system used to represent the motion of an element. Concentrated Forces. When a concentrated force is applied at a point P of an element that is not a basic point (See Figure 4.10), a transformation needs to be introduced that will transform the forces at point P to the natural coordinates (basic points and unit vectors) of the element. A typical case is the gravity force, applied at the center of gravity of an element, which needs to be referred to the natural coordinate system used to represent the motion of the given element. In order to set up the required transformation, we can directly use equation (4.50) applied to point P, so that rP = C P qe (4.90) where qe represents the natural coordinates of the element. The coefficients that form the matrix CP can be obtained using equation (4.55), such that c = X P (r P ri)
1

(4.91)

4.3 External Forces

145

rj M f rk

uF f ri

Figure 4.11. Concentrated torque M acting on element (i-j).

The vector r P contains the coordinates of P in the local frame. The matrix CP in equation (4.90) acts as a coordinate transformation matrix that may also be used to transform the forces fP into equivalent forces, Qe e x, expressed in terms of the natural coordinates of the element. We use for this purpose the principle of virtual work and impose the condition that both sets of forces perform an equal amount of virtual work. Accordingly,
eT dW = drT Q ex P fP = d q e

(4.92)

Since

d rT P

= dq

eT

T C P,

equation (4.92) becomes

dW = dq eT C P f P

(4.93)

Comparing equations (4.92) and (4.93), we obtain the following equation for the force transformation: Qex = C P fP
e T

(4.94)

The potential of this concentrated external force fP, acting at point P, is defined by the expression:
q q

V =
qo

drT P f q =
qo

dq T C P f q =
qo

dq T Q

(4.95)

that is also valid for the case in which the force depends on the position. The differential term dq has been intentionally placed in the left side of the integral, because this order leads to simpler and more congruent expressions for the derivatives of V that are calculated in Chapter 6. Concentrated Torques. The case of a concentrated torque can be dealt with in a similar way as that for the concentrated force except for the addition of a preceding and important transformation. In basic statics, any torque M may be replaced by an equivalent pair of forces, f and f, of equal magnitude and opposite directions, acting on a plane perpendicular to the direction of M, and separated by a vector d. The result is: M=df.

146

4. Dynamic Analysis. Mass Matrices and External Forces

Figure 4.11 shows a torque M acting on an element with basic points i and j, which is substituted by a pair of forces f and f applied at the beginning and end of a unit vector uf with origin at ri. This unit vector is defined by: uf = and f = uf M (4.97) Note that vector uf does not belong to the natural coordinates vector qe. The forces f and f can now be treated as concentrated forces following the results of the previous section. Accordingly, the virtual work expression becomes ( rj ri) M ( rj ri) M (4.96)

d W = d q eT (Ci f C i+ u f f)

(4.98)

Therefore the equivalent generalized force with respect to the natural coordinates finally results in Q e x = (Ci C i+ u f) f
e T T

(4.99)

where C i is very simply obtained in this case because i is a basic point. It will be shown that the force f given by equation (4.97) is the force variable conjugated with the displacement variable uf. The potential of the torque M can be calculated as the sum of the potentials of forces f and f. This potential is
q

V =
qo

dri f dri + du T f f =+
qo

du T f f =
qo

du T f M uf

(4.100)

The forces f and f can be treated as concentrated forces following the results of the previous section. Accordingly, the expression for the potential becomes
q q

V =
qo

T T drT i f dri + du f f =
qo

dqe (Ci C i+ uf) f

(4.101)

4.3.2 Forces Exerted by Springs

Springs are elements capable of storing elastic potential energy and of exerting forces that are a function of their positions. In addition, springs play an important role in all but the kinematic problems. We study in this section the forces exerted by both translational and rotational springs as well as the potential energy stored in them. As in previous sections, the study will be conducted within the context of planar and three-dimensional natural coordinates. Translational Springs Between Basic Points. Consider the translational spring shown in Figure 4.12 which connects the basic points i and j. Let L ij and L 0 be the deformed and undeformed lengths of the spring, respectively. When the

4.3 External Forces

147

f j

y i f
Figure 4.12. Translational spring with mixed coordinates.

spring is stretched (or compressed), it exerts a force between the basic points in the directions of the spring that is a function of the elongation. The value of the force is given by f = k L ij Lij L 0 (4.102)

where k L ij is the stiffness of the spring which will have a constant value k, if the spring is linear. In the planar case, according to Figure 4.12, the force vector that acts on the basic points i and j can be expressed as follows: Q ix cos y Q iy sin y Q= =f = f cos y Q jx L ij sin y Q jy and in the three-dimensional case:
T Q = f L ij

xi yi xj yj

xj yj xi yi

(4.103)

x ix j

y i y j

zizj

x jx i

y j y i

zjzi

(4.104)

Both f and L ij are functions of the natural coordinates. If the spring is linear equation (4.104) becomes Q = k 1 Lo L ij
T

x ix j

y i y j

zizj

x jx i

y j y i

zjzi

(4.105)

It is important at times to evaluate the potential energy stored in spring elements (See Chapter 6). The potential energy stored by a translational spring is equal to the integral of the force times the differential extension of the spring, between the non-deformed state and the final deformed configuration. Considering the coordinates of points i and j:

148

4. Dynamic Analysis. Mass Matrices and External Forces

u2 k 2 l 2 i 1 u1 1 j
Figure 4.13. Translational spring between any two points.

v2

v1

LLo

V=
o

dq T Q =

LLo

f dl
o

(4.106)

When the spring is linear, the integration of equation (4.106) yields V = 1 k Lij L o 2
2

=1 k 2

x jx i

+ y j y i

+ zjzi

Lo

(4.107)

Translational Springs With Relative Coordinates. Equations (4.103) and (4.104) define the spring forces in terms of the natural coordinates. However, this formulation can be greatly simplified if the relative distance s between the points i and j is introduced as a new dependent (mixed) coordinate through the following constraint condition: xi xj
2

+ yi y j

+ zi zj

s2 = 0

(4.108)

The force f can be directly entered into the formulation as the conjugate variable of distance s with a value: f = k(s) s s 0 (4.109)

where so= Lo, and s is equal to the deformed length. Proceeding with this mixed type of coordinate representation, the formulation of the forcing terms becomes much simpler at the expense of increasing the number of dependent variables by one variable.

4.3 External Forces

149

If the dependent coordinate s is used, the potential energy of the spring is simply given by
ss o ss o

V=
o

f ds =
o

k s sso ds

(4.110)

Translational Springs Between Any Two Points. When the origin and end points of the spring are not basic points but any two points corresponding to two different elements, equations (4.104) and (4.107) cannot be used. As shown in Figure 4.13, it is necessary to construct the position vectors of the origin and end of the spring, starting from the local coordinates of these points at the local coordinates frames attached to the elements to which they belong. In the case of Figure 4.13, we can write: r 1 = r i + A 1 (r 1 ri) r 2 = r k + A 2 (r 2 rk) (4.111) (4.112)

where A1 and A2 are rotation matrices that depend on the natural coordinates of the points and vectors of each element. In computing these rotation matrices, we can express the local and global coordinates of the basic points and vectors that belong to a rigid body as the columns of a (3x3) matrix X as follows: X rirj u v = A X = A rirj u v Hence, the rotation matrix A can be found as A =X
1 1 X = rirj u v X

(4.113)

(4.114)

Consequently, the rotation matrices corresponding to the rigid bodies to which points 1 and 2 belong will be defined by A 1 = X 1 X 1 = r i r j u 1 v 1 A 2 = X 2 X 2 = rk rl u 2 v 2 Using the result of equation (4.90), we can write r 1 = C 1 qe 1 r 2 = C 2 qe 2 (4.117) (4.118)
1 1

X1

(4.115) (4.116)

X2

Matrices C 1 and C 2 , defined by equation (4.50), are constant matrices that permit finding the global coordinates of points r1 and r2 in terms of the natural coordinates of the elements to which they belong. Using matrices C 1 and C 2, we can obtain two expressions similar to (4.94): Q1 = C 1 f1 Q2 = C 2 f2
e T e T

(4.119) (4.120)

150

4. Dynamic Analysis. Mass Matrices and External Forces

y i
Figure 4.14. Rotational spring with mixed coordinates.

e where f2 = f1. The generalized forces Qe 1 and Q 2 are conjugated with the virtual displacements dr1 and dr2, that can be defined as functions of dq by equations (4.117) and (4.118). It is possible to write an expression for the potential energy analogous to equation (4.106):

V=

k L 1 Lo . C 1 q1 C 2 q2 C 1 q 1 + C 2 q2 L qo
T T

C 1 d q1 C 2 d q2

(4.121)

and expanding the product of vectors, V =

T k L 1 L o . qT 1 q2 L qo q

C 1 C 1
T C 2 T T T C2

C1 0 0 C2
T

d q1 = d q2
q1 q2 (4.122)

k ( L ) 1 L o d q1 d q2 L qo
T

C 1 C 1 C 1 C 2
T C 2

C1

T C2

C2

The length L that is the distance between points r1 and r2 can be computed by using the formulae (4.117) and (4.118) to solve for the coordinates of r1 and r2; and then to find the distance directly. Rotational Springs. A rotational spring exerts, between the elements to which it is connected, a torque about the common articulation of both elements; that is a function of the relative angle twisted between them (Figure 4.14). Angles of more than 360o are possible; therefore, it is necessary to take into account not only the angle between both elements as calculated with the scalar or cross products of vectors, but also the number of complete turns that the rotational spring may have gone through.

4.3 External Forces

151

f1 k y1 y y2 i f2 j' u j f2

k' f1

Figure 4.15. Special case of rotational springs.

Consider the planar rotational spring shown in Figure 4.14 in which y 0 is the angle corresponding to the non-deformed position of the spring. If the relative angular position of the elements is considered as a new mixed variable (y+2np), the torque exerted on both elements will be directly given by M = k(f) y + 2np y 0 (4.123) Note that the new dependent coordinate y may be introduced by either one or a combination of the following constraint conditions: r ij r jk = Lij L jk cos y rij rjk = Lij L jk sin y In addition, the potential energy is given by
y +2 n p y o y +2 n p y o

(4.124) (4.125)

V=
yo

M y dy =
yo

k y yyo dy

(4.126)

Since y is one of the dependent coordinates of the system, no additional transformation is required. In the three-dimensional case, the situation is a little more complicated, particularly in the case in which the points i, j, and k are not in a plane perpendicular to the axis of the revolute joint. Here, the angle y in the joint is not the angle formed by the segments (ij) and (ik), but the angle determined by two straight lines normal to the axis of the pair. This can be seen in Figure 4.15. It is assumed in this figure that the axis of the pair is determined by the unit vector u. A similar and simpler formulation for the forces is obtained when the pair is determined by two basic points.

152

4. Dynamic Analysis. Mass Matrices and External Forces

If the angle is introduced as a dependent coordinate, the expression for the potential energy is immediate and responds to the same equation as that of the planar case (equation (4.126)). When the spring is linear, V = 1 k y + 2np y o 2 where k is the spring constant.
2

(4.127)

4.3.3 Forces Induced by Known Acceleration Fields

The simplest case of forces induced by known acceleration fields are gravitational forces. The gravitational force acting on an element is simply the product of its mass m times the gravitational acceleration g, acting on the center of gravity rG: f =m g and the potential of this force can be expressed as V = m rT Gg (4.129) (4.128)

A matrix C G similar to that of expression (4.50) can be constructed for the center of gravity to express its coordinates in terms of natural coordinates qe of the element. Consequently the potential becomes V = m q eT C G g
T

(4.130)

Another important case, is the one in which the external forces originate from a known accelerations field. This situation arises when the fixed element is moving in a prescribed mode or when the entire multibody system is subjected to a rotation. Let vo and W be the velocity of the origin and the angular velocity vector, respectively, of the reference coordinate frame whose motion is known. The acceleration of this system is defined by v o , W , and W which are known. Using the principles of relative motion (Greenwood (1988)), the motion with respect to a moving reference frame can be studied as if it were an absolute motion, introducing as known external forces all the inertia forces corresponding to the motion of the frame. Assuming a standard element with two basic points i and j and two non-coplanar unit vectors u and v, these accelerations are: ri = W ri + W W ri + v o (4.131) rj = W rj + W W rj + v o u=W u + W W u v=W v + W W v As a consequence of these accelerations, the vector of inertia forces on the element is, (4.132) (4.133) (4.134) Qe in acting

4.3 External Forces

153

ri Q in = M
e e

rj u v

=M q

(4.135)

where Me is the mass matrix developed in Section 4.2.2. The potential of these forces, which is position dependent, will be
q

V =
qo

dq

eT

e Qin

=
qo

d q eT M q

(4.136)

One will find these expressions useful for carrying out the dynamic analysis of a multibody system evolving in a known field of centrifugal forces.

References
Bastero, J.M. and Casellas, J., Curso de Mecnica, Eunsa, (1976). Goldstein, H., Classical Mechanics, 2nd Edition, Addison-Wesley, London, (1980). Greenwood, T. J., Principles of Dynamics, 2nd Edition, Prentice Hall, (1988). Hamilton, W.R., "On the General Method in Dynamics", P h i l o s o p h i c a l Transactions of the Royal Society of London, pp. 247-308, (1834). Huston, R.L., Multibody Dynamics, Butterworth-Heinemann, (1990). Kane, T.R. and Levinson, D.A., Dynamics: Theory and Applications, McGraw-Hill, (1985). Reddy J.N., Energy and Variational Methods in Applied Mechanics, Wiley Interscience, (1984). Schiehlen, W.O., "Dynamics of Complex Multibody Systems", SM Archives, Vol. 9, pp. 159-195, (1984).

Problems
4/1 Using equation (4.48), find the inertia matrix (with respect to points i and j) of the 2-D element shown in the figure, in the following cases: a) The element has a concentrated unit mass located at i. b) The element has a concentrated unit mass located at j. c) The element has a concentrated unit mass located at point (0,1) d) The element has a concentrated unit mass located at point (0,1) e) The element consists of a disk with its center at i, having unit radius, and uniformly distributed unit mass. 4/2 Using the results of Problem 4/1 and admitting the possibility of inertia matrices corresponding to negative masses, to eliminate mass from a real rigid body,

154

4.Dynamic Analysis. Mass Matrices and External Forces find the concentrated masses at points (0,0), (1,0), (0,1), and (0,1), so that the resulting mass matrix is: 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0

i 1
Figure P4/1. 4/3

x
i

Figure P4/3.

Find by integration of
L

Mqm r(s) ds L 0 the inertia matrix M of an homogeneous 3-D bar and show that M is independent of the position of the bar.

y z 1 2 4 y z x 3
Figure P4/4. Figure P4/5.

j 1 i 3 2

4/4

The inertia properties of a 3-D rigid body depend on ten parameters including the mass m , the position of center of gravity rG, and the inertia tensor IG, defined on a moving frame attached to the body. Show by intuitive reasoning (without al-

Problems

155

gebraic calculations) that it is possible to find a (12x12) constant inertia matrix in the global frame. Use as acceleration variables the Cartesian accelerations of points 1, 2, 3, and 4. 4/5 A bar 1 connects two rigid bodies 2 and 3 through two spherical joints i and j. The motion of bodies 2 and 3 is known, so the motion of the bar is known except for the rotation around the direction i-j. There is a friction torque of constant magnitude T and direction opposed to the relative angular velocity applied at the joints i and j. Find the directions of torques T i and T j that guarantee the complete equilibrium of the bar.